Model-based whole-brain effective connectivity to study distributed cognition in health and disease

Neuroimaging techniques are now widely used to study human cognition. The functional associations between brain areas have become a standard proxy to describe how cognitive processes are distributed across the brain network. Among the many analysis tools available, dynamic models of brain activity have been developed to overcome the limitations of original connectivity measures such as functional connectivity. This goes in line with the many efforts devoted to the assessment of directional interactions between brain areas from the observed neuroimaging activity. This opinion article provides an overview of our model-based whole-brain effective connectivity to analyze fMRI data, while discussing the pros and cons of our approach with respect to other established approaches. Our framework relies on the multivariate Ornstein-Uhlenbeck (MOU) process and is thus referred to as MOU-EC. Once tuned, the model provides a directed connectivity estimate that reflects the dynamical state of BOLD activity, which can be used to explore cognition. We illustrate this approach using two applications on task-evoked fMRI data. First, as a connectivity measure, MOU-EC can be used to extract biomarkers for task-specific brain coordination, understood as the patterns of areas exchanging information. The multivariate nature of connectivity measures raises several challenges for whole-brain analysis, for which machine-learning tools present some advantages over statistical testing. Second, we show how to interpret changes in MOU-EC connections in a collective and model-based manner, bridging with network analysis. Our framework provides a comprehensive set of tools that open exciting perspectives to study distributed cognition, as well as neuropathologies.


A Movie dataset
The experimental procedure was approved by the Ethics Committee of the Chieti University and the participants signed a written informed consent. We use the data from 22 right-handed young and healthy volunteers (15 females, 20-31 years old) that had recordings for both a resting state with eyes opened and a natural viewing condition (30 minutes of the movie 'The Good, the Bad and the Ugly').
The fMRI time resolution (TR) is 2 seconds and each session of 10 minutes corresponds to 300 frames (or volumes). The preprocessing performed using SPM8 involves the coregistration of structural and functional MRI, spatial normalization in the MNI coordinates and motion correction. The data are parcellated in N = 66 ROIs (Hagmann et al., 2008). The ROIs are listed in Table 1, indexed in the order of Figs. 5 and 6 in the main text. The generic structural matrix is obtained from averaging individual matrices obtained using the diffusion spectrum imaging. A threshold was then applied to retain the 27% largest connections. Inter-hemispheric connections between homotopic ROIs were added, increasing the density to 28%. The same generic SC matrix is used to determine the MOU-EC topology for all subjects.
We refer the readers to previous publications (Hlinka et al., 2011;Mantini et al., 2012; for further details.

B Functional connectivity measures
For each session of T = 300 time points (2 for rest and 3 for movie), we denote the BOLD time series by s t i for each region 1 ≤ i ≤ N with time indexed by 1 ≤ t ≤ T . The mean signal is s i = 1 T t s t i for all i. Following (Gilson et al., 2016), the empirical FC consists of two BOLD covariance matrices (see the blue and green matrices in the top box of Fig. S1), without and with time lag: Pearson correlation can be calculated from the covariances in Eq. (1): (2) C Multivariate Ornstein-Uhlenbeck (MOU) process to model whole-brain dynamics The network dynamics are determined by two sets of parameters: Functional data Whole-brain parcellation Figure S1: Dynamic model and optimization procedure to capture BOLD dynamics. From the BOLD signals of the parcellation (top left in top box), two empirical FC matrices (blue and green matrices) are calculated, both without and with time lag as depicted by the blue and green curve. The dynamic network model describes the activity of the ROIs (only two are represented by the grey circles) and comprises a parameter for the spontaneous activity of each ROI (Σ in the bottom box) and for each directed connection between ROIs (effective connectivity, MOU-EC). The topology of existing connections is determined by the structural connectivity (black matrix in top box). The optimization of the model is similar to a gradient descent during which the model FC matrices are evaluated at each step and compared to their empirical counterparts. From the differences between the corresponding matrices an update of the MOU-EC and Σ parameters is calculated and the procedure is repeated until reaching a minimum for the model error (light gray-blue curves in the left central panel between the boxes), which corresponds to a high value for the Pearson correlation between the matrix elements of the model and empirical FCs (black curve). Before the optimization, the time constant τ x of the BOLD autocovariances (see log-plot) is evaluated to calibrate the nodal dynamics in the model. Details can be found in (Gilson et al., 2016).
• the time constant τ x is an abstraction of the • the network effective connectivity (MOU-EC) between these ROIs embodied by the matrix C, whose topological skeleton is determined by structural data; • the local variability embodied in the matrix Σ inputed individually to each of the N = 66 ROIs.
The activity variables obey are described by a multivariate Ornstein-Uhlenbeck process. Each activity variable x i of node i decays exponentially according to the time constant τ x and evolves depending on the activity of other populations: Compared to our previous work (Gilson et al., 2016, we use a different convention for the weights C ji from ROI j to ROI i, in line with graph theory. In matrix notation it now reads Here the fluctuating inputs dB i are independent and correspond to a diagonal covariance matrix Σ, as represented in purple by the vector of variances Σ ii in the bottom box of Fig. S1. In the model, all variables x i have zero mean. Their spatiotemporal covariances are denoted by Q 1 ij , where τ ∈ {0, 1} is the time lag, and correspond to the the blue and green matrices in the bottom box of Fig. S1. They can be calculated by solving the consistency equations: The first equation is the continuous Lyapunov equation, which can be solved using the Bartels-Stewart algorithm in scipy (??, sci). Here J is the Jacobian of the dynamical system and depends on the time constant τ x and the network MOU-EC: J ji = −δ ii τx + C ji , where δ ji is the Kronecker delta and the superscript † denotes the matrix transpose. In the second equation e denotes the matrix exponential. Ideally, the model should be extended to incorporate subcortical areas and the relevance of input cross-correlations inputs should be evaluated for all ROI pairs .

D Parameter estimation procedures D.1 Lyapunov optimization or natural gradient descent
We tune the model such that its covariance matrices Q 0 and Q 1 reproduce the empirical FC0 and FC1 matrices in Fig. ??B, denoted here by Q 0 and Q 1 , respectively. The uniqueness of this maximum-likelihood estimation follows from the bijective mapping from the model parameters τ x , C and Σ to the FC pair (Q 0 , Q 1 ). The iterative optimization procedure for C is similar to the original version (Gilson et al., 2016), which can be related to the concept of "natural" gradient descent (Amari, 1998) that takes the non-linearity of the mapping between J and the matrix pair Q 0 and Q 1 in the second line of Eq. (5). Note that the parameters τ x and Σ also follow a gradient descent now.
For each session (subject and condition), we calculate the time constant τ ac associated with the exponential decay of the autocovariance averaged over all ROIs using time lags 0 and 1 TR: where a(v i |u) is the slope of the linear regression of v i = [log( Q 0 ii ), log( Q 1 ii )] by u = [0, 1], and τ TR is the value of the TR in seconds. Note that one can also use more time lags -e.g. up to 2 TRs-to assess the stability of the estimated τ ac , as was previously analyzed for resting-state data (Gilson et al., 2016). The model is initialized with τ x = τ ac and no connectivity C = 0, as well as unit variances without covariances (Σ ij = δ ij ). At each step, the Jacobian J is straightforwardly calculated from the current values of τ x and C. Using the current Σ, the model FC matrices Q 0 and Q 1 are then computed from the consistency equations, using the Bartels-Stewart algorithm to solve the Lyapunov equation.
The difference matrices ∆Q 0 = Q 0 − Q 0 and ∆Q 1 = Q 1 − Q 1 determine the model error where each term -for FC0 and FC1 -is the matrix distance between the model and data covariances, normalized by the latter.
The following Jacobian update can be applied to decrease the model error E at each optimization step similar to a gradient descent: It turns out that a modified update is more robust to empirical noise in practice which would corresponds to a proxy error based on the matrix logarithm: See (Insabato et al., 2018) for a comparison of optimization methods. From the Jacobian update ∆J, we obtain the connectivity update: for existing connections only; other weights are forced at 0. We also impose non-negativity for the MOU-EC values during the optimization.
To take properly the effect of cross-correlated inputs into account, the Σ update has been adjusted from the heuristic update in the first study (Gilson et al., 2016) to a gradient descent : As with weights for non-existing connections, only diagonal elements of Σ -and possible crosscorrelated inputs -are tuned; others are kept equal to 0 at all times. Last, to compensate for the increase of recurrent feedback due to the updated C, one can also tune the model time constant τ x as where λ max is the maximum negative real part of the eigenvalues of the Jacobian J. The rationale is to avoid an explosion of the network dynamics (when λ max → 0) while letting the model connectivity C develop a non-trivial structure to reproduce FC.
Repeating the parameter updates, the best fit corresponds to the minimum of the model error E.

D.2 Heuristic optimization
Instead of the derivation of the Jacobian in Eq. (10) that takes into account the nonlinearity of the mappings in Eq. (5), one can use a greedy algorithm to optimize the weights to fit an objective measure on the data. Here we update the weights for the model to reproduce the correlation matrix K in Eq. (2): Note that only correlations ∆K 0 ij for existing connections i → j and j → i are taken into account in this update.

E Dynamic communicability and flow for network analysis
Following our previous study , we firstly define dynamic communicability to characterize the network interactions due to the MOU-EC connectivity C, ignoring the input properties Σ. Our definition is adapted to study complex networks associated with realistic (stable) dynamics where time has a natural and concrete meaning. In comparison, a previous version of communicability for graphs (Estrada and Hatano, 2008) relied on abstract dynamics. The basis of our framework is the network response over time, or Green function, which is the basis of the concept of dynamic communicability, which focuses on the temporal evolution of such interactions. Although we focus on the MOU process here, our framework can be easily adapted to distinct local dynamics for which the Green function is known. In addition to the MOU-EC matrix C, the MOU dynamics is determined by the input properties, so we use the dynamic flow in the main text ( Figs. 6 and 8).

E.1 Definitions
Formally, dynamic communicability is the "deformation" of the Green function e Jt of the MOU process due to the presence of the (weighted and directed) matrix C, as compared to the Green function e J 0 t corresponding to the Jacobian with leakage only and no connectivity, J 0 ij = −δ ij /τ . It is defined as the family of time-dependent matrices depicted in Fig. 6A: The scaling factor ||J 0 || −1 = || t≥0 e J 0 t dt|| where || · || is the L1-norm for matrices (i.e., sum of elements in absolute value) is used for normalization purpose . Recall that t ≥ 0 here is the time for the propagation of activity in the network, referred to as 'impulse-response time' in the figures.
To incorporate the effect of local spontaneous activity or excitability (inputs in the model), we define the dynamic flow that fully characterizes the complex network dynamics . The input statistics of interest for a stable MOU process correspond to the input (co)variance matrix Σ, which are independent parameters from the MOU-EC matrix C. This is represented by the purple arrow in the left diagram of Fig. 6A, indicating that the fluctuation amplitude is individual for each ROI. The Σ matrix may be non-diagonal when ROIs experience cross-correlated noise , as represented by the purple dashed arrows. The dynamic flow describes the propagation of local fluctuating activity over time via the recurrent connectivity and is defined by the where √ Σ is the real symmetric "square-root" matrix of the input covariance matrix, satisfying Σ = √ Σ √ Σ † . Dynamic communicability is thus a particular case of the flow for homogeneous input statistics. From the time-dependent matrices F(t), we define the total flow that sums all interactions Total communicability for graphs has been used to define a version of centrality (Benzi and Klymko, 2013). Here the proximity between ROIs correspond to how much flow they exchange. We also define the diversity (akin to heterogeneity) among the ROI interactions in the time-dependent matrices F(t), which can be seen as a proxy for their homogenization over time: defined as a coefficient of variation where µ {i,j} and σ {i,j} are the mean and standard deviation over the matrix elements indexed by (i, j). We also define the input and output flows for each node:

E.2 Community detection
To detect communities from F(t) in Fig. 8 in the main text, we rely on Newman's greedy algorithm for modularity (Newman, 2006) that was originally designed for weight-based communities in a graph. Adapting it here to the flow matrix F(t) at a given time t, we seek for communities where ROIs have strong bidirectional flow interactions. In the same manner as with weighted modularity, we calculate a null model for MOU-EC: Note that we preserve the empty diagonal. The resulting matrix contains from the expected weight for each connection, given the observed input strengths c in and output strengths c out ; S C is the total sum of the weights in C. Then we caclulate F null (t) using Eq. (17) with C null instead of C. Starting from a partition where each ROI is a singleton community, the algorithm iteratively aggregates ROIs to form a partition of K communities denoted by S k that maximizes the following quality function: At each step of the greedy algorithm, the merging of two of the current communities that maximizes the increase of Φ is performed. Note that communicability-based communities can be defined similarly using C(t) and the corresponding null model C null (t).