Model of interconnected cortical areas

The node activity is governed by the following coupled ordinary differential stochastic equations: (3) Each experiences an exponential decay with time constant τ_x and is excited by other whose contributions are scaled by the recurrent weights C_ij. Fluctuations are generated by the Gaussian white noise that has variance (σ_i)²; formally, is a Wiener process scaled by a factor σ_i. We forbid self-connections, i.e., C_ii = 0, although this is not crucial. In practice and without loss of generality, the input e_i is homogeneous for all nodes.

Following previous work [16, 27, 33], we derive the well-known self-consistency equations for the mean and covariances Q^τ in Eq (1), assuming stationarity of the process. When the system has a stable fixed point, it is given by the zeros of the linear matrix system: (4) where x^t indicates the activity vector whose entries are for 1 ≤ i ≤ N; e is the homogeneous vector of inputs e_i. We have denoted the Jacobian of Eq (3) by (5) where the Kronecker delta δ_ij = 1 if i = j and 0 otherwise. The system has a single fixed point provided J is invertible.

Using Ito’s formula for the derivation of with respect to t, we obtain the Lyapunov equation for the steady-state of the second-order fluctuations with zero time shift: (6) The noise matrix Σ is diagonal with terms and † is the matrix transpose operator.

As done for the variable t, the derivation of with respect to τ yields an equation satisfied by the time-shifted covariance Q^τ: (7) where expm is the matrix exponential. When J is dominated by the diagonal elements, the autocovariances are close to an exponential decay with time constant τ_x, as illustrated in Fig 1E with a log y-axis by the straight black line that coincices with the mean in blue up to τ = 2.5 s.

Evaluating covariances from known network parameters

For a model with given connectivity matrix C and intrinsic noise matrix Σ, the theoretical zero-time-shift covariance matrix Q⁰ is the unique solution of Eq (6) and can be evaluated using the long-known Bartels-Stewart algorithm [16]. Then, the model Q^τ can be calculated from Q⁰ and J using Eq (7) for any given τ.

Direct estimation of network parameters from empirical covariances

The inverse problem consists in finding an estimate of the connectivity matrix C such that the model reproduces the two objective matrices and for a given τ > 0, e.g., given by empirical observations. The relation in Eq (7) allows for the calculation of J from and an occurrence of for a given τ > 0: (8) where logm is the matrix logarithm [34]. In our case, the non-diagonal elements of the Jacobian directly give the connectivity weight C_ij = J_ij. Note that the diagonal terms of J can provide an estimation of τ_x, so far as we assume C_ii = 0. Then, the noise matrix can be evaluated from the Lyapunov equation Eq (6): (9) Due to the matrix logarithm, this method is very sensitive to the noise in and : it can give non-real numbers for J, hence C. For this reason, we develop the optimization method in the next section.

Lyapunov optimization (LO) to estimate network connectivity from empirical covariances

As will be shown later in Results, the direct method does not work well with noisy empirical . We thus propose an alternative estimation method, where C is tuned iteratively such that the model reproduces the covariance observables, as illustrated in Fig 1H. Considering the noise matrix Σ to be known for now, we optimize C in order to reduce the model error (13) for a pair of objective covariance matrices and with a given τ > 0; in Results, it is referred to as τ_est. Being the sum of two matrix distances, the Lyapunov function V is positive definite and becomes zero only when the two covariance matrices are equal to the objective counterparts. Eq (8) ensures the unicity of the connectivity C for a given pair of covariances matrices, when the solution exists.

Starting from zero weights initially, each optimization step aims to reduce V. To do so, we calculate the Jacobian J in Eq (5) for the current connectivity C, then the model matrices Q⁰ and Q^τ using Eqs (6) and (7) as we know Σ. We want to update C to obtain the following desired changes for the model covariances Q⁰ and Q^τ: (14) where the optimization rate ϵ_C is a some small positive number. To evaluate the Jacobian update ΔJ that corresponds to ΔQ⁰ and ΔQ^τ, we consider the equivalent of Eq (8) for the theoretical covariance matrices, namely J = logm(X)^†/τ with X = (Q⁰)⁻¹ Q^τ. We perform an implicit differentiation with respect to J and X, yielding (15) Finally, the desired connectivity update is simply for all i ≠ j (16) and zero for the diagonal elements.

To obtain the rhs of the first line in Eq (15), we have assumed that matrices X and ΔX commute. For the LO procedure, this implies the existence of a path of matrix increments corresponding to ΔJ such that the commutation requirement is satisfied at each step, in order to reach the global optimum. Assuming X to be arbitrary, the subspace of commuting matrices in which ΔX is constrained is a linear subspace of dimension n(n−1)/2 with n the number of nodes, to be compared to the dimension n² of the space of X. Although this subspace is not dense, it is sufficiently rich to hope that suitable paths exist and the optimum can be approached; this will be verified using numerical simulation.

This optimization process may also not reach the global optimum when and do not correspond to a real solution for C, such as noisy empirical covariances. We will thus verify the estimation performance numerically, using the known theoretical covariance matrices obtained from Eqs (6) and (7) to check for a given connectivity matrix C.

Tuning of input noise variances

So far, we have considered Σ to be known and fixed during the optimization. Now we extend the theory to optimize Σ in parallel with C. We focus on the case where Σ is diagonal with elements Σ_ii = (σ_i)² that are the variances of the noise experienced by each network node. The key to adjust each Σ_ii is to reproduce the corresponding diagonal element . To do so, the Σ_ii are simply modified at each optimization step based on the difference between the variance of the i^th node of the model and its objective counterpart: (20) where the update rate is ϵ_Σ.

The implications of tuning the noise matrix Σ will be explored in a dedicated section of Results. Intuitively, it can be seen from Eq (6) that Σ affects the model covariance Q⁰, hence Q^τ via Eq (7); see Methods for reference. Therefore, assuming an erroneous Σ for the model implies incorrect desired updates ΔQ⁰ and ΔQ^τ, leading to a wrong update ΔC.

Heuristic optimization procedure for network connectivity

We also compare the LO method with a heuristic optimization, where existing weights are increased or decreased depending on the Q^τ difference between the model and objective for the corresponding matrix element for a given τ (21) This method used previously [12] is adapted here to spatio-temporal covariances.

Network configuration and model parameters used in numerical simulation

Unless stated otherwise, we use the parameters in Table 2 and Euler’s approach to simulate Eq (3) numerically. The code is written in python using the packages numpy and scipy, which incorporate a version of the Bartels-Stewart algorithm. The topology for the cluster-hub networks as in Fig 1A consists of two groups of 30% and 60% of the N nodes. Each group is connected recurrently with probability p_con, with a weight randomly chosen in [10, 100]% of c_max. We exclude self connections, although this is not crucial. The two groups are connected via hubs with a higher probability equal to 1.3p_con and a weight in the same range as other connections. Hubs do not have connections between each other. For random networks, all pairs of nodes have the probability p_con to be connected, with a random weight in [10, 100]% of c_max. Simulations of duration T with identical parameters but distinct initial conditions are repeated n_sim times in order to evaluate the empirical Q.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Table of simulation, network and optimization parameters for the ND artificial network model.

https://doi.org/10.1371/journal.pcbi.1004762.t002

We also compare our ND model with the DMF model [27] that simulates a network of neural population with AMPA and NMDA conductances, whose synaptic activity determines the BOLD signal via a hemodynamic nonlinear filter [29]. The equations regulating the synaptic activity variable S_i with 1 ≤ i ≤ N for the DMF are (22) The HRF model relies on the synaptic activity S_i to calculate the intermediate variables s_i, f_i, ν_i and q_i for each area, which give the BOLD signal B_i: (23) All parameters are recapitulated in Table 3.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Parameters for the dynamic mean-field (DMF) model and hemodynamic response function (HRF).

https://doi.org/10.1371/journal.pcbi.1004762.t003

Experimental data set

Resting-state BOLD signal time series and dwMRI data were acquired for 25 healthy subjects aged between 18 and 39 years (mean 26.7), including 12 females and 13 males. A detailed description of the data acquisition can be found in [35], which we briefly summarize here. Subjects were asked to stay awake and keep their eyes closed, during which joint MRI-EEG recordings were performed using a 3 Tesla Siemens Trim Trio MR scanner and a 12-channel Siemens head coil [36, 37]. Acquired time series are spatially aggregated for voxels corresponding to the same area according to the parcellation. Only fMRI data are used in the present study: the BOLD signals were recorded for about 20 minutes (661 time points taken every 2 s). The BOLD signals are high-passed filtered above 0.01 Hz. Then the FC matrices correspond to covariances between all 68 regions, which are calculated using Eq (2) taking each of the 25 subjects as a separate session. Time shifts are multiples of the temporal resolution after preprocessing (2 s). For the most part, we use for FC the average covariance matrix over all subjects; a comparison with optimization based on FC for individual subjects is also presented in ‘Robustness of LO applied to experimental data. For each participant, dwMRI was used to evaluate the white-matter intracortical connectivity. Tractography is performed using corrections for motion, eddy currents, crossing pathways and the size of regions. Processing steps for MRI data include 1) preprocessing of T1-weighted scans, cortical reconstruction, tessellation and parcellation, 2) transformation of anatomical masks to diffusion space, 3) processing of diffusion data, 4) transformation of anatomical masks to fMRI space, 5) Processing of fMRI data. The dwMRI matrix used in the application to experimental data is the average over all subjects, which is thresholded to determine weights to tune in the optimization.

Ethics statement

The paper makes a secondary use of experimental data already published [35]. The experiments have been approved by the Ethic Committee of Charite University (Berlin).

Noise-diffusion (ND) with non-trivial network topology

The network model consists of interconnected neural populations that experience intrinsic noise and excite each other. Details about the mathematical formalism and model parameters are provided in Methods. For artificial networks, we mainly consider the topology in Fig 1A, where two groups of interconnected populations are linked via hubs. An example of the connectivity matrix is displayed in Fig 1B, where the matrix element indexed by ij corresponds to the connection from area j to i. Although hubs (cyan) are connected to both groups and not between each other (as indicated red arrow), they seem to make one with the large group (light green) in Q⁰ and Q^τ (see the blue arrow in Fig 1D). In contrast, the small group (dark green in Fig 1A and 1B) does not appear clearly as it exhibits rather low covariances. The reason behind this choice is to check whether the method can recover the correct topology from the Q^0/τ values that significantly deviate from the underlying connectivity strengths in C.

Our estimation method uses the spatio-temporal information in the covariance of the entire network to recover the directed weights in C, implicitly relying on the asymmetric matrix for τ > 0. In general, Q^τ is not symmetric for τ > 0 for an arbitrary C: Fig 1F shows that the asymmetry of Q^τ linearly scales on average with that of C, with a slope that increases with τ. We measure the asymmetry using the following index (24) for a matrix M. It ranges from 0 for symmetric matrices to 1 for antisymmetric matrices. Large weights with no or weak reciprocal connection mainly contribute to the index. The relationship in Fig 1F can be understood looking at the expression for Q^τ in Eq (7), which amounts to Q⁰ multiplied by the matrix exponential expm(J^† τ). The Jacobian J^† is the transpose of the Jacobian matrix defined in Eq (5) and has the same off-diagonal elements as C, so asym_J = asym_C. For small τ, Q^τ can be approximated by Q⁰ + Q⁰ J^† τ, which explains why asym_Q^τ is smaller than asym_C. The variability of the Q^τ asymmetry compared to the small slopes of the curves makes it difficult to infer the asymmetry of C. This further motivates our method to estimate accurately C from Q.

We focus on the case where the recurrent weights are sufficiently large to generate a strong network effect: the mapping between C and Q is not trivially linear and many large values in Q⁰ correspond to nodes that are only weakly or not connected. Here, hubs in the network in Fig 1A are not connected, but have strong variances in Fig 1G and covariances between each other due to their numerous connections.

Test of Lyapunov optimization (LO) on artificial networks

The optimization procedure is summarized in Fig 1H: C is iteratively modified such that the model covariance matrices Q⁰ and Q^τ converge towards the objective covariance matrices and , for a chosen τ = τ_est > 0. In the first place we assume Σ to be known. Starting from a initial matrix C with all weights equal to zero, LO calculate a C update at each step such that it reduces the Lyapunov function V in Eq (13); V is the sum of the matrix distance between Q⁰ and on the one hand, and that between Q^τ and on the other hand. To measure the goodness of fit of C or Q, we define the model ‘error’ as the normalized distance for a matrix M and its objective M^obj: (25) which are involved in V.

Firstly we verify that LO recovers the correct connectivity for theoretical objectives calculated with Eqs (6) and (7) with the original C and Σ. A typical example for the evolution of the C and Q errors is illustrated in Fig 2A for τ_est = τ_x = 1 s. While the Q^0/τ errors exhibit a plateau, the C error keeps on decreasing. Fig 2B shows the faster convergence for Q than C for three stages of the optimization. The residual C error is very small; it may remain non-zero when the commutativity assumption necessary to derive Eq (15) is not satisfied.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Iterative optimization to recover C from empirical Q.

A Estimation based on theoretical covariances. Evolution of Q⁰ and Q^τ errors between the model and theoretical objectives (blue curves), as well as C error between the original and model connectivity (red). The theoretical objectives are calculated using Eqs (6) and (7). Errors are calculated using the normalized matrix distance in Eq (25). The optimization relies on Eq (15) with time shift τ_est = 1 s. B Details of the match of the model and objective matrices for Q^τ (blue dots) and C (red crosses) at three stages of the optimization in A. The black line indicates a perfect match. C Similar plot to A with empirical Q as objectives, evaluated using Eq (2) for 50 simulations of the same network as in A. D Similar plot to B for the optimization in C.

https://doi.org/10.1371/journal.pcbi.1004762.g002

For the same original C, we now simulate the network and calculate the empirical objectives using Eq (2) for 50 independent repetitions of 300 s each with timestep 0.05 s (i.e., 3 × 10⁵ data points). In Fig 2C, the Q^0/τ errors first drop, then rebound and stabilize. Interestingly, the C error continues to decrease even though the Q error worsens. Because of the noise in the empirical , there is no real solution C to Eq (7) and the residual error does not vanish. In Fig 2D, the model C gives a better fit for the original C after 5000 optimization steps than for the minimum for the Q error. LO also allows for the use of constraints on the weights in C. Non-negativity is imposed here, although it does not affect the results significantly.

Matching time shift is crucial for accurate connectivity estimation.

The estimation performance depends on the choice for τ_est related to the empirical objectives: Fig 3A shows for several τ_est the ‘best’ C, corresponding to the minimum of the Q error as in Fig 2C during the optimization. From now on, the Q error denotes the mean normalized distance for both Q matrices: . For τ_est = 0.1 and 1 s, the global characteristics of the original C topology (left panel) are recovered, namely the two interconnected groups linked via hubs that are not connected with each other. Although weights are weaker than in the original C, the detail of the connectivity is recovered for τ_est = 1 s. In contrast, the estimated C for τ_est = 0.1 s is quasi symmetric. Last, τ_est = 5 s gives a poor estimation where the cluster-hub topology is hardly detected.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Influence of choice for τ_est on accuracy of recovered C.

A Best estimated matrices C for τ_est = 0, 0.1, 1 and 5 s. For τ_est = 1, the C matrix at the end of the optimization (step 10000) is also shown. B Evolution of Q errors for the optimizations with various τ_est in A. bf C C errors corresponding to the best matrices in A (red solid curve), as well as last matrices of the optimization (red dashed curve). The green curve indicates the error of the best C with the symmetric part of the original C. The blue curve indicates the error between the empirical and the corresponding theoretical matrix. D Pearson correlation coefficient between the best and original C (red) as well as reconstructed and objective Q^0/τ. between and each corresponding theoretical Q^τ matrix for the same τ = τ_est in C. E Comparison between LO with empirical calculated with Eq (2) (red solid curve) and LO with theoretical obtained from Eq (7) (red dashed curve). The blue curve indicates the inaccuracy between the empirical and theoretical Q^τ matrices. F Same as B for the Q error with theoretical objectives .

https://doi.org/10.1371/journal.pcbi.1004762.g003

The left panel of Fig 3B shows the evolution of the Q errors for the optimizations in A. For small or large τ_est > 0 compared to τ_x = 1 s, the Q error often exhibits a large rebound similar to an explosion (dashed curves). LO being less stable for small and large τ_est > 0, it gives a poor estimate for C. In the case τ_est = 0, the optimization using Eq (19) is rather stable; however, the resulting estimate is closer to the symmetric part of the original C, not the directed matrix. The rate of change for C is displayed in Fig 3B (right): during the rebound or plateau of the Q error, the C matrix stabilizes for τ_est = 0 and 1 s (solid curves).

Fig 3C recapitulates our first main theoretical finding: τ_est needs to be matched with τ_x = 1 s in order to obtain an accurate estimation (see C errors in red). The dashed red curve corresponds to the C for the minimum of the rate of change (red curves in Fig 3B), yielding a better C recovery than the solid red curve that corresponds to the minimum of the Q error (blue curves). The green curve indicates the matrix distance to the symmetric part of the original C, namely 0.5(C+C^†). It shows that LO for small τ_est > 0 recovers the undirected connectivity like τ_est = 0, rather than the correct directed one. Finally, we observe that the best Q error corresponding to the minimum over the optimization (blue curve) is not a good indicator for the C error (red curve): it does not provide useful insight for choosing τ_est. The performance is further illustrated in Fig 3D by the Pearson correlation coefficient between the best and empirical C (red curve) and Q. The best C for τ_est = 1 s corresponds to values larger than 0.8, which means that the weight ranking of the original C is well captured.

In Fig 3E, we compare LO for empirical objectives (red solid curve) with LO for theoretical objectives (red dashed curve). For τ_est > 0, even theoretical objectives give a large C error, meaning that the poor performance comes from the procedure itself. This can be intuited from Eq (15), which relies on a “contrast” between and : for small τ_est, these two matrices are too similar and the C update is subject to numerical imprecision. This is confirmed by the evolution of the Q error in Fig 3F. For large τ_est, however, theoretical objectives yield a very good estimation. This means that the poor performance for large τ_est comes from the inaccuracy in the empirical , as represented by the blue solid curve in Fig 3E.

Need for estimating intrinsic noise in addition to connectivity.

Now we extend the LO procedure to adjust the variance of the noise individually inputted to the network nodes: each diagonal element Σ_ii is increased or decreased such that the node variance of the model converges toward its objective , see Eq (20) in Methods. To verify the estimation accuracy, we consider a target network with Σ_ii elements uniformly distributed between in [0.1, 0.6] and perform LO for the theoretical objectives and with τ = 1 s obtained using Eqs (6) and (7). The match between the model and objectives for both C and Σ is shown in Fig 4A, corresponding to a quasi perfect estimation of the connectivity as in Fig 2B for the tuning of C only.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Need for tuning intrinsic noise Σ for each node in addition to C.

A Match between the recovered and the original C (red crosses) and Σ (purple triangles) for LO based on Eqs (15) and (20) with theoretical objectives and τ_est = 1 s. The original network has inhomogeneous noise Σ_ii. B Comparison of Q error, C error and Pearson correlation coefficient for non-zero weights of the original C for four optimizations: ‘opt’ where both C and Σ are tuned as in A; ‘low’ for the tuning of C with a fixed Σ with the minimum value of the original Σ; likewise, ‘med’ and ‘high’ with a fixed Σ with the mean and maximum value of the original Σ. C Match between the original and recovered C for the optimization with fixed homogeneous Σ for three distinct values: the minimum, mean and maximum of the values in the original Σ. They are to be compared with the left panel in A.

https://doi.org/10.1371/journal.pcbi.1004762.g004

To assess the importance of tuning Σ, we perform three optimizations using a fixed homogeneous Σ for three levels of noise, namely the minimum, mean and maximum of the original Σ. As illustrated in Fig 4B, the Q error and, more importantly, the C error are larger for those cases (circle, square and triangle) than that with the tuning of Σ (crosses). Moreover, the values for the non-zero weights in the original C are estimated less accurately, as measured by the Pearson correlation coefficient. Fig 4C shows for each case the corresponding matches between the original and estimated C matrices. This is the second main theoretical finding: the optimization of Σ is as important as choosing τ_est to obtain an accurate estimate for C.

In practice, we find that a fast tuning for Σ works best in order to keep the optimization stable, i.e., avoid early “explosion” of the Q error. The presented simulations corresponds to ϵ_Σ = 10⁻¹ in Eq (20), to be compared with ϵ_C = 2 × 10⁻⁴ to calculate ΔQ in Eq (14) for the C update, cf. Methods.

Benchmark of LO procedure and comparison to the direct and heuristic methods.

Fig 5 combines the results found in Figs 3 and 4 to verify the robustness of LO when varying the network parameters. Fig 5A compares the C error for two network topologies: the lighter red curve corresponds to the cluster-hub type as in Fig 1A and the darker curve to randomly connected networks. The need for matching τ with τ_x is confirmed in the left panel for both types. The Σ error in purple is close to perfect for all τ_est. The error bars correspond to the variability over 40 networks with randomly chosen parameters: 25 ≤ n_sim ≤ 50 simulations to calculate the empirical ; network size 50 ≤ N ≤ 100; probability of connection p_con ∈ [10%, 30%]; maximum recurrent weight c_max ∈ [1.2, 2] × 1/N/p_con; input strength e ∈ [0.8, 1.2]; and average noise level σ ∈ [0.3, 0.6], with individual values shuffled ±50%; see also Table 2. No significant difference is observed there. On average over all 80 networks, LO with τ_est = 0.5 and 1 s give the best performance. On a technical ground, following the observations in Fig 3B and 3C, here the best C is chosen after the minimum for the Q error, provided the rate of change for C decreases (i.e., C keeps stabilizing) and the Q error remains smaller than 1.5 times its minimum.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Robustness of the estimation procedure of C and Σ on artificial networks.

A Left: Comparison of the C error (red curves) and Σ error (purple curve) for optimization based on several τ_est. For the C error, the cluster-hub topology in Fig 1A corresponds to lighter red and the random connectivity to darker red. The plot is similar to Fig 3C; the error bars correspond to the standard deviation over 40 network configurations. Right: Similar plot for the error of asym_C for the two network topologies. B Pearson correlation coefficient between the matrix elements of the whole recovered and original C matrices in red; idem but limited to non-zero elements of the original C in yellow; and between the recovered C and the symmetric part of the original C in green. The plot concerns 50 networks of either topology for which the asymmetry of original C is larger than 0.5; LO is performed using τ_est = 1 s. C C and Q errors (in red and blue, respectively) plotted against several parameters for the network and LO: the number of simulations used to calculate the empirical , the network size N and the asymmetry of the original C. Each cross/dot corresponds to one of the 80 simulations of either topology in A for τ_est = 1 s. D Comparison of the estimation performance of LO wit the direct and heuristic methods for the C error (left), the Σ error (middle) and asymmetry error (right). For each case, the left boxplot correspond to 40 cluster-hub networks and the right one to 40 randomly connected networks; the optimizations are performed with τ_est = 1 s. Contrary to others, the asymmetry error is not normalized.

https://doi.org/10.1371/journal.pcbi.1004762.g005

The right panel of Fig 1A shows that the asymmetry of the original C is correctly estimated for matched τ and τ_x. This means that directed connections are well estimated. Fig 5B further illustrates the quality of the estimation for the 50 networks with the strongest asymmetries: for small τ_est ≥ 0, the estimated C is closer to the symmetric part of the original C (largest Pearson correlations in green). For τ_est = 1 s, the C recovers well the directed connections in the original C, as shown by the Pearson correlations in red (all connections) and yellow (non-zero connections in the original C).

Fig 5C shows how LO performs depending on the amount of observed data and network parameters. As expected from previous work [25], the optimization performs better in Fig 5C when a larger number n_sim of simulations of the network activity is used to compute and , see Eq (2). Connectivity in larger networks is more difficult to estimate, but the asymmetry of the original C does not affect much the performance.

Finally, we compare our method with the direct and heuristic methods corresponding to Eqs (8) and (21), respectively. For each case in Fig 5D, the left and right box plots correspond to the optimization with with τ_est = 1 s for 40 cluster-hub and 40 randomly connected networks as in Fig 5A, respectively. While the Σ error is very small for all methods, the C error is much smaller for LO than the others. The direct method gives large complex values in the Jacobian obtained from Eq (8) due to the “noise” in empirical , which severely impairs the estimation. The heuristic method also yields a poor C error, although the Q error is very low: it does not capture the network effect correctly, namely the non-linear mapping between C and Q due to the recurrent feedback. Last, the asymmetry of the original C is best estimated by LO.

Hemodynamics and temporal information in BOLD covariances

Our method relies on information in time-shifted covariances to infer the connectivity of the underlying network. In the context of cortical models of BOLD activity, EC classically denotes the connections between neural populations, whose activity is transformed into BOLD signals via a hemodynamic response function (HRF). Before applying our method to fMRI experimental data, we need to test whether BOLD time series as obtained from the HRF convey information in their time-shifted measures about EC.

Therefore we compare our noise-diffusion (ND) network with the dynamic mean-field (DMF) model based on AMPA and NMDA dynamics [27] equipped with the Balloon/Windkessel model for the HRF [29]. As illustrated in Fig 6A, the synaptic input activity of the DMF (u_i in Methods) is transformed by the HRF to generate the BOLD signal. We consider the two topologies in Fig 6B: a randomly connected network and a network whose connections correspond to the SC obtained from dwMRI (see Methods). The SC connectivity has 32% density corresponding to the largest dwMRI values, as will be used in the next sections for the ND model applied to fMRI data. The model parameters are taken from previous publications [27, 29] and recapitulated in Table 3. In Fig 6C, we simulate 50 configurations of the DMF+HRF model with the two topologies (‘DMF/rnd’ and ‘DMF/SC’, respectively). For each, we choose weights randomly and calculate the empirical autocovariances of the BOLD signals. We find significant values up to time shifts equal to 3 s, which then drop. Interestingly, the slope in the log plot for the smaller time shifts is close to that found in experimental data (black line) up to 2 s of time shift. Then we estimate the similarity between the neural covariances at a short timescale and the FC given by the BOLD covariances, using the Pearson correlation coefficient between the corresponding matrix elements. As shown by the Pearson correlations in Fig 6D (left panel), the match between the neural and BOLD covariances is good, especially when EC is constrained by SC. This means that, although the HRF involves nonlinearities, the DMF model preserves the ranking of the spatio-temporal covariances is preserved between neural and BOLD activity, as illustrated in Fig 6E. In particular, the time-shifted BOLD covariance matrices are asymmetric.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Information conveyed by spatio-temporal FC and hemodynamic response function (HRF).

A Schematic representation of the properties of the two dynamic cortical model: ND network whose activity directly models the BOLD signals; dynamic mean-field (DMF) model whose synaptic activity is processed by the HRF to obtain the BOLD signal. B Two types of neural connectivity used with the DMF: random and SC matrix obtained from dwMRI. In the matrix, darker pixels indicate a higher probability of existing fibers between cortical areas. C Autocovariances of the two models for 50 different simulations. The y-axis has a log-scale. For each simulation, the curves are centered vertically with respect of the mean over all nodes to focus on the slope. The black line represents the exponential decay fitted on the mean experimental BOLD with time constant 5.3 s. D Left: Similarity between neural and BOLD covariances (excluding variances) for the two considered DMF models, as measured by the Pearson correlation coefficient. Right: Performance of EC estimation from BOLD FC for the DMF models. The performance is measured by the Pearson correlation coefficient. For each simulation, the objective is the average BOLD FC taken from 50 simulations of the same network. E Example of a typical mapping between neural and BOLD covariances for a DMF/SC network. The covariances have been rescaled. F Variability of the EC estimated from individual FCs. The grey dots represent the match of four EC each estimated for a single simulation of 300 s. The red dots correspond to the average over 50 estimated EC for 50 simulations of the same DMF/SC network. G Uncertainty of the estimated EC as a function of the estimated weight for each neural connection. The y-axis is the standard deviation over the 50 optimizations in F divided by the mean.

https://doi.org/10.1371/journal.pcbi.1004762.g006

Now we use the LO procedure to estimate the original EC from the BOLD FC with τ_est = 1 s. The right panel in Fig 6D indicates the Pearson correlation between the original and estimated connectivity. For the DMF model, the agreement corresponds to 0.6 and 0.65 for the random and SC topologies, respectively. This means a fair recovery of the ranking between the original EC weights. For the DMF/SC networks where the existing connections are “known” (32% density), LO only tunes the corresponding weights in C.

Finally, we find that the variability of the estimated ECs decreases with the strength of the neural connection. This suggests that strong estimated EC values can be trusted based on a threshold on the variability of the estimations (standard deviation divided by the mean). Together, these results support the applicability of our method to fMRI data and directly estimate the EC strengths from the average FC.

Application to estimation of EC from fMRI data during rest

Now we use the ND model to reproduce the cortical FC obtained from fMRI recordings during rest. In contrast to previous studies [12, 27], FC involves both zero-time-shift and time-shifted covariances of the BOLD signals for the 68 cortical regions here. Following the results in Fig 6, we denote by EC the connectivity C in our model, even though the network activity directly models the BOLD signals without a HRF. The data set corresponds to 25 subjects aged from 20 to 39 years. The empirical covariances are calculated using Eq (2) for τ that are multiples of 2 s, which is the temporal resolution of the BOLD time series. Fig 7A displays the objective FC matrices and with τ = 4 s, which are averages over all subjects. The autocorrelograms of all 68 nodes are displayed in the inset of Fig 7B. The main panel shows the same curves with a log-scale for the y-axis: they seem quasi straight lines indicating exponential decays, as observed for artificial networks in Fig 1E. Therefore, we estimate from the mean BOLD autocovariance over all nodes (red curve of Fig 7B) the time constant τ_x = 5.3 s, which we use to calibrate the model. We use a single τ_x for all nodes, because the distribution of individual time constants over all nodes is narrow, as shown in Fig 7C.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Fitting fMRI data to infer cortico-cortical connectivity.

A Empirical FC, namely and , averaged over 25 subjects for BOLD signal recorded during rest. B Autocovariances (grey curves, with a log y-axis) as a function of the time shift τ. From the mean curve (red) the time constant for the ND model is estimated: τ_x = 5.3 s. The inset represents the autocovariances with a standard y-axis. C Variability of similar decay time constants to τ_x in B for individual cortical area. D SC mask of weights to tune (black pixels), corresponding to the strongest connections of the dwMRI matrix in Fig 6B thresholded to obtain 32% density. E Evolution of the Q error during the optimization. The black line indicates the rate of change of C. F EC corresponding to the estimated C matrix. G Match between empirical and reconstructed Q matrices for the optimization with τ_est = 4 s while the connections for the mask in D. For both Q⁰ and Q^τ, the blue dots correspond to tuned connections, the black crosses to non-diagonal untuned connections and the cyan dots to diagonal connections. The Pearson correlation coefficient between matrix elements are given above each plot. H Weak correlation between the symmetric part of the estimated C and , indicating a strong network effect.

https://doi.org/10.1371/journal.pcbi.1004762.g007

The direct method in Eq (8) does not work here: the reason is that the matrix logarithm gives complex values with large imaginary parts. This motivated the development of the LO procedure. Moreover, our approach allows for the use of constraints on the weights in C: only a subset of all possible cortico-cortical connections can be tuned. The mask in Fig 7D represents the anatomical cortical connectome obtained from thresholding the dwMRI data averaged over all subjects: connections with large SC values are optimized, other weights are kept equal to zero. The density of existing connections is 32% in Fig 7D.

The blue curve in Fig 7E represents the evolution of the Q error during LO with τ_est = 4 s: it first decreases and then “explodes”. This instability corresponds to a transition to excessively high activity due to too strong feedback. The black curve show the normalized rate of change for C, whose minimum is close to that of the Q error. For the experimental data, the minimum of the black curve comes before the minimum of the Q error, but the rate of change of C remains close to its minimum until the explosion. In the following, we choose the estimated EC in Fig 7F and reconstructed FC as the Q and C for the minimum of the Q error. The match between the model and empirical Q matrices is illustrated in Fig 7G. Importantly, the match is similar for the Q corresponding to tuned connections (blue dots) and absent connections (black crosses). The fit of variances corresponding to the diagonal elements of Q (cyan dots) is very good.

Fig 7H indicates the poor match between the empirical FC and the symmetric part of estimated C: large values in arise from strong network effect even though the EC weight is low.

Robustness of LO applied to experimental data.

Before analyzing the properties of the estimated EC, we check the consistency of the C and Σ estimated by LO. As shown in Fig 8A, the best C matrices for τ_est ≥ 2 s describe a similar network topology, i.e., approximately the same ranking between weights (the black line indicates a perfect match). In contrast, the best C obtained from Q⁰ only does not match the best C for Q⁰ and Q⁴, as illustrated in Fig 8B. The estimated intrinsic noise Σ in Fig 8C is similar for all τ_est ≥ 0.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Robustness of estimated C and Σ from FC with different time shifts τ_est.

A Match between C for τ_est = 4 s (x-axis) and C for τ_est = 2, 6 and 8 s (y-axis). Each plotted point corresponds to an EC connection and the black line indicates a perfect match. B Same as B for the EC estimated for τ_est = 0 on the y-axis. C Similar plot to A for the diagonal elements of the estimated noise matrix Σ. D Similar plot to B for the EC estimated using the direct method. E Agreement as measured by the Pearson correlation coefficient between the pairs of EC and Σ plotted in B and D for the LO (red), heuristic (purple) and direct (cyan) methods. The green error bar on the left represents the same between EC with τ_est = 0 and the four τ_est > 0. F Comparison of the models tuned with LO using different time shifts. For each model C obtained with τ_est from 0 to 6 s, the curve represents the Q^τ error for all τ_mod on the x-axis.

https://doi.org/10.1371/journal.pcbi.1004762.g008

The direct method gives inconsistent C for the different τ_est > 0, as illustrated in Fig 8D to be compared with LO in A. Fig 8E measures the similarity between these matrices using the Pearson correlation coefficients between the weights corresponding to the mask in Fig 7D. LO (red) provides very consistent C for the different τ_est, the heuristic method performs worse (purple) and the direct method (cyan) very poorly. The Σ consistency, however, is good for all methods. These findings are in line with the conclusions were found for artificial networks in Fig 5.

Each curve in Fig 8F represents a model whose C has been estimated using LO with only (green curve); and with τ_est = 2, 4 and 6 ms (blue to cyan). For each model, the Q^τ error between the model and empirical matrices is displayed for τ = τ_mod on the x-axis. The errors for the two darker blue curves are smaller than the green one for τ_mod ≥ 2 s, although the fit for τ_mod = 0 is slighly worse. Our optimization procedure thus improves the whole spatio-temporal Q fit, from to , as compared to LO based on only. Directed connectivity is synonymous with capturing the information in the time-shifted with τ > 0. In the following, we take as reference for EC the C corresponding to τ_est = 4 s and 32% density in Fig 7F (and darkest blue curve in Fig 8F).

To further check the consistency of the estimated C, we repeat the optimization procedures with various SC masks, similar to that in Fig 7D. By moving the threshold on the dwMRI values, we obtain various densities for the mask ranging from 25% to almost 100% in Fig 9A. Although the Q error decreases for denser C, the estimated C becomes more similar to , and the mapping between C and Q more linear. Importantly, the estimated topology of C estimated by LO for the 4 densities is consistent, as measured by the Pearson correlation in red in Fig 9B. In comparison, the heuristic method is less consistent (purple). Likewise, Fig 9C shows that the asymmetry of the C matrices is consistent across SC masks and τ_est for LO (red), but less for the heuristic method (purple).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Robustness of estimated EC with respect to the EC density and individual FC.

A Fit performance as a function of the density of EC. The blue curve indicates the Q error for optimizations with several masks for C (x-axis); the error bar correspond to the 4 τ_est > 0. The red curve indicates the similarity between the estimated C and the empirical . B Similarity of estimated C measured by the Pearson correlation for the different masks in A and τ_est = 4 s. C Input-output asymmetry of the estimated C obtained for various optimization masks (error bars) and τ_est (x-axis). The red curve corresponds to LO and the purple curve to the heuristic method. The blue curve indicates the asymmetry of the empirical FCτ. D Q Error for LO based on individual FCs with τ_est = 2 and 4 s with 32% density. The horizontal dashed line corresponds to the Q error in Fig 7. E Similarity between the individual ECs and the EC obtained from the mean FC in Fig 7F. Similarity is measured by the Pearson correlation between the non-zero EC weights for the 32% density mask. F EC asymmetry of the individual ECs. The horizontal dashed line indicated the asymmetry of the EC in Fig 7F. G Agreement between the mean EC over individual FCs and the EC estimated from the mean FC.

https://doi.org/10.1371/journal.pcbi.1004762.g009

Finally, we perform LO for each of the 25 individual empirical FC. The Q errors are plotted in Fig 9D for τ_est = 2 and 4 s. The errors are similar for the two cases and are larger than the Q error for the average FC (dashed horizontal line). Fig 9E shows the good agreement between the individual ECs and the EC obtained from the mean FC in Fig 7F, with Pearson correlation coefficients around 0.7 for both τ_est = 2 and 4 s. As shown in Fig 9F, the asymmetry of the individual ECs is homogeneous around 0.4, which is slightly higher than the asymmetry of the EC from the mean FC that is equal to 0.35 (and indicated by the horizontal dashed line). Last, Fig 9G shows the excellent agreement between the mean of individual ECs and the EC from the mean FC. This confirms the consistency of the results obtained from our LO procedure.

Interpretation of estimated cortical EC

Fig 10A shows the very weak match between the C and the SC corresponding to the dwMRI data averaged over all subjects, with a Pearson correlation coefficient equal to 0.06. The estimated EC is thus structurally very different from a scaled version of SC. As dwMRI reflects the density of cortico-cortical white-matter fibers, this suggests that the efficacies of these connections are determined by other factors than their size, such as types of neurotransmitters, concentration of synaptic receptors and excitability of cortical areas.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Interpretation of estimated cortical interactions and intrinsic noise.

A Weak match between EC and SC values. Each red cross corresponds to a connection. B Plot of the sum of incoming C weights versus the sum of outgoing C weights (x-axis and y-axis, resp.) for each cortical area (red cross). C 3D mapped representations on the cortex of the variances of the empirical BOLD signal, sums of incoming weights in EC, and intrinsic noise (related to spontaneous noisy activity). Red corresponds to larger values and blue to small values.

https://doi.org/10.1371/journal.pcbi.1004762.g010

The estimated EC matrix in Fig 7F is not symmetric, meaning that cortical interactions are not reciprocal. This information is important as dwMRI data estimate the density of axonal fibers in the white matter, but do not recover the direction of those fibers. SC is by construction quasi symmetric, as is the case in Fig 7D. The matrix asymmetry for C relates to the reciprocity of intracortical connections and can be seen in the difference between incoming and outgoing strengths in Fig 10B. No area has both large incoming and outgoing weights, meaning that hubs act either as receptors or feeders.

In Fig 10C, the results are mapped onto the cortical surface. Cortical feeders with the largest outgoing weights are the left and right fusiform, middle temporal and superior temporal gyri, as well as the pre- and postcentral gyri in the left hemisphere. Cortical receivers with the largest incoming weights are the left and right precuneus, lateral occipital and superior parietal gyri, as well as the left isthmus of the cingulate gyrus. We also find that the following areas exhibit the largest values for Σ, synonymous with strong intrinsic variability: both left and right lingual gyri, pericalcarine cortices and frontal poles, as well as right cuneus and transverse temporal gyrus, and left pars orbitalis. This suggests the propagation of spontaneous activity, mainly from visual cortices and the prefrontal area.

Inference of directed connectivity from observed network activity

It is known that zero-time-shift covariances are not sufficient to retrieve directed connectivity, but only its symmetric part [16, 25, 28]. Information-based methods able to estimate directionality such as likelihood maximization [13], Granger causality [14, 23] and transfer entropy [24] also use temporal information of the observed activity. In minimizing the matrix distance between the model Q^0/τ and objective , instead of considering connections independently, the LO procedure captures network effects due to the recurrent feedback. Here we do not perform a stochastic gradient descent using many samples of the observed activity, but a deterministic optimization based on the averaged over the whole observation period (or several simulation sessions). It follows that the optimization is quick: a few minutes for 10⁴ optimization steps with a network of 50 nodes and a given τ_est on a recent desktop computer.

To obtain the best performance, we have shown that the time shift τ_est corresponding to used in LO and the time constant of the dynamical system τ_x should be matched. As shown in Fig 3E, poor estimation for large τ_est arises from the inaccuracy in empirical ; for small τ_est, LO itself is unstable (see dashed curve in Fig 3F). Fig 5C shows that the method performance is not strongly affected by the network topology or the connectivity asymmetry, but worsens with the network size and becomes better with more observations used to calculate the empirical objective . This is in line with previous results [25]. In addition to C, the intrinsic noise Σ received by the network nodes must be tuned to obtain a correct C estimation (Fig 4). Here we use a heuristic optimization for a diagonal Σ in Eq (20); the present framework should be extended to take into account correlated noise instead of white noise.

The direct method in Eq (8) to estimate C has been used previously with statistical tests to estimate the existence of connections from observed activity [16, 26]. As shown in Fig 5D, it does not work well for the level of noise in empirical observations considered here. This motivated the development of our LO procedure in Methods. Here we focus on the estimation of connection weights, i.e., their ranking; the detection of connections using statistical methods can be based on the estimated C from LO, but this is left for further work. Nevertheless, detection should be based on an as-good-as-possible estimated ranking of connections weights, which can be measured by the Pearson correlation. Our method for Ornstein-Uhlenbeck processes also bears similarities with MVAR [3], but we enforce additional constraints on the connectivity.

Recent studies [12, 28] have also used greedy algorithms to optimize symmetric C relying on using more elaborate network models. Those procedures update C step-by-step according to various measures such as the Pearson correlation between all matrix elements of zero-time-shift correlations. Here we have transposed this heuristic method to reproduce in addition to , see Eq (21). Although the resulting Q^0/τ fit is close to perfect, the C estimation remains poor in Fig 5D. Taking the network effect correctly into account via LO is important to recover the original C, as compared to tuning connections individually based on the corresponding Q^0/τ value.

More generally, the problem with inferring C lies in the definition of observables or objective functions to constrain models without ambiguity: to a set of network parameters should correspond only a single value of the observable (here a pair + ). In [25], the minimization of matrix L1 norm for sparse networks was used to reduce the indetermination in using only. Beyond our results based on noise-diffusion processes, we expect that directed C can be recovered for other network models such as Hawkes processes (a.k.a. Poisson neurons) or binary neurons using second-order moments with non-zero time shift (here covariances). This is supported by recent results that demonstrate how the covariance structures are formally related across these neural models [38].

The present framework appears well suited to model activity as continuous signals; for spike trains generated by networks of GLM or Poisson neurons with simple temporal filters, it remains to be seen whether our fast tuning procedure can be adapted. We have only considered resting-state activity, but the procedure may also be extended to the case of multiple stimulation-response pairings. In particular, the external input e may be adjusted [16], in addition to C and Σ.

Whole-cortex dynamic model fitted to fMRI data and cortical interactions

The goal of our model-based approach is to reproduce the resting-state FC obtained from fMRI. Although the ND model is not new, we propose a novel ‘LO’ procedure to tune the model parameters with suitable observables, namely time-shifted covariance matrices as FC. The estimated connectivity of the ND model relates to the intracortical EC, whose properties can then be analyzed [1]. For instance, EC can be searched for hubs, communities and similar features [30]. Here we have focused on hubs (Fig 10) and our results suggest that activity propagates from the visual, auditory and prefrontal areas. We find among the listening hubs the precuneus and superior parietal gyrus that belong to the default-mode network. These networks are usually found in resting-state FC and our results shed a new light on the architecture that shapes the activity propagation between them.

The ND model was previously used together with hemodynamics in order to reproduce similar fMRI data [12, 27]. As was demonstrated for that study for zero time shifts, the BOLD covariances with time shifts of the order of seconds also convey information about the interactions that determine the neural dynamics (Fig 6). This is in line with recent results showing that the BOLD time series convey information about cognitive processing for similar time lags [39, 40]. In our model, the EC asymmetry generates lags in covariances. Moreover, the LO procedure can estimate the EC ranking this neural connectivity from the BOLD covariances provided the ranking is preserved between the neural and the BOLD covariances for those more elaborate models involving a HRF (Fig 6). This supports the application of the ND model without HRF to model directly the BOLD data. In particular, the ND model seems well adapted to fMRI data at the considered parcellation of about 100 nodes: the BOLD autocovariances are close to exponential decays (i.e., straight lines for the log y-axis), as shown by the comparison between the ND model in Fig 1E and our experimental data for time shifts up to 8 s in Fig 7B. The corresponding time constants are rather homogeneous over all regions in Fig 7C, so we use a single τ_x = 5.3 s to calibrate the ND model, which corresponds to typical values for HRF [41]. In contrast, the autocovariances for the DMF+HRF in Fig 6C show a drop for τ ≥ 4 s that is not observed in the experimental data. We keep in mind that the precise relationship between the fluctuations of the BOLD signal and neural activity is still under debate [41–43]. Numerous studies such as those about repetition priming and suppression [44] have shown how changes in fMRI signals reflects those in neural activity, such as synchronization of the latter at a much shorter timescale. Moreover, fMRI is one of the few non-invasive methods to evaluate processing in the human brain (of course in vivo) and has many clinical applications [45], irrespective of its precise link to neural activity. This further supports efforts to develop generative models of BOLD and methods to better interpret these fMRI data.

An underlying assumption of our approach is that the individual variances of BOLD signals are meaningful [46]. This motivates the use of covariances to fit our model, rather than correlations that are often used [1]. The results in Fig 7 give a FC Pearson correlation coefficient (sometimes called predictive power) larger than 0.6 for both Q⁰ and Q⁴. In this sense, our study improves previous results [12, 31]. Importantly, the EC and Σ estimated by LO is surprisingly stable with respect to the choice of time shift for the considered experimental data. We find consistent results for a broad range of τ_est from 2 to 8 s, with a Pearson correlation coefficient larger than 0.9 in Fig 8. Choosing denser connectivity for EC improves the FC fit, but does not significantly change the EC ranking as shown in Fig 9A and 9B for results from 25% to 45% density. As an additional check, we have also shown that EC estimated from individual FCs coincides with EC obtained from the mean FC over all 25 subjects. In comparison, the heuristic and direct methods in Figs 8E and 9B and 9C give far more inconsistent results. Moreover, the heuristic method only uses part of the FC matrices when EC is sparse: it tends to overfit the existing connections compared to absent connections. A recent study [15] compared methods to estimate directed cortical interactions in a generative model of BOLD activity including hemodynamics: among those, Patel’s τ[47] gave the best performance, but that study was limited to sparse connectivity. Here we consider the FC for the whole cortex, which gives coupled constraints for a reasonably large network. Incorporating more areas limits the number of unknown contributions to the node activities, strengthening the estimation accuracy. This was observed for partial correlations in a similar manner for undirected connections [15].

Several previous studies used a scaled version of the SC matrix as EC [12, 28, 31]. As shown in Fig 8F, a symmetric EC does not satisfactorily fit time-shifted FCs. Furthermore, Fig 10A shows the weak match between dwMRI and our estimated EC values. This suggests that the fiber densities as measured by dwMRI may not be a good predictor for the dynamic interactions between cortical areas, but only for the skeleton of the cortical connectome. In addition, those previous studies used interhemispherical connections to improve the FC fit of the models. Here these connections also appear to be fairly strong in the estimated EC in Fig 7F, as shown by the two secondary diagonals in the top-left and bottom-right quadrant of the matrix. In any case, large covariances are obtained even for unconnected areas because of the strong overall feedback in the network, in line with results supporting that the cortex is in a state close to criticality [27]. As the dynamics of each cortical region has no nonlinearity in our model, large EC values are the only origin for the strong network effect, which induces the non-trivial mapping between EC and FC in Fig 7H.

As shown in Fig 1F for artificial networks, the C asymmetry is reflected in Q^τ for τ > 0. In our experimental data, the corresponding FCτ (blue line in Fig 9C) is not symmetric and relates to propagating activity [39]. In the ND network, the noise received by each area is embodied by Σ, which is then shaped by EC to generate FC. Here the LO procedure extracts the spatio-temporal FC information to estimate the EC asymmetry. The information about asymmetry in intracortical connections is thus important and complements the SC obtained from dwMRI, which is symmetric. Compared to previous studies that investigate the directionality of cortico-cortical connections [2, 4, 7], the novelty is that we estimate this property at the scale of the whole cortex. As all nodes have the same activation function, EC values indicate the relative interaction strengths between areas. What matters in our EC analysis is the ranking of EC weights: the difference between a low EC value and an absent connection is not so important here. Here we have estimated the asymmetry of cortico-cortical EC to be equal to 0.35 in Fig 8C. This is larger than macaque’s COCOMAC asymmetry that gives 0.14 for the same index [48, 49], but note that nonlinearity in the model dynamics would affect the precise EC values; what is important here is that EC is significantly asymmetric. This translates to an imbalance between incoming and outgoing EC weights. Fig 10B suggests that hubs are either feeders or receivers, but do not have both strong incoming or outgoing connections. As a conclusion, our work sets a suitable ground to study both local and global properties of the whole cortical connectivity, which will give insight in the underlying neural processing.