## Abstract

Measuring functional connectivity from fMRI recordings is important in understanding processing in cortical networks. However, because the brain's connection pattern is complex, currently used methods are prone to producing false functional connections. We introduce differential covariance analysis, a new method that uses derivatives of the signal for estimating functional connectivity. We generated neural activities from dynamical causal modeling and a neural network of Hodgkin-Huxley neurons and then converted them to hemodynamic signals using the forward balloon model. The simulated fMRI signals, together with the ground-truth connectivity pattern, were used to benchmark our method with other commonly used methods. Differential covariance achieved better results in complex network simulations. This new method opens an alternative way to estimate functional connectivity.

## 1 Introduction

Analysis of fMRI recordings to extract functional connectivity aims to characterize a pattern of statistical associations among time series measured from different brain regions (Reid et al., 2019). It has been used quite successfully in identifying brain network interactions at rest and during certain behaviors (Raichle, 2015). Moreover, interesting geometric properties of brain networks can be studied with this method (Achard, Salvador, Whitcher, Suckling, & Bullmore, 2006).

Many statistical methods have been proposed for estimating functional connectivity, including covariance-based methods, lag-based methods (Seth, 2010; Geweke, 1984), and Bayes NET methods (Mumford & Ramsey, 2014). However, without knowing the ground truth, it has been difficult to assess the validity of these methods. To compare these methods, we simulated BOLD signals from network models with known patterns of connectivity (Friston, Harrison, & Penny, 2003). The simulated neural signals were then transformed to BOLD sgnals with an fMRI forward model using the nonlinear balloon model (Buxton, Wong, & Frank, 1998). In a previous study (Smith et al., 2011), many well-known statistical methods were benchmarked with simulated fMRI data sets and concluded that covariance-based methods performed best.

As Stevenson, Rebesco, Miller, and Körding (2008) reviewed, a fundamental problem for covariance-based methods is that two correlated nodes do not necessarily have a direct physical connection. This could be problematic when one wants to compare the functional connectivity with the anatomical wiring of the brain such as fiber tracts or neural projections. We introduce a new method to reduce these false connections and provide estimates that are closer to ground truth in simulated fMRI data. Our new method is a novel form of functional connectivity with weighted strength that extracts information about the directionality of causal influence. This method, differential covariance, performs better than other covariance-based methods on synthetic data by taking advantage of differential signals.

This article is organized as follows. In section 2 we explain how synthetic BOLD signals were generated and introduce our new method. We also explain how results from different estimators were quantified. In section 3, we compare the performance of our new method with the performance of other methods. In section 4, we discuss the advantages and generalizability of the differential covariance.

## 2 Methods

This section describes in detail the methods of generating synthetic BOLD signals, estimating functional connectivity and performance quantification.

### 2.1 BOLD Signal Simulations

#### 2.1.1 Dynamic Causal Modeling

Parameter . | Description . | Value . |
---|---|---|

$\kappa $ | Rate of signal decay | 0.65 per s |

$\gamma $ | Rate of flow-dependent elimination | 0.41 per s |

$\tau $ | Hemodynamic transit time | 0.98 per s |

$\alpha $ | Grubbs exponent | 0.32 |

$Eo$ | Resting oxygen extraction fraction | 0.5 |

$\epsilon $ | Ratio of intra- and extravascular signal | 0.5 |

$\u03d1o$ | Frequency offset at the outer surface of the magnetized | 40.3 per s |

vessel for fully deoxygenated blood | ||

TE | Time of echo | 40 ms |

Parameter . | Description . | Value . |
---|---|---|

$\kappa $ | Rate of signal decay | 0.65 per s |

$\gamma $ | Rate of flow-dependent elimination | 0.41 per s |

$\tau $ | Hemodynamic transit time | 0.98 per s |

$\alpha $ | Grubbs exponent | 0.32 |

$Eo$ | Resting oxygen extraction fraction | 0.5 |

$\epsilon $ | Ratio of intra- and extravascular signal | 0.5 |

$\u03d1o$ | Frequency offset at the outer surface of the magnetized | 40.3 per s |

vessel for fully deoxygenated blood | ||

TE | Time of echo | 40 ms |

#### 2.1.2 Benchmark Data Set from Previous Work

Based on DCM, using different connection patterns, a benchmark data set of 28 simulations was previously created (Smith et al., 2011). We slightly modified the simulation process to evaluate the performance of our method and previous methods over this popular benchmark data set.

First, the original binary external stimulus was replaced by a random gaussian input because our method cannot deal with sharp state transitions. However, our method is able to give a correct estimation in a network with gradual state transition, as in the case of the thalamocortical model simulations (see section 2.1.3).

Second, we reduced the sampling interval (TR) from 3 seconds to 1 second. This is a reasonable modification because an increasing number of fMRI studies (such as the Human Connectome Project) are adopting a densely sampled protocol with TR less than or equal to 1 second. Third, because the new method needs more data to provide a good estimation, simulation duration was increased from 600 seconds to 6000 seconds for most of the simulations. Fourth, the thermal noise added to the BOLD signal was changed from 1% to 0.1%. All other procedures in the original paper (Smith et al., 2011) were followed. Twenty-eight different connection patterns were first simulated using DCM and passed through the forward balloon model, as shown in equation 2.2.

#### 2.1.3 Thalamocortical Model

To further validate our method, a more realistic Hodgkin-Huxley-based ionic neural network model was used in this article to generate neural signals, and the neural signals were then passed through the forward balloon model (see equation 2.2) to generate the BOLD signals. The thalamocortical model used in this study was based on several previous studies (Bazhenov, Timofeev, Steriade, & Sejnowski, 2002; Chen, Chauvette, Skorheim, Timofeev, & Bazhenov, 2012; Bonjean et al., 2011). It was originally used to model spindle and slow wave activity. The thalamocortical model was structured as a one-dimensional, multilayer array of cells. It consisted of 50 cortical pyramidal (PY) neurons, 10 cortical inhibitory (IN) neurons, 10 thalamic relay (TC) neurons, and 10 reticular (RE) neurons. The connection between the 50 PY neurons is shown in Figure 6A. For the rest of the connection types, a neuron connects to all target neurons within a predefined radius, as described in Bazhenov et al. (2002). The details of this model can be found in our previous publication (Lin, Das, Krishnan, Bazhenov, & Sejnowski, 2017).

For each simulation, we simulated the network for 4000 seconds with 600 seconds at each sleep stage (Awake$\u2192$Stage 2$\u2192$Slow Wave$\u2192$REM$\u2192$Stage 2) and 200 seconds of transition period between stages. The simulated neural signals were transferred to BOLD signals and sampled at 10 Hz.

### 2.2 New Method for Functional Connectivity Estimation

In our new method, we first build a backward model to reconstruct the neural signals from the BOLD signals. Then we apply our differential covariance method to the reconstructed neural signals to estimate the functional connectivity.

#### 2.2.1 Backward Model

#### 2.2.2 Differential Covariance

In our previous work (Lin et al., 2017), we demonstrated that when derivative signals are used, the differential covariance method can reduce the false functional connections in a simulated network of neuronal activities. Here we briefly explain this method.

Note that different from our previous work (Lin et al., 2017), the differential covariance is computed using the reconstructed neural signal, $z=q1dVdt+q0V$, instead of the neural signal $V$, which is not recoverable using the balloon model.

#### 2.2.3 Partial Differential Covariance

As Stevenson et al. (2008) previously noted, one problem of the covariance method is the propagation of covariance. Here we designed a customized partial covariance algorithm to reduce this type of error in our differential covariance method. We use $\Delta Pi,j$ to denote the partial covariance between $dzi(t)$ and $zj(t)$.

#### 2.2.4 Sparse Latent Regularization

Finally, we applied the sparse latent regularization method to our estimation (Chandrasekaran, Sanghavi, Parrilo, & Willsky, 2011; Yatsenko et al., 2015). As explained in section A.4, during recording, there are observable nodes and latent nodes in a network. If the connections between the observable nodes are sparse and the number of latent nodes is small, this method can separate the covariance matrix into these two parts, and one can regard the sparse matrix as the intrinsic connections between the observable nodes.

### 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 2). The last one indicates the method's ability to correctly estimate the true positive connections against all possible interference.

## 3 Results

### 3.1 False Connections in Covariance-Based Methods

It is a well-known problem that the commonly used covariance-based methods produce systematic false connections (Stevenson et al., 2008). Shown in Figure 3A is the ground truth of the connection in our DCM model (nodes 1 to 50 are the observable nodes). Figure 3B shows the results from the covariance method, Figure 3C is the precision matrix, and Figure 3D is the ICOV (sparsely regularized precision matrix). Compared to the ground truth, all of these methods produce extra false connections. Here we briefly review the cause of these three types of false connections. The detailed mechanisms have been well reviewed previously (Stevenson et al., 2008).

We define the diagonal strip of connections in the ground truth (Figure 3A) as the $+$6 to $+$10 diagonal lines, because they are 6 to 10 steps away from the diagonal line of the matrix.

Shown in Figure 2A, confounder errors are produced because two nodes receive the same input from another node. The same input that passes into the two nodes generates positive correlation between the two nodes. However, there is no physiological connection between these two nodes. For example, as shown in Figure 3A, node 1 projects to nodes 6 and 7 ($A16,A17>0$), but nodes 6 and 7 are not connected. However, because they share the same input from node 1, in Figure 3B, there are positive connections on coordinates (6,7) and (7,6). Because this connection pattern is also applied to other nodes (node 2 projects to nodes 7 and 8, node 3 projects to nodes 8 and so forth), there are false connections on the $\xb11$ to $\xb14$ diagonal lines of Figure 3B.

Shown in Figure 2B, the chain errors are due to the propagation of the covariance. Because node $A$ connects to node $B$ and node $B$ connects to node $C$, the covariance method presents covariance between nodes $A$ and $C$, which do not have a physical connection. This phenomenon is shown in Figure 3B as the extra diagonal strips.

Shown in Figure 2C, the latent confounder errors are produced due to the common input passed into two nodes. However, in this case, they are from the unobserved latent inputs. For this DCM model, it is due to the inputs from the 10 latent nodes (nodes 51 to 60), as shown in Figure 3A. Because the latent nodes have broad connections to the observed nodes, they introduce a square box shape covariance pattern into Figure 3B. And this type of false connection remains in the precision matrix and ICOV estimations (see Figures 3C and 3D) as the yellow color false connections on the bottom-left corner and the up-right corner.

### 3.2 Estimation from Differential Covariance Method

In this section, we explain the estimation results from the differential covariance method. Comparing the ground-truth connections in Figure 4A with our final estimation in Figure 4D, we see that our method essentially transformed the connections in the ground truth into a map of sources and sinks in a network. An excitatory connection, $i\u2192j$, in our estimations has negative value for $\Delta Sij$ and positive value for $\Delta Sji$, which means the current is coming out of the source $i$ and goes into the sink $j$. We note that there is another ambiguous case: an inhibitory connection $j\u2192i$, which produces the same results in our estimations. Our method cannot differentiate these two cases; instead, they indicate sources and sinks in a network.

#### 3.2.1 The Differential Covariance Method Reduces Confounder Errors

#### 3.2.2 The Partial Differential Covariance Reduces Chain Errors

#### 3.2.3 The Sparse $+$ Latent Regularization Reduces Latent Confounder Errors

Third, we use the sparse latent regularization to remove the covariance introduced by the latent variables. As already noted, when the observable nodes' connections are sparse and the number of latent inputs is small, covariance introduced by the latent inputs can be separated. As shown in Figure 4D, the external covariance in the background of Figure 4C is reduced, while the true positive connections and the directionality of the connection are maintained. Each method's ability to reduce this type of error is quantified in Figure 5C, and each method's overall performance at reducing all errors is quantified in Figure 5D.

### 3.3 Thalamocortical Model's Results

In addition to the DCM simulations, we also benchmarked various methods with a more realistic Hodgkin-Huxley-based model.

The covariance-based methods are applied to both the raw BOLD signals and the reconstructed neural signals. The differential covariance-based methods are applied to both the raw BOLD signals and the reconstructed neural signals as well. Compared to the ground truth in Figure 6A, covariance-based estimations failed to show the true connections. Also, confounder errors are noticeable on the $\xb1$1 diagonal line of Figures 6C and 6D. We show later that applying the covariance-based methods to the reconstructed neural signals does provide better results; however, covariance-based methods still suffer from the three types of false connections.

However, by applying differential covariance to the reconstructed neural signals, we see in Figure 7B that $\Delta C$ provides directionality information of the connection. However, $\Delta C$ estimation is contaminated by connections introduced from the propagation of the signals and the latent inputs. Confounder errors are greatly reduced. In Figure 7C, the partial covariance reduces chain errors and produces $\Delta P$. Finally, in Figure 7D, the sparse $+$ latent regularization is applied to reduce latent confounder errors.

In Figure 8, each estimator's performance in reducing different types of false positive connections is quantified. As explained above, while covariance-based methods perform well at chain errors and latent confounder errors, a differential covariance method performs better at reducing confounder errors. Thus, it achieves better performance overall.

### 3.4 The Backward Model and the Differential Covariance Are Both Needed to Produce Correct Estimations

We have explained our newly developed backward model and the differential covariance method. Now we show that both steps are necessary to produce good estimation.

The purpose of the backward model is to restore the dynamic form of neural signals, such that it is suitable for the differential covariance estimation. As shown in Figure 9, if we directly apply the differential covariance method to the raw BOLD signals, the estimator's performance is poor. The performance quantified using c-sensitivity was shown in Figure 8. The backward model is a better procedure (see appendix E) for reconstructing a neural signal compared to blind deconvolution (Wu et al., 2013), which is commonly used to deconvolute resting-state fMRI signals. In addition, the backward model is robust to parameter jittering because the same set of backward parameters could still reconstruct the signals simulated using forward balloon model with a wide range of parameters (see appendix F).

The differential covariance method is also needed to reduce the false connections. In Figure 10, estimations from the covariance-based methods applied to the reconstructed neural signals failed to reduce the false connections mentioned. Confounder errors that are on the $\xb11$ and $\xb12$ diagonal lines are apparent in Figures 10C and 10D. This comparison can also be seen from the quantifications on Figures 5 and 8. For both simulations, even using the same backward-processed signals, the differential covariance method performs better at reducing confounder errors than the covariance-based methods. While the covariance-based methods perform well at reducing chain errors and latent confounder errors, the differential covariance method's performance is better overall. Therefore, the differential covariance method is needed to make better estimations after the backward reconstruction step.

### 3.5 Benchmark Using Previous Study's Data Set

To show that our new method is generalizable outside of the two example simulations just given, we benchmarked various methods using a previously generated data set containing 28 BOLD signal simulations (Smith et al., 2011). Top performers mentioned in the previous study (Smith et al., 2011) (ICOV, Precision matrix, Bayes net, and Patel's $\kappa $) are included as benchmarks. Shown in Figure 11A, for the thalamocortical ionic model, we see that the $\Delta S$ outperforms others. This is mainly because of the method's better performance at reducing confounder errors.

For the 28 simulations from the previous study (Smith et al., 2011), the differential covariance method's average performance is comparable to other top methods identified previously (Figure 11B). Performance on individual data sets is provided in section A.5. Because most of these 28 simulations only have five nodes and the connection matrix is not sparse, we dropped the sparse $+$ latent regularization step in our method, except for simulation 4, the only one that has a sparse connection matrix (50 nodes, about 70 connections).

## 4 Discussion

### 4.1 Main Results

One goal of this article was to reconstruct neural signals from BOLD signals. We used synthetic BOLD data from detailed neural models to benchmark a novel differential covariance method to estimate the underlying functional connectivity. Although the relationship between fMRI and neuronal activity is a well-studied topic (for reviews, see Heeger and Ress, 2002 and Logothetis, 2002), we found that a custom-designed backward model was needed for our differential covariance method to efficiently and validly recover the underlying ground truth (Figure 9). This new method not only reduced estimation errors made by previous methods but also provided information about the directionality of the connections.

For the DCM data set from a previous study (Smith et al., 2011), differential covariances had similar performance to the previously benchmarked best methods (covariance-based methods). Each step of the differential covariance method is specialized at reducing one type of false connection. However, since the presence of these three types of false connections are minimal in this data set our new method does not outperform the covariance-based methods. Further, because the size of the network is small, the connectivity matrix is not sparse, so we cannot take advantage of the sparse $+$ latent regularization step in our method. We believe this will not be an issue when applying differential covariance method to real-world problems because the brain network has a much higher dimensionality and is sparser.

### 4.2 Caveats and Future Directions

We made several key assumptions that should be further tested and tuned in the future. Our approach depends on the neural signals reconstructed from the fMRI transfer function. Even though the fMRI transfer function we used (the balloon model) is based on the fundamental vesicular dynamics, whether high-quality neural signals can be reconstructed from the noisy real fMRI instruments should be further studied. The backward reconstruction process also had lower noise tolerance (see appendix D). In addition, the sparse and low-rank assumptions were critical to the cases we simulated. These assumptions should be validated in applications to fMRI recordings.

The next step is to evaluate the performance of differential covariance analysis on fMRI recordings and compare the results with other methods for measuring functional connectivity. Preliminary results from analyzing fMRI resting state data from humans and mice with this new method are promising and will be reported elsewhere. Estimating functional connectivity is of great importance in understanding the dynamical patterns of brain activity and differential covariance analysis of fMRI recordings offers an new way to identify causal interactions between brain region.

## Appendix A: Previous Methods

In this appendix we review several top-performing methods benchmarked in Smith et al. (2011).

### A.1 Covariance-Based Methods

The well-used covariance-based methods are used as a benchmark in this article because a prior work (Smith et al., 2011) has shown its superior performance by applying to simulated BOLD signal.

#### A.1.1 Covariance Method

#### A.1.2 Precision Matrix

We begin by considering a pair of variables $(x,y)$ and remove the correlation in them introduced from a control variable $z$.

So $pxy$ and $COVrx,ry$ differ by a ratio of $-\sigma zz|COVxyz|$.

#### A.1.3 Sparsely Regularized Precision Matrix

We used the implementation (https://www.cs.ubc.ca/∼schmidtm/Software/L1precision.html) mentioned in the previous work (Smith et al., 2011) to benchmark our new method.

### A.2 Bayes-Net Method

A variety of the Bayes net estimation methods are implemented in the Tetrad toolbox (http://www.phil.cmu.edu/tetrad/). Following the convention of the prior work (Smith et al., 2011), we tested CCD, CPC, FCI, PC, and GES. PC (“Peter and Clark”) is a causal inference algorithm used to estimate the causal relationship in a directed acyclic graph (DAG; Meek, 1995). CPC (Conservative PC) was based on the PC algorithm, but produces fewer false connections (Ramsey, Zhang, & Spirtes, 2012). GES (Greedy Equivalence Search) is another implementation based on the PC algorithm, but it uses a score system to estimate the causal connections (Chickering, 2003; Ramsey et al., 2010). The FCI (Fast Causal Inference) algorithm is a different approach, which takes into account the presence of latent confounder and selection bias (Zhang, 2008). CCD (Cyclic Causal Discovery), unlike the other algorithms mentioned, allows cycles in the network for estimation (Richardson & Spirtes, 1996). The code used in this article is based on the Tetrad Toolbox and further adapted by Smith et al. (2011) for fMRI data.

### A.3 Patel's $\kappa $ Method

Patel's $\kappa $ measurement assesses the relationship between pairs of distinct brain regions by comparing expected joint and marginal probabilities of elevated activity of voxel pairs through a Bayesian paradigm (Patel, Bowman, & Rilling, 2006). For this article we used the implementation generously provided by Smith et al. (2011).

### A.4 Sparse Latent Regularization

Prior studies (Banerjee et al., 2006; Friedman et al., 2008) have shown that regularizations can provide better estimation if the ground-truth connection matrix has a known structure (e.g., sparse). For all data tested in this article, 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.

### A.5 Methods' Performance for Individual Data Sets

In Figure 12, the performance of the different methods mentioned above is quantified using c-sensitivity on individual data sets simulated in Smith et al. (2011). For comparison, the methods' performance on thalamocortical model and their average performance across all 28 data sets are shown in the first two columns.

## Appendix B: Balloon Model Parameters

## Appendix C: Backward Model Parameters

Parameter . | Description . | Value . |
---|---|---|

$\kappa $ | Rate of signal decay | 0.65 per s |

$\gamma $ | Rate of flow-dependent elimination | 0.41 per s |

$\tau $ | Hemodynamic transit time | 0.98 per s |

$\alpha $ | Grubbs exponent | 0.32 |

$Eo$ | Resting oxygen extraction fraction | 0.5 |

$\epsilon $ | Ratio of intra- and extravascular signal | 0.5 |

$\u03d1o$ | Frequency offset at the outer surface of the magnetized | 40.3 per s |

vessel for fully deoxygenated blood | ||

TE | Time of echo | 40 ms |

Parameter . | Description . | Value . |
---|---|---|

$\kappa $ | Rate of signal decay | 0.65 per s |

$\gamma $ | Rate of flow-dependent elimination | 0.41 per s |

$\tau $ | Hemodynamic transit time | 0.98 per s |

$\alpha $ | Grubbs exponent | 0.32 |

$Eo$ | Resting oxygen extraction fraction | 0.5 |

$\epsilon $ | Ratio of intra- and extravascular signal | 0.5 |

$\u03d1o$ | Frequency offset at the outer surface of the magnetized | 40.3 per s |

vessel for fully deoxygenated blood | ||

TE | Time of echo | 40 ms |

## Appendix D: Noise Tolerance

## Appendix E: Benchmark with Blind Deconvolution

Blind deconvolution aims to reduce the confounding effects of the hemodynamic response function (HRF) on resting-state fMRI signals. For task-related fMRI, neural population dynamics can be captured by modeling signal dynamics with explicit exogenous inputs, which is absent during resting-state fMRI recordings. Blind deconvolution considers resting-state fMRI as spontaneous event related and then estimates region-specific HRF and performs deconvolution on the signal. For details, see Wu et al. (2013).

## Appendix F: Robustness of the Backward Model

## Acknowledgments

This research was supported by the Office of Naval Research (N00014-16-1-2829) and NIH/NIBIB (R01EB026899).

## References

## Author notes

T. Lin and Y. Chen contributed equally to this work.