Asymmetric high-order anatomical brain connectivity sculpts effective connectivity

Multimodal neuroimaging analyses are expected to improve our understanding of structure-function relationships in the brain (Toga et al., 2006; Honey et al., 2010; Sporns, 2014); drawing on measures of structural, functional, and effective brain connectivity (Sporns et al., 2000; Park & Friston, 2013). However, relating symmetric and static structural connectivity derived from diffusion magnetic resonance imaging (dMRI) to time-varying and context-sensitive functional dynamics (recorded by functional magnetic resonance imaging, fMRI, electroencephalography, EEG, or magnetoencephalography, MEG) remains an unresolved technical and conceptual challenge (Honey et al., 2009; Stephan et al., 2009; Pineda-Pardo et al., 2014; Uludag & Roebroeck, 2014). White matter (WM) pathways are sufficient for communication between brain regions, but functional brain dynamics can also be mediated through polysynaptic connections (Figure 1). Indeed, previous studies suggested the direct structural pathways inferred using dMRI account for only about 55% of measured resting-state functional connectivity patterns (Koch et al., 2002; Honey et al., 2009; Deligianni et al., 2011; Becker et al., 2016).

Current measures of anatomical and of resting-state functional connectivity are symmetric in the sense that they do not enable an assessment of whether one orientation of a pathway may be more prominent than the inverse (Friston, 2011). In contrast, models of effective connectivity such as dynamic causal models (DCMs) indicate the weights of specific directions of interaction (Friston et al., 2003), and recent data across species suggest that information about directed, asymmetric connectivity may more appropriately reflect brain architecture (Kale et al., 2018; Avena-Koenigsberger et al., 2019; Seguin et al., 2019).

Previous work has analysed the relationships between indirect anatomical connectivity and resting-state functional connectivity (Honey et al., 2007; Deligianni et al., 2011; Abdelnour et al., 2014; Becker et al., 2016; Meier et al., 2016; Bettinardi et al., 2017; Liang & Wang, 2017; Abdelnour et al., 2018). Recent graph-theoretic research has demonstrated that conventional, symmetric measures of brain WM architecture contain information on the differential efficiency of afferent and efferent network communication (Avena-Koenigsberger et al., 2019; Seguin et al., 2019). Furthermore, asymmetries in predicted communication efficiency were found to reflect neurobiological concepts of functional hierarchy and were correlated with directionality in resting-state effective connectivity analysed using spectral dynamic causal modelling (Seguin et al., 2019). Thus far, formal integration of effective with anatomical connectivity has only been implemented for direct and symmetric measures of structural connectivity (Stephan et al., 2009; Sokolov et al., 2018; Sokolov et al., 2019). The primary motivation for this study was to develop an integrative approach simulating symmetric and asymmetric high-order (polysynaptic) structural connectivity and using the outcomes to constrain models of task-related effective connectivity.

This central aim inspired the use of the graph Laplacian (GL: see Materials and Methods) to compute polysynaptic symmetric and asymmetric structural connectivity. The GL is a construct from spectral graph theory and represents the difference between the adjacency (indicating which network nodes are interconnected) and degree matrices (indicating the number of nodes connected with each node). The GL can be used to simulate the diffusion of a conserved quantity of particles over the network (Figure 2 and Supporting Information Video S1; Biggs, 1993). Crucially, the GL approach allows introducing asymmetry (see Materials and Methods and Supporting Information Figure S1): weighting (normalising) the structural adjacency matrix to the in-degree implies that each target node has a fixed capacity to be influenced by other nodes, and its relative sensitivity is determined by the probability of receiving inputs. Conversely, if we normalise to the out-degree, we assume each node has a fixed capacity to influence other nodes, and the relative influence is proportional to efferent particle diffusion.

In our implementation, the GL matrix L is exponentiated and raised to the order τ, exp(L)^τ. At the start of the diffusion process, each node is equipped with a large number of “spikes.” Each order τ represents a time step of the diffusion process, at which the spikes are distributed to other nodes at a rate that is proportional to the connection strengths. Each increment of τ thus indicates an extra path between any two given regions (nodes), via τ − 1 intermediate nodes (e.g., a third-order connection between two nodes means they are connected via two other nodes; Figure 1B). In a continuous time interpretation of this process, we can equate connectivity with the number of spikes accumulated as time progresses. This interpretation is effectively a constrained diffusion process, where the diffusion coefficients are determined by the GL and, ultimately, every node is connected to every other node through multiple paths. This equilibrium distribution is the principle eigenmode of the matrix exponential of the GL (Figure 2H), revealing which nodes are more strongly involved in the diffusion or propagation process.

Diffusion simulation approaches capture the propensity of information or particle distribution along all possible paths in a network, and thus approximate network communicability (Estrada & Hatano, 2008; Crofts & Higham, 2009). Network communicability has already been used to characterise brain networks in normalcy and pathology, and in different species (Crofts & Higham, 2009; Andreotti et al., 2014; Grayson et al., 2016; Shine et al., 2018). In addition to routing efficiency representing shortest paths and thus easy and speedy communication between network nodes, taking into account recurrent neuronal message passing over multiple paths may afford more optimal approximations of brain dynamics (Bullmore & Sporns, 2012; Avena-Koenigsberger et al., 2017). Using the matrix exponential of the GL allows a path length-based correction of network communicability (Estrada & Hatano, 2008; Crofts & Higham, 2009). Previous applications of the GL in neuroscience have suggested it as a promising tool for modelling neurotransmitter diffusion in the synaptic junction (Barreda & Zhou, 2011), simulating the spread of neurodegeneration (Raj et al., 2012) and comparing resting-state functional with structural connectivity (Abdelnour et al., 2014, 2018).

The present study asks whether measures of simulated symmetric and asymmetric anatomical network diffusion may usefully inform the effective connectivity that underwrites causal and asymmetric interactions among distributed neuronal populations. We demonstrate the approach in the context of fMRI responses of a brain network to emotional body language, using probabilistic tractography on high angular resolution diffusion imaging (HARDI) data from the same cohort of normal individuals.

We used fMRI and HARDI data from 17 right-handed, male normal subjects (mean age 27.9 years) from a study on emotional body language processing. The cohort overlapped with that analysed in previous research (Sokolov et al., 2012, 2014, 2018, 2019). The study was approved by the Ethics Committee of the University of Tübingen Medical School, Germany. Participants provided informed written consent and were financially compensated.

A 3T scanner (TimTrio, Siemens Medical Solutions, Erlangen, Germany; 12 channel head coil) was used for recording of three-dimensional T1-weighted structural MRI (magnetisation-prepared rapid gradient echo, MPRAGE; 176 sagittal slices, TR = 2,300 ms, TE = 2.92 ms, TI = 1,100 ms, voxel size = 1 × 1 × 1 mm³), a field map for inhomogeneity correction, HARDI data (two sessions with 64 diffusion gradient directions per subject; b-value = 2,600 s/mm², 54 axial slices, TR = 7,800 ms, TE = 108 ms, slice thickness = 2.5 mm, matrix size = 88 × 88, field of view = 216 mm) and functional echo-planar imaging (EPI; 171 volumes, 56 axial slices, TR = 4,000 ms, TE = 35 ms, in-plane resolution 2 × 2 mm², slice thickness = 2 mm, 1 mm gap).

Participants viewed animations of an arm represented by bright dots placed on the head and main upper limb joints, facing to the right and knocking on an invisible door with different emotional content (happy, angry, neutral; Pollick et al., 2001; Sokolov et al., 2011). In an event-related design, the participants had to indicate which emotion was expressed by button press (button assignment counterbalanced between participants). Stimulus duration was 1,000 ms, and each stimulus category (emotion) was presented 30 times throughout the experiment. To optimise event-related response function estimation, we applied jittering of stimulus onset intervals (between 4,000 and 8,000 ms in steps of 500 ms) and stimulus order pseudo-randomisation.

Preprocessing of fMRI data was performed using Statistical Parametric Mapping software (SPM12, Wellcome Centre for Human Neuroimaging, Institute of Neurology, UCL, http://www.fil.ion.ucl.ac.uk/spm) and included slice timing correction, realignment, unwarping, image co-registration, segmentation-based normalisation, and smoothing. HARDI data preprocessing with the FMRIB’s Diffusion Toolbox within the FMRIB Software Library (FSL5, Oxford Centre for Functional MRI of the Brain, UK, http://www.fmrib.ox.ac.uk/fsl) consisted of brain extraction (Smith, 2002), motion and eddy current correction, followed by co-registration with the anatomical reference image and normalisation to Montreal Neurological Institute (MNI) space using the FMRIB Linear Image Registration Tool (FLIRT; Jenkinson et al., 2002). Gradient directions were adjusted according to the FLIRT parameters.

Analysis of fMRI data was conducted by first specifying a general linear model (GLM). Trials with correctly classified emotional expression of point-light knocking (happy, angry, neutral) were assigned distinct regressors, and regressors of no interest were modelled for trials with incorrect classification (e.g., neutral stimulus classified as happy), trials with missing responses, six head motion parameters, and time series from WM and cerebrospinal fluid. The regressors were convolved with the haemodynamic response function. High-pass filtering was performed (cutoff 1/256 Hz), and serial autocorrelations were accounted for by a first-order autoregressive process (coefficient of 0.2) plus white noise model. The GLM was applied to individual preprocessed EPI data, and the contrasts happy versus neutral, angry versus neutral, and neutral versus emotional knocking were specified. Individual contrast images were submitted to second-level random effects analyses, and regional activations (at a p < 0.05 family-wise error corrected voxel-wise threshold for multiple comparisons) were identified using the automated anatomical labelling in SPM (Tzourio-Mazoyer et al., 2002) and the NeuroSynth.org database (http://neurosynth.org; Yarkoni et al., 2011).

A one-state, bilinear, and deterministic DCM with mean-centred inputs and reciprocal extrinsic connections between all nodes (full model) was created for each subject. This DCM included seven regions showing differential activation with respect to correctly classified emotional expressions of point-light knocking and three regions for activation versus baseline (Table 1). For each region, time series were extracted as the first eigenvariate of all activated voxels within a sphere with a radius of 8 mm, centred on each individual maximum (p < 0.05, uncorrected). The individual maxima were found within 7 mm of the group activation coordinate in every subject. The time series extraction was adjusted to remove effects that were not related to the task such as motion. According to previous data on the architecture of the brain network for body motion processing (Sokolov et al., 2018), driving input was specified on the left middle temporal cortex, right fusiform gyrus, and right superior temporal sulcus. Modulating input of different emotional content was expressed in the DCM B-matrix (see Supporting Information Methods) by modelling the influence of the corresponding regressors for happy, neutral, and angry stimuli on all extrinsic connections in the network, as well as on intrinsic coupling within the seven nodes showing differential activation depending on emotional content (Table 1).

Individual preprocessed and normalised HARDI data were submitted to Bayesian Estimation of Diffusion Parameters Obtained using Sampling Techniques with modelling of Crossing Fibres (BEDPOSTX; Behrens et al., 2007) in FSL to obtain diffusion parameters for each voxel. Probabilistic tractography with crossing fibres (PROBTRACKX; step length = 0.5 mm, number of steps = 2,000, number of pathways = 5,000, curvature threshold = 0.2, modified Euler integration; Behrens et al., 2007) was performed for each DCM node as a seed, and the other DCM nodes as classification targets. The fibre pathway outputs were visually controlled for plausibility. Structural connection strength between a seed region i and a classification target region j was obtained by averaging the number of streamlines connecting every voxel in i to one voxel of j, across both directions of tractography. Further averaging across all subjects provided a symmetric group structural adjacency matrix (Figure 2A). We eschewed thresholding and considered weighted adjacency matrices. Each element of the group adjacency matrix Z was normalised to represent direct structural connection strength or probability φ, relative to the greatest connection strength within the matrix. The between-region elements of matrix Z were used to inform models of effective connectivity by direct structural connectivity, and for GL-based simulation of network diffusion to obtain measures of indirect structural connectivity.

We used the GL matrix L to construct a connectivity operator simulating diffusion of a conserved quantity (heat, spikes) along direct and indirect pathways between the network nodes of the structural adjacency matrix Z. As per definition, each column of the GL matrix L expressing probabilities φ of extrinsic structural connections has to sum to zero, which consequently applies to any linear mixture of the columns of L.

After defining prior beliefs about the effective connectivity parameters, the estimation of DCMs affords posterior estimates of the parameters as well as the evidence for the respective model (Friston et al., 2003, Supporting Information Methods). The priors for extrinsic (off-diagonal; between-region) connections in dynamic causal modelling form a multivariate normal distribution, defined by a vector of expectations and a prior covariance matrix Σ_y. By default, the prior expectation is zero and the variance is equal across all extrinsic connections. The greater the prior variance, the further the connectivity parameters can deviate from their prior expectation of zero.

The priors also contribute to the calculation of the model evidence—the quantity used to compare DCMs—which is the trade-off between model accuracy and complexity. In this context, complexity is defined as the discrepancy between prior assumptions and posterior estimates, where greater complexity decreases model evidence. Optimising the priors according to measures of structural connectivity may therefore increase model evidence, through a reduction of model complexity (Stephan et al., 2009; Sokolov et al., 2019).

However, the precise relationship between structure and function is unknown and likely varies for different networks. Previous research provided support for the intuition that the strength of direct structural connections relates to the prior effective connectivity in a positive monotonic fashion (Stephan et al., 2009; Sokolov et al., 2018, 2019). This means that for lower structural connection strengths, the prior variance shrinks to a small value, precluding strong effective connectivity. Conversely, for greater structural connection strengths, our prior belief that the effective connectivity is close to zero can be relaxed by increasing the prior variance.

Crucially, both model spaces (for direct and indirect structural connectivity) include flat mappings (i.e., δ = 0 ), where structural constraints do not matter and Σ_{y red} is the same across all extrinsic connections, thus representing an intrinsic control (null hypothesis) for the assumption that structural constraints usefully shape effective connectivity.

As we analysed second-level measures of structural connectivity, we used parametric empirical Bayes (PEB) (Friston et al., 2015; Supporting Information Methods) to make inferences on effective connectivity at the group level. PEB is a hierarchical model, in which the average group connectivity acts as an empirical prior on individual connectivity (Friston et al., 2016). PEB estimation thus represents an iterative process between individual and group effective connectivity. Therefore, PEB properly partitions within- and between-subject random effects (Zeidman et al., 2019) and is robust to local minima problems (Friston et al., 2016). The individual DCMs were estimated using the PEB scheme with a prior variance of 0.5 for all extrinsic connections. Subsequently, the prior PEB variance of each extrinsic connection was adapted (reduced) using Equations 5 and 6 for measures of direct and indirect structural connection strengths, respectively.

The search for the models with the greatest evidence was afforded by Bayesian model reduction (BMR) (Friston et al., 2016; Supporting Information Methods). This recently introduced statistical device enables analytical evaluation of large model spaces in a matter of seconds, based on the estimation of a single, so-called “full” model. In contrast, the use of conventional dynamic causal modelling would have required separate estimation of each of the 1,941 alternative models (estimated processing duration: 1,300 days). By comparing the log evidences of the different models, we assessed whether effective connectivity was better explained by (1) indirect as opposed to direct measures of structural connectivity, (2) a particular order of the connectivity operator Γ^τ, and (3) a particular normalisation scheme of the structural adjacency matrix underlying Γ. Very strong evidence that one model provides a better account for the observed data than another is concluded from a relative log-model evidence of three (Penny, 2012), corresponding to a posterior probability of 95% or above.

We assessed the value of simulated indirect (high-order) anatomical connectivity afforded by network diffusion under the GL for sculpting effective connectivity, relative to models informed by direct structural connectivity and to DCMs without anatomical information.

For indirect structural connectivity, we used a sigmoid function (Equation 6) with the hyperparameters δ (slope; range 0–0.25 in seven equal steps) and Σ_{y min} (lower boundary on prior variance; range 0.0156–0.0625 in seven equal steps) to map the log-transformed group connectivity values provided at eight different, equally distributed orders τ of the diffusion process (from 1 to 64) onto prior variance of second-level effective connectivity.

A grid search over the 1,536 candidate models resulting from the diffusion process under three normalisation schemes (symmetric, weighted to the out-degree, weighted to the in-degree) using BMR (overall computation time 5.98 seconds) indicated the best constraints on effective connectivity were provided by the largest order (τ = 64) of a GL normalised to the in-degree. The structure-function mapping at this order was governed by the hyperparameters δ = 0.15 and Σ_{y min} = 0.05 (Figure 3).

This model clearly outperformed models informed by direct anatomical connectivity (log-evidence difference 33.13 in favour of indirect structural connectivity) and those without anatomical information (log-evidence difference 35.3 in favour of indirect structural connectivity). Very strong evidence that one model provides a better account for the observed data than another is concluded for a log-evidence difference of three or above (Penny, 2012).

As shown in Figure 3B, the evidence for DCMs of effective connectivity informed by indirect structural connectivity priors increased with progression of the particle diffusion process simulated by the GL, saturating at orders above τ = 50, corresponding to the principle eigenmode of the GL. This suggests the structural connectivity that matters for dynamical coupling and effective connectivity is best conceived in terms of reciprocal message passing over long (polysynaptic) paths, or periods of time.

Bayesian model comparison (Penny, 2012) across the three normalisation schemes provided consistently very strong evidence in favour of normalisation to the in-degree (log-evidence difference between this and the next probable normalisation scheme: 3.06; Figure 4). This result implied that effective connectivity is best predicted by the relative sensitivity of nodes to incoming information or afference and confirmed our hypothesis that introduction of asymmetry in the diffusion process may offer a more plausible characterisation of asymmetric brain dynamics.

When examining how the principle eigenmode of the GL normalised to the in-degree related to the functional afference of each node, one can see that the four nodes receiving the greatest input during the GL diffusion process (midcingular cortex (MCC), superior temporal sulcus (STS), middle temporal cortex (MTC) and insula (INS)) were also those with the highest functional in-degree based on effective connectivity (Figure 4). This further speaks to the utility and construct validity of the GL approach to inform effective by indirect structural connectivity.

Random permutations (n = 256) of the network nodes in the adjacency matrix (thereby preserving the distribution of edge weights) were used to assess how often a random structural adjacency matrix would afford a greater log-model evidence than reported above (after optimising the normalisation, order, and sigmoid hyperparameters). This permutation testing suggested the improvement in evidence afforded by applying the GL to the actual structural adjacency matrix was significant in a classical sense (p = 0.01), with respect to a null distribution of largest log-model evidences (Figure 5).

This study makes several contributions to the understanding of structure-function relationships in the brain. Based on previous research (Kale et al., 2018; Avena-Koenigsberger et al., 2019; Seguin et al., 2019), we hypothesised that asymmetric polysynaptic anatomical connectivity would better predict the directed causal dynamics between neuronal populations (effective connectivity) than conventional (i.e., symmetric and monosynaptic) information on WM pathways. The introduction of the GL allowed us to parameterise a diffusion process on the structural graph, providing symmetrically and asymmetrically weighted adjacency matrices of increasing order. Of note, other methods for modelling network communication based on structural connectomes can inherently inform on asymmetry. Such approaches include navigation (Seguin et al., 2018), search information (Goni et al., 2014), linear transmission models (Misic et al., 2015), and diffusion efficiency (Goni et al., 2013). The novel approach presented here enables hypotheses to be tested about the mapping from indirect anatomical connectivity to effective connectivity via a (variational) Bayesian framework. Using a dataset with fMRI and HARDI measures, we found that high-order structural connectivity simulated using the GL greatly improved the evidence for DCMs of effective connectivity, compared to DCMs informed by direct structural connectivity and anatomically uninformed models. Most importantly, input sensitivity during the diffusion process best predicted effective connectivity.

We introduced a computationally efficient approach to map indirect anatomical to effective connectivity, using hierarchical PEB models and BMR for DCMs (Friston et al., 2016; Zeidman et al., 2019). This procedure is designed to account inherently for possible variations in network architecture, normalisation scheme, value of anatomical information and mapping between structural and effective connectivity, for any given study and context. By definition, effective connectivity is determined by the context, such as the specific experimental task or cognitive set (Friston et al., 2003). For this reason, we would not expect a universally optimal set of priors on effective connectivity that could explain cognition per se. We sought to provide an efficient method for finding the best effective connectivity priors for any specific context, informed by indirect anatomical connectivity. An interesting future question will be whether the utility of the GL approach to assess indirect structural connectivity generalises to models of effective connectivity for resting-state data. The value of understanding asymmetric coupling between functionally related regions at rest, through combination of neuronal and observational models such as those used by dynamic causal modelling, has become increasingly recognised (Friston et al., 2014; Razi et al., 2015, 2017). Recent work has already linked asymmetries in indirect anatomical connectivity to resting-state effective connectivity (Seguin et al., 2019).

Using high-quality dMRI along with task-related and resting-state fMRI datasets such as from the Human Connectome Project (Van Essen et al., 2013) may contribute to further test and refine the outlined approach and conclusions, complementing previous research on these datasets (Seguin et al., 2019). Furthermore, it will be of interest to perform large-scale analyses integrating the structural connectome with whole-brain effective connectivity using the recently introduced regressive dynamic causal modelling (Frassle et al., 2018). Such approaches may also help to inform generative models of how WM pathways give rise to brain dynamics (Robinson, 2012; Deco et al., 2013; Sanz Leon et al., 2013; Melozzi et al., 2017; Messe et al., 2018), and to better understand how network dynamics may shape cognitive function and behaviour (Aertsen et al., 1994; Gerraty et al., 2018; Sokolov et al., 2018). Relating the diffusion properties of brain networks to their functions in extension of previous approaches based on direct anatomical connectivity (Deco et al., 2012; Senden et al., 2012; Hermundstad et al., 2014) may further improve our conceptualisation of distributed information processing in the brain.

Endowing the GL diffusion process with asymmetry clearly outperformed a symmetric diffusion process in sculpting effective connectivity, suggesting that ensemble dynamics in the brain may be shaped by the sensitivity of regions to their distributed input. These findings agree with and extend recent work on resting-state functional connectivity in macaques and humans showing that synchrony between nodes does not only depend on their direct or second-order connectivity, but also the similarity of afferents they receive from the entire network and other adjacent network characteristics (Adachi et al., 2012; Goni et al., 2014; Bettinardi et al., 2017). Consequently, densely connected regions may not necessarily be in a best position to influence or to be influenced by other regions (Avena-Koenigsberger et al., 2017). This follows from the fact that being locked into dense subgraphs may preclude a more widespread sensitivity to distributed dynamics (Pillai & Jirsa, 2017). Clarifying whether and why some networks may be better characterised in terms of their nodes’ sensitivity to inputs as opposed to their capacity to influence other nodes is a promising avenue for future research on normal and altered brain network function that can be pursued formally by the procedure described here.

The results presented here agree with and extend previous research employing network communication models (Avena-Koenigsberger et al., 2019; Seguin et al., 2019). These studies inferred the ease of sending and receiving information from undirected structural connectomes. The differences between send and receive efficiencies were mapped onto functional brain topography and the outcomes of separate resting-state effective connectivity analyses. For instance, anatomical connectivity-based classification suggested that unimodal areas such as the primary visual cortex or sensorimotor cortices were predominantly sending information, whereas multimodal regions were mainly receivers (Seguin et al., 2019). The diffusion efficiency approach, described as the (inverse) mean first passage time of a Markov chain process (Goni et al., 2013; Seguin et al., 2019), is similar to in-degree normalisation. In contrast to the present work, diffusion efficiency is derived using a matrix of transition probabilities. Nonetheless, measuring diffusion efficiency in a structural adjacency matrix yielded similar results, revealing regional variability in input sensitivity and a rather uniform capacity to influence other regions (Seguin et al., 2019). Future investigations are needed to fully explore the implications of the various measures to model diffusion processes.

The other important finding was that higher orders and the principle eigenmode of the GL afforded better priors on effective connectivity than lower GL orders. This indicated that WM connections and distributed neural dynamics give rise to brain communication through recurrent neuronal message passing over multiple paths. GL eigenmodes and eigenvalues are closely related to network communicability, representing the ease of information transmission along all possible paths in a network (Estrada & Hatano, 2008; Crofts & Higham, 2009; Andreotti et al., 2014; Grayson et al., 2016; Shine et al., 2018). GL eigenmodes of structural adjacency matrices exhibit a high degree of similarity between healthy subjects, as well as consistent and meaningful alterations in developmental and virtual agenesis of the corpus callosum (Wang et al., 2017). Laplacian eigenvalue spectra have been used for cross-species comparison of anatomical networks and revealed specific characteristics of neural networks as opposed to other network classes (de Lange et al., 2014). Furthermore, the anatomical graph energy (connectedness measure representing the sum of all absolute GL eigenvalues) has been shown to be significantly lower in patients with Alzheimer’s disease than in controls (Daianu et al., 2015). A greater number of apolipoprotein E4 gene copies predicted this reduction in graph energy. The use and interpretation of metrics such as eigenmodes and eigenvalues afforded by diffusion processes simulated by the exponentiation of a GL matrix could lead towards consideration of more global network characteristics beyond the conventionally assessed single hub or subgraph properties.

In clinical research, truly integrative computational, graph theoretic, or even simple correlative analyses of multimodal connectivity remain rather sparse. However, the assessment and comparison of network communicability using the GL may be of potential relevance to clinical neuroscience. Implementation of the GL in patients to assess how local and global changes in anatomical connectivity affect functional dynamics may shed new light on pathophysiology. Other comparative network measures afforded by the GL are topological similarity, persistent homology and graph diffusion distance (Hammond et al., 2013a; Bettinardi et al., 2017; Liang & Wang, 2017). Furthermore, simulation of network diffusion by means of the GL has been used to predict neuronal spreading and resulting brain atrophy patterns in Alzheimer’s and frontotemporal dementia (Raj et al., 2012) and to infer sources of disease propagation in mild cognitive impairment and Alzheimer’s dementia (Hu et al., 2016). Ultimately, the global measures afforded by the GL may be used towards assessment of the relationships between connectivity and behaviour at the network level (Sokolov et al., 2018). As efficiency and ease-of-use are of primary significance in everyday clinical practice, this relatively straightforward and rapid approach could potentially afford useful network biomarkers in neurological and psychiatric disorders.

Dynamic causal modelling for fMRI has contributed to establishing or refining various neuroanatomical and neurobiological concepts and hypotheses in functional realms such as reading, mental imagery, memory retrieval, or body language reading (Chow et al., 2008; Sokolov et al., 2012; Dijkstra et al., 2017; Ren et al., 2018; Sokolov et al., 2018). Inclusion of electrophysiological data from EEG, MEG, or intracranial recording (David et al., 2013; Almashaikhi et al., 2014; Proix et al., 2017), which enable more detailed biophysical modelling due to their high temporal resolution, may offer additional insight. Indeed, structure-function relationships appear to depend on different timescales of functional dynamics. Changes over periods of several minutes largely reflect underlying anatomical connectivity, whereas dynamics lasting a few seconds do so to a lesser degree (Honey et al., 2007; Shen et al., 2015; Cabral et al., 2017). Interestingly, resting-state functional connectivity data derived from simultaneous EEG and fMRI outperforms fMRI data alone in modelling structural connectivity, when comparing these predictions to actual dMRI measures (Wirsich et al., 2017). Network models based on anatomical connectivity and driven by EEG source activity accurately predict individual fMRI resting-state patterns and reproduce neurophysiological phenomena and mechanisms observed with different imaging modalities (Schirner et al., 2018). GL approaches to structural connectivity have also been used to improve EEG source localisation (Hammond et al., 2013b). Another potentially exciting development could be the computational modelling of sparse impulse stimulations as input in the GL diffusion process, similar to that implemented for modelling of disease propagation sources (Hu et al., 2016). This could further approximate the neuronal model within dynamic causal modelling (Friston et al., 2003) and information flow in the brain.

In summary, we have studied symmetric and asymmetric diffusion processes based on anatomical connectivity to optimise prior constraints on models of directed effective connectivity. Bayesian model comparison indicated the best effective connectivity priors are provided by indirect, high-order structural connectivity determined by the regional sensitivity to inputs that would be seen under equilibrium states of particle diffusion within an anatomical network. This may speak to a reappraisal of how we characterise the anatomical connectome, when trying to understand asymmetric functional dynamics arising from structure of the sort measured in neuroimaging.

Supporting Information for this article is available at https://doi.org/10.1162/netn_a_00150.

Arseny A. Sokolov: Conceptualization; Formal analysis; Funding acquisition; Investigation; Methodology; Writing - Original Draft. Peter Zeidman: Methodology; Supervision; Writing - Review & Editing. Adeel Razi: Conceptualization; Methodology; Writing - Review & Editing. Michael Erb: Methodology; Supervision. Philippe Ryvlin: Supervision; Writing - Review & Editing. Marina A. Pavlova: Conceptualization; Funding acquisition; Methodology; Resources; Supervision; Writing - Review & Editing. Karl J. Friston: Conceptualization; Formal analysis; Funding acquisition; Methodology; Resources; Supervision; Writing - Review & Editing.

Arseny A. Sokolov, Baasch-Medicus Foundation. Arseny A. Sokolov, Fondation Leenaards (http://dx.doi.org/10.13039/501100006387). Arseny A. Sokolov, Schweizerische Neurologische Gesellschaft (http://dx.doi.org/10.13039/100010766). Arseny A. Sokolov, Helmut Horten Foundation. Arseny A. Sokolov, Synapsis Foundation Alzheimer Research Switzerland, Award ID: 2019-CDA03. Marina A. Pavlova, Reinhold Beitlich Stiftung (http://dx.doi.org/10.13039/501100003541). Marina A. Pavlova, BBBank Foundation. Marina A. Pavlova, Deutsche Forschungsgemeinschaft (http://dx.doi.org/10.13039/501100001659), Award ID: DFG PA 847/22-1. Karl J. Friston, Wellcome Trust (http://dx.doi.org/10.13039/100004440), Award ID: 088130/Z/09/Z.

The authors wish to thank Richard Frackowiak, Patric Hagmann, Alexander Sokolov, and Klaas-Enno Stephan for valuable discussion. Technical support was provided by Ric Davis, Jürgen Dax, Chris Freemantle, Bernd Kardatzki, Rachael Maddock and Liam Reilly, and administrative assistance by Marcia Bennett, David Blundred, Kamlyn Ramkissoon, Tracy Skinner, and Daniela Warr.

Effective connectivity:

A measure of the directed (causal) influence of one neural system over another using a model of neuronal interactions.

Dynamic causal modelling:

A Bayesian framework which is used to infer causal interactions between coupled or distributed neuronal systems (effective connectivity).

High-order structural connectivity:

Structurally unconnected regions communicate polysynaptically to engender indirect connectivity over multiple hops.

Graph Laplacian:

A matrix representation of a graph that combines node adjacency and node degree in mathematical formulation and belongs to spectral graph theory.

Spectral graph theory:

A study of the relationship between a graph and the eigenvalues and eigenvectors of its Laplacian matrix.

Adjacency matrix:

Square matrix representation of a graph which is either binary (presence or absence of connections) or weighted (showing strength of connections).

Eigenmode:

A stable state (i.e., mode) of a dynamic system in which all parts of the system oscillate at the same frequency.

Network communicability:

A graph-theoretic measure of the ease of information propagation in a network of interconnected regions.

Routing efficiency:

A measure of communication efficiency, representing the average inverse shortest path length between all pairs of nodes in a complex network.

Network diffusion:

A process simulating the propagation of heat (or information) within a network.

Probabilistic tractography:

Estimation of white matter pathway trajectories based on diffusion MRI data.

Model evidence:

The model evidence, or marginal likelihood, represents the probability of observing the measured data under a certain model and is used for Bayesian model comparisons.