Abstract

With our ability to record more neurons simultaneously, making sense of these data is a challenge. Functional connectivity is one popular way to study the relationship of multiple neural signals. Correlation-based methods are a set of currently well-used techniques for functional connectivity estimation. However, due to explaining away and unobserved common inputs (Stevenson, Rebesco, Miller, & Körding, 2008), they produce spurious connections. The general linear model (GLM), which models spike trains as Poisson processes (Okatan, Wilson, & Brown, 2005; Truccolo, Eden, Fellows, Donoghue, & Brown, 2005; Pillow et al., 2008), avoids these confounds. We develop here a new class of methods by using differential signals based on simulated intracellular voltage recordings. It is equivalent to a regularized AR(2) model. We also expand the method to simulated local field potential recordings and calcium imaging. In all of our simulated data, the differential covariance-based methods achieved performance better than or similar to the GLM method and required fewer data samples. This new class of methods provides alternative ways to analyze neural signals.

1  Introduction

Simultaneous recording of large populations of neurons is an inexorable trend in current neuroscience research (Kandel, Markram, Matthews, Yuste, & Koch, 2013). Over the last five decades, the number of simultaneously recorded neurons has doubled approximately every seven years (Stevenson & Kording, 2011). One way to make sense of these big data is to measure the functional connectivity between neurons (Friston, 2011) and link the function of the neural circuit to behavior. As previously reviewed (Stevenson, Rebesco, Miller, & Körding, 2008), correlation-based methods have been used to estimate functional connectivity for a long time. However, they are suffering from the problem of explaining away unobserved common inputs (Stevenson et al., 2008), which makes it difficult to interpret the link between the estimated correlation and the physiological network. More recently, Okatan, Wilson, and Brown (2005), Truccolo et al. (2005), and Pillow et al. (2008) applied the generalized linear model to spike train data and showed good probabilistic modeling of the data.

To overcome these issues, we developed a new class of methods that use not only the spike trains but the voltages of the neurons. They achieve better performance than the GLM method but are free from the Poisson process model and require fewer data samples. They provide directionality information about sources and sinks in a network, which is important to determine the hierarchical structure of a neural circuit.

In this article, we further show that our differential covariance method is equivalent to a regularized second-order multivariate autoregressive model. The multivariate autoregressive (MAR) model has been used to analyze neuroscience data (Friston, 2011; McIntosh & Gonzalez-Lima, 1991). In particular, this model with an arbitrary order has been discussed previously (Harrison, Penny, & Friston, 2003). Continuous AR process, which is known as the Ornstein-Uhlenbeck (OU) process (Uhlenbeck & Ornstein, 1930), has also been applied to model neuronal activity (Burkitt, 2006; Ricciardi & Sacerdote, 1979; Lánskỳ & Rospars, 1995). However, the modified AR(2) model in this has not been used previously.

All of the data that we analyze in this article are simulated data so that we can compare different methods with ground truth. We first generated data using a simple passive neuron model and provided theoretical proof for why the new methods perform better. Then we used a more realistic Hodgkin-Huxley (HH)–based thalamocortical model to simulate intracellular recordings and local field potential data. This model can successfully generate sleep patterns such as spindles and slow waves (Bazhenov, Timofeev, Steriade, & Sejnowski, 2002; Chen, Chauvette, Skorheim, Timofeev, & Bazhenov, 2012; Bonjean et al., 2011). Since the model has a cortical layer and a thalamic layer, we further assume that the neural signals in the cortical layer are visible by the recording instruments, while those from the thalamic layer are not. This is a reasonable assumption for many experiments. Since the thalamus is a deep brain structure, most experiments involve only measurements from the cortex.

Next, we simulated 1000 Hodgkin-Huxley neurons networks with 80% excitatory neurons and 20% inhibitory neurons sparsely connected. As in real experiments, we recorded simulated calcium signals from only a small percentage of the network (50 neurons) and compared the performance of different methods. In all simulations, our differential covariance-based methods achieve performance better than or similar to the GLM method, and in the LFP and calcium imaging data sets, they achieve the same performance with fewer data samples.

The article is organized as follow. In section 2, we introduce our new methods. In section 3, we show the performance of all methods and explain why our methods perform better. In section 4, we discuss the advantage and generalizability of our methods. We also propose a few improvements for the future.

2  Methods

2.1  Simulation Models Used to Benchmark the Methods

2.1.1  Passive Neuron Model

To validate and test our new methods, we first developed a passive neuron model. Because of its simplicity, we can provide some theoretical proof for why our new class of methods is better. Every neuron in this model has a passive cell body with capacitance and a leakage channel with conductance . Neurons are connected with a directional synaptic conductance ; for example, neuron receives inputs from neurons and :
formula
2.1
Here, we let , and is gaussian noise with a standard deviation of 1. The connection pattern is shown in Figure 3A. There are 60 neurons in this circuit. The first 50 neurons are connected with a pattern in which neuron projects to neurons and . To make the case more realistic, aside from these 50 neurons that are visible, we added 10 more neurons (neuron 51 to 60 in Figure 3A) that are invisible during our estimations (only the membrane voltages of the first 50 neurons are used to estimate connectivity). These 10 neurons send latent inputs to the visible neurons and introduce external correlations into the system. Therefore, we update our passive neuron's model as
formula
2.2
where . We choose the latent input's strength in the same scale as other connections in the network. We tested multiple values between [0, 10]; higher value generates more interference and therefore makes the case more difficult.

We added the latent inputs to the system because unobserved common inputs exist in real-world problems (Stevenson et al., 2008). For example, one could be using two-photon imaging to record calcium signals from the cortical circuit. The cortical circuit might receive synaptic inputs from deeper layers in the brain, such as the thalamus, which is not visible to the microscope. Each invisible neuron projects to many visible neurons, leading to common synaptic currents to the cortical circuit, and causes neurons in the cortical circuit to be highly correlated. We will discuss how to remove interference from the latent inputs.

2.1.2  Thalamocortical Model

To test and benchmark the differential covariance-based methods in a more realistic model, we simulated neural signals from a Hodgkin-Huxley–based spiking neural network. The thalamocortical model used in this study was based on several previous studies, which were used to model spindle and slow wave activity (Bazhenov et al., 2002; Chen et al., 2012; Bonjean et al., 2011). A schematic of the thalamocortical model in this work is shown in Figure 1. As shown, the thalamocortical model was structured as a one-dimensional, multilayer array of cells. The thalamocortical network consisted of 50 cortical pyramidal (PY) neurons, 10 cortical inhibitory (IN) neurons, 10 thalamic relay (TC) neurons, and 10 reticular (RE) neurons. The connections between the 50 PY neurons follow the pattern in the passive neuron model and are shown in Figure 7A. For the rest of the connection types, a neuron connects to all target neurons within the radius listed in Table 1 (Bazhenov et al., 2002). The network is driven by spontaneous oscillations. Details of the model are explained in appendix C.
Figure 1:

Network model for the thalamocortical interactions, with four layers of neurons: thalamocortical (TC), reticular nucleus (RE) of the thalamus, cortical pyramidal (PY), and inhibitory (IN) neurons. Filled circles indicate excitatory synapses; open circles indicate inhibitory synapses.

Figure 1:

Network model for the thalamocortical interactions, with four layers of neurons: thalamocortical (TC), reticular nucleus (RE) of the thalamus, cortical pyramidal (PY), and inhibitory (IN) neurons. Filled circles indicate excitatory synapses; open circles indicate inhibitory synapses.

Table 1:
Connectivity Properties.
PYTCPYRETCPYTCINPYININPYRERETCRERETC
Radius 15 
PYTCPYRETCPYTCINPYININPYRERETCRERETC
Radius 15 

For each simulation, we simulated the network for 600 secs. The data were sampled at 1000 Hz.

2.1.3  Local Field Potential (LFP) Model

To simulate local field potential (LFP) data, we expanded the previous thalamocortical model by 10 times (see Figure 2). Each connection in the previous network (see Figure 1) is transformed into a fiber bundle connecting 10 pairs of neurons in a subnetwork. For cortical pyramidal neurons (PY), we connect each neuron to its four neighboring neurons inside the subnetwork. Other settings for the neurons and the synapses are the same as the previous thalamocortical model. We simulated the network for 600 secs. The data were sampled at 1000 Hz.
Figure 2:

The LFP model was transformed from the thalamocortical model. In this figure, we used the TC PY connection in the red box as an example. Each neuron in the original thalamocortical model was expanded into a subnetwork of 10 neurons. Each connection in the original thalamocortical model was transformed into 10 parallel connections between two subnetworks. Moreover, subnetworks transformed from PY neurons have local connections. Each PY neuron in this subnetwork connects to its two nearest neighbors on each side. We put an LFP electrode at the center of each PY subnetwork. The electrode received neurons' signals inversely proportional to the distance.

Figure 2:

The LFP model was transformed from the thalamocortical model. In this figure, we used the TC PY connection in the red box as an example. Each neuron in the original thalamocortical model was expanded into a subnetwork of 10 neurons. Each connection in the original thalamocortical model was transformed into 10 parallel connections between two subnetworks. Moreover, subnetworks transformed from PY neurons have local connections. Each PY neuron in this subnetwork connects to its two nearest neighbors on each side. We put an LFP electrode at the center of each PY subnetwork. The electrode received neurons' signals inversely proportional to the distance.

We plant an LFP electrode at the center of each of the 50 PY neuron local circuits. The distance between the 10 local neurons is , and the distance between the PY subnetworks is 1 cm. The LFP is calculated according to the standard model in Bonjean et al. (2011), Destexhe (1998), and Nunez and Srinivasan (2005). The LFPs are mainly contributed by elongated dendrites of the cortical pyramidal neurons. In our model, each cortical pyramidal neuron has a 2 mm long dendrite.

For each LFP ,
formula
2.3
where the sum is taken over all excitatory cortical neurons. is the postsynaptic current of neuron , is the distance from the electrode to the center of the dendrite of neuron , and is the distance from the electrode to the soma of neuron .

2.1.4  Calcium Imaging Model

Vogelstein et al. (2009) proposed a transfer function between spike trains and calcium fluorescence signals,
formula
2.4
where is a step influx of calcium molecules at each action potential. is the number of spikes at each time step, is the saturation concentration of calcium, and is a gaussian noise with a standard deviation of 0.000003. Since our data are sampled at 1000 Hz, we can resolve every single action potential. So in our data, there are no multiple spikes at one time step. is the decay constant for calcium molecules. To maintain the information in the differential signal, instead of setting a hard cutoff value and transforming the intracellular voltages to binary spike trains, we use a sigmoid function to transform the voltages to the calcium influx activation parameter (),
formula
2.5
where is the threshold potential.

In real experiments, we can image only a small percentage of neurons in the brain. Therefore, we simulated HH-based networks of 1000 neurons and record from only 50 neurons. We used the four connection patterns (normal-1 normal-4) provided online (https://www.kaggle.com/c/connectomics) (Stetter, Battaglia, Soriano, & Geisel, 2012). Briefly, the networks have 80% excitatory neurons and 20% inhibitory neurons. The connection probability is 0.01—one neuron connects to about 10 neurons—so it is a sparsely connected network.

Similar to our thalamocortical model, we used AMPA and NMDA synapses for the excitatory synapses and GABA synapses for the inhibitory synapses. The simulations ran 600 sec and were sampled at 1000 Hz. The intracellular voltages obtained from the simulations were then transferred to calcium fluorescence signals and downsampled to 50 Hz. For each of the four networks, we conducted 25 recordings. Each recording contains calcium signals from 50 randomly selected neurons.

For accurate estimations, the differential covariance-based methods require the reconstructed action potentials from the calcium imaging. While this is an active research area and many methods have been proposed (Quan, Liu, Lv, Chen, & Zeng, 2010; Rahmati, Kirmse, Marković, Holthoff, & Kiebel, 2016), in this study, we simply reversed the transfer function. By assuming the transfer function from action potentials to calcium fluorescence signals is known, we can reconstruct the action potentials. Given as the observed fluorescence signal,
formula
2.6

2.2  Differential Covariance-Based Methods

In this section, we introduce a new class of methods to estimate the functional connectivity of neurons (code is provided online at https://github.com/tigerwlin/diffCov).

2.2.1  Step 1: Differential Covariance

The input to the method, , is an neural recording data set. is the number of neurons or channels recorded, and is the number of data samples during recordings. We compute the derivative of each time series with . Then the covariance between and is computed and denoted as , which is an matrix defined as
formula
2.7
where is the differential signal of neuron/channel , is the signal of neuron/channel , and is the sample covariance function for two time series. In appendix A, we provide a theoretical proof about why the differential covariance estimator generates fewer false connections than the covariance-based methods do.

2.2.2  Step 2: Applying Partial Covariance Method

As Stevenson et al. (2008) previously mentioned, one problem of the correlation method is the propagation of correlation. Here we designed a customized partial covariance algorithm to reduce this type of error in our methods. We use to denote the differential covariance estimation after applying the partial covariance method.

Using the derivation from section B.2, we have
formula
2.8
where is a set of all neurons/channels except :
formula
and were computed from section 2.2.1, and is the covariance matrix of set . is a flat matrix denoting the covariance of neuron and neurons in set . is the matrix dot product.

As we explain in section B.2, the partial covariance of the differential covariance is not equivalent to the inverse of the differential covariance estimation. The two are equivalent only for the covariance matrix and when the partial correlation is controlling on all variables in the observed set. In our case, the differential covariance matrix is nonsymmetric because it is the covariance between recorded signals and their differential signals. We have signals and differential signals in our observed set; however, we are controlling only on the original signals for the partial covariance algorithm. Due to these differences, we developed this customized partial covariance algorithm, equation 2.8, which performs well for neural signals in the form of equation 2.2.

2.2.3  Step 3: Sparse Latent Regularization

Finally, we applied the sparse latent regularization method to partial covariance version of the differential covariance (Chandrasekaran, Sanghavi, Parrilo, & Willsky, 2011; Yatsenko et al., 2015). As explained in section B.3, in the sparse latent regularization method, we made the assumption that there are observed neurons and unobserved common inputs in a network. If the connections between the observed neurons are sparse and the number of unobserved common inputs is small, this method can separate the covariance into two parts and the sparse matrix is the intrinsic connections between the observed neurons.

Here we define as the sparse result from the method and as the low-rank result from the method. Then we solve
formula
2.9
under the constraint that
formula
2.10
where is the L1-norm of a matrix and is the trace of a matrix. is the penalty ratio between the L1-norm of and the trace of . It was set to for all our estimations. is the partial differential covariance computed from section 2.2.2. We receive a sparse estimation, , of the connectivity.

2.2.4  Computing Derivative

In differential covariance, computing the derivative using provides better estimation results than using . To elaborate this point, we first explain our method's connection to the autoregressive (AR) method.

Detailed discussion of the multivariate autoregressive model and its estimators has been discussed in Harrison et al. (2003). In discrete space, our differential covariance estimator is similar to the mean squared error (MSE) estimator of the AR model. Following the definition in equation 2.2, we have
formula
2.11
where are neurons' membrane voltages. is the connection matrix that describes the conductance between each pair of neurons. is the gaussian noise.
For , we note here:
formula
2.12
As Harrison et al. (2003) explained, in the AR model, the MSE estimator of is . We note that the numerator of this MSE estimator is our differential covariance estimator. Therefore, in this case, the model we proposed is equivalent to a regularized AR(2) model, where the transition matrix of the second order is restricted to be an identity matrix.
In the case, we note that the differential covariance estimator is similar to an AR(1) estimator,
formula
2.13
and the MSE estimator of is . Therefore:
formula
2.14
where the numerator of this MSE estimator is our differential covariance estimator.

As explained above, using different methods to compute the derivative will produce different differential covariance estimators, and they are equivalent to estimators from different AR models for the connection matrix. In section 3, we show that the performances of these two estimators are significantly different.

2.3  Performance Quantification

The performance of each method is judged by four quantified values. The first three values indicate the method's abilities to reduce the three types of false connections (see Figure 4). The last one indicates the method's ability to correctly estimate the true positive connections against all possible interference.

We define as the ground-truth connectivity matrix, where:
formula
2.15
Then we can use a three-dimensional tensor to represent the false connections caused by common inputs. For example, neurons and receive common input from neuron :
formula
2.16
Therefore, we can compute a mask that labels all of the type 1 false connections:
formula
2.17
For the type 2 false connections (e.g., neuron projects to neuron ; then neuron projects to neuron ), the mask is defined as
formula
2.18
or, in simple matrix notation,
formula
2.19
Similar to , the false connections caused by unobserved common inputs are
formula
2.20
Finally, is defined as
formula
2.21
Given a connectivity matrix estimation result, , the four values for the performance quantification are computed as the area under the ROC curve for two sets: the true positive set and the false positive set:
formula
2.22

3  Results

3.1  False Connections in Correlation-Based Methods

When applied to neural circuits, the commonly used correlation-based methods produce systematic false connections. As shown, Figure 3A is the ground truth of the connections in our passive neuron model (neurons 1 to 50 are the observed neurons). Figure 3B is from the correlation method, Figure 3C is the precision matrix, and Figure 3D is the sparse + latent regularized precision matrix. As shown, all of these methods produce extra false connections.
Figure 3:

Ground-truth connections from a passive neuron model and estimations from correlation-based methods. (A) Ground-truth connection matrix. Neurons 1 to 50 are observed neurons. Neurons 51 to 60 are unobserved neurons that introduce common inputs. (B) Estimation from the correlation method. (C) Estimation from the precision matrix. (D) Estimation from the sparse + latent regularized precision matrix. (E) Zoom-in of panel A. (F) Zoom-in of panel B. (G) Zoom-in of panel C. (H) Zoom-in of panel D.

Figure 3:

Ground-truth connections from a passive neuron model and estimations from correlation-based methods. (A) Ground-truth connection matrix. Neurons 1 to 50 are observed neurons. Neurons 51 to 60 are unobserved neurons that introduce common inputs. (B) Estimation from the correlation method. (C) Estimation from the precision matrix. (D) Estimation from the sparse + latent regularized precision matrix. (E) Zoom-in of panel A. (F) Zoom-in of panel B. (G) Zoom-in of panel C. (H) Zoom-in of panel D.

For convenience of explanation, we define the diagonal strip of connections in the ground truth (the first 50 neurons in Figure 3A) as the 3 and 4 diagonal lines, because they are three and four steps away from the diagonal line of the matrix. As shown in Figure 3, all of these methods produce false connections on the diagonal lines. The precision matrix method (see Figure 3C) also has square box shape of false connections in the background.

3.1.1  Type 1 False Connections

The type 1 false connections shown in Figure 4A are produced because two neurons receive the same input from another neuron. The same synaptic current that passes into the two neurons generates a positive correlation between the two neurons. However, there is no physiological connection between these two neurons. In our connection pattern (see Figure 3A), we notice that two neurons next to each other receive common synaptic inputs; therefore, there are false connections on the diagonal lines of the correlation-based estimations.
Figure 4:

Illustrations of the three types of false connections in the correlation-based methods. Solid lines indicate the physical wiring between neurons, and the filled circles at the end indicate the synaptic contacts (i.e., the direction of the connections). The dotted lines are the false connections introduced by the correlation-based methods. (A) Type 1 false connections, which are due to two neurons receiving the same synaptic inputs. (B) Type 2 false connections, which are due to the propagation of correlation. (C) Type 3 false connections, which are due to unobserved common inputs.

Figure 4:

Illustrations of the three types of false connections in the correlation-based methods. Solid lines indicate the physical wiring between neurons, and the filled circles at the end indicate the synaptic contacts (i.e., the direction of the connections). The dotted lines are the false connections introduced by the correlation-based methods. (A) Type 1 false connections, which are due to two neurons receiving the same synaptic inputs. (B) Type 2 false connections, which are due to the propagation of correlation. (C) Type 3 false connections, which are due to unobserved common inputs.

3.1.2  Type 2 False Connections

The type 2 false connections shown in Figure 4B are due to the propagation of correlation. Because one neuron connects to another neuron and neuron connects to another neuron , the correlation method presents the correlation between and , which do not have a physical connection. This phenomenon is shown in Figure 3B as the extra diagonal strips. This problem, shown in Figure 3C can be greatly reduced by the precision matrix/partial covariance method.

3.1.3  Type 3 False Connections

The type 3 false connections shown in Figure 4C are also due to the common currents that pass into two neurons. However, in this case, they are from the unobserved neurons. For this particular passive neuron model, it is due to the inputs from the 10 unobserved neurons (neurons 51–60) as shown in Figure 3A. Because the latent neurons have broad connections to the observed neurons, they introduce a square box shape correlation pattern into the estimations. (See Figure 3C. Figure 3B also contains this error, but it is hard to see.) Since the latent neurons are not observable, the partial covariance method cannot be used to regress out this type of correlation. The sparse latent regularization can be applied if the sparse and low-rank assumption is valid, and the sparse + latent regularized result is shown in Figure 3D. However, even after using this regularization, the correlation-based methods still leave false connections in Figure 3D.

3.2  Estimations from Differential Covariance-Based Methods

Comparing the ground-truth connections in Figure 5A with our final estimation in Figure 5D, we see that our methods essentially transformed the connections in the ground truth into a map of sources and sinks in a network. An excitatory connection, , in our estimations has a negative value for and a positive value for , which means the current is coming out of the source and goes into the sink . We note that there is another ambiguous case, an inhibitory connection , which produces the same results in our estimations. Our methods cannot differentiate these two cases; instead, they indicate sources and sinks in a network.
Figure 5:

Differential covariance analysis of the passive neuron model. The color in panels B, C, D, F, G, and H indicate the direction of the connections. For element , a warm color indicates is the sink and is the source—that is, —and the cool color indicates is the sink and is the source—that is, . (A) Ground-truth connection matrix. (B) Estimation from the differential covariance method. (C) Estimation from the partial differential covariance method. (D) Estimation from the sparse + latent regularized partial differential covariance method. (E) Zoom-in of panel A. (F) Zoom-in of panel B. (G) Zoom-in of panel C. (H) Zoom-in of panel D.

Figure 5:

Differential covariance analysis of the passive neuron model. The color in panels B, C, D, F, G, and H indicate the direction of the connections. For element , a warm color indicates is the sink and is the source—that is, —and the cool color indicates is the sink and is the source—that is, . (A) Ground-truth connection matrix. (B) Estimation from the differential covariance method. (C) Estimation from the partial differential covariance method. (D) Estimation from the sparse + latent regularized partial differential covariance method. (E) Zoom-in of panel A. (F) Zoom-in of panel B. (G) Zoom-in of panel C. (H) Zoom-in of panel D.

3.2.1  The Differential Covariance Method Reduces Type 1 False Connections

By comparing Figure 3B with Figure 5B, we see that the type 1 false connections on the diagonal lines of Figure 3B are reduced in Figure 5B. This is reflecting the theorems we prove in appendix A, in particular theorem 8, which shows that the strength of the type 1 false connections is reduced in the differential covariance method by a factor of . Moreover, the performance of the differential covariance method on reducing type 1 false connections is quantified in Figure 6.
Figure 6:

Performance quantification (area under the ROC curve) of different methods with respect to their abilities to reduce the three types of false connections and their abilities to estimate the true positive connections using the passive neuron data set.

Figure 6:

Performance quantification (area under the ROC curve) of different methods with respect to their abilities to reduce the three types of false connections and their abilities to estimate the true positive connections using the passive neuron data set.

3.2.2  The Partial Covariance Method Reduces Type 2 False Connections

Second, we see that due to the propagation of correlation, there are extra diagonal strips in Figure 5B. These are removed in Figure 5C by applying the partial covariance method. Each estimator's performance for reducing type 2 false connections is quantified in Figure 6.

3.2.3  The Sparse + Latent Regularization Reduces Type 3 False Connections

Third, the sparse + latent regularization to remove the correlation is introduced by the latent inputs. When the observed neurons' connections are sparse and the number of unobserved common inputs is small, the covariance introduced by the unobserved common inputs can be removed. As shown in Figure 5D, the external covariance in the background of Figure 5C is removed, while the true diagonal connections and the directionality of the connections are maintained. This regularization is also effective for correlation-based methods, but type 1 false connections are maintained in the estimation even after applying this regularization (see Figure 3D). Each estimator's performance for reducing type 3 false connections is quantified in Figure 6.

3.2.4  The Differential Covariance-Based Methods Provide Directionality Information of the Connections

Using this passive neuron model, in section A.3, we provide a mathematical explanation for why the differential covariance-based methods provide directional information for the connections. Given an excitatory connection (neuron projects to neuron ), from theorem 9 in section A.3, we have
formula
3.1
We note here that there is another ambiguous setting that provides the same result, which is an inhibitory connection . Conceptually, the differential covariance indicates the current sources and sinks in a neural circuit, but the exact type of synapse is unknown.

3.2.5  Performance Quantification for the Passive Neuron Data Set

In Figure 6, we quantified the performance of the estimators for one example data set. We see that with the increase in the sample size, our differential covariance-based methods reduce the three types of false connections, while maintaining high true positive rates. Also, as we apply more advanced techniques (), the estimator's performance increases in all four panels of the quantification indices. Although the precision matrix and the sparse + latent regularization help the correlation method reduce type 2 and type 3 errors, all correlation-based methods handle the type 1 false connections poorly. We also note that the masks we used to quantify each type of false connections are not mutually exclusive (i.e., there are false connections that belong to more than one type of false connections). Therefore, in Figure 6, it seems that a regularization is reducing the multiple types of false connections. For example, the sparse + latent regularization is reducing both type 2 and type 3 false connections.

In Table 2, we provide quantified results (area under the ROC curve) for two connection patterns (cxcx34, and cxcx56789) and three conductance settings (g5, g30, and g50). We see that the key results in Figure 6 are also generalized here. By applying more advanced techniques to the original differential covariance estimator (), the performance increases with respect to the three types of error, while the true positive rate is not sacrificed. We also note that although the precision matrix and the sparse + latent regularization help the correlation method reduce type 2 and type 3 errors, all correlation-based methods handle the type 1 error poorly.

Table 2:
Performance Quantification (Area under the ROC Curve) of Different Methods with Respect to Their Abilities to Reduce the Three Types of False Connections and Their Abilities to Estimate the True Positive Connections under Five Different Passive Neuron Model Settings.
CovPrecisionPrecision+SL-reg
cxcx34 g5       
Error 1 0.6327 0.3469 0.8776 
Error 2 0.1469 0.9520 0.9915 0.3757 0.8347 1.0000 
Error 3 0.4638 0.9362 0.9797 0.6541 0.8391 0.9986 
True positive 0.7312 1.0000 1.0000 0.9677 0.9946 1.0000 
cxcx34 g30       
Error 1 0.0510 0.5816 0.9490 
Error 2 0.0056 0.8927 0.9972 0.2881 0.9548 1.0000 
Error 3 0.2164 0.9188 0.9942 0.5430 0.9662 1.0000 
True positive 0.5591 1.0000 0.9892 0.6559 1.0000 1.0000 
cxcx34 g50       
Error 1 0.2041 0.6531 
Error 2 0.7034 0.9944 0.0523 0.9054 1.0000 
Error 3 0.3179 0.8000 0.9894 0.4145 0.9309 1.0000 
True positive 0.9140 1.0000 0.9946 0.9516 0.9785 1.0000 
cxcx56789 g5       
Error 1 0.0053 0.0053 0.6895 0.6263 0.8526 
Error 2 0.1896 0.6229 0.7896 0.5240 0.7896 0.9938 
Error 3 0.3573 0.6085 0.7659 0.6957 0.7591 0.9817 
True positive 0.6884 0.9442 0.6674 0.9930 0.9605 0.9837 
cxcx56789 g50       
Error 1 0.0263 0.2816 0.6842 
Error 2 0.1083 0.5312 0.8240 0.2844 0.6990 0.9979 
Error 3 0.4256 0.4927 0.7762 0.5091 0.7116 0.9835 
True positive 0.9256 0.9116 0.6698 0.9395 0.9279 0.9419 
CovPrecisionPrecision+SL-reg
cxcx34 g5       
Error 1 0.6327 0.3469 0.8776 
Error 2 0.1469 0.9520 0.9915 0.3757 0.8347 1.0000 
Error 3 0.4638 0.9362 0.9797 0.6541 0.8391 0.9986 
True positive 0.7312 1.0000 1.0000 0.9677 0.9946 1.0000 
cxcx34 g30       
Error 1 0.0510 0.5816 0.9490 
Error 2 0.0056 0.8927 0.9972 0.2881 0.9548 1.0000 
Error 3 0.2164 0.9188 0.9942 0.5430 0.9662 1.0000 
True positive 0.5591 1.0000 0.9892 0.6559 1.0000 1.0000 
cxcx34 g50       
Error 1 0.2041 0.6531 
Error 2 0.7034 0.9944 0.0523 0.9054 1.0000 
Error 3 0.3179 0.8000 0.9894 0.4145 0.9309 1.0000 
True positive 0.9140 1.0000 0.9946 0.9516 0.9785 1.0000 
cxcx56789 g5       
Error 1 0.0053 0.0053 0.6895 0.6263 0.8526 
Error 2 0.1896 0.6229 0.7896 0.5240 0.7896 0.9938 
Error 3 0.3573 0.6085 0.7659 0.6957 0.7591 0.9817 
True positive 0.6884 0.9442 0.6674 0.9930 0.9605 0.9837 
cxcx56789 g50       
Error 1 0.0263 0.2816 0.6842 
Error 2 0.1083 0.5312 0.8240 0.2844 0.6990 0.9979 
Error 3 0.4256 0.4927 0.7762 0.5091 0.7116 0.9835 
True positive 0.9256 0.9116 0.6698 0.9395 0.9279 0.9419 

3.3  Thalamocortical Model Results

We tested the methods further in a more realistic Hodgkin-Huxley–based model. Because the synaptic conductances in the model are no longer constants but become nonlinear dynamic functions, which depend on the presynaptic voltages, the previous derivations can be considered only a first-order approximation.

The ground-truth connections between the cortical neurons are shown in Figure 7. These neurons also receive latent inputs from and send feedback currents to inhibitory neurons in the cortex (IN) and thalamic neurons (TC). For clarity of representation, these latent connections are not shown here; the detailed connections are described in section 2.
Figure 7:

Analysis of the thalamocortical model with correlation-based methods. (A) Ground-truth connections of the PY neurons in the thalamocortical model. (B) Estimation from the correlation method. (C) Estimation from the precision matrix method. (D) Estimation from the sparse + latent regularized precision matrix method. (E) Zoom-in of panel A. (F) Zoom-in of panel B. (G) Zoom-in of panel C. (H) Zoom-in of panel D.

Figure 7:

Analysis of the thalamocortical model with correlation-based methods. (A) Ground-truth connections of the PY neurons in the thalamocortical model. (B) Estimation from the correlation method. (C) Estimation from the precision matrix method. (D) Estimation from the sparse + latent regularized precision matrix method. (E) Zoom-in of panel A. (F) Zoom-in of panel B. (G) Zoom-in of panel C. (H) Zoom-in of panel D.

Similar to the passive neuron model, in Figure 7B, the correlation method still suffers from those three types of false connections. As shown, the latent inputs generate false correlations in the background, and the diagonal line false connections, which are due to the common currents, exist in all correlation-based methods (see Figures 7B7D). Comparing Figures 7C and 7D because the type 1 false connections are strong in the Hodgkin-Huxley–based model, the sparse + latent regularization removed the true connections but kept these false connections in its final estimation.

As shown in Figure 8B, differential covariance method reduces the type 1 false connections. In Figure 8C, the partial differential covariance method reduces type 2 false connections in Figure 8B (the yellow connections around the red strip in Figure 8B). Finally, in Figure 8D, the sparse latent regularization removes the external covariance in the background of Figure 8C. The current sources (positive value, red) and current sinks (negative value, blue) in the network are also indicated on our estimators.
Figure 8:

Analysis of the thalamocortical model with differential covariance-based methods. The color in panels B, C, D, F, G, and H indicates the direction of the connections. For element , a warm color indicates is the sink and is the source—. A cool color indicates is the sink and is the source—. (A) Ground-truth connection matrix. (B) Estimation from the differential covariance method. (C) Estimation from the partial differential covariance method. (D) Estimation from the sparse + latent regularized partial differential covariance method. (E) Zoom-in of panel A. (F) Zoom-in of panel B. (G) Zoom-in of panel C. (H) Zoom-in of panel D.

Figure 8:

Analysis of the thalamocortical model with differential covariance-based methods. The color in panels B, C, D, F, G, and H indicates the direction of the connections. For element , a warm color indicates is the sink and is the source—. A cool color indicates is the sink and is the source—. (A) Ground-truth connection matrix. (B) Estimation from the differential covariance method. (C) Estimation from the partial differential covariance method. (D) Estimation from the sparse + latent regularized partial differential covariance method. (E) Zoom-in of panel A. (F) Zoom-in of panel B. (G) Zoom-in of panel C. (H) Zoom-in of panel D.

In Figure 9, each estimator's performance on each type of false connections is quantified. In this example, our differential covariance-based methods achieve similar performance to the GLM method.
Figure 9:

Performance quantification (area under the ROC curve) of different methods with respect to their abilities to reduce the three types of false connections and their abilities to estimate the true positive connections using the thalamocortical data set.

Figure 9:

Performance quantification (area under the ROC curve) of different methods with respect to their abilities to reduce the three types of false connections and their abilities to estimate the true positive connections using the thalamocortical data set.

3.4  Simulated LFP Results

For population recordings, our methods have similar performance to the thalamocortical model example. While the correlation-based methods are still suffering from the problem of type 1 false connections (see Figure 10), our differential covariance-based methods can reduce all three types of false connections (see Figure 11). In Figure 12, each estimator's performance on LFP data is quantified. In this example, with sufficient data samples, our differential covariance-based methods achieve performance similar to that of the GLM method. For smaller sample sizes, our new methods perform better than the GLM method.
Figure 10:

Analysis of the simulated LFP data with correlation-based methods. (A) Ground-truth connection matrix. (B) Estimation from the correlation method. (C) -score of the correlation matrix. (D) Estimation from the precision matrix method. (E) Estimation from the sparse + latent regularized precision matrix method. (F) Zoom-in of panel A. (G) Zoom-in of panel B. (H) Zoom-in of panel C. (I) Zoom-in of panel D. (J) Zoom-in of panel E.

Figure 10:

Analysis of the simulated LFP data with correlation-based methods. (A) Ground-truth connection matrix. (B) Estimation from the correlation method. (C) -score of the correlation matrix. (D) Estimation from the precision matrix method. (E) Estimation from the sparse + latent regularized precision matrix method. (F) Zoom-in of panel A. (G) Zoom-in of panel B. (H) Zoom-in of panel C. (I) Zoom-in of panel D. (J) Zoom-in of panel E.

Figure 11:

Analysis of the simulated LFP data with differential covariance-based methods. The color in panels B, C, D, F, G, and H indicates the direction of the connections. For element , a warm color indicates is the sink and is the source—. A cool color indicates is the sink and is the source—. (A) Ground-truth connection matrix. (B) Estimation from the differential covariance method. (C) Estimation from the partial differential covariance method. (D) Estimation from the sparse + latent regularized partial differential covariance method. (E) Zoom-in of panel A. (F) Zoom-in of panel B. (G) Zoom-in of panel C. (H) Zoom-in of panel D.

Figure 11:

Analysis of the simulated LFP data with differential covariance-based methods. The color in panels B, C, D, F, G, and H indicates the direction of the connections. For element , a warm color indicates is the sink and is the source—. A cool color indicates is the sink and is the source—. (A) Ground-truth connection matrix. (B) Estimation from the differential covariance method. (C) Estimation from the partial differential covariance method. (D) Estimation from the sparse + latent regularized partial differential covariance method. (E) Zoom-in of panel A. (F) Zoom-in of panel B. (G) Zoom-in of panel C. (H) Zoom-in of panel D.

Figure 12:

Performance quantification (area under the ROC curve) of different methods with respect to their abilities to reduce the three types of false connections and their abilities to estimate the true positive connections using the simulated LFP data set.

Figure 12:

Performance quantification (area under the ROC curve) of different methods with respect to their abilities to reduce the three types of false connections and their abilities to estimate the true positive connections using the simulated LFP data set.

3.5  Simulated Calcium Imaging Results

Because current techniques allow recording of only a small percentage of neurons in the brain, we tested our methods on a calcium imaging data set of 50 neurons recorded from networks of 1000 neurons. In this example, our differential covariance-based methods (see Figure 14) match better with the ground truth than the correlation-based methods (see Figure 13).
Figure 13:

Analysis of the simulated calcium imaging data set with correlation-based methods. (A) Ground-truth connection matrix. (B) Estimation from the correlation method. (C) Estimation from the precision matrix method. (D) Estimation from the sparse + latent regularized precision matrix method. For clarity, panels B, C, and D are thresholded to show only the strongest connections, so one can compare it with the ground truth.

Figure 13:

Analysis of the simulated calcium imaging data set with correlation-based methods. (A) Ground-truth connection matrix. (B) Estimation from the correlation method. (C) Estimation from the precision matrix method. (D) Estimation from the sparse + latent regularized precision matrix method. For clarity, panels B, C, and D are thresholded to show only the strongest connections, so one can compare it with the ground truth.

Figure 14:

Analysis of the simulated calcium imaging data set with differential covariance-based methods. The color in panels B–D indicates the direction of the connections. For element , a warm color indicates is the sink and is the source—. A cool color indicates is the sink and is the source—. (A) Ground truth connection matrix. (B) Estimation from the differential covariance method. (C) Estimation from the partial differential covariance method. (D) Estimation from the sparse + latent regularized partial differential covariance method. For clarity, panels B, C, and D are thresholded to show only the strongest connections, so one can compare it with the ground truth.

Figure 14:

Analysis of the simulated calcium imaging data set with differential covariance-based methods. The color in panels B–D indicates the direction of the connections. For element , a warm color indicates is the sink and is the source—. A cool color indicates is the sink and is the source—. (A) Ground truth connection matrix. (B) Estimation from the differential covariance method. (C) Estimation from the partial differential covariance method. (D) Estimation from the sparse + latent regularized partial differential covariance method. For clarity, panels B, C, and D are thresholded to show only the strongest connections, so one can compare it with the ground truth.

We performed 25 sets of recordings with 50 neurons randomly selected in each of the four large networks and quantified the results. The markers on the plots in Figure 15 are the average area under the ROC curve values across the 100 sets, and the error bars indicate the standard deviations across these 100 sets of recordings. Our differential covariance-based methods perform better than the GLM method, and the performance differences seem to be greater in situations with fewer data samples.
Figure 15:

Performance quantification (area under the ROC curve) of different methods with respect to their abilities to reduce the three types of false connections and their abilities to estimate the true positive connections using the simulated calcium imaging data set. The error bar is the standard deviation across 100 sets of experiments. Each experiment randomly recorded 50 neurons in a large network. The markers on the plots indicate the average area under the ROC curve values across the 100 sets of experiments.

Figure 15:

Performance quantification (area under the ROC curve) of different methods with respect to their abilities to reduce the three types of false connections and their abilities to estimate the true positive connections using the simulated calcium imaging data set. The error bar is the standard deviation across 100 sets of experiments. Each experiment randomly recorded 50 neurons in a large network. The markers on the plots indicate the average area under the ROC curve values across the 100 sets of experiments.

4  Discussion

4.1  Generalizability and Applicability of the Differential Covariance-Based Methods to Real Experimental Data

Many methods have been proposed to solve the problem of reconstructing the connectivity of a neural network. Winterhalder et al. (2005) reviewed the nonparametric methods and Granger causality–based methods. Much progress has been made recently using the kinetic Ising model and Hopfield network (Huang, 2013; Dunn & Roudi, 2013; Battistin, Hertz, Tyrcha, & Roudi, 2015; Capone, Filosa, Gigante, Ricci-Tersenghi, & Del Giudice, 2015; Roudi & Hertz, 2011) with sparsity regularization (Pernice & Rotter, 2013). The GLM method (Okatan et al., 2005; Truccolo et al., 2005; Pillow et al., 2008) and the maximum entropy method (Schneidman, Berry, Segev, & Bialek, 2006) are two popular classes of methods and the main modern approaches for modeling multi-unit recordings (Roudi, Dunn, & Hertz, 2015).

In current research, people are recording more and more neurons and looking for new data analysis techniques to handle bigger data with higher dimensionality. The field is in favor of algorithms that require fewer samples and scale well with dimensionality, but without sacrificing accuracy. Also, an algorithm that is model free or makes minimum assumptions about the hidden structure of the data has the potential to be applied to multiple types of neural recordings.

The key difference between our methods and other ones is that we use the relationship between a neuron's differential voltage and its voltage rather than finding the relationship between voltages. This provides better performance because the differential voltage is a proxy for a neuron's synaptic current. Also, the relationship between a neuron's synaptic current and its input voltages is more linear, which is suitable for data analysis techniques like the covariance method. While this linear relationship holds only for our passive neuron model, we still see similar or better performance of our methods in our examples based on the Hodgkin-Huxley model, where we relaxed this assumption and allowed loops in the networks. This implies that this class of methods is still applicable even when the ion channels' conductances vary nonlinearly with the voltages, which makes the linear relationship hold only weakly.

4.2  Caveats and Future Directions

One open question for the differential covariance-based methods is how to improve the way they handle neural signals that are nonlinearly transformed from the intracellular voltages. Currently, to achieve good performance in the calcium imaging example, we need to assume we know the exact transfer function and reversely reconstruct the action potentials. We find this reverse transform method to be prone to additive gaussian noise. Further study is needed to find a better way to preprocess calcium imaging data for the differential covariance–based methods.

Throughout our simulations, the differential covariance method has better performance and needs fewer data samples in some simulations but not in others. Further investigation is needed to understand where the improvement comes from.

Our Hodgkin-Huxley simulations did not include axonal or synaptic delays, a critical feature of a real neural circuit. Unfortunately, it is nontrivial to add this feature to our Hodgkin-Huxley model. Nevertheless, we tested our methods with the passive neuron model using the same connection patterns but with random synaptic delays between neurons. In appendix D, we show that for up to a 10 ms uniformly distributed synaptic delay pattern, our methods still outperform the correlation-based methods.

Appendix A:  Differential Covariance Derivations

In this appendix we first build a simple three-neuron network to demonstrate that our differential covariance–based methods can reduce type 1 false connections. Then we develop a generalized theory showing that the type 1 false connections' strength is always lower in our differential covariance–based methods than the original correlation-based methods.

A.1  A Three-Neuron Network

Let us assume a network of three neurons, where neuron A projects to neurons B and C:
formula
A.1
Here, the cell conductance is , neuron A's synaptic connection strength to neuron B is , and neuron A's synaptic connection strength to neuron C is . are independent white gaussian noises.
From equation 18 of Fan, Shan, Yuan, and Ren (2011), we can derive the covariance matrix of this network:
formula
A.2
where
formula
A.3
is the transpose of the ground-truth connection of the network. And
formula
A.4
since each neuron receives independent noise. is an identity matrix of the size of , and is an identity matrix of the size of . is the Kronecker product, and is the column vectorization function.
Therefore, we have the covariance matrix of the network as:
formula
A.5
When computing the differential covariance, we plug in equation A.1—for example:
formula
A.6
Therefore, from equation A.5, we can compute the differential covariance as
formula
A.7

Notice that because the ratio between and is , in differential covariance, the type 1 false connection has value 0.

A.2  Type 1 False Connection's Strength Is Reduced in Differential Covariance

In this section, we propose a theory. A network consists of passive neurons in the following form,
formula
A.8
where is the set of neurons that project to neuron , is the synaptic conductance for the projection from neuron to neuron , and is a Brownian motion. Further assume that:
  • All neurons' leakage conductance and membrane capacitance are constants and the same.

  • There is no loop in the network.

  • , where is the maximum of , for

Then we prove below that:

  • For two neurons that have physical connection, their covariance is .

  • For two neurons that do not have physical connection, their covariance is .

  • For two neurons that have physical connection, their differential covariance is .

  • For two neurons that do not have physical connection, their differential covariance is .

  • The type 1 false connection's strength is reduced in differential covariance.

Lemma 1.
The asymptotic autocovariance of a neuron is
formula
A.9
Proof.
From equation 9 of Fan et al. (2011), we have
formula
A.10
where is the expectation operation. is dropped from when the meaning is unambiguous.
From theorem 4 of Fan et al. (2011), integrating by parts using Itô calculus gives
formula
A.11
Taking the expectation of both sides with equation A.10 gives
formula
A.12
When , equation A.12 becomes
formula
A.13
Lemma 2.
The asymptotic covariance between two neurons is
formula
A.14
Proof.
From equation 9 of Fan et al. (2011), we have
formula
A.15
From theorem 4 of Fan et al. (2011), integrating by parts using Itô calculus gives
formula
A.16
Taking the expectation of both sides with equation A.15 gives
formula
A.17
When , equation A.17 becomes
formula
A.18
Theorem 1.

The autocovariance of a neuron is . The covariance of two different neurons with or without physical connection is .

Proof.

We prove this by induction.

The basis. The base case contains two neurons:
formula
A.19
From lemma 1, we have
formula
A.20
Then, from lemma 2, we have
formula
A.21
From lemma 1, we have
formula
A.22

So the statement holds.

The inductive step. If the statement holds for a network of neurons, we add one more neuron to it.

Part 1. First, we prove that the covariance of any neuron with is also .

From lemma 2, we have:
formula
A.23
where are the neurons projecting to neuron and are the neurons projecting to neuron .

Note that because are neurons from the old network, is at most , and it is only when .

Now we need to prove that is also . We prove this by contradiction. Suppose that is larger than . Then similar to equation A.23, we have:
formula
A.24

Here we separate the problem into two situations: cases 1 and 2.

Case 1: Neuron i projects to neuron n Since there is no loop in the network, . Therefore, is the covariance of two neurons from the old network and is . must be larger than , such that is larger than .

Therefore, if a neuron's covariance with neuron is larger than , one of its antecedents' covariance with neuron is also larger than . Since we assume there is no loop in this network, there must be at least one antecedent (say, neuron ) whose covariance with neuron is larger than and it has no antecedent.

However, from lemma 2:
formula
A.25
Since is , is , which is smaller than . This is a contradiction. So is no larger than . Therefore, is .

Case 2: Neuron i does not project to neuron n Now, in equation A.24, it is possible that . However, in case 1, we just proved that the covariance of any neuron that projects to neuron n is . Therefore, is regardless of whether .

Then, similar to case 1, there must be an antecedent of neuron (say neuron ), whose covariance with neuron is larger than and it has no antecedent. Then from equation A.25, we know this is a contradiction. So is .

Part 2: Then for the autocovariance of neuron ,
formula
A.26

As we already proved that any neuron's covariance with neuron is , the dominant term in equation A.26 is . Therefore, the autocovariance of neuron is also .

Theorem 2.

The covariance of two neurons that are physically connected is . The covariance of two neurons that are not physically connected is .

Proof.
From lemma 2, we have:
formula
A.27

If neuron and neuron are physically connected, let's say ; then . Thus, one of the for is . Therefore, is . Since there is no loop in the network, , so is . Therefore, is .

If neuron and neuron are not physically connected, we have , so is . And , so is . Therefore, is .

Lemma 3.
The differential covariance of two neurons,
formula
A.28
Proof.
From lemma 2, we have
formula
A.29

Theorem 3.

The differential covariance of two neurons that are physically connected is .

Proof.
Assume two neurons have physical connection as . The differential covariance of them is
formula
A.30

From theorem 4, we know is , and

  • If is projecting to , is .

  • If is not projecting to , is .

Therefore, the dominant term is , and is .

From lemma 5,
formula
A.31
Therefore, is .
Theorem 4.

The differential covariance of two neurons that are not physically connected is .

Proof.
First, we define the antecedents of neurons and . Shown in Figure 16, is the set of common antecedents of neuron and neuron . is the set of antecedents of . is the set of exclusive antecedents of neuron . is the set of exclusive antecedents of neuron .
Figure 16:

The network used in the proof of theorem 7.

Figure 16:

The network used in the proof of theorem 7.

From lemma 2, we have, for any neuron ,
formula
A.32
formula
A.33
For simplicity, we define
formula
A.34
From lemma 2, we also have
formula
A.35
For simplicity, we define
formula
A.36
Plug in equations A.32 and A.33 to A.35, and we have
formula
A.37
Now we look at the differential covariance between neuron and neuron :
formula
A.38
Plug in equations A.33 and A.37 and we have:
formula
A.39

Note:

  • There is no physical connection between a neuron in and neuron ; otherwise, this neuron belongs to . Therefore, from theorem 4, is .

  • There is no physical connection between a neuron in and neuron ; otherwise, this neuron belongs to . Therefore, from theorem 4, is .

  • There could be physical connections between neurons in and , so is .

  • There could be physical connections between neurons in and , so is .

  • There could be physical connections between neurons in and neuron , so is .

  • There could be physical connections between neurons in and neuron , so is .

Therefore,
formula
A.40

From lemma 5, we know that is also .

Theorem 5.

Type 1 false connection's strength is reduced in differential covariance.

Proof.

From theorems 3 and 4, we know that in the correlation method, the strength of a nonphysical connection is times that of a physical connection .

From theorems 6 and 7, we know that in the differential covariance method, the strength of a nonphysical connection is times that of a physical connection .

Because , the relative strength of the nonphysical connections is reduced in the differential covariance method.

A.3  Directionality Information in Differential Covariance

Theorem 6.

If neuron projects to neuron with an excitatory connection, and .

Proof.
Given the model above, similar to theorem 7, from lemma 2, we have, for any neuron :
formula
A.41
formula
A.42
For simplicity, we define
formula
A.43
From lemma 2, we also have
formula
A.44
For simplicity, we define
formula
A.45
Plug in , , , , , and :
formula
A.46
Now we look at the differential covariance between neuron and neuron :
formula
A.47

Note that in theorem 7, we already proved the scale of . Also:

  • There are physical connections between neurons in and neuron , so is .

  • From lemma 1, the autocovariance of neuron is . So is .

Therefore, is the dominant term in . Since , for excitatory connection , .

From lemma 5,
formula
A.48
Therefore, .
Corollary 1.

If neuron projects to neuron with an inhibitory connection, and .

Proof.

The proof is similar to theorem 9. Again, we know is the dominant term in . Since , for an inhibitory connection , .

From lemma 5,
formula
A.49
Therefore, .

Appendix B:  Benchmarked Methods

We compared our methods to a few popular methods.

B.1  Covariance Method

The covariance matrix is defined as
formula
B.1
where and are two variables and and are their population mean.

B.2  Precision Matrix

The precision matrix is the inverse of the covariance matrix:
formula
B.2
It can be considered one kind of partial correlation. Here we briefly review this derivation because we use it to develop our new method. The derivation is based on and adapted from Cox and Wermuth (1996).

We begin by considering a pair of variables and remove the correlation in them introduced from a control variable .

First, we define the covariance matrix as
formula
B.3
By solving the linear regression problem,
formula
B.4
we have
formula
B.5
Then we define the residual of as
formula
B.6
Therefore, the covariance of is
formula
B.7
If we define the precision matrix as
formula
B.8
using Cramer's rule, we have
formula
B.9
Therefore,
formula
B.10

So and differ by a ratio of .

B.3  Sparse Latent Regularization

Prior studies (Banerjee, Ghaoui, d'Aspremont, & Natsoulis, 2006; Friedman, Hastie, & Tibshirani, 2008) have shown that regularizations can provide better a estimation if the ground-truth connection matrix has a known structure (e.g., sparse). For all data tested in this letter, the sparse latent regularization (Yatsenko et al., 2015) worked best. For a fair comparison, we applied the sparse latent regularization to both the precision matrix method and our differential covariance method.

In the original sparse latent regularization method, people made the assumption that a larger precision matrix is the joint distribution of the observed neurons and latent neurons (Yatsenko et al., 2015):
formula
where corresponds to the observable neurons. If we can measure only the observable neurons, the partial correlation computed from the observed neural signals is
formula
B.11
because the invisible latent neurons as shown in equation 2.2 introduce correlations into the measurable system. We denote this correlation introduced from the latent inputs as
formula
B.12
If we can make the assumption that the connections between the visible neurons are sparse ( is sparse and the number of latent neurons is much smaller than the number of visible neurons, that is, ), then prior work (Chandrasekaran et al., 2011) has shown that if is known, is sparse enough and L's rank is low enough (within the bound defined in Chandrasekaran et al., 2011); then the solution of
formula
B.13
is uniquely defined and can be solved by the following convex optimization problem,
formula
B.14
under the constraint that
formula
B.15
Here, is the L1-norm of a matrix, and is the trace of a matrix. is the penalty ratio between the L1-norm of and the trace of L and is set to for all our estimations.
However, this method is used to regularize precision matrix. For our differential covariance estimation, we need to make small changes to the derivation. Note that if we assume that the neural signals of the latent neurons are known and let be the indexes of these latent neurons, from section 2.2.2,
formula
B.16
removes the terms in equation 2.2.
Even if is unknown,
formula
is low rank because it is bounded by the dimensionality of , which is . And is the internal connection between the visible neurons, which should be a sparse matrix. Therefore, letting
formula
B.17
we can use the original sparse + latent method to solve for . In this article, we used the inexact robust PCA algorithm (http://perception.csl.illinois.edu/matrix-rank/sample_code.html) to solve this problem (Lin, Liu, & Su, 2011).

B.4  The Generalized Linear Model Method

As summarized by Roudi et al. (2015), GLMs assume that every neuron spikes at a time-varying rate that depends on earlier spikes (both those of other neurons and its own) and on external covariates (such as a stimulus or other quantities measured in the experiment). As they explained, the influence of earlier spikes on the firing probability at a given time is assumed to depend on the time since they occurred. For each prepostsynaptic pair , , it is described by a function of this time lag (Roudi et al., 2015). In this article, we average this temporal dependent function over time to obtain the functional connectivity estimation of this method.

The spike trains used for the GLM method were computed using the spike detection algorithm from Quiroga, Nadasdy, and Ben-Shaul (2004) (https://github.com/csn-le/wave_clus). The article's default parameter set is used except that the maximum detection threshold is increased. The code is provided in our github repository (https://github.com/tigerwlin/diffCov/blob/master/spDetect.m) for reference. The GLM code was obtained from Pillow et al. (2008) (http://pillowlab.princeton.edu/code_GLM.html) and applied on the spike trains.

Appendix C:  Details about the Thalamocortical Model

C.1  Intrinsic Currents

For the thalamocortical model, a conductance-based formulation was used for all neurons. The cortical neuron consisted of two compartments, dendritic and axo-somatic compartments, similar to previous studies (Bazhenov et al., 2002; Chen et al., 2012; Bonjean et al., 2011) and is described by the following equations,
formula
C.1
where the subscripts and correspond to axo-somatic and dendritic compartment, I is the Cl leak currents, I is fast Na channels, I is persistent sodium current, I is fast delayed rectifier K current, I is slow voltage-dependent noninactivating K current, I is slow Ca-dependent K current, I is high-threshold Ca current, I is hyperpolarization-activated depolarizing current, and I is the sum of synaptic currents to the neuron. All intrinsic currents were of the form , where is the conductance, is the voltage of the corresponding compartment, and is the reversal potential. The detailed descriptions of individual currents are provided in previous publications (Bazhenov et al., 2002; Chen et al., 2012). The conductances of the leak currents were 0.007 mS/cm for I and 0.023 mS/cm for I. The maximal conductances for different currents were I: 2.0 mS/cm; I: 0.8 mS/cm; I: 0.012 mS/cm; I: 0.015 mS/cm; I: 0.012 mS/cm; I: 3000 mS/cm; I: 200 mS/cm; and I: 15 mS/cm. C was .
The following describes the IN neurons:
formula
C.2
The conductances for leak currents for IN neurons were 0.034 mS/cm for I and 0.006 mS/cm for I. Maximal conductances for other currents were I : 0.8 mS/cm