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.

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, (C_ij) the structural connectivity between each pair of regions i and j. J_i represents the local feedback inhibitory synaptic coupling, I₀ the external input, w_(E,I) the population external input scaling weight, w_p the local excitatory recurrence, and J_N 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 J_i, and (iii) recurrent excitation, modulated by the local excitatory recurrence w_p and the local excitatory synaptic coupling J_N. We investigated simulated data based on a range of values for G, J_N, J_i, and w_p.

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.

In order to reduce the parameter space, we performed an initial exploration of parameters to identify a range for J_N, J_i, and w_p 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, J_i, and w_p 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.

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.

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² 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² 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).

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.

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.

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 J_N. 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 J_N, 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 J_N across all regions (mean r < −0.16). There were no significant associations between J_i, representing inhibitory synaptic coupling, or w_p, 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 J_i and w_p do not.

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.

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²(sFC to FCV) = 0.44, R²(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²(sFC to sFC) = 0.72, R²(TC to FCV) = 0.68). Neither of the models performed particularly well in predicting the fit according to NC (R² < 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² 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 J_N reached a moderate value (R²(TC to J_N) = 0.12). The sFC-based prediction model performed worse than a constant model for all parameter targets.

For those prediction models with a high R² 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 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 (σ²(Δr) = 0.001), while the TC fit varied most strongly when the coupling was perturbed in the left pars triangularis and left pars orbitalis (σ²(Δr) > 0.008). Over all regions, the fit generally decreased the further G was altered from the benchmark (Supplementary Figure 3).

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.