Computational modelling reveals neurobiological contributions to static and dynamic functional connectivity patterns
Linnea Hoheisel, Hannah Hacker, Gereon R Fink, Silvia Daun, Joseph Kambeitz
TL;DR
This study uses computational models to show how brain connectivity patterns are shaped by neurobiological factors and specific brain regions.
Contribution
The study links neurobiological parameters to both static and dynamic functional connectivity using computational modeling in healthy individuals.
Findings
Model parameters like G and JN showed strong correlations with static and dynamic FC features.
TC fit was most influenced by coupling changes in the left pars triangularis, while sFC fit was affected by the left paracentral gyrus.
sFC and dFC features were predictive of their respective model fits with high accuracy (R2 > 0.5).
Abstract
Functional connectivity (FC) is a widely used indicator of brain function in health and disease, yet its neurobiological underpinnings still need to be firmly established. Recent advances in computational modelling allow us to investigate the relationship of both static FC (sFC) and dynamic FC (dFC) with neurobiology non-invasively. In this study, we modelled the brain activity of 200 healthy individuals based on empirical resting-state functional magnetic resonance imaging (fMRI) and diffusion tensor imaging (DTI) data. Simulations were conducted using a group-averaged structural connectome and four parameters guiding global integration and local excitation-inhibition balance: (i) G, a global coupling scaling parameter; (ii) Ji, an inhibitory coupling parameter; (iii) JN, the excitatory NMDA synaptic coupling parameter; and (iv) wp, the excitatory population recurrence weight. For…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
FIGURE 1
FIGURE 2
FIGURE 3
FIGURE 4
FIGURE 5Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsFunctional Brain Connectivity Studies · Neural dynamics and brain function · Complex Network Analysis Techniques
Introduction
Resting-state functional connectivity (FC), defined as the correlation of activity between different brain regions at rest, is a key signature of brain functioning (Finn et al., 2015; Ooi et al., 2022; van den Heuvel and Hulshoff Pol, 2010). It is related to cognitive ability and personality (Kong et al., 2019; Sripada et al., 2020) as well as the brain’s health status and individual symptoms (Diener et al., 2012; Gallo et al., 2023; Mehta et al., 2021; Parkes et al., 2020; Wang et al., 2020). Recent research has indicated that FC is not static but changes over time, exhibiting dynamic characteristics related to behavior and clinical features (Hutchison et al., 2013; Patanaik et al., 2018; Hoheisel et al., 2024a). However, the biological processes underlying both static and dynamic FC remain unclear. Novel computational approaches allow us to investigate these mechanisms by simulating brain activity using empirical neuroimaging data (Wang et al., 2024), revealing new insights into communication between brain regions and its disturbances in brain disorders (Popovych et al., 2018). As previous models primarily focused on replicating static FC, they were not suitable for the investigation of dynamic aspects of FC, which are crucial for understanding brain connectivity and related disorders. In the present work, we utilize brain network models to elucidate the neurobiological underpinnings of FC. We incorporate both static and dynamic FC into model fitting, providing a comprehensive picture of the relationship between structural connectivity (SC), neurophysiology, and FC with the aim of improving our understanding of how neurobiological processes contribute to individual differences in brain network architecture and function.
Individual brain FC patterns are influenced by several mechanisms. The role of SC which represents the pattern of anatomical links between brain regions, in static FC origin has long been established (Honey et al., 2009; Sporns et al., 2000). However, SC does not account for the entirety of static FC variability between subjects, and cannot account for patterns of dynamic changes in FC (Liégeois et al., 2020). Mechanisms of neurotransmission and neuromodulation also regulate communication between regions (Brink et al., 2016; Hansen et al., 2022b) through factors such as dopaminergic and serotonergic signaling (Klaassens et al., 2015), and neuroreceptor expression patterns (Hansen et al., 2022a). In addition, neurobiological characteristics of brain areas which determine regional activity, especially the balance between excitatory and inhibitory populations, impact FC (Gu et al., 2019; Kapogiannis et al., 2013; Levar et al., 2019). Excitatory and inhibitory dynamics within brain regions, which are understood to be governed by excitatory and inhibitory synaptic coupling and excitatory self-excitation, shape regional activity as well as long-range integration (Deco et al., 2014). In line with this, disruptions of N-methyl-D-aspartate (NMDA)-related excitatory signaling have been suggested as a possible mechanism leading to increased connectivity across the brain (Anticevic et al., 2015; Driesen et al., 2013) and increased network flexibility in dynamic FC (Braun et al., 2016).
Brain network modelling allows for a holistic analysis of SC and FC by inferring neurobiological processes from data acquired via magnetic resonance imaging (MRI) (Breakspear, 2017; Schirner et al., 2018; Shine et al., 2021). These models describe regional brain activity using a set of equations with variable parameters governing excitation-inhibition balance and the integration of long-range input from other regions, while empirical SC determines the strength and timing of signals exchanged between regions (Deco et al., 2014, 2017; Sanz-Leon et al., 2015). By determining which optimal parameters best replicate the empirical brain activity of each individual, we can explore systems underlying individual variability in brain activity and common processes shaping network dynamics in health and disease (Klein et al., 2021; Zimmermann et al., 2018; Popovych et al., 2018).
While global neurobiological characteristics shape individual FC patterns, brain disorders are often related to changes in the properties of only a few areas (Fornito et al., 2015). By exploring which regions play an essential part in maintaining healthy FC, we can discover the neurobiological underpinnings of brain network architecture, as well as potential mechanisms of brain disorders. Previous research suggests that brain networks rely on a few highly connected regions (hub nodes) that significantly affect overall FC (Bullmore and Sporns, 2009; Fransson and Thompson, 2020; Zamani Esfahlani et al., 2020). The strength of these hub nodes is linked to gene expression profiles (Vértes et al., 2016), and their dysfunction is implicated in several brain disorders (Crossley et al., 2014; Royer et al., 2022; Rubinov and Bullmore, 2013). Brain network modelling allows for the precise manipulations of regional neurobiological processes to explore how changes in specific areas affect both static and dynamic FC (Aerts et al., 2016; Alstott et al., 2009).
Traditionally, brain network modelling studies have attempted to replicate the empirical static FC of an individual (Deco et al., 2014; Domhof et al., 2021). While static FC can robustly identify individuals, it cannot capture a wealth of information linked to the temporal evolution of the connectivity pattern (Calhoun et al., 2014; Chiang et al., 2018; Hutchison et al., 2013). In addition to the structure of network connectivity, dFC captures the evolution of connectivity over time, potentially providing evidence of dynamic physiological and cognitive processes. dFC complements and in some cases exceeds sFC as an indicator of behavioral variability, and allows for the consideration of time-varying mechanisms that are averaged out in traditional sFC analysis (Lurie et al., 2020). Research has shown that dynamic network characteristics of FC also vary between individuals (Chen et al., 2016; Davison et al., 2016) and exhibit diagnosis- and symptom-related alterations in brain activity of patients that are not apparent from static FC alone (Pang et al., 2022; White and Calhoun, 2019). These findings suggest that dynamic FC should be considered in investigations of brain alterations.
Here, we present findings from a computational modelling investigation considering dynamic FC. We determined the optimal values of parameters representing excitation-inhibition balance and global integration that best reproduced empirical static and dynamic FC, as well as the optimal correlation between empirical and simulated static and dynamic FC that could be achieved for each individual. We analysed the associations of fits and parameters with static and dynamic FC features across the brain in order to discover which patterns of FC determine model fits and how neurobiological parameters contribute to individual variability in regional static and dynamic FC. In addition, we investigated the role of different brain regions in generating brain network dynamics by systematically altering the coupling in each region and measuring the resulting changes in static and dynamic FC. Our findings provide new evidence for the role of the interaction of structural connectivity and neurophysiological mechanisms modulating communication between brain regions in the origin of static and dynamic FC. By elucidating this pathway, these findings highlight processes that, when disturbed, might contribute to impairments in brain disorders.
Materials and methods
Data acquisition and preprocessing
We modelled individual brain activity using empirical structural and functional MRI data from 200 healthy, not related subjects from the Human Connectome Project (HCP) S1200 release (Van Essen et al., 2013). The HCP study was approved by the institutional review board of Washington University and written informed consent was given by all participants. Empirical SCs and parcellated rs-fMRI time courses were generated by Domhof et al. (2021, 2022b,2022a). We used the Desikan-Killiany cortical parcellation (Desikan et al., 2006) for our analyses, which delineates 70 regions of interest based on anatomical structures.
Computational modelling
Whole-brain network models describe the brain as a set of regions joined by large-scale connections representing white matter tracts (Sanz-Leon et al., 2015). Each brain area is described as a network of neuronal populations, with a set of parameters regulating the balance between them. Information is passed from each region to regions with which it exhibits empirical connections, scaled by the relevant tract weights. Tract lengths determine the delay with which the information arrives at the target region. In this study, we used the reduced Wong-Wang model with excitatory and inhibitory components (Deco et al., 2014) to simulate the regional activity, which was then converted to a simulated fMRI signal. This dynamic mean field model describes the change of the average firing rates and synaptic activities of a population of excitatory and a population of inhibitory neurons with a set of coupled non-linear differential equations:
r_i_^(E,I)^ represents the firing rate in excitatory (E) or inhibitory (I) population of node I, S_i_^(E,I)^ the average synaptic gating variable, and I_i_^(E,I)^ the input current. G denotes the global coupling, (Cij) the structural connectivity between each pair of regions i and j. Ji represents the local feedback inhibitory synaptic coupling, I0 the external input, w(E,I) the population external input scaling weight, wp the local excitatory recurrence, and JN the excitatory synaptic coupling.
Three forms of input drive the regional activity at each point in time: (i) long-range excitatory inputs from other regions, modulated by the structural connectivity, the global coupling parameter G, and the local excitatory synaptic coupling J_N, (ii) regional inhibitory currents, modulated by the local feedback inhibitory synaptic coupling Ji, and (iii) recurrent excitation, modulated by the local excitatory recurrence wp_ and the local excitatory synaptic coupling J_N. We investigated simulated data based on a range of values for G, JN, Ji, and wp_.
A graphical representation of the modelling approach is shown in Figure 1. We performed all simulations using an average empirical SC, represented by the mean white matter tract weights and lengths across all subjects in the HCP sample. We transformed the two matrices by removing the median and scaling to the interquartile range to obtain a scaling robust to outliers, and rounded to even numbers. We relied on the implementation of the model in the Virtual Brain toolbox (Sanz Leon et al., 2013) to perform the simulations. All simulation settings can be found in Supplementary Table 1.
Computational modelling workflow. We produced a range of simulated fMRI time courses based on empirical structural connectivity (SC) and regional models of activity. Each regional model describes the interaction between an excitatory and an inhibitory population and outside input. The excitatory population receives input from all other regions in the brain scaled by the empirical SC and the global coupling parameter G and the excitatory synaptic coupling parameter JN. It transfers information back to itself in the form of recurrent excitation, scaled by the excitatory recurrence parameter wp, and forward to the inhibitory population, scaled by JN. The inhibitory population returns an inhibitory current, scaled by the inhibitory synaptic coupling parameter Ji. The resulting simulated brain activity is translated to a simulated fMRI time course, which is compared with the empirical fMRI time course of each individual using the static functional connectivity (FC) and three dynamic FC metrics, the FC variance, the temporal correlation, and the node cohesion to assess the fit between empirical and simulated data. This process is repeated for a range of values for the four parameters G, JN, Ji, wp, to identify the combination of parameters producing the optimal fit, which is then analysed further.
In order to reduce the parameter space, we performed an initial exploration of parameters to identify a range for J_N, Ji, and wp_ at which the model exhibited multistable behavior for 101 coupling values equally distributed between 0 and 10. In this interval, the model oscillates between two stable states, representing biologically plausible activity. The procedure for this analysis is outlined in the Supplementary Information. We then generated simulations for each combination of 10 values for J_N, Ji, and wp_ in the multistability range, amounting to 1000 simulations per coupling value, for a total of 101000 simulations. We obtained a simulated blood oxygen level-dependent (BOLD) fMRI signal in each region by transforming the excitatory synaptic activity using the Balloon-Windkessel hemodynamic model (Buxton et al., 1998; Friston et al., 2003; Mandeville et al., 1999). We then determined the optimal parameter combination by comparing the FC of these simulated time courses with the empirical FC of each subject.
Calculation of FC measures
We evaluated the models based on both static and dynamic FC. The static FC (sFC) denotes the correlation in the activity of each pair of regions over the course of the scan, providing a metric for the average of the connectivity over time. We selected three established metrics which capture important and distinct aspects of dynamic FC (Barber et al., 2021; Sizemore and Bassett, 2018), the FC variance (FCV), the temporal correlation (TC), and the node cohesion (NC). FCV is the change in this correlation over time. This measure captures differences in the temporal variability between connections. TC is the consistency of a region’s neighbors from one time point to the next, indicating whether changes in connectivity are abrupt or more gradual. NC is the number of times each pair of regions switches communities together. This measure reveals whether nodes are likely to be involved in the same networks and processes. The static FC was determined by calculating the Pearson correlation between the time courses of each pair of regions, yielding a 70 × 70 matrix for each subject. We used Fisher’s z-transformation to normalize these matrices. In order to compute the dynamic FC parameters, we first calculated the FC in windows of approximately 60 s length, overlapping by 2 s and convolved with a Gaussian kernel of σ = 6 s, producing a 70 × 70 × 420 matrix for each subject (Leonardi and Van De Ville, 2015). We then calculated the variance in the connectivity of each pair of regions over time. In addition, we transformed this time course of FC matrices into a dynamic graph by binarizing it, keeping only the top 10% of connectivities. From this graph, we first derived the temporal correlation, resulting in a vector of 70 elements. Then, we detected communities of recurrently connected nodes using the tnetwork python library (Cazabet et al., 2021, 2023) and computed the node cohesion, yielding again a 70 × 70 matrix. The HCP sample contains two resting-state fMRI scans, one recorded using left-right and the other right-left phase encoding, in order to enable researchers to reduce phase encoding-related artefacts by averaging (Van Essen et al., 2013). We computed the dynamic measures separately for each subject, and used the mean of the two matrices or vectors for further analysis.
We calculated the Pearson correlation of these metrics between the simulated and the empirical data. For the measures represented as matrices, we compared only the upper triangulars of the simulated and empirical FC. In order to produce models that reliably reproduced both static and dynamic FC, we identified the model that provided an optimal fit according to both the sFC and the TC using the l2-norm as a global criterion, here referred to as the “combined” metric. In addition, we evaluated the models that produced the highest correlation on each of the individual metrics for comparison. For each of these five optimal models of each subject, we computed the fit across all four of the measures, as well as the corresponding parameter set. We used the optimal model according to the combined metric for all further analyses. We repeated the calculation of optimal parameters and fits and all follow-up analyses in the rs-fMRI data from the second scanning session available in the HCP data set to validate our findings.
Correlation of fits and parameters with brain-derived features
Based on the optimal models, we considered the relationships between connectivity features and parameters as well as fits. This analysis allows us to discover whether model fits are determined by FC patterns and whether individual variability in FC is related to differences in model parameters. In order to investigate the association between fits and optimal parameters and static and dynamic FC, we computed Pearson correlations of each of the four parameters and the fit according to each of the four measures with a range of features over all 200 subjects of the HCP data. Those features included (a) the sFC between each of the 2415 pairs of regions and (b) the TC in each of the 70 regions. We performed permutation testing to estimate the significance of the resulting values. A null distribution for each correlation was produced by randomly permuting the parameter or fit variable over the subjects 100000 times and recomputing the correlation. All p-values were adjusted for multiple comparisons using the false discovery rate (FDR) unless otherwise specified. To enhance visualizations, we additionally summarized sFC and TC correlations over seven canonical resting-state networks (Yeo et al., 2011).
To add to this analysis of univariate relationships between model outcomes and individual connectivities, we also attempted to detect multivariate patterns of connectivity features across regions associated with model outcomes. We produced 16 ridge regression models that aimed to predict each parameter or fit value. These models considered either the sFC of the 2415 pairs of regions or the TC in each of the 70 regions as features. For the sFC-based model, we reduced the number of features by performing a principal component analysis (PCA) (Maćkiewicz and Ratajczak, 1993), and selected the minimal number of components that explained more than 90 % of the variance in the data. Additionally, we transformed the smaller number of features for the TC-based model using a PCA, keeping all components. We employed a nested cross-validation (CV) approach to determine generalization performance (Varma and Simon, 2006). The inner CV cycle determined the optimal alpha value from a range of 100 values logarithmically distributed between 0.01 and 100. In contrast, the outer CV cycle estimated the mean R^2^ score of the prediction on previously unseen data. Each level of the CV consisted of 5 folds with 2 permutations, with the PCA transformations estimated individually for each training and test set. In order to determine which regions contributed most strongly to the prediction, we computed the feature importance of each connectivity (Breiman, 2001). We re-trained the model on the entire data set, and calculated the importance of each feature based on the mean drop in the R^2^ score over 1000 permutations of that feature. For the sFC-based models, we computed the mean of the importance of the connections of each region to obtain a single importance value for each region. To enhance visualizations, we additionally summarized feature importance over resting-state networks (RSNs).
Region-wise perturbation of parameters
We used a perturbation approach to investigate how FC responds to changes in regional coupling. We systematically altered G in each region and analysed which regional coupling values most strongly impacted FC, and which were particularly important for maintaining biologically accurate FC. We determined a default set of parameters that produce the optimal fit on the combined metric between the simulated data and an empirical sFC matrix and TC vector averaged over the 200 subjects in our sample. We then produced a set of simulations using these default parameters, varying G in each region using 101 values uniformly distributed from 0 to 10. For each of the resulting time series, we computed the sFC, and obtained the global efficiency for each perturbation-derived and the default simulated sFC. The efficiency represents the average inverse shortest path length between each pair of regions, and measures how efficiently information is exchanged across brain networks (Latora and Marchiori, 2001). We determined the difference in the sFC and TC fit of simulated and averaged empirical data and the difference in the efficiency between the result obtained by simulating with the default parameter values and with the perturbation in each region. In order to determine why certain regions might contribute more strongly to global FC than others, we investigated the relationship between perturbation-related differences and graph properties, dynamic properties and gene expression across the brain. We determined correlations between the mean and variance of the sFC fit difference, TC fit difference, and efficiency difference and three sets of additional features: (i) graph metrics including centrality and degree measures of the empirical SC weights and lengths, and the average empirical FC; (ii) dynamic metrics of empirical regional BOLD time courses (Lubba et al., 2019; Supplementary Table 2); and (iii) expression data of 15653 genes extracted from the Allen brain atlas (Hawrylycz et al., 2012; Markello et al., 2021). We determined significance levels for each correlation based on nulls generated by permuting the feature importance matrices while preserving their spatial autocorrelation using spin-testing (Alexander-Bloch et al., 2018). These steps were performed with the neuromaps python toolbox (Markello et al., 2022). To facilitate the interpretation of the differential gene expression, we used the MetaScape platform (Zhou et al., 2019) to functionally enrich the genes whose expression correlated significantly (p_uncorrected_ < 0.05) with the perturbation-induced differences. The gene lists relating to each of the six difference metrics were compared to gene sets included in each Gene Ontology (GO) (GO Consortium, 2004) biological process term, and significantly overrepresented gene sets were extracted. Redundant terms were collapsed into one representative term via clustering.
Results
Simulation outcomes
While we were able to achieve high fits between empirical and simulated sFC (r = 0.389, σ = 0.037) and TC (r = 0.388, σ = 0.046), the FC variance (r = 0.082, σ = 0.032) and NC (r = 0.196, σ = 0.026) could not be replicated well (Figure 2). Simulated sFC matrices exhibited strong overlap with empirical sFC matrices (Supplementary Figure 1). Using the combined metric of sFC and TC for optimisation, sFC fits averaged 0.339 (σ = 0.044) and TC fits 0.378 (σ = 0.044). sFC fit and FCV fit were significantly higher in female than male participants (p(sFC) < 0.001, p(FCV) = 0.003), while NC fit was significantly lower in the 22–25 than in the 26–30 age group (p = 0.039). There were no other significant differences in any fits or parameters between any groups.
Optimal fits and parameters determined by comparing simulated and empirical FC. (A) Correlations between simulated and empirical data of optimal solutions based on each of the five metrics considered. (B) Parameters and fits split by covariate values. Stars indicate significant differences (p < 0.05).
Association of optimal parameters and fits with FC
The optimal values for the model parameters G, representing global coupling, and J_N, representing excitatory synaptic coupling, obtained for each subject when considering both static and dynamic FC showed a significant correlation with FC between and TC within several regions (Figure 3i). The sFC of multiple regions showed a strong positive correlation with the coupling parameter G but a negative correlation with the excitatory synaptic coupling JN. G exhibited particularly high correlations with connections within the dorsal attention network (DAN) (mean r = 0.22), as well as between the somatomotor network (SMN) and the DAN (mean r = 0.12) and visual network (VN) (mean r = 0.13). For JN, the strongest associations were present in connections within the VN (mean r = −0.13) and those between the DMN and the SMN and VN (mean r = −0.06). The TC correlated positively with G (mean r > 0.19), particularly in the DMN (mean r = 0.10) and DAN (mean r = 0.11), and negatively with JN_ across all regions (mean r < −0.16). There were no significant associations between J_i, representing inhibitory synaptic coupling, or wp, representing excitatory recurrence, and sFC or TC in any regions. These findings indicate that G and J_N contribute individual variations in specific static and dynamic FC patterns, while Ji_ and w_p_ do not.
Association between model outcomes and static and dynamic FC features. (A) Correlation between TC and (i) optimal model parameters and (ii) optimal fits. Only significant correlations (pFDR < 0.05) are displayed. (B) Mean correlation between TC and (i) parameters and (ii) fits in each resting-state network. Non-significant correlations (pFDR > = 0.05) were set to 0. (C) Correlation between sFC and (i) parameters, and (ii) fits. Only the strongest 30% of significant correlations (pFDR < 0.05) are displayed. (D) Mean correlation between sFC and (i) parameters and (ii) fits within and between resting-state networks. Non-significant correlations (pFDR > = 0.05) were set to 0. J_i, w_p and NC fit did not exhibit significant correlations with the FC of any region. DAN, dorsal attention network; DMN, default mode network; FPN, frontoparietal network; LN, limbic network; SMN, somatomotor network; VAN, ventral attention network; VN, visual network.
The fit between simulated and empirical data according to the four measures was also strongly connected to some TC and sFC features (Figure 3ii). Subjects with strong sFC also showed high correspondence between empirical and simulated data when considering either sFC or FCV as a target measure. Functional connections between the FPN and the SMN were mainly associated with a strong fit according to sFC (mean r = 0.26). Functional connections between the default mode network (DMN) and the SMN were associated with a strong fit according to FC variance (mean r = 0.56). TC fit exhibited negative correlations with most sFC features, involving the DAN, SMN and VN most strongly (mean r < −0.13). The sFC fit (mean r > 0.15) and especially the FC variance fit (mean r > 0.47) exhibited strong positive, the TC fit (mean r < −0.06) negative associations with the TC in most regions. The NC fit did not correlate significantly with sFC or TC in any regions, showing that sFC fit, FC variance fit, and TC fit, but not NC fit, are determined by individual static and dynamic FC patterns.
Prediction of optimal fits and parameters from FC
Multivariate analysis showed that temporal correlation and static connectivity features could predict some model fits reasonably well. Both models were strongly predictive of the correlation between empirical and simulated FC variance (Figure 4A, R^2^(sFC to FCV) = 0.44, R^2^(TC to FCV) = 0.53). The fit according to sFC and TC could be predicted with a high score by the model using sFC and TC as features respectively (R^2^(sFC to sFC) = 0.72, R^2^(TC to FCV) = 0.68). Neither of the models performed particularly well in predicting the fit according to NC (R^2^ < 0). The model parameters could be best predicted using the TC in each region as features. For J_i_ and w_p, the achieved scores were close to 0, the R^2^ value expected from a constant model which predicts the mean of the target vector for each subject regardless of feature input. The score for G was somewhat below that, while the score for JN_ reached a moderate value (R^2^(TC to J_N_) = 0.12). The sFC-based prediction model performed worse than a constant model for all parameter targets.
Prediction of fits and parameters from static FC and TC. (A) Mean R2 score achieved by models predicting each fit and parameter value based on either the sFC matrix or the TC vector in the outer CV. Error bars represent the standard error of the mean across folds and permutations in the CV. An R2 of 0 indicates that the model performed as well as a constant model which always predicts the average target value regardless of features. (B) Feature importance of TC values of each region. Only the top 30% most important regions are displayed. (C) Mean feature importance of TC values in each RSN. (D) Feature importance of sFC values of each pair of regions. Only the top 5% most important connections are displayed. (E) Feature importance of sFC values within and between RSNs. Only predictions with an R2 score above 0.25 are featured in panels (B–E). DAN, dorsal attention network; DMN, default mode network; FPN, frontoparietal network; LN, limbic network; SMN, somatomotor network; VAN, ventral attention network; VN, visual network.
For those prediction models with a high R^2^ score, feature importance analysis showed connectivities that were particularly predictive of the targets (Figures 4B–E). The fit measured by the FC variance was associated most strongly with the TC in the left temporal pole and the left isthmus cingulate, as well as a distributed network of static connections. Regions within the VN, SMN and between the VN and DAN were particularly important for the sFC fit. In contrast, the TC in the LN and VAN, especially in the left medial orbitofrontal gyrus, left pars opercularis, and right fusiform gyrus, had a strong effect on the TC fit.
Perturbation-induced alterations in FC
Perturbation analysis revealed that perturbing the coupling in most regions led to a decrease in fits and efficiency (Figure 5A). However, some regions had a stronger influence on the overall network connectivity than others, while altering the coupling in some regions slightly increased fits and efficiency. The alteration in the efficiency and the sFC fit over all the selected G values was particularly pronounced in the left and right paracentral, the right pre- and postcentral and the right transverse temporal gyrus (Δr > 0.03, Δeff > 0.01). The sFC fit could be improved by perturbations in the right medial orbitofrontal gyrus (Δr = −0.01). The effect on the TC fit was much stronger than on the sFC fit, with perturbation in the left pars triangularis in particular causing a reduction (Δr = 0.24), and in the left inferior temporal gyrus an improvement in fit (Δr = −0.01). In some regions, the alterations in fits and efficiency varied strongly depending on the extent of the perturbation. A change in the coupling of the left paracentral gyrus produced a particularly high variance in sFC fit and efficiency differences (σ^2^(Δr) = 0.001), while the TC fit varied most strongly when the coupling was perturbed in the left pars triangularis and left pars orbitalis (σ^2^(Δr) > 0.008). Over all regions, the fit generally decreased the further G was altered from the benchmark (Supplementary Figure 3).
Effects of regional coupling perturbation. (A) Alterations in (i) sFC fit, (ii) TC fit, and (iii) network efficiency produced by varying the coupling parameter G in each region. Node colors represent the mean alteration. Node sizes represent the variance of the alteration over all G values used. (B) Correlation of alteration mean and variance of perturbation-induced differences with a range of graph metrics. Orange bars indicate that features did not exhibit significant correlations, red bars show that they correlated significantly before FDR correction (puncorrected < 0.05). (C) Top 10 biological processes involving genes which correlated significantly (p < 0.05) in their expression with mean and variance of perturbation-induced differences. bc, betweenness centrality; cc, closeness centrality.
All six outcome metrics of the perturbation analysis, the mean and variance of the sFC fit difference, the TC fit difference and the efficiency difference, correlated highly with graph parameters of empirical structural and functional connectivity (Figure 5B). The association with the betweenness centrality and closeness centrality of the empirical SC weights were strongest for the means and variances of the sFC fit difference and the efficiency difference (r > 0.44), with the variance of the efficiency differences also correlating strongly with the degree of the SC weights (r = 0.51). The mean of the TC fit difference correlated most strongly with the closeness centrality of the SC weights (r = 0.35) and the degree of the SC lengths (r = 0.38), while the variance correlated most strongly with the degree (r = −0.43) and closeness centrality of the empirical sFC (r = −0.42). Regions which were connected more strongly within the network also led to greater changes in global static and dynamic FC if their coupling was altered. Correlations with dynamic parameters of the regional time courses were weaker and generally did not survive correction for multiple comparisons (Supplementary Figure 4). All six metrics exhibited strong correlations with the expressions of several genes. Enrichment analysis revealed several biological processes involving relevant genes (Figure 5C). Similar pathways were associated with the mean of sFC fit, TC fit and efficiency differences, as well as the variance of sFC fit and efficiency differences, and included brain-related processes such as nervous system development, modulation of chemical synapse transmission, and behaviour. The variance of the TC fit differences was only associated with a subset of the pathways. This structure of association indicates that not all biological processes related to those regional couplings which most impact static FC were associated with the same variability in dynamic FC, but they were associated with regions vital to maintaining biologically accurate global FC.
Discussion
In this study, we investigated individual neurobiology through brain models which simulate static and dynamic functional connectivity. We then analysed the relationship between optimal model parameters and regional static and dynamic FC features, and determined the effect of each region on overall network connectivity.
We achieved reasonably high maximal fits between simulated and empirical sFC as well as TC for each individual. Given that models optimized for sFC did not manage to reliably replicate empirical TC and vice versa, it was advisable to optimize models for both metrics. When we selected optimal individual models based on a combination of both factors rather than each factor individually, the obtained fits were on par with those reported previously (Klein et al., 2021; Zimmermann et al., 2018). Our models could not replicate empirical NC or FC variance well. The fits of the optimal models were strongly correlated with a set of connectivity features. sFC and FC variance fit showed significant positive correlations with the connectivity in multiple regions, indicating that the models are better able to replicate the sFC of strong and stable connections. This is likely due to the time-independent influence of the underlying SC. TC fit exhibited predominantly negative correlations, showing that the approach favors more random dynamic fluctuations. This effect is presumably caused by the noise added to the regional activity, and suggests that the models can better reproduce arbitrary changes in activity over time, rather than the slow fluctuation between distinct states observed in empirical data (Allen et al., 2014; Vidaurre et al., 2017). sFC features could predict FC variance and sFC fit; in contrast, TC features achieved high prediction scores for FC variance and TC fit, with a small number of features proving to be particularly predictive. These findings show that the modelling framework can better replicate certain connectivity signatures, particularly stronger sFCs within and between the DAN, VAN and LN, and weaker TCs across the brain, particularly in the DAN and VN. Patterns of sFC in SMN, VN and VN-DAN connections, as well as in the TC of LN and VAN, specifically in the left medial orbitofrontal gyrus, left pars opercularis, and right fusiform gyrus, contributed most strongly to fit predictions. The relationship between FC patterns and fit should be considered in future studies, as it provides a potential source of bias if the connectivity in these regions is differentially distributed between groups or inconsistent across sites. In addition, we identified some significant age and sex effects on model fits, potentially due to an overrepresentation of FC patterns that can be simulated well in some demographic groups. While we did not find sex or age differences in the resulting optimal parameters, future studies should investigate this discrepancy to ensure model-derived findings are generalizable across populations.
Several of the model parameters exhibited related signatures in the dynamic and static FC. G was generally positively correlated with connectivities, likely because a higher G leads to a higher impact of signals from other areas on regional activity, resulting in increased integration (Sanz-Leon et al., 2015). Connections that were significantly associated with G were concentrated within and between the DAN, SMN and VN. Evidence suggests that multiple neurotransmitters modulate FC in distinct ways. While the impact of dopamine is strongest in the DMN, the VAN (Conio et al., 2020), and the SMN, serotonin affects FC in the DMN, SMN, the FPN and the auditory network (Klaassens et al., 2015). G integrates the effect of disparate systems, suggesting that the regions in which the connectivity correlates with G will see a higher change in connectivity if the overall coupling is altered. Since G scales the combined input from other areas into regional activity, those connections that are strongly linked to G are likely between two regions that receive similar inputs. J_N_ correlated negatively with the connectivity of regions across the brain, particularly in the VN and in regions connecting the VN and DMN to other networks, although we did not successfully validate these associations in the data from the second scanning session (Supplementary Figure 2). This finding matches previous reports linking disruption in NMDA signaling, via administration of NMDA receptor antagonist ketamine, to brain hyperconnectivity (Anticevic et al., 2015; Driesen et al., 2013) as well as increased dynamic network flexibility (Braun et al., 2016). G and J_N_ were strongly associated with TC in many regions, indicating that differential moderation of long-range connections and excitatory input is a major driver of TC variability. G correlated positively with regions across the brain, showing that an increase in coupling also leads to increased temporal stability of FC. J_N_, on the other hand, correlated negatively with virtually all regions, indicating that increased excitatory signaling might lead to more random FC dynamics. The associations of neurophysiological characteristics with static and dynamic FC highlights the promise of the combined computational modelling approach for the investigation of biological mechanisms governing communication between brain regions. Further research should determine whether model parameters can capture pathological alterations in patients with brain disorders, and might serve as biomarkers or highlight therapeutic targets.
Perturbation analysis revealed that the coupling in a few regions had a disproportionate effect on overall network connectivity. Altering the coupling in some regions led to particularly large variations in sFC and TC fit and efficiency, suggesting that the coupling level in these regions has an outsized impact on brain connectivity. Other regions exhibited a common shift towards lower fits for most levels of perturbation. These regions appear to be particularly biologically relevant, as the optimal global parameter converged to their regional optimal parameter when fitting the model to empirical data. On the other hand, the few regions in which perturbation led to a mean increase in fit appear to be less relevant to shaping biologically accurate connectivity patterns but could be useful for further improving model accuracy. The fact that the regions which particularly affected sFC fit and efficiency differed from those affecting TC fit indicates that static and dynamic connectivity are governed by distinct neurobiological systems. While the regions central to sFC and efficiency, specifically the left and right paracentral, and the right pre- and postcentral as well as the right transverse temporal gyrus, govern the average functional network architecture, the regions relevant for TC, particularly the left pars triangularis, appear to be vital to the gradual transition between short-term network configurations.
Regions in which perturbation had particularly strong effects exhibited high centrality in the SC, showing that FC is driven to a large extent by regions which facilitate many structural connections. In addition, the means of the perturbation-related differences correlated highly with regional centrality in the functional connectivity, indicating that regions which constitute hubs in the FC are particularly important for maintaining biologically accurate connectivity. These findings support the evidence from lesion modelling studies, which suggests that disturbances in hub nodes (Achard et al., 2006; Aerts et al., 2016), particularly those along the cortical midline (Alstott et al., 2009), have the strongest effect on FC. The role of the regions we identified, particularly the pre-, post- and paracentral gyri, in maintaining normal brain connectivity further explains the observation that these regions are particularly well protected against stroke (Thirugnanachandran et al., 2024).
The expression of genes associated with a number of biological processes was related to the strength of the perturbation effect. The relevant mechanisms included some specific to the central nervous system, such as nervous system development and behavior. Other, more general processes, such as membrane organization and phospholipid metabolism, might nonetheless influence brain connectivity, as they are involved in myelination (Raasakka et al., 2017), a key factor in preserving connections between brain regions (Huntenburg et al., 2017; Vandewouw et al., 2021). In addition, the impact of regional perturbations also correlated strongly with the expression of genes relating to the modulation of chemical synapse signaling, indicating that regions which affect FC most strongly exhibit more neuromodulatory activity. This further supports the hypothesis that neuromodulation and neurotransmission are essential contributors to the establishment and maintenance of global FC.
While the majority of our findings could be replicated in a second set of resting-state fMRI data of the same participants (Supplementary Figure 2), some limitations should be considered in the interpretation of these results. Firstly, we modelled individual brain activity based on an averaged empirical SC for all subjects. While this provided a clearer insight into the relationship between FC and model parameters, employing individual structural connectomes could help elucidate individual structure-function relationships. The use of functionally informed SCs has been shown to further improve model fits (Manos et al., 2023). Further, fitting model parameters individually for each region could contribute to an even more detailed picture of individual neurobiology. Given that we found only a subset of regional connectivity features reflected individual variability in model parameters, analysing these relationships in patients with brain disorders might provide new information on brain alterations. In addition, we focused our analysis on a single atlas. As research has shown that the choice of atlas can impact the quality of models (Domhof et al., 2021), future studies should attempt to validate our results using a different atlas.
In this study, we showed that modelling the static and dynamic architecture of FC allows us to investigate neurobiological correlates of brain network dynamics. We found that differences in long-range inputs drive individual variability in some aspects of global dynamic FC, while variability in static FC is shaped by a combination of regional parameters. In addition, we identified that the coupling in a subset of frontal regions has a major impact on global network connectivity. A future investigation of the relationships between neurobiology and dynamic FC in brain disorders might reveal new insights into the origins of brain abnormalities in patients.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Achard S.Salvador R.Whitcher B.Suckling J.Bullmore E. (2006). A resilient, low-frequency, small-world human brain functional network with highly connected association cortical hubs. J. Neurosci. 26 63–72. 10.1523/JNEUROSCI.3874-05.2006 16399673 PMC 6674299 · doi ↗ · pubmed ↗
- 2Aerts H.Fias W.Caeyenberghs K.Marinazzo D. (2016). Brain networks under attack: Robustness properties and the impact of lesions. Brain 139 3063–3083. 10.1093/brain/aww 194 27497487 · doi ↗ · pubmed ↗
- 3Alexander-Bloch A. F.Shou H.Liu S.Satterthwaite T. D.Glahn D. C.Shinohara R. T. (2018). On testing for spatial correspondence between maps of human brain structure and function. Neuroimage 178 540–551. 10.1016/j.neuroimage.2018.05.070 29860082 PMC 6095687 · doi ↗ · pubmed ↗
- 4Allen E. A.Damaraju E.Plis S. M.Erhardt E. B.Eichele T.Calhoun V. D. (2014). Tracking whole-brain connectivity dynamics in the resting state. Cereb. Cortex 24 663–676. 10.1093/cercor/bhs 352 23146964 PMC 3920766 · doi ↗ · pubmed ↗
- 5Alstott J.Breakspear M.Hagmann P.Cammoun L.Sporns O. (2009). Modeling the impact of lesions in the human brain. P Lo S Comput. Biol. 5:e 1000408. 10.1371/journal.pcbi.1000408 19521503 PMC 2688028 · doi ↗ · pubmed ↗
- 6Anticevic A.Corlett P. R.Cole M. W.Savic A.Gancsos M.Tang Y. (2015). N-Methyl-D-aspartate receptor antagonist effects on prefrontal cortical connectivity better model early than chronic schizophrenia. Biol. Psychiatry 77 569–580. 10.1016/j.biopsych.2014.07.022 25281999 · doi ↗ · pubmed ↗
- 7Barber A. D.Hegarty C. E.Lindquist M.Karlsgodt K. H. (2021). Heritability of functional connectivity in resting state: Assessment of the dynamic mean, dynamic variance, and static connectivity across networks. Cereb. Cortex 31 2834–2844. 10.1093/cercor/bhaa 391 33429433 PMC 8325018 · doi ↗ · pubmed ↗
- 8Braun U.Schäfer A.Bassett D. S.Rausch F.Schweiger J. I.Bilek E. (2016). Dynamic brain network reconfiguration as a potential schizophrenia genetic risk mechanism modulated by NMDA receptor function. Proc. Natl. Acad. Sci. U.S.A. 113 12568–12573. 10.1073/pnas.1608819113 27791105 PMC 5098640 · doi ↗ · pubmed ↗