Abstract

Although intracellular electrophysi-ology is a powerful tool for studying synaptic integration and plasticity, it is limited by the small number of neurons that can be recorded simultaneously in Vitro and by the technical difficulty of intracellular recording in vivo. One way around these difficulties may be to use large-scale extracellular recording of spike trains and apply statistical methods to model and infer functional connections between neurons. These techniques have the potential to reveal large-scale connectivity structure based on the spike timing alone. However, the interpretation of functional connectivity is often approximate, since only a small fraction of pre-synaptic inputs are typically observed. Here we use in Vitro current injection in layer 2/3 pyramidal neurons to validate methods for inferring functional connectivity in a setting where input to the neuron is controlled. In experiments with partially-defined input, we inject a single simulated input with known amplitude on a background of fluctuating noise. In a fully-defined input paradigm, we than control the synaptic weights and timing of many simulated presynaptic neurons. By analyzing the firing of neurons in response to these artificial inputs, we ask 1) How does functional connectivity inferred from spikes relate to simulated synaptic input? and 2) What are the limitations of connectivity inference? We find that individual current-based synaptic inputs are detectable over a broad range of amplitudes and conditions. Detectability depends on input amplitude and output firing rate, and excitatory inputs are detected more readily than inhibitory. Moreover, as we model increasing numbers of presynaptic inputs, we are able to estimate connection strengths more accurately and detect the presence of connections more quickly. These results illustrate the possibilities and outline the limits of inferring synaptic input from spikes.

Author Summary

Synapses play a central role in neural information processing — weighting individual inputs in different ways allows neurons to perform a range of computations, and the changing of synaptic weights over time allows learning and recovery from injury. Intracellular recordings provide the most detailed view of the properties and dynamics of individual synapses, but studying many synapses simultaneously during natural behavior is not feasible with current methods. In contrast, extracellular recordings allow many neurons to be observed simultaneously, but the details of their synaptic interactions have to be inferred from spiking alone. By modeling how spikes from one neuron, statistically, affect the spiking of another neuron, statistical inference methods can reveal “functional” connections between neurons. Here we examine these methods using neuronal spiking evoked by intracellular injection of a defined artificial current that simulates input from a single presynaptic neuron or a large population of presynaptic neurons. We study how well functional connectivity methods are able to reconstruct the simulated inputs, and assess the validity and limitations of functional connectivity inference. We find that, with a sufficient amount of data, accurate inference is often possible, and can become more accurate as more of the presynaptic inputs are observed.

Introduction

Changes in these transformations underlie learning, memory, and recovery from injury [4,5]. Tools for identifying synaptic weights and tracking their changes, thus, play a key role in understanding neural information processing. Traditionally, synaptic integration and plasticity are studied using intracellular recordings in vitro, where synaptic weights can be directly measured as the amplitude of postsynaptic potentials or currents. Although there are singular studies employing simultaneous intracellular recordings from several neurons in vivo [6—8], recording intracellularly from connected neurons in vivo is technically prohibitive. On the other hand, methods for recording extracellular spike trains are advancing at a rapid pace [9,10] and allowing the simultaneous recording of hundreds of neurons. Estimation of synaptic interactions from extracellularly recorded spike trains requires development of sensitive data analysis tools. Although strong synapses are usually readily detectable using cross-correlation analysis [11—17], where they appear as asymmetric, short latency peaks on cross-correlograms [18,19], in general, it is difficult to link the statistical relationships between spike trains to specific synaptic processes [20,21]. Here we provide empirical tests of statistical tools for such analysis using in vitro current injection where the true synaptic input is known.

Although cor-relational methods have long been used to study interactions between pairs of neurons [18,19], recent work has shown that statistical inference methods may be able to substantially improve our ability to detect neuronal connectivity and predict neural activity [24—26]. These model-based methods [22,27,28] are important in removing the confounds that occur with simultaneous recordings [20,29] and have revealed highly structured functional interactions, that accurately reflect the known circuit architecture, in the retina [30] and invertebrate systems [31]. However, it has proven difficult to relate functional connectivity reconstructed from spikes to the known anatomy and physiology of cortical connectivity [26,32—34]. Sparse sampling of neurons and large electrode spacing may contribute somewhat to the difficulty in interpreting the results of functional connectivity analyses of cortical circuits, but it is also unclear what information these inference methods can provide about actual synaptic inputs and what limitations there are to the use of these methods in general.

We use in vitro intracellular recordings from layer 2/3 pyramidal cells in slices from rat neocortex as they respond to simulated, current-based presynaptic input. The fully-defined input is composed of excitatory and inhibitory postsynaptic currents produced by firing of large number of simulated presyn-aptic neurons. Since we know both the spike timing of the input presynaptic neurons and spike timing of the postsynaptic neuron, we can examine the limits of functional connectivity inference. We ask how well synaptic inputs of different amplitudes can be detected, how much data is necessary to reconstruct the amplitudes of excitatory and inhibitory synaptic inputs, and how precisely synaptic weights can be estimated from spikes alone. We briefly examine the feasibility of tracking changes in synaptic weight over time. Finally, we examine how accurately the firing of the postsynaptic neuron can be predicted from presynaptic spiking, and to what extent knowledge of multiple presynaptic inputs improves the accuracy of spike prediction.

Results

Here we examine the relationship between simulated synaptic input and functional connections estimated from spikes using in vitro current injection experiments.

Detection of artificial EPSCs immersed in fluctuating noise

We injected into layer 2/3 pyramidal cells current consisting of three components: simulated, artificial excitatory postsynaptic current (aEPSC) from a single presynaptic neuron, fluctuating noise with standard deviation 0, and a DC offset (Fig. 1A). We adjusted the gain of the injected fluctuating current to produce membrane potential fluctuations with ~ 15—20 mV peak to peak amplitude and DC current to achieve postsynaptic spiking ~ 5Hz. After this adjustment, the overall standard deviation of the scaled fluctuating current awas between 70 and 110pA. On top of this fluctuating noise, we injected aEPSCs produced by the single presynaptic input and varied the amplitude of aEPSCs in different realizations. We then analyzed the postsynaptic responses aiming 1) to examine whether input can be detected based on spikes alone, 2) to quantify how much data is necessary to detect a synaptic input of a given strength, 3) to quantify how much data is necessary to detect changes in input strength, and 4) to determine how accurately such pairwise models describe and predict spiking of the postsynaptic neuron.

1C). Here we use a common model for estimating functional connectivity from point-process observations: a generalized linear model (GLM). We assume that postsynaptic spiking is generated by a Poisson process with a rate determined by a baseline firing rate, the recent history of the neuron’s firing, as well as input produced by presynaptic spikes (see Methods for details).

Specifically we compare a model that predicts postsynaptic spikes based only on the recorded neuron’s spike history (Model 1) with a model that uses both spike history and presynaptic input (Model 2). These models accurately capture two different aspects of postsynaptic spiking: the spike history term captures the fact that immediately after an action potential, the probability of spike generation decreases, while the coupling terms capture the variable (excitatory, in these experiments) effect of the presyn-aptic input (Fig. 2A). In Model 2 the estimated presynaptic inputs correlate well with the actual amplitude of EPSCs (Fig. 2A and 2E).

The model does not aim to distinguish between the underlying processes that determine spike timing. For instance, in vivo results show that, although absolute refractoriness lasts only few milliseconds, generation of an action potential can influence the spike threshold for up to 1s [35,36]. Similarly, although postsynaptic potentials (PSPs) persist only over 10’s of ms, the effect of an EPSP on postsynaptic spike timing may last substantially longer [36,37]. Here, the estimated post-spike history and coupling effect in the models persist over 50—200ms. Both the history-only Model 1 and Model 2 with constant coupling capture statistics of the spike trains and predict spikes with reasonable accuracy. Model 2 provides an increasingly accurate prediction as the input amplitude increases (Fig. 2B). Using ROC analysis (1ms timescale, 10-fold cross-validated), we find that spikes are predicted with an area-under-the-curve (AUC) of 0.71:0.01 for Model 1, while Model 2 yields 0.72:0.01, 0.74:0.01, and 0.77:0.01 as the input amplitude increases from 0.5 a, to 1.00 a, to 1.5 a. The relatively high accuracy of history-only Model 1 predictions reflects the fact that these neurons show substantial regularity in their ISIs (average CV = 0.54:0.02). However, Model 2 provides better spike prediction (AUC) in all cases (paired t-tests, p = 0.01, p<10_5, p<10_3). Model 1 accounts for auto-correlations fairly well, but the spike history alone cannot account for the cross-correlations between pre and postsynaptic spikes (Fig. 2C, black lines). Model 2, using the spike history and coupling parameters together, captures both the auto and cross-correla-tions present in the observed spiking (Fig. 2C, red lines). Consistent with results of cross-corre-lation and area-under-the-curve analysis, the model with coupling provides a better fit to the data than the spike-history alone model when the amplitude of added synaptic input is 0.5 a or larger. For a fixed amount of data (recording length 200s), modeling a synaptic input with amplitude < 0.50 does not improve prediction of the postsynaptic spikes over a spike history-only model (Fig. 2D). Nevertheless, even for these weaker inputs (0.1 a, 0.2 a, 0.25 a, and 0.3 a) the A M1: No Coupling Model 2: Constant Coupling B ROC Post-Spike Gain Post-Spike Gain Coupling Gain Gain Time [ms Time [ms

n

Observed 8 Model 1 Model 2

U

Statistical inference of synaptic connections of different strength from spike trains. A) Model fits for a typical neuron (top), and results for N = 10 cells (bottom). Model 1 (M1, N0 Coupling) describes only the effects of prior spikes in the postsynaptic cell on the generation of future spikes (Post-Spike Gain). Model 2 (M2, Constant Coupling) describes both the post-spike history effects (Post-Spike Gain) as well as the coupling effects following a pre-synaptic spike (Coupling Gain). Gray lines show parameters estimated from bootstrap samples; solid colored lines show their averages. B) Model accuracy: Receiver operating characteristic (ROC) curves for the example cell and area under the curve (AUC) for all cells from A. Curves show the cross-validated false positive rate (FPR) vs true positive rate (TPR) for spike detection in 1ms bins. Error bars denote standard deviation across cells. Note that the accuracy of Model 2, with constant-gain coupling, (M2) increases with increasing input amplitude. C) Cross and auto-correlations forthe typical cell with the model predictions. D) Detectability of synaptic connections from spike trains: Dependence of the log likelihood ratio between Models M1 and M2 on the input amplitude. Each dot represents data for one cell and one input amplitude. Black dots denote trials where M2 significantly outperforms M1 (Chi-squared test, p<0.05). E) The true amplitude of aEPSCs vs the estimated amplitude (coupling coefficient) using Model 2 with all data (length varies across cells 216—720s). Colors denote different cells recorded with inputs of different amplitudes, and error bars denote bootstrap standard error. F) Dependence of the log likelihood ratios of models M1 and M2 relative to a homogeneous Poisson process on the length of data used for analysis. lnput amplitude of 10, data for N = 10 cells. Error bands denote standard deviation across cells. “Detection time” is defined as minimal data-length at which Model 2 with constant-gain coupling starts to outperform the Model 1 with no connection (marked by the red point). Also note that for <50s of data the models tend to over-fit (the accuracy on training data is higher than the accuracy on test data). G) Detection time as a function of input amplitude. The dependence is well-approximated by the black curve c/x2. Gray dots denote recording length for those cells where no change was detected. If it exists the detection time must be longer than the recording length, as indicated by the arrow.

2E). Here we use the kernel mean over the first 25ms to summarize the coupling strength in Model 2 (see Methods).

With only a short recording of pre and postsynaptic spikes it is difficult to determine if two neurons are “connected”. Both Model 1 with spike history-only and Model 2 with coupling tend to over-fit data from recordings shorter than ~ 5—10s (LLR<0 when compared to a homogeneous Poisson model which describes only the baseline firing; Fig. 2F). As the amount of available data increases over-fitting is reduced, and, with sufficient data, Model 2 is more accurate than Model 1 if there is truly a synaptic connection between the neurons. We define detection time as the length of data at which Model 2 with coupling provides a better fit than Model 1 with history only (Fig. 2F, red point). For the data in Fig. 2F , where the pre and postsynaptic neurons each have 5Hz firing rates and the amplitude of the synaptic connection is lo, this crossing point occurs around 20s. By varying the input amplitude we find that the recording time needed to detect an input of amplitude x falls off as approximately c / x2 (Fig. 2G), with c = 16.0:0.1s.

To study the detectability of a change in connection strength, we constructed artificial data sets from the recorded data, in which the strength of the synaptic connection between two neurons changes at time point t, from amplitude a1, during time period from 0 to t, to amplitude a2, for time t to 2t (Fig. 3A). As in previous analysis, we use the likelihood ratio to determine whether the synaptic weight has changed. We compare Model 2, used previously, with a single synaptic input of a constant amplitude throughout the recording time (0 to 2t) to a new Model 3 which allows the synaptic input to differ before and after time t (Fig. 3B).

3C, left). However, only Model 3 can accurately track changes of firing statistics produced by changes of synaptic strength and accurately describe the cross-correlation before and after the change-point t (Fig. 3C, right). As in the case of detecting a connection, the parameters of Model 3 accurately reproduce the amplitude changes in the presynaptic input (R2 = 0.96:0.01, data not shown). Detectability of synaptic weight changes in long recordings is comparable to the detectability of connections of constant strength. In 200s long recordings, a likelihood ratio test between

3D, compare to Fig. 2D). Again, the ability to detect a change strongly depends on the recording length. With less data available, detecting changes in synaptic strength becomes more difficult than simply detecting a connection. The recording length necessary to detect a change x 2 Au of connection strength after a known change-point falls off as approximately c / x2 With c = 21.5:0.1s (Fig. 3E).

Detection of connectivity in fully-defined input setting

1B). We used a presynaptic population consisting of the equal number of excitatory and inhibitory neurons, with log-normal distribution of synaptic amplitudes (same distribution, positive weights for excitatory, negative weights for inhibitory), and PSC kernels consisted of the same difference of two exponentials with time constants of 0.5ms and 5ms. The average input had an amplitude of 0.15 0 (corresponding to ~ 15pA, depending on o), and membrane potential responses to injection of this current again mimicked the statistics of membrane potential fluctuations in vivo with amplitudes of 15—20mV [8,38,39]. This paradigm for injection of fully-defined current allows us to examine the detectability of excitatory and inhibitory inputs of multiple amplitudes using the same recording [40]. Using this paradigm we examined 5 additional cells, 3 of which were driven to fire at ~ 5Hz, and the other 2 were each driven at several different rates: ~ 1Hz, ~ 5H2, and ~ 10Hz (produced by varying the DC offset). We then analyze the postsynaptic responses using the same pairwise model comparison techniques as we used for analyzing single

input experiments.

High-throughput experiments with injection of the fully-defined input allowed us to reveal further features of input detection. The detectability is influenced not only by amplitude, but also by the sign of the input (excitatory vs inhibitory) and the postsynaptic spike rates (Fig. 4A). Inhibitory inputs are detected less readily than excitatory. They have less impact on the postsynaptic firing, and thus are less accurate in predicting output spikes compared to excitatory inputs of the same magnitude (the log likelihood ratios comparing Model 2 with coupling to Model 1 with spike-history alone are 58i2% smaller for inhibitory inputs). As in single-input experiments, detection time decreases with amplitude approximately as c / x2. With these data we find c = 111, 32, and 18s for detection of excitatory inputs at 1Hz, 5Hz, and 10Hz output rates and c = 182, 46, and 29s for detection of inhibitory inputs at these rates. Across postsynaptic firing rates r, these times are well approximated by c/ rx2 (Fig. 4B). On average, inhibitory inputs required 32:3% more data for detection than excitatory inputs. Although detection time likely depends also on presynaptic firing rates as well as time course of PSCs (not just their amplitude), here, for making the comparison clear, pre-synaptic rates for all inputs were held at 5Hz and PSC kernels had the same time course, differing only in amplitude.

The coupling coefficients accurately reconstruct both excitatory and inhibitory input amplitudes over a broad range, and this reconstruction becomes more accurate with higher postsynaptic firing rates (Fig. 4C). The pairwise cross-correlations (0—25ms following presynaptic spikes) also reflect the input amplitudes fairly well but are much more nonlinear (Fig. 4D). Uncorrected cross-correlations tend to underestimate the relative amplitude of strong inhibitory inputs, but overestimate the relative amplitude of strong excitatory inputs.

Model 2, with post-spike history and coupling to a single input, can be extended to capture multiple presynaptic inputs simply by adding extra coupling parameters for each additional input. However, modelling all 1024 inputs in detail is computationally prohibitive. Here we consider two alternatives: 1) a model with 64 coupling terms where inputs of the same amplitude are grouped together, and 2) a bilinear model where the 1024 inputs are restricted to have the same shape. Since the presynaptic inputs in our experiments were specifically grouped together and all have the same PSC shape, these alternative models both provide accurate descriptions of the data (see Methods). It is important to note, though, that these models are simplifications. In typical multi-electrode spike recordings grouping information would not be available, and PSCs at different synaptic connections would have different shapes.

The accuracy of the reconstruction for both the grouped and bilinear models increases as more inputs are included in the model (Fig. 5). With 5% of the inputs observed, the coupling coefficients estimated by the GLM (average kernel from 0—25ms) are already well correlated with the underlying input amplitudes (R2 = 0.87 for the grouped model and R2 = 0.69 for the bilinear model). With all inputs observed this correlation increases to R2 = 0.94 and R2 = 0.82 for the grouped and bilinear models, respectively. Since the bilinear model has many more parameters, it is not unsurprising that there is more uncertainty in the parameter estimates given the same amount of data (200s in this case). However, in both the grouped and bilinear models, including multiple inputs provides more accurate reconstruction than pairwise model fit to the same data (Fig. 5B).

Previous methods for detecting synaptic inputs, by using pairwise spike statistics alone, do not take advantage of this additional predictability. For instance, in the traditional, nonparametric tests the cross-correlation between spiking of two neurons is compared to the cross-correlation obtained when spike timings of the presynaptic neuron are shuffled (Fig. 6A). Using bootstrapping we can create a distribution of the total spike count in a cross-correlation over a window following presynaptic spikes (Fig. 6B). Comparison of the distributions obtained using the observed vs shuffled spike trains allows us to test whether an input has a statistically significant effect on the firing of the postsynaptic neuron (see Methods). As with the model-based methods, we can vary the amount of data (recording length) used for computing this statistic, and determine a minimal recording time necessary for the detection that, as expected, depends strongly on the amplitude of an input (Fig. 6C). Using data from experiments with fully-defined input we find that, even after optimizing the detection window for the nonparametric test (Fig. 6D), the pairwise, model-based method yields slight but significant reductions in detection time (to 94:2% on average compared to the nonparametric test with 15ms time window, paired two-sided sign test p = 0.005). A grouped model that includes all inputs yields additional reductions in detection time (to 86i4%) compared to the nonparametric approach (paired two-sided sign test p<10_5). Using all inputs in the grouped model also reduces detection times relative to the pairwise, model-based tests, where spiking of only one presynaptic neuron is included in the model, while all other inputs are treated as noise (to 94i2%, paired two sided sign test p = 0.04).

Making this interval too long or too short will reduce the sensitivity of the test (Fig. 6D). However, the model-based methods are more robust to these differences, since the coupling kernel can flexibly account for changes in postsynaptic spike statistics produced by inputs with different time-courses. This advantage may be especially useful for reconstructing connectivity from in vivo data, where the shape of PSCs varies broadly across connections. Indeed, previous work has shown how model-based methods can produce a more accurate picture of connectivity compared to descriptive (cross-correlation) methods [22]. These results obtained in a fully-defined experimental setting illustrate how model-based methods can be applied to improve input detection over descriptive methods, as well.

Prediction of spikes

The pairwise analyses above indicate that the spike times of individual presynaptic inputs are useful in predicting postsynaptic spikes. Since inputs are independent, observing more presynaptic neurons should improve the overall model accuracy. Because spike times of all N = 1024 excitatory and inhibitory presynaptic neurons composing the total input to the postsynaptic neuron (Fig. 7A) are defined, we can analyze how increasing the fraction of observed inputs improves the prediction of postsynaptic spikes.

This term describes the strong influence that generation of an action potential has on the probability of future spikes. These dynamics are intrinsic to the spike generation mechanism and are present no matter how many inputs are included in the model.

Fig. 7B shows the time course of changes of gain that are due to the post-spike history term and due to the influence of N = 4 inputs included in the model for a typical neuron fit using the grouped model (middle plot). A linear combination of these few sparse, excitatory (spikes indicated by small bars above the trace) and inhibitory (spikes below the trace) inputs has only weak influence on the spiking, and can account for only a small portion of the variations in spike timing. However, as the number of inputs included in the model grows, their total contribution becomes comparable to that of the post-spike history term (Fig. 7B, bottom). Although the model used here does not explicitly describe the underlying fluctuations in membrane potential that result from the current injection, the contribution of the cumulative coupling terms of all inputs (N = 1024, red trace in Fig. 7B, bottom) tends to be correlated with the total injected current (Fig. 7B, bottom, gray trace, R = 0.54 for this example).

Fig. 7C shows a combined influence of the post-spike history term and an increasing number of inputs (N = 4; 256; 512; 1024) on the probability of postsynaptic spiking. With few inputs, the post-synaptic firing is dominated by post-spike effects and timing of postsynaptic spikes is predicted relatively inaccurately (AUC = 0.78 for this example with 4 inputs). With an increasing fraction of inputs included in the model, their contribution to spike prediction increases, and timing of individual spikes is predicted with progressively increasing accuracy (AUC = 0.82, 0.86, and 0.98 for 256, 512, and 1024 inputs, respectively). When all inputs are observed, the model spike probability tends to be sharply peaked around actual spike times (vertical lines in Fig. 7C), and the coupling effects contribute substantially to the rate variation (81% of the variance in the log rate is due to coupling).

When only few synaptic inputs are included in the grouped model the post-spike history accounts for nearly all of the variability in the firing rate (Fig. 8A). As more inputs are modeled the contribution of coupling becomes progressively stronger, and, in the grouped model, with inclusion of ~ 70—90% of the inputs it becomes as important as the intrinsic dynamics of spike generation described by the post-spike history (Fig. 8A). Note that the effect of individual inputs remains weak (Fig. 8A, bottom), but the total contribution from all inputs becomes substantial.

8B, black curve). Since the total current is simply the sum of all PSCs produced by all presynaptic inputs, we can also extract the “partial” current that is the current contributed only by a fraction of the observed neurons. This partial current is fairly consistently correlated with the coupling term (Fig. 8B, red curve). The reduced correlation of both the total and the partial current with the coupling term when the fraction of observed inputs approaches 100% appears to be due to the differences in the timescales of the coupling terms under the grouped model. The coupling term tends to follow the injected current with a slightly longer timescale as more inputs are grouped together (Fig. 8B, inset).

Log likelihood ratios (relative to a homogeneous Poisson model) increase monotonically with the increasing fraction of observed inputs (Fig. 8C). Note that, unlike in vivo data where the LLR saturates as more neurons are observed [26], here the model accuracy shows no signs of saturation. This is likely due to the fact that the synaptic inputs simulated here are independent, while firing of presynaptic neurons in vivo can be much more correlated. Using ROC analysis we find that cross-validated area-under-the-curve increases from 0.72:0.01 when 5% of the inputs are observed to 0.98:0.01 in the model that includes all (100%) inputs (Fig. 8D). Assessing spike prediction accuracy with single trial data is not always intuitive. Here, an area-under-the-curve of 0.98 signifies that a randomly chosen time bin containing a spike in the recorded response will have a higher p(spike) than a randomly chosen non-spike bin 98% of the time.

9). Here we simulate a neuron receiving the same presynaptic inputs on each trial and use the simulated post-spike history to generate model responses. Simulating multiple trials, we find that there is substantial trial-to-trial variability when only a small fraction (5%) of the inputs is used to generate the model response (Fig. 9, top). With this small number of inputs, the average response only weakly matches the peri-stimulus time histogram (PSTH) of observed spike responses (R = 0.09:0.02, PSTHs estimated by a maximum likelihood, fixed bandwidth kernel density estimator). When all inputs are modeled, trial to trial variability is substantially reduced (Fig. 9, middle), and the simulated PSTH more accurately matches the observed responses to repeated trials (R = 0.73:0.06). Similarly, the bilinear model reproduces the PSTH with R = 0.10:0.02 and R = 0.55:0.04 When 5% and 100% of the inputs are observed, respectively. Even when all inputs are modeled, the GLMs do not capture the full precision of the data (Fig. 9, right). These results suggest that, although single trial spike prediction with the GLM is quite accurate, more precise models, with additional nonlinearities [41,42] or explicit estimates of the underlying membrane potential [43], might be necessary to provide a full account of the transformations that occur as fluctuating current input leads to spiking output.

Current-based vs conductance-based synaptic input

I n vivo, nonlinear dendritic integration, probabilistic release, and plasticity Will all introduce certain variability in transmission at a synapse, making it more difficult to detect the connection and to estimate its strength. It is also important to note that here we use current-based synaptic input, While at real synapses currents are generated by changes of conductance. This causes the PSC amplitude to vary as a function of the membrane potential. To test whether the nature of current-based and conductance-based inputs could have a large effect on their detectability and estimation of their strength, we used post-synaptic spikes generated by two model neurons with either current-based or conductance-based synapses, and fit pairwise models to these new data.

Using the same presynaptic spike times and weights delivered to the observed neurons, we then optimize the parameters of these models to match both the observed membrane potential and spike timing (Fig. 10A, see Methods). After optimization, the subthreshold fluctuations of the models reproduce the observed fluctuations to 4.1mV RMSE — with similar accuracy between current-based and conductance-based models. Overall spike statistics (post-synaptic inter-spike interval distributions and cross-correlations) were well-matched between the models and data (81 Fig. ), and the spike coincidence factor [45] with A = 4ms averaged F = 0.35 for the current-based models and with the conductance-based models performing slightly worse at F = 0.31 (82 Fig. ).

10B-D). As with the experimental data, we find that, in both the current-based and conductance-based models, higher amplitude synapses are easier to detect and there is an asymmetry between excitatory and inhibitory inputs with inhibitory inputs being more difficult to detect than excitatory inputs of the same magnitude. As before, the postsynaptic firing rate has a large effect on detectability and estimation accuracy with higher rates resulting in faster detection of inputs. Estimation accuracy for the current-based model, in this case, was R2 = 0.85 for 1Hz output, 0.93 for 5Hz output, and 0.96 for 10Hz output (for clarity only 5Hz condition is shown in Fig. 10C). The conductance-based model yields similar results with R2 = 0.84, 0.93, and 0.96 for the three output rates. Finally, we find that detection times decrease as ~ c/x2 with increasing input amplitudes, and are shorter for higher post-synaptic firing rates. Together, these results show that the form of synaptic input — conductances vs currents — does not substantially alter our results once spike statistics are matched.

Discussion

Our results show that 1) Simulated synaptic inputs with amplitudes >0.25onoise can be reliably inferred from several minutes of spiking activity of pairs of neurons. Stronger connections can be detected even faster, but the amount of data required for detecting weaker connections increases rapidly. Thus, in typical experiments only a subset of connections can be detected, with a low amplitude limit depending on recording time and firing rate. 2) The inferred coupling strength is strongly correlated with the true amplitudes of synaptic inputs. 3) Excitatory inputs are more readily detected than inhibitory. 4) Detectability of changes of synaptic weights follows same rules and has same limitations as detection of individual synaptic connections. 5) Observing and modeling larger fractions of the presynaptic input improves model accuracy and reduces the necessary recording time for detecting inputs. Thus, model-based inference of functional connectivity allows accurate detection of individual inputs, reconstruction of synaptic amplitudes, and prediction of postsynaptic spiking.

These recordings provide a link between properties of synaptic connections typically measured in intracellular experiments in vitro, such as PSC amplitude, and estimates of connectivity made from in vivo extracellular recordings, such as inferred functional connectivity and coupling strength. We have used a real spike generation mechanism of an actual neuron to study spike responses to artificial PSCs with realistic time courses, distribution of weights and firing statistics of presynaptic neurons [46,47], consistent with experimental observations. Injected currents induced membrane potential fluctuations typical for in vivo activity [8,48,49]. However, we have also made some simplifying assumptions and neglected to model many known effects, such as nonlinear dendritic integration [50], short-term plasticity [51], possible input correlations [52] , probabilistic release, channel noise [53], and the fact that actual synapses are conductance-based [54,55]. Since these effects can have a substantial impact on spike statistics, future experimental work and statistical modeling may yield improvements over the current methods. These additional experimental parameters will dramatically increase the complexity of future experimental results and their interpretation. Here we have neglected these effects in order to simplify our analyses and determine the basic constraints of what can be inferred about simulated synaptic connectivity from spiking of neuronal ensembles. Although these experimental simplif1cations may limit the extent to which our results can be directly generalized to extracellular recordings, the trends described here provide concrete links between synaptic currents and the detectability of synaptic inputs. Applying model-based methods to spike trains generated by model neurons fit to our observed data, we found that our results are robust to both the form of input (conductance-based vs current-based synapses) and the spike generation mechanism (integrate-and-fire vs biophysical in real neurons).

The model-based methods used here are based on a tractable class of rate models and fail to describe some fundamental aspects of spike statistics. Here we used a common GLM approach where neurons are assumed to emit spikes according to a Poisson random variable with a rate that is determined by recent pre and postsynaptic spiking. While these models can describe the auto and cross-correlations present in the data, neurons often have more complex nonlinear behavior [56,57] and generate spikes much more reliably than Poisson processes [58]. In general, rate models may fail to capture some of the important nonlinear dynamics of real neural systems [59]. However, several statistical models have been developed that explicitly aim to describe the underlying membrane potential dynamics [60,61] and tend to yield more accurate spike prediction. Although inference from spike trains is more difficult in these models, they are able to model the high reliability of spikes in many experimental settings [43] and provide detailed insight into the membrane dynamics and changing conductances underlying the observed spiking [62]. Including coupling between neurons in these voltage-based models, as latent synaptic conductances, may lead to improved detection, estimation, and tracking of inputs from spikes.

Descriptive statistics of functional connectivity [18,19,63,64] and recent extensions [65] are widely applied tools for understanding interactions between neurons. However, as many previous studies have noted, correlations that are not due to direct monosynaptic connections can confound estimates of connectivity and compromise the interpretation of cross-correlograms [20,29]. Common-input, the dynamic effects of post-spike history on future spikes, poly-synaptic effects, and changing synaptic weights will all affect basic pairwise spike statistics. One key advantage of statistical inference model-based methods, over descriptive techniques, is that the many factors contributing to postsynaptic spiking can be modeled simultaneously and, to some extent, disentangled.

We demonstrated that detecting a synaptic weight change that might occur after a specific experimental intervention, such as tetanization, requires only moderately more data than detecting the presence of a connection. However, statistical methods that can infer synaptic changes from extracellularly recorded spike trains could also be applied to situations where change-points are unknown or not defined, for instance, when plasticity is a result of ongoing natural activity. New statistical methods are beginning to address these issues by explicitly modeling ongoing plasticity [66,67] , considering state-changes with unknown change points [68] , and providing more flexible descriptions of time-varying neural dynamics [69—72].

We find that detection time for an input of amplitude x is approximately proportional to 1 / x2 and also depends on firing rates. Importantly, inhibitory inputs are more difficult to detect than excitatory inputs of the same amplitude, and inhibition required ~ 30% more data for detection with pairwise tests. We further show that when an increasing fraction of presyn-aptic neurons can be included in the model, spike prediction accuracy increases and less data is necessary for detection of individual connections. Here we found that with ~ 15min of data and the pre and postsynaptic cells firing at ~ 5Hz, synapses could only be reliably detected when they had amplitudes > 0.250. This corresponds to PSCs amplitudes of 18—27 pA, which is in the range of average EPSCs amplitudes in L2/ 3 cells (15—25 pA, [73]). These results suggest that, to fully map functional connectivity estimated from spikes to actual physiological connectivity, longer recording lengths may be necessary to detect weak connections. At the same time, pairs of neurons with strong synapses or high firing rates may be evident even in short in vivo recordings.

Here we consider a limiting case where the input is fully controlled during the experiment. We find that a basic rate model (GLM) can provide surprisingly accurate spike predictions and approximately capture the true input current. These findings suggest that models of functional connectivity may begin to provide insight into actual synaptic effects as modern techniques allow recording of simultaneous spiking from increasing numbers of neurons in vivo.

Here we focused on the detectability and identification of individual inputs from simultaneously recorded spike trains. Our results highlight some of the important scientific possibilities offered by the use of statistical methods for understanding large-scale interactions between neurons. Modeling the interplay between post-spike history effects and presynaptic input will be essential to describe biophysical diversity within [77,78] and across [79,80] cell types. Extending both the statistical and experimental tools used here will provide a deeper understanding of the transformations of synaptic inputs into a neuron’s spiking output.

Methods

Ethics statement All animal use procedures were approved by the institutional animal care and use committee at the University of Connecticut, and conform to the principles outlined in the Guide for the Care and Use of Laboratory Animals (National Institutes of Health publication no. 86-23, revised 1985).

Electrophysiology and current injection

1) Partially-defined current was produced by a single simulated, excitatory presynaptic input immersed in fluctuating noise. Several amplitudes of artificial EPSCs were used in these experiments. 2) Fully-defined current was produced by the firing of a large population of simulated presynaptic excitatory and inhibitory neurons whose postsynaptic currents (PSCs) sum to mimic fluctuating, naturalistic input.

Experiment 1.

Partially-defined input: Artificial EPSCs immersed in fluctuating noise

Seven cells were injected with (0.5, 1.0, 1.5 0 inputs), one cell with (0.2, 0.3, 0.5, and 1 0 inputs), and two cells with (0.1, 0.25, 0.5, 1.0 0 inputs). This approach was similar to that used in previous work [2,3]. Note that this approach based on injection of current through the intracellular electrode ignores the transformation of synaptic conductances into postsynaptic currents, as well as effects of dendritic integration of synaptic inputs. However, the injected fluctuating current reproduces well the membrane potential fluctuations recorded in the soma of neocortical neurons in vivo [8,48,49]. Presynaptic spike timing for the artificial EPSCs was generated by a gamma renewal process (shape k = 2, scale 6 = 2.5). This corresponds to a neuron firing at a rate of 5Hz with more regularity (CV % 0.7) than a Poisson process. Spikes were then convolved with a synaptic integration kernel generated by a difference of exponentials — rise time of 1ms and decay time of 10ms — to generate an EPSC trace. This kernel was previously used in [3]. I n vivo, there is substantial variation in spike statistics [80,81] and dendritic integration [82,83]. Although a complete exploration of the parameter space is beyond the scope of these experiments, the values used here were chosen to be consistent with in vivo and in vitro observations. Current was injected in episodes of 46s, with intervals of 40—90s between the episodes. We recorded from 10 cells with 10—60 episodes recorded from each cell.

That is, if episodes are recorded with aEPSC amplitudes 0.1, 0.2, 0.5, 0.2 0, we assume that we can examine an amplitude change of 0.4 o by concatenating files 1 and 3 and that concatenating files 2 and 4 is equivalent to a single longer recording with 0.2 0 input. In all cases, the intervals between episodes, where there was no postsynaptic spiking, as well as the first and last 3s of each record were excluded from the analysis. In this first set of experiments we do not consider the effects of IPSCs.

Experiment 2.

Fully-defined input produced by a population of spiking neurons

Several factors determine the statistics of the injected current: 1) the statistics of presynaptic spiking, 2) the amplitudes of synaptic “weights”, and 3) the time-course of the PSCs. Here we make several simplifying assumptions. We assume that population of pre-synaptic neurons consists of an identical pool of gamma renewal processes, 50% of which are excitatory and 50% of which are inhibitory, and we assume a log-normal distribution of synaptic weights based on observations from paired in vitro cortical recordings [46]. Finally, we generated the PSCs as a difference of exponentials with a rise time of 0.5ms and decay time of 5ms.

Here we use N = 1024 with amplitudes drawn pseudo-randomly from a discrete approximation to the log-normal distribution described in [46] (1,4 = 0.702, 0 = 0.9355). Absolute PSC amplitudes are constrained to 32 discrete, log-spaced values to simplify the analysis, and the proportion of amplitudes in each bin is exact to maintain excitatory-inhib-itory balance. Excitatory and inhibitory inputs are assumed to have the same distribution, differing only by the sign. We use k = 2, 6 = 2.5 for the gamma renewal processes and resample short (<10ms) ISIs to avoid strongly overlapping PSCs.

Moreover, this basic population spiking model does not capture all of the details of cortical circuitry. There are known differences in the spike statistics and time course of PSCs for different cortical cell types [80]. However, this model maintains balanced, approximately OU statistics and allows us to examine the detectability of a broad range of synaptic weights. The simulation parameters here were chosen with the goal of making inputs physiologically realistic, but structured enough to analyze statistically. Using these parameter settings we recorded an additional 5 cells where the DC component was tuned to maintain a desired firing rate. Three cells were recorded around ~ 5 Hz, while two cells were driven at multiple firing rates: ~ 1 Hz, ~ 5 Hz and ~ 10 Hz.

Inferring functional connectivity from spikes

For detection of synaptic inputs from spiking, we use a basic statistical model that describes the probability of postsynaptic spikes as a function of the presynaptic spike timing — a generalized linear model (GLM) with Poisson observations. For a single presynaptic input we assume

This rate is determined by a linear combination of a baseline firing rate 190, the influence of spiking history, that is the previous spikes of the post-synaptic neuron, parameterized by bpost, and a term that describes the effect of the input from presynaptic spikes, parameterized by bpre. A set of basis functions fk provides a smooth expansion of each of these effects. The parameters of this model can then be readily estimated by maximizing the log-likelihood [84] In the results presented here we often summarize the effect of coupling with a single “cou-K pling weight”, defined as the average kernel over the first 25ms: > k=1 This basic GLM can then be extended to include multiple inputs by using In this case, with O(kN) parameters, rather than straightforward maximum likelihood estimation, regularization is important to prevent over-fitting. Here we use L1-regularization and find the maximum a posteriori (MAP) parameters estimates The regularization hyper-parameter is optimized by maximizing cross-validated log-likeli-hood, and no regularization is performed on the baseline or post-spike history terms. Here we use 10-fold cross-validation both to optimize the hyper-parameter and evaluate model accuracy (see below).

This assumption is useful for analysis of in vivo data, since only a small number of potential presynaptic inputs are likely to actually be connected to the postsynaptic neuron, so such sparse solutions are desirable. Other regularizers, such as spike-slab, group-L1, or even the log-normal prior used to generate the input weights here, may produce more accurate estimates. However, the fact that the log-posterior is concave with L1-regu-larization provides the additional advantage of a unique, global optimum. For both MLE and MAP (L1 regularized) optimization we use scaled sub-gradient maximization using the L1 General toolbox [85] . Even with regularization, fitting a model with 1024 inputs is computationally prohibitive. Here we use two alternative simplifications: 1) a grouped GLM where inputs of the same weight are grouped into a single covariate, and 2) a bilinear extension of the GLM Where the inputs are constrained to take the same shape.

Namely, groups of presynaptic neurons have identical amplitudes distributed according to a discrete approximation of the log-normal distribution [46]. We can thus collect inputs into groups to simplify the analysis from Eq. (3) — instead of N = 1024 simulated presynaptic input neurons, we can study an equivalent system with only N = G inputs, Where G is the number of amplitude “groups” (G = 64 used in our simulations) by combining spike trains Fitting a GLM for this reduced set of inputs as is much more computationally tractable than modeling all individual inputs when N > G. Here we make use of this group simplification to examine how detectability and overall model accuracy are affected by the fraction of observed inputs. In this case we randomly sample a fraction p of the total inputs N and include these “observed” presynaptic spikes in a grouped-GLM as described above.

The results from the grouped model likely underestimate parameter uncertainty and overestimate spike prediction accuracy compared to the full estimation problem. However, accurately estimating connectivity using models with O(N) parameters is not a trivial problem for N> 1000; structured regularization techniques [86] and detailed model comparison [34,87] will be essential as the number of connections in these types of models increases. Here as a point of comparison for the grouped model we also use a bilinear model that models individual inputs with separate parameters, but restricts all inputs to have the same shape. Here, rather than K parameters per input, there is only a single “weight” parameter w,per input. In this case we fit the model by coordinate ascent — alternating between maximizing the posterior with respect to the input weights while holding the shape constant and maximizing the likelihood (no regularization) with respect to the shape while holding the weights constant. This approach tends to work well in practice [66,88] , and effectively approximates the full model, which would have O(Nk) parameters, with a model with O(N+k) parameters.

. .fK} to be a set of gamma distributions 1‘" exp(—t/ m) / Z where n = {0. . .K — 1}5 and Z denotes normalization by the maximum. For the post-spike history term and coupling terms during pairwise analyses we use m = 25ms, s = 1.2, and K = 8. In the full model, this same K = 8 basis is used for the post-spike history, but for simplicity coupling terms use a reduced basis with m = 45ms, s = 1.5, and K = 4. These basis sets allow us to fit a variety of band-limited, smooth functions on a range ~ 0—500ms with faster variation near Oms. Similar bases have been used in previous work [30,86]. Together with the baseline parameter, these settings give 9 total parameters for Model 1, 17 for Model 2, 25 for Model 3, 265 for the full grouped model with all inputs, and 1037 parameters for the full bilinear model with all inputs.

Quantifying accuracy and detecting functional connections

Given this model framework, we can now examine whether the effect of presynaptic spiking on postsynaptic spiking is detectable. In general, if we have two models H1 and Hz with Poisson observations the log likelihood ratio is given by where the two models have conditional intensities defined by 11 and 12 (log base 2 is used LLR (H 1,H2) * log2 when reporting bits). Importantly, the log likelihood ratio quantifies the relative accuracy of the two models. For instance, when H2 is a homogeneous Poisson model that only describes the mean firing rate, the log likelihood ratio quantifies how much more accurately spikes are predicted by the model H1 over just predicting the mean.

A second approach is to use an explicit likelihood-ratio test (without cross-validation). This test makes use of the fact that, for nested models, log likelihood ratios are approximately [2 distributed with dfz — dfl degrees of freedom when df]L and dfz denote the number of free parameters in H1 and H2, respectively. We then consider the effect of the presynaptic input to be statistically significant when the ZZ-test gives p < 0.05. Here we base “detection time” on the minimum amount of recorded data (averaged over blocks length T) needed to satisfy the cross-validated LLR criterion.

Especially for models with many parameters (i.e. the bilinear model with many pre-synaptic inputs), cross-validation is essential to avoid over-fitting. Even though estimating confidence intervals for test-set log-likelihood is problematic [89], defining “detectability” qualitatively, based on generalization performance, is more appealing in larger models. Here we use a convention of reporting and plotting the cross-validated LLR even for simple, pairwise models. In evaluating statistical significance, we limit our analysis to the pairwise models and use the un-cross-validated likelihood-ratio test.

Although the models used here assume Poisson observations, potentially allowing >1 spike per bin, we can also perform ROC (receiver operating characteristic) analysis by limiting our predictions to no-spike p(n(t) = 0|l(t)) vs spike (1 — p(n(t) = 0|l(t))) classification. A bin size of 1ms is used for all analysis, such that, in practice, more than one postsynaptic spike per bin never occurred.

Detection of weight changes

The model, tests for detectability, and accuracy analyses described above are generic tools for model-based estimation of spike train statistics and have been applied to a variety of in vivo eXtracellular data [22] , in which the underlying synaptic weights are unknown. The data collected here, with the controlled synaptic input, offer several unique opportunities to study how functional connectivity estimation is related to synaptic input.

We assume that responses to injected current containing PSCs of different amplitudes can be concatenated in order to study putative changes in synaptic strength. As with the detection of connections, here we define a weight change as detectable if modeling the effects 191”ng and bpre¢>g around a known change-point t’ results in a cross-validated log likelihood ratio >0 when compared to a model with a single coupling effect bpre for all observations.

Nonparametric testing on functional connections

However, one traditional method for determining whether a functional connection is present is to examine the cross-correlations between two spike trains. Using bootstrapping we can generate a distribution that reflects the expected rate change in the postsynaptic firing following a presynaptic spike, as well as the uncertainty in those rate changes. We then construct a null distribution by measuring the cross-correlations in data where the presynaptic ISIs have been shuffled. By taking differences between samples of the observed and shuffled distributions we can then estimate a p-value: the probability of observing, by chance, a rate change at least as extreme as the one actually observed. To compare this nonparametric test with the model-based framework, we define a test statistic based on the average cross-correlation in a range of lags 0 g T g 1’ ms following the presynaptic spikes After computing this statistic for bootstrap samples of the observed data robs and data with shuffled ISIs rshufl; we can test against the null-hypothesis that there was no difference between observed and shuffled statistics. The null hypothesis is rejected at a confidence level a, when P(T'0bs — rshufir> 0) > 1 — 06/ 2 01‘ 2 P(T'obs — rshufic> 0) < 06/ 2. Here we use I’ = 15ms with 1ms bins and a = 0.05. To compute detection times we use 4096 total samples for each distribution with recording length T, aggregated over 16 randomly selected time blocks.

Comparing current and conductance-based inputs with simulated neurons

Our in vitro data were collected by injecting current into L2/ 3 pyramidal cells, where the current is generated by a sum of fixed post-synaptic currents from simulated pre-synaptic neurons. Since real synaptic currents are generated by changing conductances, it is possible that our detection results overestimated the detectability and estimation accuracy of synaptic weights from spikes. To examine this potential confound we simulate an adaptive exponential integrate-and-fire model neuron in two situations: receiving either current-based or conductance-based synaptic inputs. The adaptive integrate-and-fire model takes the general form Where the dynamics of the membrane potential V depend on the capacitance C, leak conductance gL, resting potential EL, an adaptation variable w, DC current input IO, and fluctuating synaptic currents Isyn(t). The behavior of the membrane near the spike threshold VT is deter-V— VT AT

mined by the exponential nonlinearitygLATexp< ), and the adaptation variable has its own

As with other integrate-and-fire models, spikes occur when the membrane potential crosses a threshold VT, after which the potential is immediately reset to Vreset. After spiking the adaptation term is also updated instantaneously with w <— w + 19. Here the exponential nonlinearity and adaptation term allow us to better describe the post-spike and near-threshold behavior of the observed neurons.

For current-based inputs we know Isyn(t) exactly. Given the simulated presynaptic spike times ni(t), PSC kernel p, and weights for each synapse v1- these currents were generated as To modify the synaptic inputs for changing conductances we split the inputs into excitatory and inhibitory contributions with two separate reversal potentials ry and inhibitory contributions.

For the current-based models this formulation leaves the parameters {C,gL,VT, Vreset,AT,TW,a,b} to be optimized. For the conductance-based models, since we want the average PSC for each presynaptic input to match the original inputs, we introduce an additional constraint that the excitatory and inhibitory input be balanced: 611 < V(t) — VI > = — < V(t) — VE > aE. Since Isyn(t) has 0 mean, we approximate the expectations using the equilibrium voltage far from threshold: Veq = VL + Io / gL. Given this constraint we have the parameters {C,gL,VT,Vre$et, AT,TW,a,b,aE} for the conductance-based model. For both the current-based and conductance-based inputs we then optimize model parameters to match the observed voltage fluctuations and spike timing using derivative-free search (Nelder-Mead) with random restarts. Specifically we aim to minimize a cost function that mixes voltage and spike timing terms Here the first term allows us to minimize the discrepancy between the observed and simulated membrane potential (neglecting voltage changes within 5ms of observed post-synaptic spikes N). The second term, the van Rossum distance [90], accounts for discrepancy between observed and simulated spike times, and here we use TR 2 2 / f Where f denotes the post-synap-tic firing rate of the observed neuron. Together these terms allow us to closely approximate the current-injection experiments and simulate model neurons receiving similar current-based or

Supporting Information

Inter-spike interval distributions and cross-correlograms for synaptic inputs of different strength in a neuron recorded in a slice (“Observed”) along with results from two adaptive-exponential integrate-and-fire simulated neurons, receiving inputs either Via current-based or conductance-based synapses. In all three cases, the same presynaptic spike times and weights were used. Model neurons were optimized to fit the observed voltage fluctuations and spike timing.

Parameters and error metrics for the current-based and conductance-based adaptive-exponential integrate-and-flre models fit to two different observed neurons (Layer 2/3 pyramidal cells). Box plots denote median, interquartile range and 1.5x interquartile range for 60 ($2905) or 61 ($2906) experimental trials with post-synaptic firing rates of le, 5Hz, and 10Hz (combined). Outliers not shown for clarity. R denotes membrane resistance (1 / gL), 1‘ denotes the membrane time constant, DT determines the strength of the exponential nonlinearity near threshold, While a, b, and TW determine the dynamics of the adaptation variable. For the conductance-based models the scaling factor A (aE in the text) converts the presynaptic conductances to currents.

Author Contributions

Performed the experiments: MV VI. Analyzed the data: IHS. Wrote the paper: MV IHS.