In this work, we investigate the network reconstruction of the conductance-based integrate-and-fire (I&F) neuronal network. Experimentally, it has been shown that the I&F models can capture linear subthreshold properties as well as firing statistics of a real neuron (Carandini et al., 1996; Rauch et al., 2003; Burkitt, 2006). Theoretically, the conductance-based I&F neuronal models have been widely applied in large-scale neuronal network modeling and simulations to investigate information processing in various brain areas (Somers et al., 1995; Troyer et al., 1998; McLaughlin et al., 2000; Tao et al., 2004; Cai et al., 2005; Rangan et al., 2005; Zhou et al., 2013a). The I&F neuronal dynamics is governed by the following equation,

where Vⁱ is the membrane potential for the ith neuron in the network, GEi and GIi are its excitatory and inhibitory conductances, respectively. G_L is the leakage conductance, ϵ_L is the resting potential, ϵ_E and ϵ_I are the excitatory and inhibitory reversal potentials, respectively. When the membrane potential of a neuron is below the threshold V_th, the neuronal dynamics is described by Equation (1). When the membrane potential of a neuron reaches the threshold V_th, it is reset to the resting potential ϵ_L for a refractory period τ_ref, and evolves again following Equation (1). An action potential of a real neuron is signified by the event of the threshold crossing and voltage resetting in the I&F model. This threshold-reset event is referred to as a firing event (spike) of a neuron. The time evolution of conductances can be explicitly expressed as GEi=∑j≠i∑ksij+α(σGE, σHE, t-Tj,k)+f∑lα(σGE, σHE,t-Ti,lP) and GIi=∑j≠i∑ksij-α(σGI, σHI, t-Tj,k), where T_j,k is the firing time of the kth spike of the jth neuron, Ti,lP is the arrival time of the lth spike of the external Poisson input to the ith neuron with strength f and rate μ. The alpha function α(σd, σr, t)=σrσdσr-σd[exp(-t/σr)-exp(-t/σd)]Θ(t), where Θ(t) is the Heaviside function, has rising time constant σ_r and decay time constant σ_d. The nonnegative values sij+ and sij- are the excitatory and inhibitory coupling strengths from neuron j to neuron i, respectively. For the sake of convenience, we unify the notation of excitatory coupling strength sij+ and the inhibitory coupling strength sij- into s_ij, which is defined as sij=sij+ if neuron j is excitatory and sij=-sij- if neuron j is inhibitory. Therefore, a positive s_ij indicates an excitatory coupling whereas a negative s_ij an inhibitory coupling.

In our work, we use dimensionless unit for membrane potentials, in particular, V_th = 1, ϵ_L = 0, ϵE=143, ϵI=-23. The corresponding unscaled physiological values are V_th = −55 mV, ϵ_L = −70 mV, ϵ_E = 0 mV, ϵ_I = −80 mV (McLaughlin et al., 2000; Cai et al., 2005; Rangan et al., 2005; Zhou et al., 2013a). Time constants retain their dimension with the unit ms. We set τ_ref = 2 ms and the conductance time constants σ_GE = 2 ms, σ_HE = 0.5 ms, σ_GI = 5 ms, σ_HI = 0.8 ms. Conductance has the unit ms⁻¹. In our simulation, GL=0.05ms-1, which corresponds to the physiological membrane time constant 20 ms. The numerical method we use to evolve the system (Equation 1) is a fourth order Runge-Kutta method with spike-spike corrections (Rangan and Cai, 2007; Zhou et al., 2009, 2010).

It can be seen intuitively from Equation (2) that αijl encodes the information of how the history of spike train of neuron j contributes to the prediction of the subthreshold voltage of neuron i. However, it remains an important issue of how αijl is related to the underlying neuronal dynamics, in particular, the synaptic coupling structure s_ij. In the following, we investigate the conductance-based I&F dynamics (Equation 1) to address this issue.

We have established above the linear relation between M_ij and the coupling strength s_ij. Next, we investigate the issue of whether this linear relation is invariant across different dynamical regimes of the neuronal network. Specifically, we study the behavior of the proportionality constants B_E and B_I in Equation (4) for different network input rate μ and strength f. Figure 3A (or Figure 3D) displays the result of our scanning of M_ij as a function of f and μf for a fixed synaptic coupling s₂₁ = 0.01 (or s₂₁ = −0.01), s₁₂ = 0 for a unidirectional two-neuron network. Note that μf is the mean input current of the Poisson input. In our parameter scan, the Poisson input strength f ranges from 0.0016 to 0.1, the Poisson input current μf ranges from 0.012 to 0.062, thus the range of the Poisson input rate μ is from ~0.1 ms⁻¹ (100 Hz) to ~40 ms⁻¹ (40 kHz) correspondingly. The network dynamics driven by the above Poisson inputs covers a wide variety of dynamical regimes from the highly fluctuating regime to the mean driven regime as firing rate changes from ~1 to 150 Hz.

As observed from Figure 3, for the uncoupled direction, M₁₂ always remains nearly 0 across different dynamical regimes; whilst for the excitatory (inhibitory) coupling, M₂₁ is nearly a positive (negative) constant across different dynamical regimes. Therefore, we can conclude that, over broad dynamical regimes, the coupled direction can be well distinguished from the uncoupled direction. It is important to emphasize that B_E (or B_I) is approximately a constant over a wide range of dynamical regimes for the excitatory (or inhibitory) coupling between neurons. As will be seen below, B_E (or B_I) remains the same for the inference of synaptic coupling between different pairs of neurons in a neuronal network. If the proportionality constant B_E (or B_I) is known, we can recover the excitatory (or inhibitory) coupling strength s_ij by s_ij = M_ij/B_E (or s_ij = −M_ij/B_I). In the following, we show that the coupling strengths between neurons obtained from different dynamical regimes can be used to verify the consistency of prediction of the network structure.

For the I&F dynamics (Equation 1), it can be clearly seen that from Figures 3C,F that the time lag l_ij ≡ 2 (i.e., Mij≡αij(2)), which corresponds to the peak of the response kernel αijl, for both excitatorily and inhibitorily coupled directions across different dynamical regimes. This is consistent with what is presented in Figures 2B,D, where M_ij attains the peak (absolute) value of αijl at l = 2 for both excitatory and inhibitory couplings. Therefore, in the following discussion, we will set l_ij ≡ 2 and H0 will be rejected when |M~ij|>F-1(1-r/2).

We have presented above the case of two-neuron I&F networks to illustrate our method. We now turn to the question of how well our STR method can perform for a multi-neuron network with complex connectivity motifs and different neuron types. Here, we apply our STR method to an I&F network of five neurons consisting of three excitatory neurons and two inhibitory neurons with various coupling strengths (see Figure 4A). Using 20 s of the voltage and spike train time series, we can successfully reconstruct the topological connectivity of the network (see Figure 4B). By the property of inference invariance, using B_E = 0.32 and B_I = −0.15 obtained from the two-neuron network in Figure 1, we can even recover excitatory and inhibitory coupling types and strengths as shown in Figure 4B. The corresponding 99% confidence intervals for the values of couplings are also displayed in Figure 4B. Note that the true values of the coupling strength indicated in Figure 4A all fall within the 99% confidence intervals. In this network there is a type of motif of neuron X presynaptically coupled to neuron Y, Y in turn to neuron Z, i.e., X→Y→Z (for example, in Figure 4A, 1 → 2 → 3, 3 → 4 → 5, 4 → 5 → 1 and 4 → 5 → 2 are of this motif). In this motif, X does not influence Z directly, and it influences Z only indirectly through Y. As demonstrated in Figure 4, our STR method is able to differentiate the direct influence from the indirect influence and reconstruct the true synaptic connectivity. We emphasize that it takes only 20 s of data for our STR method to provide a rather precise reconstruction of synaptic connectivity with an accurately recovered connectivity strength s_ij.

We next illustrate our STR method can be successfully used to recover a neuronal network of large size. We apply our STR method to 100-neuron network (80 excitatory neurons and 20 inhibitory neurons) with random topological connectivity of different coupling strengths. The absolute coupling strength |s_ij| is generated from the uniform distribution U(0, 0.01) to describe the diversity of coupling strength. Figure 5 displays the results for a sparse network and a dense network with connection probability of 15 and 70%, respectively. The sparse network exhibits asynchronous firing activity as in Figure 5A whereas the dense network exhibits nearly synchronous firing activity as shown in Figure 5D. The effectiveness of our STR reconstruction is demonstrated in Figures 5B,C,E,F, respectively, for both networks.

In Figures 5B,E, each connected pair of neurons is represented by a dot which describes the relation between M_ij and s_ij. The dots are tightly concentrated around the straight lines of the linear relations (Equation 4) between M_ij and s_ij. It should be stressed that, as with the case of the five-neuron network above, the linear relation between M_ij and s_ij with B_E and B_I obtained from the two-neuron network persists and is robustly preserved over different pairs of neurons with different coupling strengths for a large neuronal network. Therefore, the value of M_ij can be used to predict both the coupling type and the coupling strength from neuron j to neuron i. The width of the spread of dots around the straight lines is determined by the mean value of θ_ij (standard deviation of M_ij) averaged over all directions, which can be used to quantify the uncertainty of reconstruction. In both Figures 5B,E, the mean values of θ_ij are ~ 9 × 10⁻⁵. As will be discussed below, the spread width can be reduced by increasing the data length in the STR analysis.

We further note that there are some coupled directions incorrectly predicted as uncoupled (green dots in Figures 5B,E). This arises from the fact that, when the coupling is sufficiently weak, the strength of M_ij can be comparable to or even smaller than its standard deviation, making H0 unlikely to be rejected when a finite length of data is used. On the other hand, for a fixed length of time series, when a coupling is sufficiently strong, the value of M_ij can be much larger than its standard deviation, yielding a correct rejection of H0. To quantify to which extent a weak coupling can still be successfully predicted, we define excitatory and inhibitory critical values SEc and SIc as follows. First, we reconstruct the topological connectivity of the network with significance level r (we usually set r = 0.01) using time series of certain fixed duration. Then we compare the predicted and the true connectivity and locate the values of SEc and SIc so that, for all couplings satisfying sij>SEc or sij<SIc, 99% of them can be correctly predicted. Clearly, the closer the critical value is to 0, the weaker a coupling that we can correctly predict through our STR, hence the more accurate the network reconstruction. In Figures 5B,E, which use 100 s of voltage and spike train times series, the critical values for our STR reconstruction are SIc~-0.002 for inhibitory couplings and SEc~0.0008 for excitatory couplings.

Figures 5C,F displays the histogram of the predicted M_ij for the uncoupled pairs. As expected, the predicted M_ij for the uncoupled directions distributes around 0 and ~99% of the uncoupled directions are correctly predicted, which conforms with the significance level r = 0.01. The histogram of M_ij shown in Figures 5C,F is not necessarily Gaussian because each M_ij is from a Gaussian distribution with a different variance.

In Figure 5D, the neuronal network exhibits a nearly synchronous global oscillation of ~35 Hz and the firing pattern of each individual neuron is relatively irregular due to the stochastic Poisson input. In such case, our STR method can still reconstruct the underlying synaptic coupling. Note that the network dynamics shown in Figure 5D is not dominated by inhibition and the firing rate of each neuron is close to the global oscillation frequency. These dynamical features are different from those of the synchronous irregular regime as described in Brunel (2000) where the dynamics are dominated by inhibition and the firing rate of each neuron is much lower than the global oscillation frequency. Furthermore, the dense connectivity also gives rise to a strong effect of indirect influence (i.e., influence mediated by other neurons) between neurons, making it difficult to uncover their true couplings. However, as demonstrated in Figure 5, our STR network reconstruction method performs surprisingly well even for nearly synchronous, densely connected I&F neuronal networks.

Through further examination of many networks with connection probability from 5 to 70%, we can conclude that our STR is able to reconstruct accurately the connectivity of both sparse and dense networks in either asynchronous or nearly synchronous dynamical regime.

As is demonstrated above that our STR method can reconstruct accurately the I&F neuronal synaptic connectivity using time series of 20 ~ 100 s, we further address issues of how one can improve the prediction accuracy. Intuitively, a longer data length should reduce the statistical error and thus improve the network reconstruction accuracy. To determine how long one needs for time series in order to achieve a desired accuracy, we investigate the relation between the STR reconstruction accuracy and the data length T. In Figures 6A,B, we present numerical results of two types of quantifications of the network reconstruction accuracy.

First, we consider θ_ij, which represents the statistical uncertainty of M_ij. Given a coupled direction with coupling strength s_ij (without loss of generality, we consider s_ij > 0), the mean value of M_ij is approximately B_Es_ij. It can be easily seen that, the smaller θ_ij, the more likely the hypothesis M_ij = 0 can be correctly rejected. The confidence interval of the predicted coupling strength is determined by θ_ij as Mij/BE±F-1(1/2+c/2)θij/BE, where c is the confidence level, F⁻¹(·) is the inverse function of F(x)=12π∫-∞xexp(-t22)dt. Clearly, given c and B_E, the smaller θ_ij is, the more accurate the predicted coupling strength M_ij/B_E. To represent the overall accuracy of our synaptic connectivity reconstruction, in Figure 6A, we plot the mean value of θ_ij averaged over all directions, which controls the width of the spread of dots around the straight lines as in Figures 5B,E, as a function of data length T. We can see that it follows a T-12 scaling as expected from the central limit theorem.

As a second quantification, we use the excitatory and inhibitory critical values SEc and SIc of coupling as the indicator of the accuracy of the network reconstruction. In Figure 6B, both SEc and SIc possess approximately a power law decay as the data length increases. Their decay rate is determined empirically as T^−0.65, which is slightly faster than T-12. As discussed above, the smaller the critical (absolute) values SEc and SIc, the more accurate the weak coupling strengths can be predicted, thus the better network reconstruction.

Based on the above observations, we can conclude that the accuracy of the STR reconstruction follows a scaling of at least T-12 as the data length T increases. Therefore, if one desires to achieve a factor of Q (Q > 1) improvement in the accuracy of STR reconstruction (i.e., a 1/Q times smaller confidence interval, or to detect a 1/Q times weaker coupling strength), data of Q² times longer length is required.

The accuracy of reconstruction can be further improved through reducing the sampling interval length τ, i.e., increasing the sampling rate. Intuitively, as one samples at a higher rate, i.e., a smaller τ, more information of detailed dynamics is preserved, hence more accurate the network reconstruction. This is confirmed in Figure 6C. Note that a change of the τ value leads to a change of the values of B_E and B_I in the linear relation (Equation 4) between M_ij and s_ij. The underlying reason is as follows. As τ vanishes, with the self-prediction residual ΔVti decreasing for the smooth dynamics, the reduction of the magnitude of its response kernel αijl follows the same τ scaling. Therefore, B_E and B_I should be recalibrated when a different τ is used in our STR reconstruction. Furthermore, the firing rate of presynaptic neurons could also influence the accuracy of synaptic connectivity reconstruction. The covariance of αijl is approximately inversely proportional to the covariance of the presynaptic neurons spike trains. Therefore, if the firing rate of presynaptic neurons is too low, one may need to use longer time series to obtain a high signal-to-noise ratio for an accurate recovery of synaptic connectivity.

Under the situation where voltage and spike train information of all neurons in the I&F network can be obtained, we have demonstrated the efficiency of our STR method for connectivity reconstruction. However, to apply the STR method to physiological experiments, one may encounter certain constraints from the recording techniques. In experiments, the voltage trace of very few neurons can be acquired through intracellular recording since it is rather difficult to perform intracellular recording on a large set of neurons simultaneously. Nevertheless, we can consider a setup in which the voltage trace of only a target neuron is recorded and the spike trains of many other neurons in the network are obtained through other means, say, calcium imaging (Stosiek et al., 2003; Grewe and Helmchen, 2009; Grewe et al., 2010) or Multielectrode array (MEA) (Litke et al., 2004; Field et al., 2010; Shimono and Beggs, 2015). Applying the STR method to this type of data, we are able to reconstruct the synaptic connectivity from the recorded neuronal population to the target neuron. An example is illustrated in Figure 7. For the same network of 100 neurons as in Figure 5A, we choose Neuron 13 as the target neuron and use our STR method to recover the network couplings to Neuron 13. Note that, we can think of this example as sampling 100 active neurons from a large cortical circuit and the effect of all the unrecorded neurons on each recorded neuron is modeled by an external Poisson input with strength f = 0.012 and rate μ = 1 ms⁻¹. In general, if spike trains of a much larger neuronal population can be recorded, we could recover these corresponding synaptic connections on the dendritic tree of the target neuron. Here, 100 s of data is used with significance level r = 0.01 for reconstruction. Incidentally, we can choose any neuron in the network as a target neuron and obtain a reconstruction similar to what is shown in Figure 7. Comparing Figures 7A,B, we can observe that, except for one very weakly coupled direction (46 → 13), the other couplings to the target neuron are successfully recovered. Displayed in Figures 7C is a detailed comparison between the predicted and the underlying true coupling strengths. Evidently, the underlying true coupling strength falls within the 99% confidence interval of the predicted coupling strength.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer GM and handling Editor declared their shared affiliation.