Research ArticleCOGNITIVE NEUROSCIENCE

Dynamical consequences of regional heterogeneity in the brain’s transcriptional landscape

See allHide authors and affiliations

Science Advances  14 Jul 2021:
Vol. 7, no. 29, eabf4752
DOI: 10.1126/sciadv.abf4752

Abstract

Brain regions vary in their molecular and cellular composition, but how this heterogeneity shapes neuronal dynamics is unclear. Here, we investigate the dynamical consequences of regional heterogeneity using a biophysical model of whole-brain functional magnetic resonance imaging (MRI) dynamics in humans. We show that models in which transcriptional variations in excitatory and inhibitory receptor (E:I) gene expression constrain regional heterogeneity more accurately reproduce the spatiotemporal structure of empirical functional connectivity estimates than do models constrained by global gene expression profiles or MRI-derived estimates of myeloarchitecture. We further show that regional transcriptional heterogeneity is essential for yielding both ignition-like dynamics, which are thought to support conscious processing, and a wide variance of regional-activity time scales, which supports a broad dynamical range. We thus identify a key role for E:I heterogeneity in generating complex neuronal dynamics and demonstrate the viability of using transcriptomic data to constrain models of large-scale brain function.

INTRODUCTION

Brain dynamics are often described as complex, displaying properties that are interposed between order and disorder (1, 2). These complex dynamics arise from two cardinal factors: (i) the properties of local population activity within each brain region and (ii) the mutual influences that these populations exert on each other, as mediated by their interconnecting anatomical projections.

Biophysically plausible nonlinear models of brain dynamics have been used to clearly demonstrate the importance of anatomical connectivity in supporting complex patterns of network activity. For example, when simulated on human interregional structural connectivity (SC) matrices (connectomes) generated with diffusion magnetic resonance imaging (MRI), these models can reproduce empirically observed patterns of interregional functional connectivity (FC) [reviewed in (2)]. In macaque, a biophysical model simulated on a connectome derived from viral tact tracing indicated that both feedback projections and a heterogeneous distribution of connection strengths across projections are required for the emergence of ignition-like dynamics (3), which are thought to be a necessary condition for conscious processing (4).

By comparison, variations of local circuitry within each region have received less attention. Precise coordination of firing across spatially distributed neuronal ensembles depends on a delicate balance between excitatory and inhibitory activity. The first wave of biophysical models of large-scale brain activity treated all local population dynamics as homogeneous, driven by subpopulations of inhibitory and excitatory neurons with uniform properties across all brain regions (57). However, heterogeneity in regional cytoarchitecture has been noted for more than a century, as famously mapped by Brodmann (8), and subsequent work has provided mounting evidence for large-scale cortical gradients in gene expression, cellular composition, connectivity, and function (9, 10). Thus, quantitative variations of local circuit properties across the cortex also play an important role in shaping complex neural dynamics, likely through their influence on the local ratio, or balance, of excitatory and inhibitory cell activity (E:I balance) (9).

Ideally, one could model regional heterogeneity using empirically observed estimates of areal variations in excitatory and inhibitory cell counts or a related measure of cell type–specific activity. No such data are available on a large scale for primate brains, but some early modeling work has used alternative methods to constrain biophysical models. Deco et al. (11) showed that heterogeneity in local feedback inhibition, in which E:I balances were algorithmically adjusted to achieve a uniform firing rate of 3 Hz across all regions [in line with experimental data; (1214)], yields networks with more stable dynamics and improves model fits to static measures of FC when compared with a strictly homogeneous model. Using data in macaque, Chaudhuri et al. (15) scaled the excitatory input strength of each of 29 cortical areas according to its assumed hierarchical position, as determined by laminar patterning of interregional projections (16). They found that this heterogeneity was essential for generating dynamics characterized by a hierarchy of regional-activity time scales, in which areas with a higher hierarchical position show a prolonged response to a simulated sensory input. This result mirrors previous single-unit findings in macaque cortex (17) and aligns with other evidence for regional differences in activity time scales (1820). Two subsequent studies modeling human functional MRI (fMRI) found that, when compared to a regionally homogeneous model, scaling a region’s excitatory activity according to its assumed hierarchical position improves fits to measures of edge-level and node-level FC, more faithfully reproduces empirically observed regional variations in activity time scales, and better explains FC disturbances in people with schizophrenia (21, 22). A complementary study used a stochastic optimization procedure to invert local dynamical parameters such as recurrent synaptic connectivity strengths to optimally fit fMRI data (23). Regional variations in the estimated local dynamical parameters showed hierarchical ordering and correlated with variations in myelin and layer-specific neuronal density, presumed cognitive functions (derived from meta-analysis of task-based fMRI studies), and a principal gradient of resting-state FC (24).

These initial models of regional heterogeneity demonstrate the critical role of local variations in E:I balance for generating realistic patterns of brain-wide dynamics. However, in the absence of a gold-standard biological constraint, there has been considerable variability in how heterogeneity has been instantiated within a given model. There is also substantial variability in the dynamical properties that have been modeled. For example, modeling work in macaques shows that it is, in principle, possible to find a parameter regime that results in ignition-driven and hierarchical dynamics (3, 15), but the parameters were not fitted to empirical data, making it difficult to evaluate the model’s capacity to reproduce other key features, such as realistic patterns of static or dynamic FC. Similarly, studies in humans have modeled specific outcome measures, such as static or dynamic FC at the level of individual connections (6, 7, 22, 25), the mean static FC of each region (22), or time scale hierarchies (3, 22), meaning that it is unclear whether the working point for a model fitted to one property can also capture other properties. Developing a coherent approach that can account for each of these diverse features is an essential step toward developing a unified model that can parsimoniously explain diverse empirical properties of complex whole-brain dynamics.

Here, we seek to address this aim by introducing a model that leverages brain-wide transcriptional data (26) to constrain regional heterogeneity and tune each region’s dynamics according to regionally specific measures of inhibitory and excitatory receptor gene expression. We evaluate the performance of this model with respect to some of the major outcome measures evaluated thus far in the literature; namely, edge-level static FC, node-level static FC, FC dynamics (FCD), ignition capacity, and time scale hierarchies. We show that our model more faithfully reproduces empirical properties of static and dynamic FC than other models in which heterogeneity is imposed using plausible biological constraints, such as regional variations in the ratio of T1-to-T2–weighted MRI signal (T1w:T2w), which is thought to index intracortical myelin content (27) and has been proposed as a marker of a cortical region’s hierarchical position (28), and regional loadings on the dominant mode of gene expression (22, 28). We then show that transcriptomically constrained models show the greatest capacity for ignition-like dynamics and display a broad range of regional-activity time scales. Our results indicate that regional variations in the transcriptional activity of inhibitory and excitatory receptor genes provide a viable means for constraining biophysical models of large-scale neural dynamics and, when coupled with an empirically derived connectivity matrix, can offer a parsimonious account of diverse properties of large-scale human brain activity.

RESULTS

Overview

Our overall analysis strategy is outlined schematically in Fig. 1, and the details of our approach are provided in Methods. Briefly, we use diffusion MRI to generate a group-level representative structural connectome for a sample of 293 healthy individuals, representing SC between each pair of 68 regions defined according to an extensively used neuroanatomical atlas (Fig. 1A). We also estimate, for a partially overlapping sample of 389 individuals, various empirical FC properties using task-free resting-state fMRI to probe spontaneous blood oxygenation level–dependent (BOLD) dynamics (Fig. 1B).

Fig. 1 Workflow for model fitting and evaluation.

We use empirical diffusion MRI and fMRI data to generate sample-averaged SC (A) and FC (B) matrices, respectively. The SC matrix, and different forms of regional heterogeneity (C; for full maps, see fig. S1), are used to respectively define the coupling between, and modulate the gain of, neuronal populations within a large-scale balanced excitation-inhibition (BEI) dynamical model (D). For each model, the global coupling parameter G is tuned to optimize the fit to three empirical FC properties (E): static edge-level FC (left), static, node-averaged FC (middle), and FCD (FC; right). Gain parameters B and Z are then tuned to further improve fits to these properties. (F) After finding the optimal working point for each model, we simulate a focal perturbation of activity in occipital cortex and measure two dynamical properties of the network: ignition and activity time scales. (G) To measure ignition, we characterize the response of each area to focal stimulation as the product of two quantities: the maximum firing rate achieved at the highest stimulation intensity, rmax, and the concavity of the regional response function (i.e., the maximal second derivative), cmax. The ignition capacity of the model is taken as the average of these regional values. (H) To measure the decay time scale of each region, we fit an exponential curve to each region’s firing rate as it returns to baseline after termination of the maximal stimulation intensity (H). We also measure intrinsic time scales during a continuous white noise stimulation paradigm (see Methods).

We incorporate different forms of regional heterogeneity (Fig. 1C and fig. S1) into a whole-brain model in which the local neuronal dynamics of each brain region evolve according to a dynamic mean field reduction analytically derived from the collective behavior of empirically validated integrate-and-fire (80% excitatory and 20% inhibitory) neurons (13, 14). Regional dynamics in the model are driven by a pool of excitatory neurons, a pool of inhibitory neurons, the interactions between them, and the net output of other areas, as mediated by the SC matrix (Fig. 1D). In this model, local feedback inhibition is adjusted separately for each region to ensure a consistent firing rate of 3 Hz across regions (11). We thus refer to this model as the balanced excitation-inhibition (BEI) model. Although the local adjustments in this model introduce some degree of regional heterogeneity, the firing rates are constrained to be uniform across regions, so we consider this BEI model as a homogeneous benchmark against which we evaluate more sophisticated heterogeneous models that allow intrinsic dynamical properties to vary across regions [see (22) for a similar approach]. More specifically, the activity levels of excitatory and inhibitory populations in each region, i, of the BEI model are given by independent sigmoid functions, regulated by a single gain parameter Mi. In our benchmark BEI model, this parameter is set to a fixed value for all regions, imposing uniform excitability across brain regions (Fig. 1C).

We compare the performance of this homogeneous BEI model to three different heterogeneous models. In the first, we use an approach similar to (22) and adjust Mi according to variations in the regional mean T1w:T2w, which has been proposed as a proxy for cortical hierarchy (28). In the second heterogeneous model, we adjust Mi according to the first principal component (PC1) of transcriptional activity for 1926 brain-specific genes (28) that survived our quality control criteria, as quantified in the Allen Human Brain Atlas (AHBA) (26, 29, 30). The first PC of these genes correlates with the T1w:T2w ratio and other measures of cortical hierarchy (28). In the final model, we use a more hypothesis-driven approach and adjust Mi according to variations in the ratio of regional AHBA expression values for genes specifically coding for AMPA, N-methyl-d-aspartate (NMDA), and γ-aminobutyric acid (GABA) receptors (see Methods). We refer to this model as the E:I model. Cortical surface renderings displaying regional variations in T1w:T2, PC1, and E:I values are shown in Fig. 1C and fig S1.

In all three heterogeneous models, we assume a linear scaling between the regional biological measures of heterogeneity, Ri, and the effective gain within a region (22), given by Mi = 1 + B + ZRi, where the two unknown free parameters correspond to a bias, B, and a scaling factor, Z (see Methods). Studying how these two free parameters affect the global dynamics of the model allows us to investigate the role of regional heterogeneity.

For the homogeneous BEI and three (T1w:T2w, PC1, and E:I) heterogeneous models, we assume that all diffusion MRI–reconstructed streamline fibers have the same conductivity, and thus, the coupling between different brain areas is scaled by a single global parameter, G. We first tune the G parameter of the BEI model to adjust the strength of effective coupling in the model and identify the brain’s dynamic working point by fitting the model to three empirical FC properties (Fig. 1E): (i) the Pearson correlation between model and empirical estimates of static (i.e., time-averaged) FC estimated across all pairs of brain regions (more specifically, the correlation between the values in the upper triangles of the model and empirical FC matrices); (ii) the Pearson correlation between model and empirical node-level estimates of average FC (22); and (iii) similarity in FCD, estimated as the distance between model and empirical distributions of frame-by-frame FC properties (as quantified using the Kolmogorov-Smirnov distance DKS) (25).

After finding the working point of each model, we evaluate its dynamical properties in two ways. First, we consider the ignition capacity of the model, quantified by examining how activity propagates through the network following a simulated focal perturbation (Fig. 1, F and G). Second, we examine regional heterogeneity in intrinsic time scales in two ways, following (15), by quantifying the signal decay rate of regional activity following a simulated focal perturbation (Fig. 1, F and H) and by studying regional signal autocorrelations during continuous, focal white noise input. Our approach thus allows us to determine the performance of each model in capturing diverse dynamical phenomena at a single working point.

Fitting the homogeneous model

We first evaluate the performance of the homogeneous BEI model in reproducing empirical properties of resting-state FC data. We determine and fix the global coupling parameter G with all regional gain parameters Mi = 1 (i.e., heterogeneity is not considered). As done previously (11, 31), we examine how well the model fits, as a function of G, three different properties of empirical resting-state fMRI recordings: edge-level static FC, FCD, and node-level FC (see Fig. 1B and Methods for further details). The results of this analysis are shown in Fig. 2 (A to C). In all cases, we consider group-averaged values for both empirical data (across participants) and the model (across an equivalent number of simulated trials with the same duration as the human experiments; see Methods). Our results align with previous investigations to clearly show (25) that fitting FCD, which captures the rich spatiotemporal structure of the fMRI data, is a stronger constraint on the model. More specifically, where static edge and node FC fits are consistently high across a broad range of G, FCD yields a clear global optimum at G = 2.1 (Fig. 2A).

Fig. 2 Model optimization and evaluation.

The homogeneous BEI model is first fitted to the empirical data. (A) The correlation between model and empirical static edge-level and node-level FC, as well as the DKS between empirical and model FCD distributions, as a function of the global coupling parameter G. Only the fit to FCD shows a clear optimum at G = 2.1 (note that lower FCD fit values indicate better fit). (B) and (C) show the FC matrix, node-level average FC values, FCD matrices, and FCD distributions for the empirical data (B) and model (C). Using G = 2.1, we introduce regional heterogeneity (Fig. 1C) and further optimize the two gain parameters B and Z. (D) to (F) show the parameter landscapes for edge-level FC, node-level FC, and FCD, respectively, for the E:I model. Below each two-dimensional (2D) space, we show the 1D graph of the level heterogeneity scaling (Z) for the optimal bias parameter found as the one corresponding to the optimum FCD fitting (see Methods). This procedure is repeated for each heterogeneous model to find its optimal working point. Lower values indicate better fit in (F). (G) to (I) compare the performance of the homogeneous BEI model and heterogeneous T1w:T2w, PC1, E:I, and null expression models according to the edge-level FC, node-level FC, and FCD benchmarks, respectively. Distributions are for values obtained over multiple runs of the optimum models (see Methods). We harmonized the fit statistics so that higher numbers indicate better fit (i.e., for the case of FCD, where a low DKS indicates better fit, we show 1 − DKS). The E:I model shows the best performance in all cases.

Introducing regional heterogeneity

We now study how regional heterogeneity affects the fitting of static edge-level FC, node-level FC, and FCD. Spatial maps for each form of biological heterogeneity used in our modeling are shown in Fig. 1C and fig. S1. The T1w:T2w and PC1 maps are strongly correlated (ρ = 0.82, Pspatial < 0.001), as shown previously (28), whereas the PC1 and T1w:T2w maps show weaker correlations with the E:I map (ρ = 0.40, Pspatial = 0.005 and ρ = 0.39, Pspatial < 0.011, respectively). This result indicates that each spatial map introduces a different form of biological heterogeneity to the benchmark BEI model. We introduce this heterogeneity by modulating the regional gain functions Mi, at the optimal working point of the homogeneous BEI model (G = 2.1), through the bias and scaling parameters introduced above, denoted B and Z, respectively. Critically, we exhaustively search a broad range of values for B and Z to comprehensively map the parameter landscape and select the corresponding optimal working point for the E:I model (B = −0.3; Z = 1.8), the PC1 model (B = −0.75; Z = 1.6), and the T1w:T2 model (B = −0.7; Z = 1.4), as shown in Fig. 2 (D to F) for the E:I model. As can be seen in this figure, there is some degeneracy in the parameter landscapes but there is also a clear regime in which all three empirical properties are fitted well by the model, particularly for Z ≈ 1.7.

For these optimal values, we simulate each dynamical model 1000 times to account for the inherent stochasticity of the models and compute the respective measures of model fit. Figure 2 (G to I) shows the distributions of fit statistics across runs for the homogeneous and three heterogeneous models. In addition, we show results for a null ensemble of models in which the regional E:I receptor gene expression values were spatially rotated (32, 33) to generate surrogates with the same spatial autocorrelation as the empirical data. Across all three benchmark properties to which the data were fitted—static FC, GBC, and FCD—the heterogeneous models perform better than the homogeneous model (8.64 < z < 38.70, all Pbonf < 0.05). We also find a consistent gradient of performance across all benchmarks, with the E:I model performing best, followed by the PC1 and T1w:T2w models, and the homogeneous model showing the poorest performance. For each benchmark metric, the performance of the E:I model was significantly better than all other models (33.12 < z < 38.70, all Pbonf < 0.05). Out-of-sample model fit statistics yielded similar findings (fig. S2), supporting the generalizability of our results. Nonetheless, differences in model fit statistics were small. For instance, in Fig. 2G, the median static edge-level FC correlation between model and data was r = 0.71 for the E:I model, r = 0.69 and r = 0.70 for the T1w:T2w and PC1 models, respectively, and r = 0.67 for the BEI model. For static node-level FC, the correlations ranged between r = 0.91 (E:I model) and r = 0.87 (BEI model) (Fig. 2H). Similarly, the fits to FCD ranged between DKS = 0.08 for the E:I model and DKS = 0.11 for the BEI model (Fig. 2I). Despite these small differences in model fits, we show in the following sections marked differences between models in other dynamical properties.

Ignition capacity

We now evaluate additional dynamical properties of each model that cannot be directly evaluated with respect to empirical fMRI data but that are thought to be essential features of complex neural dynamics. We first consider ignition capacity. Ignition refers to the capacity of a sensory input to trigger self-supporting, reverberating, temporary, metastable, and distributed network activity and is thought to be a necessary condition for conscious perception of the stimulus (4, 34). This process is nonlinearly related to stimulus intensity, such that evoked activity is confined to local sensory areas for low levels of stimulation (which supports subliminal perception) and then leads to explosive activation of distributed cortical systems once an appropriate threshold is reached.

Figure 1 (F and G) shows the strategy that we follow to quantify ignition capacity. For a given specific working point of the whole-brain model, we compute how the brain broadcasts information after artificially stimulating the lateral occipital cortex region of the Desikan-Killiany atlas (35) bilaterally, which contains the foveal representation of V1. We measure the evoked responses at the level of population firing rates rather than simulated BOLD signal changes to have direct access to the millisecond time scale. To quantify the effect of occipital stimulation on activity in each of the other 66 brain regions, we plot, for each nonstimulated region, how its population firing rate changes as a function of occipital stimulation intensity. Two quantities are relevant here: (i) the maximum firing rate achieved at the highest stimulation intensity, rmax; and (ii) the speed with which the firing rate increases beyond a given intensity threshold, cmax, which we quantify as the concavity of the regional response function (i.e., the maximal second derivative). We then summarize the ignition capacity of each nonstimulated region i as I(i) = cmax(i) × rmax(i) and estimate the global ignition capacity of the brain as the mean across all regions, I=1N1I(i).

Figure 3A shows the response time courses of all left hemisphere regions of the E:I model to the maximum stimulus intensity. Strong responses were evoked throughout the brain, including prefrontal cortex. Figure 3B shows cortical surface renderings of regional I(i) values obtained following occipital stimulation over 30 runs of the homogeneous BEI and heterogeneous E:I models. The occipital stimulation in the E:I model triggers more widespread ignition-like propagation than in the BEI model. This observation is confirmed quantitatively in Fig. 3C. Specifically, the E:I model shows higher ignition capacity than all other models (20 < z < 39, all Pbonf < 0.05), suggesting that this model is able to rapidly propagate focal perturbations with higher fidelity than in the other models. Notably, the spatial null models showed higher ignition capacity than the PC1, T1w:T2w, and BEI models (17.4 < z < 33.5, all Pbonf < 0.05), suggesting that generic spatial gradients with the same autocorrelation as the E:I maps can yield strong ignition-like dynamics on average, although considerable variability was observed around the median (Fig. 3C). The PC1 model showed significantly greater ignition capacity than the both the T1w:T2w and BEI models (z = 32.5, P = < 0.001 and z = 29.2, P < 0.001, respectively), which did not differ from each other (z = 1.22, P = 0.22). A separate evaluation of ignition-like dynamics following stimulation of the postcentral gyrus, which contains primary somatosensory cortex, also revealed greater ignition capacity for the E:I compared with the BEI model, although the differences were not as pronounced as in the visual system (fig. S3). Together, these results indicate that ignition-like propagation following focal visual stimulation is most potent when regional heterogeneity is constrained by the transcriptional activity of E:I receptor genes.

Fig. 3 Evaluating the ignition capacity of each dynamical model.

(A) Time course of the impulse stimulation (top left) used to study ignition-like dynamics and decay time scales, followed by regional responses to the maximum stimulus intensity for all left hemisphere regions in the E:I model. Bottom right shows each region on a cortical surface map with the colors corresponding to the response traces and the stimulation site indicated by the yellow marker. (B) Cortical surface renderings of average regional ignition capacity values, measured in arbitrary units (a.u.), obtained after lateral occipital cortex stimulation in the E:I model (top) and BEI model (bottom). The E:I model shows widespread ignition-like propagation across distributed areas extending to prefrontal cortex, whereas the BEI model shows negligible ignition of nonvisual areas. (C) Distributions of global ignition capacity values obtained for 1000 different runs of the optimal E:I, PC1, T1w:T2w, BEI, and spatial null models. Ignition capacity is highest in the E:I model.

Regional activity time scales

We characterize temporal hierarchies in the models using two approaches, following Chaudhuri et al. (15). First, we examine the decay rate of each region’s activity following impulse stimulation of lateral occipital cortex using the same virtual stimulation approach as in our analysis of ignition capacity (Fig. 1H). Specifically, we apply a virtual impulse stimulation and quantify the temporal decay of activity in each nonstimulated target area by fitting an exponential curve to its firing rate as it returns to spontaneous conditions after termination of the maximal stimulation intensity (Fig. 1H). The value of the exponent can be used to index the decay time scale of each brain region and offers one way to examine temporal hierarchies in the brain (see Methods) (15). These fitted exponential curves are shown for all regions in the E:I model in Fig. 4A. A clear trend is evident for rapid decay in visual regions and more prolonged responses in prefrontal areas. Supporting the validity of this putative temporal hierarchy, we find a significant positive correlation between regional decay values and an independent measure of cortical hierarchy, T1w:T2w (ρ = 0.44, Pspatial < 0.001; Fig. 4B), such that higher T1w:T2w is associated with more rapid decay time scales (higher decay exponents).

Fig. 4 Time scales of regional dynamics.

(A) Normalized decay curves for all left hemisphere brain regions showing the rate at which regional activity returns to baseline following an impulse stimulation to lateral occipital cortex in the E:I model. Decay rates are slower in prefrontal compared to visual areas. Colors correspond to the regions shown in the inset. (B) Association between regional decay rates (in s−1) in the E:I model and empirical T1w:T2w values. (C) Intrinsic activity time scales (in milliseconds), quantified as ACFs of each brain region obtained following an input of continuous white noise to lateral occipital cortex. Prefrontal areas again show prolonged time scales. (D) Variability of intrinsic regional time scales (interquartile range) obtained for 1000 runs of each model.

We next estimate the intrinsic time scale of each region’s activity using the autocorrelation of its activity time series simulated during continuous white noise stimulation of lateral occipital cortex. This approach is more directly comparable to the experimental paradigms previously used to characterize intrinsic time scales of neuronal activity (17, 36, 37). Figure 4C shows the fitted autocorrelation functions (ACFs) for every brain region for the E:I model (see Methods). Once again, a clear trend is evident for more rapid time scales in visual areas and slower time scales in prefrontal regions. We define the richness of regional time scales in the model as the interquartile range of decay rates across target areas, such that high values reflect a greater diversity of time scales across the brain, which implies a broader dynamic range for information processing. We term this measure regional decay variability. Figure 4D indicates that the PC1 model shows the highest decay variability, followed by the E:I and spatial null models, all of which show greater decay variability than the homogeneous BEI model (14.2 < z < 38.2, all Pbonf < 0.05). The decay variability of the BEI and T1w:T2w models is comparable (z = 0.25, P = 0.80). Both the E:I and PC1 models show greater decay variability than the spatial null model (z = 16.45, P < 0.001 and z = 38.54, P < 0.001, respectively), indicating that empirical transcriptomic maps of heterogeneity confer a wider dynamic range than generic spatial gradients with similar autocorrelation.

DISCUSSION

It has long been recognized that cortical areas vary along multiple dimensions, including cytoarchitecture, myeloarchitecture, chemoarchitecture, and interregional connectivity (810, 38). The specific dynamical consequences of each of these variations remain a mystery, but they are likely to have an ultimate influence on the local ratio of excitatory to inhibitory cell activity, resulting in a variable E:I balance across different brain regions. Incorporating such heterogeneity into large-scale models of neuronal dynamics has proven challenging, but the recent availability of diverse data that can be used to index different forms of cellular, molecular, and anatomical variability across brain areas has provided a new means for imposing local variations in these models, thus offering an opportunity to understand in detail the dynamical effects of different forms of regional heterogeneity.

To this end, we compared models in which regional heterogeneity was constrained by transcriptional measures of E:I receptor gene expression, the dominant mode of brain-specific gene expression, and T1w:T2w, which is sensitive to variations in myeloarchitecture and closely tracks regional position in the cortical hierarchy (28). We show that the transcriptional E:I model more faithfully reproduces empirical properties of static and dynamic FC while also showing the highest ignition capacity and a broad range of regional-activity time scales. Our findings thus highlight a central role for regional variations of E:I balance in generating ignition-like dynamics and time scale hierarchies and suggest that transcriptional atlas data provide a viable means for constraining heterogeneity in large-scale models of brain activity.

Evaluating different forms of regional heterogeneity

The starting point for our analysis was an evaluation of how each model captured diverse empirical properties of FC. Early work focused on reproducing the correlation between model and empirical static FC matrices at the edge-level (11, 31). Subsequent studies introduced additional constraints, such as evaluating the fit to dynamic FC properties (25) or regional variations in average FC (22). To our knowledge, no prior study has simultaneously evaluated all three properties. Our analysis shows that the edge- and node-level measures of static FC offer loose constraints for model optimization, showing comparably high fit statistics across a broad range of values of the global coupling parameter G. In contrast, fits to FCD show a clear optimum (at G = 2.1), mirroring similar results reported previously (25). Fitting models to both static and dynamic properties of FC is thus important for identifying an appropriate working point for each model.

Across all three FC properties, heterogeneous models provide a better match to the data than the nominally homogeneous BEI model. This result suggests that regional heterogeneity plays a role in shaping the empirical FC architecture of spontaneous BOLD dynamics. Critically, the E:I model was the best-performing model across all three benchmarks, indicating that constraining regional heterogeneity by transcriptional markers of E:I balance yields a more faithful replication of empirical FC. However, the differences in fit statistics between models are small, suggesting that these statistics only have a limited capacity to tease apart dynamical differences between the models.

The model-based evaluations of ignition capacity and decay variability provide a more complete picture. In both cases, differences between the homogeneous BEI model and heterogeneous models, particularly those constrained by transcriptomic measures, were apparent. The median global ignition capacity across runs was 112% larger for the E:I model compared with the BEI model, while the median intrinsic time scale variability was 34% larger for the E:I model. The median time scale variability of the PC1 model was particularly high, being 81% larger than the E:I model, which was the second most variable. These results indicate that both ignition-like dynamics and time scale hierarchies depend on some degree of regional heterogeneity and that biologically plausible forms of transcriptional heterogeneity increase these properties beyond the expectations of a generic spatial gradient (as represented by the spatially constrained null ensemble).

The E:I model showed the greatest ignition capacity, indicating that regional differences in E:I receptor gene expression may shape local dynamics in a way that promotes rapid and high-fidelity communication within the visual system. The E:I model also showed a wide range of intrinsic time scales, although the PC1 model showed the highest variability, implying that the dominant mode of gene expression in the brain, which correlates with cortical hierarchical position (28), may facilitate a broad dynamic range of information processing. However, whether increased time scale variability in a model provides a more faithful representation of neural dynamics remains unclear. Further work could aim to validate model results either against invasive recordings, which allow precise estimation of decay time scales to focal perturbations, or through specifically designed fMRI experiments (39, 40). Critically, we evaluated the dynamical properties of each model at their empirically optimal working points, meaning that the models show these properties in a parameter regime that closely matches actual data. This is an important constraint, given that prior modeling work evaluating ignition dynamics and time scale hierarchies in macaque did not fit empirical FC data (3, 15). It is thus difficult to ascertain the biological plausibility of the dynamical regime that was studied. The E:I model considered here not only showed the best fit to empirical data but also showed the highest ignition capacity and a broad range of decay time scales, suggesting that it offers a parsimonious account of these different dynamical phenomena.

How does heterogeneity shape neuronal dynamics?

In our models, we introduced heterogeneity by modifying the regional excitability of local population activity, as determined by each region’s gain response function. Different approaches have been used to incorporate regional heterogeneity in past work (15, 22, 23, 41, 42). For example, Wang et al. (23) used a “bottom-up” approach in which large-scale model inversion was used to infer regional variations in recurrent excitation strengths. The inferred parameters suggested that recurrent excitation is stronger in sensory and weaker in association cortices. This result contradicts the “top-down” approaches used by Chaudhuri et al. (15) and Demirtaş et al. (22), where regional heterogeneity was constrained using an independent biological property of hypothetical relevance. More specifically, Chaudhuri et al. (15) used a hierarchical ordering of cortical areas in macaque based on laminar patterning of interregional connectivity to linearly scale the excitatory input strength of each population. Demirtaş et al. (25) used T1w:T2w to linearly scale the strengths of local excitatory-to-excitatory and excitatory-to-inhibitory populations when modeling cortical BOLD FC in humans, given evidence that T1w:T2w tracks hierarchical position and is negatively correlated with pyramidal spine density (15, 28, 43). In both cases, recurrent excitation was assumed to increase from sensory to association areas. Our analysis also used a top-down approach, relying on transcriptomic estimates of inhibitory and excitatory receptor expression to modify regional gain functions. Our resulting E:I ratio map showed only a moderate positive correlation with T1w:T2w (ρ = 0.40), indicating that it quantifies a distinct aspect of regional microcircuit heterogeneity. Our evaluation of model fits to empirical FC properties indicates that the specific aspect of heterogeneity captured by the transcriptomic E:I ratio plays a more central role in driving realistic brain-wide dynamics when compared to the global mode of gene expression and T1w:T2w.

A comprehensive comparison between these different approaches for incorporating heterogeneity is complicated by the different methods used for defining heterogeneity (top-down versus bottom-up), different constraints on heterogeneity (T1w:T2w versus transcriptomic data), and different ways of modeling heterogeneity (modifying recurrent connection strengths versus regional gain functions). We focused on modulating population gain response functions because local variations in E:I balance will affect the net excitability of the population, which is captured by the gain function parameter Mi. We thus assume that changes in regional gain are the final common effect of different ways in which heterogeneity might influence a specific population or interaction between populations. Critically, regional variations of Mi are modulated by just two terms—the bias B and the scaling factor Z—which leads to a marked improvement in computational efficiency when compared with more highly parameterized models [e.g., (22, 23)] that require the use of heuristic optimization algorithms. These algorithms can shift some parameter values from empirically constrained ranges and do not guarantee a globally optimum solution. The further refinement of optimization algorithms capable of robustly and effectively searching high-dimensional parameter spaces will allow a more comprehensive comparison of different heterogeneous models.

Implications for understanding ignition-like dynamics

Ignition-like dynamics are thought to be a necessary condition for conscious experience (4, 34). As the intensity of a focal sensory stimulus increases, the ignition model predicts that a threshold is reached beyond which widespread activation of other areas is rapidly triggered. Stimuli that do not reach this threshold may only be perceived subconsciously.

Our results indicate that regional heterogeneity of the E:I ratio plays an important role in driving ignition-like dynamics. To the extent that our transcriptional E:I measures can be taken as a proxy for receptor abundance (see below), this result suggests that regional variations specifically in E:I receptor activity may play an especially critical role in supporting rapid, high-fidelity broadcasting of sensory stimuli throughout the network.

Recent work using a dynamical model simulated on the macaque connectome also showed evidence of ignition-like dynamics (3). Critically, the authors showed that removal of feedback projections and randomization of the connectivity weights of the SC matrix were each sufficient to abolish this dynamical behavior, suggesting that feedback connectivity and heterogeneity in SC connection strength give rise to ignition-like activity. In our analysis, the SC matrix is derived from diffusion MRI data and thus cannot distinguish feedforward from feedback connections, although we did incorporate empirical variability in connection strengths. Thus, the fact that we observe ignition-like dynamics with heterogeneous but not homogeneous models suggests that feedback connectivity and heterogeneous edge weights alone are not sufficient for ignition. Instead, heterogeneity of regional population dynamics, and of E:I balance in particular, may represent a minimal additional requirement for the emergence of ignition-like propagation throughout the network.

While this propagation was most pronounced following visual stimulation, we also found some evidence for greater ignition capacity in the E:I model compared with the homogeneous BEI model following stimulation of somatosensory cortex. This result supports the generality of our findings, although the differences between the models were less pronounced than those observed in the visual system. This effect could be driven by variations in the connectivity of the stimulation sites, the specificity with which these sites are parcellated, intrinsic differences in the ignition capacity of visual and somatosensory systems, or the possibility that regional E:I variations may be well placed to support ignition-like propagation through the visual processing hierarchy, which is consistent with anterior-to-posterior gradients in excitatory and inhibitory cell density (44, 45). The phenomenon of ignition has been most extensively described in the visual system [e.g., (46)]. Studying the properties of ignition-like propagation in different sensory systems will be an important avenue of further work (47). Moreover, while our approach offers one way to quantify ignition-like propagation throughout the brain, it will be important to develop more refined quantitative estimates that can be directly related to physiologically validated markers of ignition-driven changes in regional activity.

Implications for understanding time scale hierarchies

We assessed time scale hierarchies by examining regional dynamics following both impulse and continuous stimulation of visual cortex. These analyses converged to identify prolonged activity time scales in prefrontal and association areas compared to primary regions. Comparison of regional intrinsic time scale variability across models indicated that the E:I, PC1, and spatial null models showed a wider dynamic range than the T1w:T2w and homogeneous BEI models, with the E:I and PC1 models both showing greater variability than the spatial null model. Thus, spatial maps with an autocorrelation that matches the gene expression data can enhance the dynamic range of the biophysical model, but the largest gains come when using empirical maps of transcriptional heterogeneity. The relatively wide range of the PC1 model implies that the dominant mode of gene expression, which correlates with T1w:T2w and is related to cortical hierarchy (28), may play a particularly important role in shaping time scale differences between regions.

Past work has shown that several other properties correlate with various measures of regional-activity time scales. In a comprehensive survey of 6390 different regional time series properties estimated for fMRI recordings in the mouse brain, Sethi et al. (19) found that time scale–related measures, such as relative high-frequency signal power, correlate with structural measures of in-degree, such that the time series of regions with many incoming connections show greater low-frequency power. Similar results were later reported in human (18). Other work in humans has also shown that regional variations in autocorrelation and other time series properties track the cortical hierarchy (20, 48). These additional measures could be used to provide additional constraints on model fitting, although the temporal resolution of fMRI may be too limited to provide a detailed understanding of regional variations in activity time scales.

Limitations

We relied on a single coarse scale parcellation of the brain to allow exhaustive parameter sweeps over several different models. This choice was essential to reduce computational burden, and it ensures that our results are directly comparable with other studies using the same parcellation [e.g., (11)]. Nonetheless, parcellation type and scale can have an effect on empirical connectome properties (49, 50), and further evaluation of the degree to which our results generalize to other parcellations is required.

Our E:I model used regional measures of excitatory and inhibitory receptor gene expression as a proxy for the relative abundance of each receptor type in each region, which, in turn, is assumed to be related to the relative activity levels of inhibitory and excitatory cells. Gene expression in the AHBA was assayed using microarray, which quantifies mRNA abundance. The correspondence between mRNA levels and protein abundance can be complex and may vary across different tissues or brain regions (51). Recent work has used autoradiographic data to derive maps of regional E:I (52), which provides a more direct estimate of protein abundance. Qualitatively, the maps resulting from these two methods show some convergence, although some discrepancies are also apparent. Direct comparison between these maps is difficult because of the different parcellations used, the incomplete anatomical coverage and laminar segregation of the autoradiographic data, and different methods for estimating the E:I ratio. More specifically, whereas Goulas et al. (52) used data from 15 different receptors, including those for neuromodulatory systems such as dopamine, serotonin, and acetylcholine [derived from (53)], we relied on a more conservative list of genes directly related to receptors with a clear role in excitation and inhibition. Future work can use the analysis framework developed here to directly compare the efficacy of different estimates of receptor abundance and E:I ratios in improving model performance. These efforts will also benefit from the increasing availability of single-cell RNA sequencing data, which can provide more precise estimates of regional cell type abundances (54). Critically, any imprecision in our estimates should limit the ability of the E:I model to yield physiologically plausible dynamics. Our results may therefore represent a lower bound on the accuracy that such a model may achieve.

In summary, we show that regional heterogeneity improves model fits to empirical FC properties, particularly when models are constrained by transcriptomic estimates of E:I receptor gene expression. Critically, at the same working point, this model also shows the highest ignition capacity and a broad range of activity time scales, suggesting that regional variations in the activity levels of excitatory and inhibitory receptors may play an important role in the emergence of these dynamical properties. Our analysis indicates that recently constructed transcriptional atlas data provide a fruitful and biologically principled means for imposing regional heterogeneity in large-scale dynamical models.

METHODS

Empirical MRI data

Image acquisition. Eyes-closed resting-state fMRI data were acquired in 389 healthy individuals (170 males; mean age, 24 years; SD, 4.2 years) using a Siemens Skyra 3 Tesla scanner with a 32-channel head coil located at Monash Biomedical Imaging in Clayton, Victoria, Australia using the following parameters: time to repetition (TR), 754 ms; time to echo (TE), 21 ms; flip angle, 50o; multiband acceleration factor, 3; field of view (FOV), 190 mm; and voxel size, 3-mm3 isotropic. A total of 620 volumes were acquired. T1-weighted structural images were also acquired using 1-mm3 isotropic voxels; TR, 2300 ms; TE, 2.07 ms; and FOV, 256 mm by 256 mm. Diffusion MRI data were available for a subset of 289 individuals, acquired using the following parameters: 2.5-mm3 voxel size; TR, 8800 ms; TE, 110 ms; FOV of 240 mm by 240 mm, and 60 directions with b = 3000 s/mm2 and seven b = 0 volumes. A single b = 0 s/mm2 volume was obtained with the reverse-phase encoding for use in distortion correction. Further details can be found in (55, 56). All participants were recruited and data required in accordance with local ethics committee guidelines.

Diffusion MRI processing and connectome construction. The diffusion data were processed using MRtrix3 version 3.0 (www.mrtrix.org/) and FSL (FMRIB Software Library) version 5.0.11 (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki), as detailed by Oldham et al. (56). Briefly, the diffusion images for each individual were corrected for eddy-induced current distortions, susceptibility-induced distortions, intervolume head motion, outliers in the diffusion signal (57), within-volume motion (58), and B1 field inhomogeneities (59). Tractography was conducted using the fiber orientation distributions (iFOD2) algorithm, as implemented in MRtrix3 (60), which uses fiber orientation distributions estimated for each voxel using constrained spherical deconvolution, which can improve the reconstruction of tracts in highly curved and crossing fiber regions (61, 62). Streamline seeds were preferentially selected from areas where streamline density was underestimated with respect to fiber density estimates from the diffusion model (63). We used anatomically constrained tractography to further improve the biological accuracy of streamlines (64). To create a SC matrix, streamlines were assigned to each of the closest regions in the parcellation within a 5-mm radius of the streamline endpoints (65), yielding an undirected 82 × 82 connectivity matrix. A group-representative connectome was generated by selecting 30% of edges showing the lowest coefficient of variation (based on the streamline count) across participants (66). Threshold selection is a challenging problem in human connectomics (67), and the correct connection density depends on numerous factors, such as the sensitivity/specificity of the tractography method used, the spatial resolution of the parcellation, and the size of the human brain (relative to other species) (68). An important goal for future work will be to understand how different parcellation and thresholding strategies affect model fits.

fMRI processing and FC estimation. The fMRI data were processed as detailed in (55). Briefly, the functional images were initially preprocessed in FSL FEAT following a basic pipeline, which included removal of the first four volumes, rigid-body head motion correction, 3-mm spatial smoothing to improve signal-to-noise ratio, and a high-pass temporal filter of 75 s to remove slow drifts (69). The data were then denoised using FSL-FIX, an independent component analysis–based denoising method that uses an automated classifier to identify noise-related components for removal from the data (70). The algorithm was trained on a manually labeled held-out set of 25 individuals scanned with identical imaging parameters. Time courses for noise-labeled components, along with 24 head motion parameters (6 rigid-body parameters, their backward derivatives, and squared values of the 12 regressors), were then removed from the voxelwise fMRI time series using ordinary least squares regression.

Denoised functional data were spatially normalized to the International Consortium for Brain Mapping 152 template in Montreal Neurological Institute (MNI) space using Advanced Normalization Tools (version 2.2.0) (71), via a three-step method: (i) registration of the mean realigned functional scan to the skull-stripped high-resolution anatomical scan via rigid-body registration; (ii) spatial normalization of the anatomical scan to the MNI template via a nonlinear registration; and (iii) normalization of functional scan to the MNI template using a single transformation matrix that concatenates the transforms generated in steps i and ii. Mean time series for each parcellated region were then extracted, and Pearson correlations between each pair of regional time series were used to estimate interregional FC matrices. FCD were also calculated for the empirical data, as outlined below.

Cross-validation sample. To test the generalizability of our model fit estimates, we examined the out-of-sample performance of our models in fMRI data from an independent sample of 1003 individuals drawn from the Human Connectome Project (HCP). The HCP website (www.humanconnectome.org/) provides the full details of the participants, acquisition protocol, and preprocessing of the resting state fMRI data, which were processed using the minimal preprocessing pipeline and ICA-FIX for denoising [see (72)]. Regional time series and FC estimates were generated using the same parcellation.

The BEI dynamic mean field model

We simulated spontaneous neuronal activity using a whole-brain model that couples the local population dynamics of different brain areas through an anatomical interregional connectivity matrix derived from diffusion-weighted MRI data. The local population dynamics of each area are described by a dynamic mean field model, as proposed in (11). The model expresses the activity of large ensembles of interconnected excitatory and inhibitory spiking neurons as a reduced set of dynamical equations describing the population rate activity of coupled excitatory (E) and inhibitory (I) pools of neurons, following the original derivation of Wong and Wang (73). At the local neuronal population level, the excitatory synaptic currents, I(E), are mediated by NMDA receptors and the inhibitory currents, I(I), are mediated by GABAA receptors. Within each brain area, the E and I neuronal pools are mutually connected. The inter-area coupling between local neuronal population dynamics corresponding to two different brain areas i and j is mediated only at the E-to-E level and defined by the scaled SC matrix C. The elements of this matrix, Cij, define the interregional anatomical connectivity with connection weights estimated using tractography of diffusion-weighted MRI data, as described above.

The balanced dynamic mean-field whole-brain model is mathematically described by the following system of coupled differential equationsIi(E)=Iiext+WEI0+w+JNMDASi(E)+GJNMDAjCijSj(E)JiSi(I)(1)Ii(I)=WII0+JNMDASi(E)Si(I)(2)ri(E)=H(E)(Ii(E))=aEIi(E)bE1exp(dE(aEIi(E)bE))(3)ri(I)=H(I)(Ii(I))=aIIi(I)bI1exp(dI(aIIi(I)bI))(4)dSi(E)(t)dt=Si(E)τE+(1Si(E))γri(E)+σvi(t)(5)dSi(I)(t)dt=Si(I)τI+ri(I)+σvi(t)(6)

In the last equations, for each excitatory (E) or inhibitory (I) pool of neurons in each brain area i, ri(E,I) (in hertz) indicates the firing rate, Ii(E,I) (in nanoamperes) denotes the total input current, and Si(E,I) represents the synaptic gating variable. The total input current received by the E and I pools is transformed into firing rates, ri(E,I), by H(E,I), which corresponds to the sigmoidal neuronal response functions defined in (74). In the excitatory (E) and inhibitory (I) neuronal response functions, aE = 310 nC−1 and aI = 615 nC−1 define the gain factors determining the slope of H, bE = 0.403 nA and bI = 0.288 nA are the threshold currents above which the firing rates increase linearly with the input currents, and dE = 0.16 and dI = 0.087 are constants determining the shape of the curvature of H around b, respectively. The synaptic gating variable of excitatory pools, Si(E), is mediated by NMDA receptors with a decay time constant τNMDA = 0.1 s and γ = 0.641, whereas the average synaptic gating in inhibitory pools is mediated by GABA receptors with a decay time constant τGABA = 0.01 s. The local external input impinging in each population is I0 = 0.382 nA and is weighted with WE = 1 for the excitatory pools and with WI = 0.7 for the inhibitory pools. The parameter Iiext captures external stimulation to the excitatory population. In all cases, we set Iiext = 0 for all regions, except when assessing the ignition capacity and response decay profiles of the model, in which case we set Iiext>0, as detailed below. The local recurrent excitation in each excitatory pool is weighted by w+ = 1.4. All excitatory synaptic couplings are weighted by JNMDA = 0.15 nA. In Eqs. 5 and 6, υn denotes uncorrelated standard Gaussian noise with an amplitude of σ = 0.01 nA. All parameters were set as in (11). We have shown previously that the synaptic gating variable of NMDA receptors has a much longer decay time constant (100 ms) than the AMPA receptors, which means that the dynamics of the NMDA gating variable dominate the time evolution of the system, while the AMPA synaptic variable instantaneously reaches its steady state. Our research has shown that we can leave out the contribution of the AMPA receptors and that this greatly reduced approximation does not compromise performance in terms of fitting the empirical data (31). Note that the balanced model maintains the spontaneous rate activity of the local population at the same low level across the whole brain, irrespective of the anatomical coupling. This makes the model even more biophysically realistic, as shown in detail in (11), and improves the model’s dynamical properties and fits to empirical data and improves the encoding capabilities.

One method for introducing heterogeneity to this model is to define the parameters such that each isolated brain area shows asynchronous spontaneous activity with a low firing rate (r(E) ~ 3 Hz and r(I) ~ 9 Hz), as shown in electrophysiological experiments (1214). In this regime, the whole-brain model is balanced. The balance is obtained by adjusting each local feedback inhibition weight Ji for each brain area i such that the firing rate of the excitatory pools ri(E) remains clamped at 3 Hz even when receiving excitatory input from connected areas. This regulation for achieving balance is known as feedback inhibition control (FIC) and is obtained numerically by the method described in (11). This work demonstrated that FIC improves the fitting of empirical resting-state FC and, more crucially, achieves more realistic evoked activity. We refer to this model as the BEI model. Despite the regional heterogeneity introduced by local tuning of the FIC parameter, we use this model as our homogeneous benchmark for comparing to other models, because the FIC adjustments ensure uniform firing rates across all regions and all regional gain parameters are fixed (as detailed below).

The BEI model is fitted to empirical fMRI data by finding the optimal global coupling factor G, which uniformly scales all inter-area E-to-E connections. The parameter G is the only control parameter that is adjusted to move the homogeneous model to its optimal working point, where the simulated activity maximally fits the empirically observed resting-state FC architecture. In particular, we fit static edge-level FC, node-level FC, and FCD, by choosing the coupling factor G such that the correlation between the empirical and simulated FC is high and the distance between the data and model FCD distributions is minimal (i.e., the minimal DKS between empirical and simulated FCD; see below).

Simulations were run for a range of G between 0 and 3 with an increment of 0.01 and with a time step of 0.1 ms. For each value of G, we averaged 389 simulations of 464.464 s each, which were then put through a BOLD forward model (see below) to emulate the empirical resting-state data, which comprised hemodynamic recordings in 616 volumes acquired in each of 389 human participants. We optimized the coupling factor G before introducing regional heterogeneity to the model using the gain and scale parameters (i.e., B = 0 and Z = 0 during this fitting procedure; see below).

Simulating BOLD signals

Regional BOLD signals are simulated in the model using the generalized hemodynamic model in (75). We calculate the BOLD signal in each brain area i from the simulated firing rate of the excitatory pools ri(E). In this hemodynamic model, an increase in the firing rate causes an increase in a vasodilatory signal, si, that is subject to autoregulatory feedback. Blood inflow, fi, responds in proportion to this signal, inducing changes in blood volume, vi, and deoxyhemoglobin content, qi. The differential equations describing the hemodynamic model coupling these biophysical variables aredsi(t)dt=0.5riE+3kisiγi(fi1)(7)dfi(t)dt=si(8)τidvi(t)dt=fivi1/α(9)τidqi(t)dt=fi(1(1ρi)1/fi)ρiqivi1/αvi(10)

In these equations, ρ denotes the resting oxygen extraction fraction, τ is a time constant, and α represents the resistance of the veins. To compute in each area i, the BOLD signal Bi, we calculate a volume-weighted sum of extra- and intravascular signals, which comprises a static nonlinear function of volume, vi, and deoxyhemoglobin content, qi, and is expressed as followsBn=V0[(k1(1qi)+k2(1qivi)+k3(1vi)](11)

We adopted the biophysical parameters recommended in (75). To focus on the most functionally relevant frequency range in resting-state conditions, we band-pass–filtered both empirical and simulated BOLD signals (0.08 > f > 0.008 Hz).

Introducing regional heterogeneity

One of the key drivers of local population activity in the dynamical model is the balance between excitatory and inhibitory activity. E:I balances across brain regions are thought to play a critical role in synchronized dynamics, and their disruptions have been implicated in diverse diseases (76, 77). In the absence of accurate density estimates for each cell type in different regions of the human brain, regional variations in E:I balance could, in principle, be parameterized in many different ways. Here, we use data on regional expression levels of excitatory and inhibitory receptors to constrain our model of E:I heterogeneity across the brain. To this end, we use data from the AHBA, which comprises microarray data quantifying the transcriptional activity of >20,000 genes in >4000 different tissue samples distributed throughout the brain, taken from six postmortem samples (26, 30).

The AHBA data were processed according to the pipeline developed by Arnatkevičiūtė et al. (29). Briefly, probe-to-gene annotations were updated using the Re-Annotator toolbox (78), resulting in the selection of 45,821 probes corresponding to a total of 20,232 genes. Tissue samples derived from brainstem and cerebellum were removed. Then, the probes that did not exceed background noise in more than 50% of remaining samples were removed, excluding 13,844 probes corresponding to 4486 genes. A representative probe for each gene was selected on the basis of the highest mean intensity across regions (79). To increase the anatomical accuracy of sample-to-region matching, the samples were first divided into four separate groups on the basis of their location—hemisphere (left/right) and structure assignment (cortex/subcortex)—and then assigned to the corresponding parcellation-defined regions of interest separately for each of the four groups. Donor-specific gray matter parcellations were generated, and samples located within 2 mm of the parcellation voxels were mapped to the closest voxel assigning almost 90% of all tissue samples. We retained only data for the left cortex, as anatomical coverage in the AHBA is more complete for this hemisphere. Last, gene expression measures within a given brain were normalized using a scaled robust sigmoid normalization for every gene across samples. Normalized expression measures in samples assigned to the same region were averaged and aggregated into a region × gene matrix consisting of expression measures for 15,745 genes over 34 left hemisphere regions. Given the lack of notable interhemispheric differences in regional gene expression identified in the AHBA (26), the left hemisphere values were reflected on to the right hemisphere to enable whole-brain model simulations.

To capture regional heterogeneity in E:I balance, we extracted expression measures specifically for genes coding for the excitatory AMPA and NMDA receptors and inhibitory GABAA receptor isoforms and subunits. We excluded genes encoding GluN3A and GluN3B subunits because of their low expression levels in adult cortex (80, 81). More specifically, the genes selected for AMPA were GRIA1, GRIA2, GRIA3, and GRIA4; the genes for NMDA were GRIN1, GRIN2A, GRIN2B, and GRIN2C; and the genes for GABAA were GABRA1, GABRA2, GABRA3, GABRA4, GABRA5, GABRB1, GABRB2, GABRB3, GABRG1, GABRG2, and GABRG3. To increase the accuracy of gene expression measures, probes representing the expression of individual genes were selected on the basis of the highest correlation to RNA sequencing data in two of the six donor brains (79).

In our implementation, we quantify regional E:I balance simply by summing all expression values for AMPA and NMDA, summing all values for GABA, dividing the former by the latter, and then normalizing the values to the unit interval. Specifically, for each brain region i, we estimate E:Ii=E:IiE:IE:IminE:ImaxE:ImaxE:Imin+1, where E:Ii = (∑ NMDA + ∑ AMPA)/∑ GABA and E:Imax and E:Imin the respective maximum and minimum E:I ratios observed across regions. We use the same normalization for all other sources of heterogeneity. We note that estimating E:I balance from gene expression data is challenging given that each receptor is a heteromer, consisting of differing subunit compositions (AMPAR: heterotetrameric assembly from 4 different subunits; NMDAR: heterotetrameric assembly from 7 different subunits; and GABA-AR: homo- and heteropentameric assembly from 19 subunits) (80, 8284). The exact combination of these subunits can alter receptor properties, including biophysical, pharmacological, and signaling attributes. For example, the GluN2 subunit is essential for determining the gating properties of NMDA receptors, with different subunit compositions affecting excitatory postsynaptic decay times (84). Here, we have adopted a pragmatic approach and chosen genes encoding the major subunit types present in the human adult cortex for each receptor, assuming that the relative expression levels provide an adequate estimate of the E:I ratio within a given cortical region.

We used the transcriptomically constrained local E:I ratio Ri to modulate the gain of the local neuronal response functions of the corresponding excitatory and inhibitory pools of each brain region, i, following a number of experimental and theoretical studies (85, 86). Specifically, we consider that Ri modulates the gain of the neuronal response function H(E,I) in each brain area according to the modified Eqs. 3 and 4ri(E)=H(E)(IiE)=Mi(aEIi(E)bE)1exp(dEMi(aEIi(E)bE))(12)ri(I)=H(I)(IiI)=Mi(aIIi(I)bI)1exp(dIMi(aIIi(I)bI))(13)

Similar to (22), we assume that the heterogeneity factors (in our case, the transcriptomically defined local excitation/inhibition ratios Ri) are linearly transformed by two unknown free parameters corresponding to a bias B and a scaling factor Z. Thus, in Eqs. 12 and 13, Mi = 1 + B + ZRi. In other words, the effect of the transcriptomically constrained E:I ratio Ri is to modulate locally the slope of the corresponding excitatory and inhibitory neuronal response functions, H. Note that the slope of H is determined by M using the input-output function of Abbott and Chance (74). We studied exhaustively the two free parameters B and Z regulating the influence of transcriptional heterogeneity. Homogeneous conditions correspond to cases where Z = 0. Note that the balancing FIC procedure is applied after fixing the gain scaling, i.e., for a specific bias B and a specific scaling factor Z. Also note that although we constrained the model using a transcriptional index of E:I balance, we used this constraint to modify the gain of each population model. We did this under the assumption that the population gain function, capturing net regional excitability, will be sensitive to regional variations in E:I balance and to limit the number of fitted parameters in the model. More specifically, we fitted only two parameters—B and Z—to model regional heterogeneity, whereas other approaches have fitted four or more (22, 23). We adopted this strategy to enable a more comprehensive search of parameter space, which is important to identify a unified working point for all outcome properties of the model. A more direct modulation of E:I weights might yield different results.

We additionally evaluated two other forms of biologically constrained heterogeneity: the PC1 and T1w:T2w models. The PC1 model was also informed by AHBA data. In this case, we identified the dominant mode of spatial variation in gene expression as the first principal component (PC1) of 1926 genes that overlapped between our final quality-controlled gene list and those previously deemed to be brain specific in (28), where it was shown that this mode correlates with T1w:T2w. Regional loadings on this PC were mapped to a positive range before they were incorporated into the model.

The T1w:T2w model incorporated heterogeneity by directly relying on empirical estimates of the T1w:T2w ratio, which is sensitive to myelin content (27). Prior work has shown that this measure can also be used to map cortical hierarchies (28) and that large-scale heterogeneous biophysical models constrained by T1w:T2w more accurately reproduce empirical FC properties than homogeneous models (22). Regional T1w:T2w values were obtained from the HCP dataset, averaged over 1200 participants following 2 mm FHWM Gaussian smoothing in grayordinates space (https://balsa.wustl.edu/study/show/RVVG), and then parcellated using the Desikan-Killany atlas (35).

Last, to further ensure that any results obtained for the E:I model cannot be explained by a generic spatial gradient, we also evaluated model performance after permuting the assignment of E:I values to regions. A total of 10,000 permutations were performed by spatially rotating the expression values using the approach described in (87). Critically, this approach preserves the spatial autocorrelation of the expression values, which is essential for ensuring that any resulting effects are not driven by low-order spatial gradients [see also (32)]. Thus, this method allows us to distinguish the effects of actual E:I heterogeneity from generic spatial gradients with the same autocorrelation structure [see also (32)]. In all cases, we normalized the corresponding heterogeneity factor to the [0,1] interval and treat it similarly to the ratio Ri.

Evaluating model performance

Models were fitted to three empirical properties of the FC data, as detailed in the following. To compare model performance, we simulated 1000 runs of each model at its own optimal parameters to obtain a distribution of fit statistics that captures the inherent stochasticity of the model. To evaluate out-of-sample performance, we repeated the same analyses using the same optimal model working points in an independent sample of 1003 individuals drawn from the HCP (fig. S2).

Static edge-level FC. The static edge-level FC is defined as the N × N matrix of BOLD signal correlations between brain areas computed over the entire recording period (see Fig. 1B). We computed the empirical FC for each human participant and for each simulated trial (the total number of trials matched the number of participants). The group-averaged empirical and simulated FC matrices were compared by computing the Pearson correlation between their upper triangular elements (given that the FC matrices are symmetric).

Static node-level FC. Node-level FC, or node strength, is another static spatial measure that characterizes the average FC strength for each area (see Fig. 1E) (68). It has also been called global brain connectivity in previous work (22). Thus, node-level FC strength is defined asFCSi=1Nj=1NFCij(14)

We computed the empirical strength vectors for each human participant and for each simulated trial (total number of trials matched the number of participants). The node FC fit is quantified by computing the Pearson correlation between the group-averaged empirical and simulated strength vectors.

Functional connectivity dynamics. The spatiotemporal dynamics of resting-state activity are captured by the statistics of how FC computed over sliding windows evolves across time (see Fig. 1E). Here, we computed the FC over a sliding window of 80 TRs (corresponding approximately to 1 min) with incremental shifts of 18 TRs. To quantify the spatiotemporal statistics of time-evolving FC, we calculated a time-versus-time matrix, called FCD (25). This FCD matrix is defined so that each entry, FCD(tx, ty), corresponds to the correlation between the FC centered at times tx and the FC centered at ty. The statistical distribution of FCD values captures the spatiotemporal architecture of the recording session. To quantitatively compare the spatiotemporal dynamical characteristics between empirical data and model simulations, we generate the distributions of the upper triangular elements of the FCD matrices over all participants and of the FCD matrices obtained from all simulated trials for a given parameter setting. The similarity between the empirical and model FCD distributions is then compared using the DKS, allowing for a meaningful evaluation of model performance in predicting the changes observed in dynamic resting-state FC. Note that FCD represents a more stringent benchmark for model performance than static FC, as it requires the model to capture the statistics of how FC varies at each edge over time.

Ignition capacity. In the literature, ignition has been described using multiple complementary perspectives as comprising three key properties: (i) a rapid, nonlinear amplification of activity once stimulus intensity has surpassed a certain threshold; (ii) widespread activation of spatially distributed systems, particularly frontoparietal areas; and (iii) sustained reverberatory activity persisting after termination of the stimulus (4, 88). Our model provides insights into the first two properties but has reduced sensitivity for capturing reverberatory activity because of its focus on noisy fluctuations and the incorporation of local feedback inhibition. In our analysis, we therefore aimed to quantify ignition-like dynamics in our model using a metric that captures elements of both nonlinear amplification and widespread activation. More specifically, we quantify the ignition capacity of each model by examining how regional dynamics respond to a simulated focal perturbation. First, for a given working point, we run the whole-brain model for 3000 ms under spontaneous conditions to eliminate transients. Second, from 3000 to 4000 ms, we apply a stimulation of strength Iiext to the bilateral lateral occipital cortex region in the Desikan-Killiany atlas (Fig. 3A), which contains the foveal representation of V1 and the confluence of other visual areas. Third, from 4000 to 7000 ms, we reset the external stimulation to zero again and let the system relax again to the spontaneous state. For each condition, i.e., for each parameter setting, and each particular stimulation Iiext, we run 30 trials and compute the average, across trials, of the elicited activity as a function of time on all brain areas, to smooth the resulting neuronal whole-brain activity.

To compare ignition capacity across the whole cortex, we fix the location i where the stimulation is applied (and also stimulate the region’s contralateral homolog). Then, we simulate the averaged neuronal response elicited in a given brain area as a consequence of the applied artificial stimulation and as a function of the strength Iiext (we modify Iiext from 0 to 0.2 in steps of 0.001). For each brain region, we perform a sigmoidal fitting of the elicited averaged rate activity in the period from 3500 to 4000 ms as a function of the applied strength Iiext in the stimulated specific site i. The ignition elicited for a given nonstimulated target brain area after the application of the artificial stimulation in the specific region i is defined by the product between the level of concavity of the neuronal response function elicited on that target area and the maximal rate elicited when the maximal strength of stimulation is applied (as fitted by the sigmoidal function). The level of concavity is defined by the maximal value of the second derivative of the fitted sigmoidal function. Figure 1 (F and G) shows graphically the key concept. The global ignition capacity is the average ignition capacity across all brain areas. The ignition measure is basically describing two effects, namely, the fidelity (captured by the maximal elicited rate in the target area) and explosiveness (captured by the level of concavity) of ignition-like activity. In secondary analyses, we used the same approach to simulate the effects of somatosensory stimulation by applying the stimulation to the postcentral gyrus region of the Desikan-Killiany atlas (fig. S3).

Activity time scales. We follow Chaudhuri et al. (15) and examine activity time scales in two ways. First, we examine the rate at which regional activity decays following the same impulse stimulation paradigm used in the analysis of ignition capacity. More specifically, we compute the temporal decay of each target brain region following simulation of the lateral occipital region in both hemispheres. We calculate the decay rate by fitting an exponential curve to the regional firing rates after stimulation, averaged over 1000 simulations (the length of total simulation as described above). The functional form of the exponential curve is riE=A(exp(Dt)+B), where riE is the firing rate (in hertz) after stimulation, t is the time after stimulation (in seconds), and the free parameters of decay rate, D, scaling constant, A, and offset, B, are fitted using nonlinear least squares.

Second, we examine intrinsic time scales by estimating each region’s ACF measured during an input of continuous white noise to the same stimulation target (lateral occipital cortex). This approach more closely mimics empirical characterizations of intrinsic time scales [e.g., (17, 36, 37)]. To this end, we provide white noise input currents with a mean of 0.356 nA and an SD of 0.05 nA, as in the work of Chaudhuri et al. (15). We set the noise levels of all other regions to be low, where the noise level in each model was chosen as to maximize the range of time scales with no stimuli present. This resulted in noise levels of 10−8,10−5,10−8,10−4, and 10−8 in the E:I, principal components analysis, T1w:T2w, BEI, and spatial null models, respectively. For each model, we calculated the average ACF across 30 simulations. To calculate the time scale estimates from these average ACFs, we truncated the ACF to be above 0.05 and fitted the curve with either a single exponential, ACF(t = lag) = exp(−t/D), or a double exponential, AC(t = lag) = A[exp(−t/D1)] + B[exp(−t/D2)] to account for the longer tails in some ACFs. In these equations, t is the lag time (ms), D, D1, and D2 are decay constants (ms−1), and A and B are weightings of the exponential functions. Following Chaudhuri et al. (15), if the sum of squared errors in the double exponential fit was eight times larger than in the single exponential, then we used the sum of the two time constants (1/D1,1/D2) from the double exponential fit, with each weighted by its amplitude (A and B, respectively); otherwise, the time scale was calculated as 1/D. These simulations and ACF calculations were repeated 1000 times for each model to build up distributions of the time scale estimates. We then characterized time hierarchy variability as the interquartile range of intrinsic ACF time scales across brain areas. A larger interquartile range implies greater variability in intrinsic time scales across brain regions, which is consistent with a broader dynamic range of activity and a more pronounced hierarchical ordering of brain areas.

Statistical analyses

Differences in model fits to empirical properties, as well as in ignition capacity and decay variability, were assessed using pairwise Wilcoxon rank sum tests, Bonferroni-corrected for 10 possible comparisons (corrected P values are denoted Pbonf). Spatial correlations between regional maps were assessed with reference to empirical null distributions generated by rotating one of the spatial maps using the same procedure as used when generating the spatial null models (89). The corresponding P values are denoted Pspatial.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/7/29/eabf4752/DC1

https://creativecommons.org/licenses/by-nc/4.0/

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: Funding: G.D. is supported by a Spanish national research project (ref. PID2019-105772GB-I00 MCIU AEI) funded by the Spanish Ministry of Science, Innovation and Universities (MCIU), State Research Agency (AEI); HBP SGA3 Human Brain Project Specific Grant Agreement 3 (grant agreement no. 945539), funded by the EU H2020 FET Flagship programme; SGR Research Support Group support (ref. 2017 SGR 1545), funded by the Catalan Agency for Management of University and Research Grants (AGAUR); Neurotwin Digital twins for model-driven non-invasive electrical brain stimulation (grant agreement ID: 101017716) funded by the EU H2020 FET Proactive programme; euSNN European School of Network Neuroscience (grant agreement ID: 860563) funded by the EU H2020 MSCA-ITN Innovative Training Networks; CECH The Emerging Human Brain Cluster (Id. 001-P-001682) within the framework of the European Research Development Fund Operational Program of Catalonia 2014-2020; Brain-Connects: Brain Connectivity during Stroke Recovery and Rehabilitation (id. 201725.33) funded by the Fundacio La Marato TV3; Corticity, FLAG–ERA JTC 2017, (ref. PCI2018-092891) funded by the Spanish Ministry of Science, Innovation and Universities (MCIU), State Research Agency (AEI). M.L.K. is supported by Center for Music in the Brain, funded by the Danish National Research Foundation (DNRF117); and Centre for Eudaimonia and Human Flourishing, funded by the Pettit and Carlsberg Foundations. AF was supported by the Sylvia and Charles Viertel Charitable Foundation and National Health and Medical Research Council (ID: 3251549). This work was supported by the computational infrastructure provided by the MASSIVE HPC facility (www.massive.org.au). Author contributions: G.D., M.L.K., K.M.A., and A.F. conceived the project and wrote the manuscript. G.D., K.M.A., M.L.K., A.A., S.O., N.C.R., and K.S. processed and analyzed the data and provided feedback on the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data and code needed to evaluate the conclusions in this paper are present in the paper and the Supplementary Materials and are available at https://github.com/KevinAquino/HNM. As per ethics guidelines, raw neuroimaging data are not available.

Stay Connected to Science Advances

Navigate This Article