Knowing brain connectivity is of great importance both in basic research and for clinical applications. We are proposing a method to infer directed connectivity from zero-lag covariances of neuronal activity recorded at multiple sites. This allows us to identify causal relations that are reflected in neuronal population activity. To derive our strategy, we assume a generic linear model of interacting continuous variables, the components of which represent the activity of local neuronal populations. The suggested method for inferring connectivity from recorded signals exploits the fact that the covariance matrix derived from the observed activity contains information about the existence, the direction and the sign of connections. Assuming a sparsely coupled network, we disambiguate the underlying causal structure via L1-minimization, which is known to prefer sparse solutions. In general, this method is suited to infer effective connectivity from resting state data of various types. We show that our method is applicable over a broad range of structural parameters regarding network size and connection probability of the network. We also explored parameters affecting its activity dynamics, like the eigenvalue spectrum. Also, based on the simulation of suitable Ornstein-Uhlenbeck processes to model BOLD dynamics, we show that with our method it is possible to estimate directed connectivity from zero-lag covariances derived from such signals. In this study, we consider measurement noise and unobserved nodes as additional confounding factors. Furthermore, we investigate the amount of data required for a reliable estimate. Additionally, we apply the proposed method on full-brain resting-state fast fMRI datasets. The resulting network exhibits a tendency for close-by areas being connected as well as inter-hemispheric connections between corresponding areas. In addition, we found that a surprisingly large fraction of more than one third of all identified connections were of inhibitory nature.
The networks of the brain are key to understanding its function and dysfunction [1]. Depending on the methods employed to assess structure and to record activity, networks may be defined at different levels of resolution. Their nodes may be individual neurons, linked by chemical or electrical synapses. Alternatively, nodes may also be conceived as populations of neurons, with links represented by the net effect of all synaptic connections that exist between two populations. In any case, this defines the structural substrate of brain connectivity, representing the physical (causal) basis of neuronal interactions. Nodes in a brain network influence each other by sending signals. For example, the activities of nodes in a network are generally not independent, and neuronal dynamics are characterized by correlations among the nodes involved in the network. This suggests an alternative perspective on active brain networks: Functional connectivity assigns a link to a pair of nodes to the degree to which their activities are correlated. It has been argued that this concept emphasizes connections that “matter”, including the possibility that the same substrate may give rise to different networks, depending on how they are used. As a consequence, functional connectivity and structural connectivity are not equivalent. A well-known phenomenon is that two nodes may be correlated, even if there is no direct anatomical link between them. For example, a shared source of input to both nodes may generate such a correlation, which does not correspond to a direct interaction between the two nodes. Apart from that, correlation is a symmetric relation between two nodes, whereas a physical connection implies a cause-effect relation that is directed. There have, in fact, been multiple attempts to overcome the shortcomings of functional connectivity, especially the lack of directed interaction. The term effective connectivity has been suggested for this [2]. The idea is to bring the networks, inferred from activity measurements, closer to structural connectivity, which can only be inferred with anatomical methods. The dichotomy between structural and functional aspects of connectivity raises the general question whether it is possible to infer brain networks from recorded activity. We are only beginning to understand the forward link between structural connectivity and functional connectivity. As a consequence, it is possible to compute correlations from connectivity in certain simplified network scenarios [3]. The correspondence between connectivity and correlation, however, is not one-to-one. Networks with different connectivity may lead to exactly the same correlations between nodes. As a consequence, the inverse problem of inferring connectivity from correlation is generally ill-defined. As we will demonstrate in this paper, additional assumptions about the connectivity can help to resolve the ambiguity. Specifically, we search for the network with the lowest number of nonzero edges (via L1-minimization) to disambiguate the problem. Structural, functional and effective connectivity are not equally well accessible. Some aspects of the anatomical structure can be assessed post mortem by invasive tracing methods, or non-invasively by Diffusion Tensor Imaging, DTI. In contrast, functional connectivity is based on statistical relationships between the activity of neuronal populations and can be easily estimated from recorded signals. For estimating effective connectivity there are methods like Dynamic Causal Modelling, DCM [4,5], Granger causality [6] and others [7–13]. Only few methods to infer effective connectivity, however, can deal with large numbers of nodes (40 or more) based on zero-lag correlation only. However, they are either limited to small networks [14], or to directed acyclic graphs [15]. Here, we are proposing a new method for the estimation of effective connectivity from population activity in the brain, especially BOLD-related signals. The new method is a variant of the procedure described in [16], based on a L1-minimization. For the method proposed here it is sufficient to use zero-lag covariances to estimate directed effective connectivity.
The main idea of our estimation method is inspired by the finding, “that the key to determining the direction of the causal relationship between X and Y lies in ‘the presence of a third variable Z that correlates with Y but not with X, ’ as in the collider X → Y ← Z…” [17,18]. Similarly, assuming a linear interaction model, the presence of a collider structure in a network (see Fig 1) produces specific entries in the corresponding inverse covariance (precision) matrix. Fig 1 shows a disconnected network in the left column, and a network which induces the same covariance matrix if all links have opposite direction in the middle column. In the latter case an estimation of the direction is impossible, because there is simply no information about it in the covariance matrix. Whenever a collider structure is present, however, the entry in the inverse covariance matrix for the two source nodes (here, 2 and 3) is non-zero. This is due to the fact that in a linear model the entry in the inverse covariance matrix depends not only on the connections of the nodes 2 and 3, but also whether these nodes have a common target. This means the presence of a collider structure allows us to disambiguate the direction of this particular connection. We consider here a scenario, where the interaction between nodes is described by a generic linear model. Assuming stationarity, let the neural activity x (t) be implicitly defined by the consistency equation x (t) = (G * x) (t) + v (t) (1) where G (t) is a matrix of causal interaction kernels and v (t) denotes fluctuating external inputs (“driving noise”). All variables are also listed in Table 1. Fourier transformation of Eq (1) yields x ^ (f) = G ^ (f) x ^ (f) + v ^ (f) and simple rearrangement leads to x ^ (f) = (1 - G ^ (f) ) - 1 v ^ (f) where x ^ denotes the Fourier transform of x. The cross spectral density of the signals is then given by C ^ (f) = (1 - G ^ (f) ) - 1 Z ^ (f) (1 - G ^ * (f) ) - 1 where Z ^ (f) is the cross-spectral density of the external inputs. It follows C ^ - 1 (f) = (1 - G ^ * (f) ) Z ^ - 1 (f) (1 - G ^ (f) ) = B * (f) B (f) (2) with B (f) = Z ^ - 1 (f) (1 - G ^ (f) ). In our model, we assume that the components of the external fluctuating input are pairwise stochastically independent for all nodes. Then, Z ^ is a diagonal matrix, and we make the additional assumption that Z ^ = 1. For the linear model considered here, there is a relation between covariance and connectivity, which can be exploited for the estimation of connectivity from correlation. In the case Z ^ = 1 it is given by C ^ - 1 (f) = (1 - G ^ * (f) ) (1 - G ^ (f) ) = 1 - G ^ (f) - G ^ * (f) + G ^ * (f) G ^ (f) where the last term contributes the information of the collider structures. If the matrix product G ^ * (f) G ^ (f) has a non-zero off-diagonal entry the corresponding nodes have outgoing connections terminating at the same node, which means these nodes form a collider. It is clear that for any unitary matrix U ∈ U (n) the product UB is still a solution of Eq (2), as U * U = 1. We will resolve this ambiguity with an L1 minimization which is known to prefer sparse solutions under certain conditions [19]. In order to find G from a given C we first fix an initial matrix B, and then search for a unitary matrix U ∈ U (n) such that ‖UB‖1 is minimal, so we are minimizing the function Γ: U (n) ⟶ R U ⟼ ∥ U B ∥ 1. (3) When estimating connectivity from simulations with known underlying network structure (ground truth), one can quantify the performance of the estimation. For measuring the accuracy of our estimation we employ three different methods. First, we use the area under the ROC-curve (AUC). The ROC (receiver operating characteristic) curve is obtained as following: For each possible parameter value (in our case the threshold for the existence of a connection), the number of true-positives (TP) and false-positives (FP) is used to calculate the true-positive rate (or recall) TP/ (TP+ FN) and the false-positive rate FP/ (FP+ TN). The ROC curve is then obtained by plotting the true-positive-rate against the false-positive rate. Secondly, we use the average precision score (PRS) which is the area under the precision-recall curve. This also includes the false-negatives (FN) (precision: TP TP + FP). If both AUC and PRS are equal to 1, the connections in the network are perfectly estimated. Sample curves are shown in Fig 2D. Thirdly, we calculate the Pearson Correlation Coefficient (PCC) which in contrast to the measures defined before also take the strength and the sign of the interactions into account. This also means that this measure is less suited to assess whether a connection exists or not. It rather measures whether the estimated connections have the same strength as the original ones. We consider all three performance measures simultaneously to establish the quality of our estimates. Seven healthy subjects underwent a 20-minute resting-state fMRI experiment on a 3 T Siemens Prisma scanner. The data was acquired using the MREG sequence [28], yielding a high temporal resolution (TR = 0. 1 s, 12000 time points) that facilitates functional connectivity analyses [29]. The other sequence parameters were TE = 36 ms, FA = 25°, 64 × 64 × 50 matrix and 3mm isotropic voxel size. Additionally, cardiac and respiratory signals were recorded with the ECG and abdominal breathing band from the scanner’s physiological monitoring unit. Motion correction was done with FSL and physiological noise correction was performed with RETROICOR [30]. Average CSF and white matter signals were regressed out, but no global signal regression was performed. Following image normalization to MNI space, voxels were parcellated according to the AAL atlas (excluding the cerebellum), and the mean activity within each atlas region was calculated. The connectivity was then estimated using zero-lag covariances of the standardized signals.
Intrinsic properties of our new estimation procedure can be identified by studying the performance of the method for perfectly estimated (noise-free) covariance matrices. This way we address properties that do not depend on any particular feature of the underlying data, and that are not due to the success of the measurement process. In particular, we show for which types of networks our estimation procedure gives good results on technical grounds, with a wide range of networks hopefully including those arising in applications. We used random Erdős-Rényi connectivity profiles for all simulations. The macro-connectivity between neuronal populations has to satisfy certain conditions in order to be tractable by our methods. Two of these conditions concern the dynamic stability of the network and the strength of the interactions. There is a trade-off between the number of physical links and the resulting strength of macro-connections, and the dynamic stability of the network. To study the performance of our method in these various regimes, we separately varied the network size N, the connection probability p, and the absolute strength of connections |J| in the connectivity matrix G, while the fraction of inhibitory couplings was kept at 50%. The spectral radius ρ of the bulk eigenvalue spectrum is approximately given by ρ 2 = J 2 p (1 - p) N. (8) The default values of the parameters used in our study were N = 100, p = 0. 1 and ρ = 0. 7, where only one of them at a time was systematically varied. Low values of the spectral radius ρ correspond to networks with weak recurrent interaction and high values to networks with strong interaction, respectively. According to the model of network interaction assumed here, the networks need to have a spectral radius ρ > 0 for network interaction to be present and ρ < 1 for the dynamics to be stable. First, our results in Fig 3A indicate that a certain minimal level of interaction is necessary to be able to estimate the connections reliably. Above a value of ρmin = 0. 2, the influence of the spectral radius on the performance of the estimation is weak, but the larger the spectral radius is the better the estimation gets. Secondly, the connection probability of the network influences the quality of the estimation. For all connection probabilities tested here the network size was kept constant at N = 100 nodes. The networks were constructed such that the strength |J| of all connections was the same and such that the spectral radius ρ was constant according to Eq (8). Fig 3B shows that the estimation works very well for sparse matrices with a connection probability in the range between 5% and 15%. For networks with higher connection probability and equally strong connections, the performance decreases as expected, due to the bias associated with L1-minimization. But even for a connection probability of p = 0. 21, a fraction of 14. 2% of the estimated connections are false negative, and 3. 3% are false positive. More than 90% of the correctly estimated connections have the correct sign. In applications, the focus of the estimation often lies on the strongest connections in the network. In networks with a background of weak connections and a sparse skeleton of stronger connections, it is possible to selectively estimate these strong links although, strictly, the assumption of a sparse network is violated. Fig 3D shows the performance of our method for such networks: the networks consist of a skeleton of strong connections with connection probability p = 10% and a connection strength derived from Eq (8) for ρ = 0. 7. Additionally, we created a second network with weaker connections for various connection probabilities q. The two networks were combined by adding the connectivity matrices. The connection strength of this weaker connections is also derived from Eq (8), with a spectral radius of the background network being 20% of the spectral radius of the skeleton network. Then the performance of the estimation is calculated with respect to the skeleton of strong connections. Thirdly, to be applicable to a broad range of data types, a method of connectivity estimation should perform stable for different network sizes N. For most common types of non-invasive recordings of population activity the number of nodes considered is in the range between 30 and 150. It is, of course, possible to consider larger networks, although the estimation becomes computationally more expensive. The runtime of the algorithm for networks with 200 nodes still in the range of seconds on a state-of-the-art desktop computer, but even networks with 1000 nodes or more are tractable. The strength of the connections |J| are set such that the spectral radius ρ of G is constant; the connection probability is constant at p = 0. 1. Fig 3C shows that our method performs better for bigger networks. We have observed that the L1 cost landscape becomes smoother for larger networks. In order to create surrogate data which fit fast fMRI data [28], we simulated interacting stochastic processes known as Ornstein-Uhlenbeck processes. In this case, the performance of the network inference depends on how well the inverse covariance matrix, which is the basis of the estimation, can be derived from the data. In addition to finite size effects, we studied the impact of observation noise on the performance, see Fig 4. We used N = 100, p = 0. 1, dt = 0. 1 s, ρ = 0. 74 and τ = 0. 1 s as default values of the parameters. Generally, it seems natural to use Welch’s method to calculate cross-spectral densities directly, and then to estimate the connectivity for each frequency band separately. For the data described here, however, we can estimate the connectivity from zero-lag sample covariances in the time domain. This is possible when the mass of the covariance function is concentrated very close around lag 0. Then lag 0 is the only one contributing to the integral of the covariance function, which corresponds to the cross-spectral density C ^ (0). As shown in Fig 4A, with noisy data the AUC is still good, but the PRS is lower than in the case, where the covariance is known without error. However, for a signal-to-noise ratio above 1 the performance improves very quickly. In the case of fMRI usually the whole brain is scanned, and there are no unobserved nodes in the network. However, for other data types (e. g. fNIRS) only parts of the brain can be observed. The question then is, whether this sub-network can nevertheless be reconstructed from the recorded signals. To model this scenario, we took simulated data and removed randomly a certain subset of components from the dataset. The interaction of the removed nodes is then not part of the covariance matrix of the reduced dataset, although the unobserved nodes of course still exert their influence on the observed ones. The performance of the estimation of the sub-network based on the reduced dataset is shown in Fig 4B. Our analysis shows very clearly that the estimation still leads to reasonable results under these conditions. In fact, we can demonstrate that we are inferring causal connections only: For unconnected observed nodes X, Y and a latent node L connected to both X and Y, our method does not erroneously indicate a link between X and Y. One key factor for a reliable estimation of the covariance matrix is the amount of data available. This depends on the length of the measurement or simulation, and on the sampling rate. Since fast fMRI time series are obtained by measuring the BOLD response as a proxy of neuronal activity, the time scale of the measured data is relatively slow compared to the time scale of the underlying neuronal activity. Fig 5 shows the performance of network inference depending on the amount of data available, and on the time-scale of the neuronal activity. Not surprisingly, the more voluminous the dataset is, the better the estimation gets. On the other hand, it shows that the estimation generally leads to better results for slower temporal dynamics. Also, for data of sufficient length with a fairly good signal-to-noise ratio, the estimation of the connectivity is possible even when only a part of the network is observed. To allow comparison of our new method with other known methods for network inference [31–33], we applied it to the NetSim dataset provided by [34]. For details on the result of this, please see Fig A in S1 Text. We estimated connectivity from seven fast fMRI datasets, for details see the methods section. The resulting networks, after a threshold of 10% was applied, consist of 810 connections for each dataset. The threshold of 10% was chosen arbitrarily. In Fig B in S1 Text we show the histogram of estimated connection strengths for all seven reconstructed networks before thresholding. The threshold is derived from the 10% strongest connections, disregarding their signs. As there is generally no full ground truth for the connectivity inferred from human fMRI recordings available [31,32,35], we cannot definitely assess the degree to which the result of our inference are correct. We can, however, establish whether they are plausible. One representative connectivity matrix is shown in Fig 6. On average, 34% of the connections were inhibitory, with negligible variability across subjects. Of all connections found, 301 (37%) were found in four subjects or more, and 4872 out of 8100 possible connections were absent in all subjects. On average, 245 of the connections were bi-directional and 565 connections were identified only for one direction. In general, close-by areas are more likely to be connected than more distant ones. This fact is (approximately) represented by a concentration of connections along secondary diagonals in the within-hemisphere blocks. Also, there are frequent inter-hemispheric connections between corresponding areas. This fact is represented by the diagonal entries in the across-hemisphere blocks. As mentioned above, different heuristics have been suggested to reconstruct networks from neuronal signals. In Fig 7 we compare the performance of the new method we propose here and the established method of Regularized Inverse Covariance [34], based on the implementation provided at https: //fsl. fmrib. ox. ac. uk/fsl/fslwiki/FSLNets. Our comparison clearly shows that our new method performs significantly better than the Regularized Inverse Covariance method, mainly, because the latter cannot establish the direction of connections. The superior performance of the new method is reflected in higher values for all three performance measures, in particular PRS and PCC. As regularization parameter required by the software toolbox, we used λ = 5. Furthermore, we applied the RIC method on all seven MREG datasets described before. A threshold was applied, such that only the 10% strongest connections are retained. To compare the outcome of both methods, we only condidered the existence of connections (binary and symmetric connectivity) and disregarded weights and directions (weighted nonsymmetric connectivity). One representative example of the comparison of both methods is shown in Fig 8. For RIC, 376 out of 810 possible connections where identified in four subjects or more out of seven, the corresponding number for our method is 392 out of 810 possible connections. If any method produced directed networks with 10% connection probability at random, this would yield an average count of less than 25 connections (3% of 810 connections) that agree for least four out of seven independently generated networks. On average, 290. 5 out of 4050 possible connections (undirected) are identified by both methods, 3530. 5 connections were found by neither of the methods. This means that both methods agree on 3821 out of 4050 connections on average. The two methods disagreed on the remaining 229 connections.
With the presented method we can estimate directed effective connectivity on a whole-brain scale. Also we are able to detect whether connections are excitatory or inhibitory. The estimation is possible based on zero-lag covariances, but can also be applied to frequency-resolved cross spectral densities.
Changes in brain connectivity are considered an important biomarker for certain brain diseases. This directly raises the question of accessibility of connectivity from measured brain signals. Here we show how directed effective connectivity can be inferred from continuous brain signals, like fMRI. The main idea is to extract the connectivity from the inverse zero-lag covariance matrix of the measured signals. This is done using L1-minimization via gradient descent algorithm on the manifold of unitary matrices. This ensures that the resulting network always fits the same covariance structure as the measured data, assuming a canonical linear model. Applying the estimation method on noise-free covariance matrices shows that the method works nicely on sparsely coupled networks with more than 40 nodes, provided network interaction is strong enough. Applying the estimation on simulated Ornstein-Uhlenbeck processes supposed to model BOLD signals demonstrates robustness against observation noise and unobserved nodes. In general, the proposed method can be applied to time-resolved covariance matrices in the frequency domain (cross-spectral densities), leading to frequency-resolved networks. We are able to demonstrate that our method leads to reliable results, if the sampled signals are long enough.
medicine and health sciences
diagnostic radiology
functional magnetic resonance imaging
neural networks
particle physics
random variables
neuroscience
covariance
colliders
magnetic resonance imaging
simulation and modeling
mathematics
algebra
network analysis
brain mapping
neuroimaging
research and analysis methods
computer and information sciences
operator theory
imaging techniques
probability theory
physics
radiology and imaging
diagnostic medicine
linear algebra
biology and life sciences
physical sciences
eigenvalues
From correlation to causation: Estimating effective connectivity from zero-lag covariances of brain signals