The impact of regional heterogeneity in whole-brain dynamics in the presence of oscillations

The collective behavior of human brain emerges from the nonlinear dynamics of millions of neurons interacting at billions of synaptic connections (Sporns et al., 2005). Nevertheless, despite this microscopic scale complexity, the brain spontaneously self-organizes into large-scale neuronal networks that give rise to a reduced number of discrete states, where the resting state in particular has been proposed to be a fundamental brain state (Fox, 2005; Gusnard & Raichle, 2001; Raichle & Mintun, 2006). Local and global recordings of the brain have been performed to understand the relation between the neuronal responsiveness and these complex spatiotemporal patterns (Destexhe & Contreras, 2006). In addition, modeling large-scale brain dynamics has been established as a crucial path to unveil mechanisms underlying the emergence of global brain behavior as the interactions of small dynamical units (Deco et al., 2009; Ghosh et al., 2008; Honey et al., 2007).

Whole-brain models based on conceptually simple local dynamical rules coupled according to empirical measurements of anatomical connectivity (Hagmann et al., 2010; Hagmann et al., 2007; Johansen-Berg & Rushworth, 2009) have been successfully implemented to explore the interplay between local dynamics, long-range structural coupling, and the formation of large-scale activity patterns in empirical data (Deco, Cabral, et al., 2017; Deco, Kringelbach, et al., 2017; Ipiña et al., 2020; Jirsa et al., 2002; Piccinini et al., 2021). The local dynamics are described by a set of equations that can be built from knowledge based on biologically realistic mechanisms underlying brain activity, or from neurobiological assumptions concerning the collective dynamics that an ensemble of small units is able to produce (Breakspear, 2017; Cofré et al., 2020; Deco et al., 2008). Several whole-brain models have been proposed based on different local dynamics: from the most abstract such as cellular automata (Haimovici et al., 2013; Tagliazucchi et al., 2016), Ising spin-models (Marinazzo et al., 2014; Ponce-Alvarez et al., 2018), or nonlinear oscillators (Deco, Kringelbach, et al., 2017; Saenger et al., 2017), up to more descriptive models such as neural masses (Breakspear et al., 2003; Honey et al., 2009) or mean-field models (Deco et al., 2014; di Volo et al., 2019; Goldman et al., 2020). However, different sources of empirical data have revealed local variations within each region (Huntenburg et al., 2018), and during the past years whole-brain models based on mean-field local models have been demonstrated to be suitable to unveil how these heterogeneities impact global brain dynamics (Deco et al., 2018; Deco, Kringelbach, et al., 2021; Demirtaş et al., 2019; Kong et al., 2021; P. Wang et al., 2019).

Specifically, the so-called dynamic mean-field whole-brain model (DMFM), as proposed in Deco et al. (2014), was implemented to demonstrate the dynamical consequences of local heterogeneity in reproducing the spatiotemporal structure of empirical data (Demirtaş et al., 2019), its ignition capacity (Deco, Kringelbach, et al., 2021), as well as the role of neuromodulation in altered states of consciousness (Deco et al., 2018). In this model the local activity is represented by a set of differential equations describing the interaction between inhibitory and excitatory pools of neurons (adapted from Wong & Wang, 2006). For each population, three variables are modeled: the synaptic current, the firing rate, and the synaptic gating, where the NMDA and gamma-aminobutyric acid (GABA)-A play a role in the excitatory and inhibitory coupling, respectively. This model is widely used to represent the asynchronous state of global behavior of large ensembles of interconnected neurons, which is guaranteed by the inclusion of a feedback inhibitory control (Deco et al., 2014). Nevertheless, how the heterogeneity shapes whole-brain activity in synchronous states, which is a ubiquitous phenomenon in brain activity, supported by the presence of oscillatory local dynamics is still unclear.

Here, we investigated the dynamical consequences of including regional heterogeneity when brain dynamics are supported by a synchronous state instead of an asynchronous state, as was reported by studies mentioned above. For this purpose, we assessed the impact of regional heterogeneities in the presence of oscillations by implementing whole-brain models capable of representing those synchronous oscillating states. Briefly, we used two models with different levels of abstraction that include a Hopf bifurcation in their local dynamical landscape, which is the simplest representation of both regimes by the same set of equations. The first model is the well-known Stuart–Landau (SL) whole-brain model comprising a set of nonlinear oscillators represented by the normal form of a Hopf bifurcation (Deco, Kringelbach, et al., 2017; Jobst et al., 2017; Saenger et al., 2017; Senden et al., 2017; Wiggins, 2003). The second model is an exact mean-field whole-brain model, which describes the mean activity of a local population of all-to-all coupled quadratic integrate-and-fire (QIF) neurons (Montbrió et al., 2015). Specifically, we fitted both whole-brain models to the large-scale neuroimaging resting-state fMRI data from healthy human participants in the Human Connectome Project (HCP; Van Essen et al., 2013). We demonstrated that both models more faithfully reproduce empirical properties when including a plausible biological source of local heterogeneity, such as regional variation in the ratio of T1- to T2-weighted (T1w/T2w) MRI signal, which is thought to index intracortical myelin content (Glasser & Van Essen, 2011). We tested that this inclusion improves different empirical properties such as the causal ignition (CI), defined as the causal effect that the dynamics of one brain region has on the whole-brain level of synchronization proposed as a quantification of the hierarchical organization of brain function (Deco, Vidaurre, et al., 2021) and static metrics as the functional global brain connectivity (GBC; Deco, Kringelbach, et al., 2021; Demirtaş et al., 2019). Finally, we applied the same framework fitting the SL whole-brain model to neuroimaging resting-state fMRI data from patients with diffuse brain volume heterogeneities (atrophy) triggered by neurodegeneration: Alzheimer’s disease (AD). We found that the inclusion of specifical functional local heterogeneity improves the capacity to reproduce empirical properties, such as functional connectivity, causal ignition, and global brain connectivity.

The dataset used for this investigation was selected from the March 2017 public data release from the Human Connectome Project (HCP) where we chose a sample of 1,000 participants in resting state. The full informed consent from all participants was obtained by the Washington University–University of Minnesota (WU-Minn HCP) Consortium, and research procedures and ethical guidelines were followed in accordance with Washington University institutional review board approval.

The study comprised 39 patients diagnosed with AD (26 female, 76.6 ± 7 years [mean ± SD]). Additionally, a group of 35 healthy controls (38 female, 69.8 ± 7.9 years) was used to compute the functional regional heterogeneity map. An expert neurologist following current criteria diagnosed patients following NINCDS-ADRDA clinical criteria for AD (Dubois et al., 2007; McKhann et al., 2011). Participants were recruited in clinical sites by a multidisciplinary team as part of an ongoing multicentric protocol taking part in the Multi-Partner Consortium to Expand Dementia Research in Latin America (ReDLat) and assessed following harmonized procedures (Ibañez, Fittipaldi, et al., 2021; Ibañez, Pina-Escudero, et al., 2021) as in previous works (Donnelly-Kehoe et al., 2019; Legaz et al., 2022; Melloni et al., 2016; Salamone et al., 2021; Sedeño et al., 2017). Clinical diagnoses were established by neurodegenerative disease experts through an extensive neurological, neuropsychiatric, and neuropsychological examination comprising semi-structured interviews and standardized cognitive assessments. All participants with neurodegenerative conditions were in early/mild stages of the disease. No participant reported a history of other neurological disorders, psychiatric conditions, primary language deficits, or substance abuse. Whole-brain gray matter was compared between AD and controls, showing an extended bilateral temporal with less extended frontoparietal atrophy in AD (Du et al., 2007; Landin-Romero et al., 2017; Pini et al., 2016). All participants provided written informed consent pursuant to the Declaration of Helsinki. The study was approved by the ethics committees of the involved institutions.

All neuroimaging data from HCP were processed using two standard cortical parcellations. For a fine-scale parcellation, we used the Glasser parcellation with 360 cortical regions (180 regions in each hemisphere; Glasser et al., 2016). For a coarser scale parcellation, we used the Desikan–Killiany parcellation (Desikan et al., 2006) with a total of 68 cortical regions (34 regions per hemisphere).

The HCP website (https://www.humanconnectome.org/) provides the full details of participants and the acquisition protocol of the data for resting state. We used one resting state fMRI acquisition of approximately 15 minutes acquired on the same day of 1003 HCP participants scanned on a 3-T connectome-Skyra scanner (Siemens).

The participants’ data were acquired in three-dimensional volumetric and 10-min-long resting-state MRI sequences. Participants were instructed not to think about anything in particular while remaining still, awake, and with eyes closed. Two independent centers recorded the data, using the parameters described below.

Center 1 (Argentina): Using a 3-T Phillips scanner with a standard head coil, we acquired whole-brain T1-rapid anatomical 3D gradient echo volumes, parallel to the plane connecting the anterior and posterior commissures, with the following parameters: repetition time (TR) = 8,300 ms; echo time (TE) = 3,800 ms; flip angle = 8°; 160 slices, matrix dimension = 224 × 224 × 160; voxel size = 1 mm × 1 mm × 1 mm. Also, functional spin echo volumes, parallel to the anterior-posterior commissures, covering the whole brain, were sequentially and ascendingly acquired with the following parameters: TR = 2,640 ms; TE = 30 ms; flip angle = 90°; 49 slices, matrix dimension = 80 × 80 × 49; voxel size in plane = 3 mm × 3 mm × 3 mm; slice thickness = 3 mm; sequence duration = 10 min; number of volumes = 220. A total of 18 AD patients and 23 controls were scanned in this center.

Center 2 (Chile): Using a 3-T Siemens Skyra scanner with a standard head coil, we acquired whole-brain T1-rapid gradient echo volumes, parallel to the plane connecting the anterior and posterior commissures, with the following parameters: TR = 2,400 ms; TE = 2,000 ms; flip angle = 8°; 192 slices, matrix dimension = 256 × 256 × 192; voxel size = 1 mm × 1 mm × 1 mm. Finally, functional EP2D-BOLD pulse sequences, parallel to the anterior-posterior commissures, covering the whole brain, were acquired sequentially intercalating pair-ascending first with the following fMRI parameters: TR = 2,660 ms; TE = 30 ms; flip angle = 90°; 46 slices, matrix dimension = 76 × 76 × 46; voxel size in plane = 3 mm × 3 mm × 3 mm; slice thickness = 3 mm; sequence duration = 10.5 min; number of volumes = 240. A total of 21 AD patients and 34 controls were scanned in this center.

The preprocessing of the HCP resting state is described in detail on the HCP website. In short, the data are preprocessed using the HCP pipeline, which is using standardized methods of FSL (FMRIB Software Library), FreeSurfer, and the Connectome Workbench software (Glasser et al., 2013; Smith et al., 2013). This preprocessing included correction for spatial and gradient distortions and head motion, intensity normalization and bias field removal, registration to the T1-weighted structural image, transformation to the 2-mm Montreal Neurological Institute (MNI) space, and using the FIX artifact removal procedure (Schröder et al., 2015; Smith et al., 2013). The head motion parameters were regressed out and structured artifacts were removed by ICA + FIX processing (independent component analysis followed by FMRIB’s ICA-based X-noiseifier (Griffanti et al., 2014; Salimi-Khorshidi et al., 2014). Preprocessed time series of all grayordinates are in HCP CIFTI grayordinate standard space and available in the surface-based CIFTI file for each participant for resting state.

We used a custom-made MATLAB script using the ft_read_cifti function (Fieldtrip toolbox; Oostenveld et al., 2011) to extract the average time series of all the grayordinates in each region of the Schaefer parcellation, which are defined in the HCP CIFTI grayordinate standard space. Furthermore, the BOLD time series were transformed to phase space by filtering the signals with the corresponding TR = 0.72 s in the range between 0.008 Hz and 0.08 Hz (Fox, 2005) and the low-pass cutoff to filter the physiological noise, which tends to dominate the higher frequencies (Cordes et al., 2001; Fox, 2005).

The first five volumes of each fMRI resting-state recording were discarded to ensure a steady state. Data Processing Assistant for Resting-State fMRI (DPARSF V2.3) was used to preprocess the images; DPARSF is an open-access toolbox that generates and implements an automatic pipeline for fMRI analysis within SPM12 and the Resting-State fMRI Data Analysis Toolkit (REST V.1.7). Preprocessing steps included slice-timing correction and realignment to the first scan of the session to correct head movement. Using least squares, we regressed out six motion parameters, as well as cerebrospinal fluid and white matter signals to attenuate the potential effects of residual movement and physiological noise. For this purpose, motion parameters were estimated during the realignment step, and cerebrospinal fluid and white matter masks were obtained from the tissue segmentation of each subject’s T1 scan in native space. None of the participants presented head movements larger than 3 mm and/or rotations higher than 3°, and no differences in head motion among groups were found. As a final step, images were normalized to common MNI space and smoothed using an 8-mm full-width-at-half-maximum isotropic Gaussian kernel. Furthermore, the BOLD time series were transformed to phase space by filtering the signals with the corresponding TR = 2.65 s (rounded for two sites) in the range between 0.008 Hz and 0.08 Hz (Fox, 2005) and the low-pass cut-off to filter the physiological noise, which tends to dominate the higher frequencies (Cordes et al., 2001; Fox, 2005).

To obtain the structural connectivity we use the Human Connectome Project (HCP) database, which contains diffusion spectrum and T2-weighted imaging data from 32 participants (acquisition parameters are described in detail on the HCP website (Setsompop et al., 2013). Briefly, the data were processed using a generalized q-sampling imaging algorithm implemented in DSI studio (https://dsi-studio.labsolver.org). Segmentation of the T2-weighted anatomical images produced a white matter mask and coregistering the images to the b0 image of the diffusion data using SPM12. In each HCP participant, 200,000 fibers were sampled within the white matter mask. Fibers were transformed into MNI space using Lead-DBS (Horn & Blankenburg, 2016). We used the standardized methods in Lead-DBS to produce the structural connectomes for both the Glasser 360 parcellation (Glasser et al., 2016) and the Desikan–Killiany 68 parcellation (Desikan et al., 2006), where the connectivity has been normalized to a maximum of 0.2. The freely available Lead-DBS software package (https://www.lead-dbs.org/) provides the preprocessing that we implemented and is described in detail in Horn et al. (2017).

The coordinated dynamics of the resting-state activity are modeled by introducing coupling between these oscillators. Several works have demonstrated that whole-brain models based on Stuart–Landau oscillators ruling the local dynamical behavior have the capability to describe different empirical observables from brain recordings when the oscillators are coupled with the structural connectivity (Deco, Kringelbach, et al., 2017; Hansen et al., 2015; Jobst et al., 2017).

The equations of the full coupled exact mean-field whole-brain model have five parameters for each population and the coupling strength parameter. Following Devalle et al. (2017), it is possible to nondimensionalize these equations so that the system can be written solely in terms of six parameters. Generally, we adopt the following notation: We use capital letters to refer to the original variables and parameters of the exact mean-field model, and lowercase letters for nondimensional quantities.

We define the causal ignition based on the previous work in which Deco and colleagues established a novel framework named normalized directed transfer entropy (NDTE) to characterize how different brain regions communicate with each other in terms of information-theoretical statistical criterion (Deco, Vidaurre, et al., 2021). This framework uses only the second-order statistics of the involved entropies, which means that instead of estimating the probabilities, the method estimates the covariance, which massively facilitates computation. Then, four computations are performed: normalization, multiple time points in the past, circular surrogates, and aggregation of p values to improve the reliability and robustness of the NDTE framework.

Finally, following the NDTE framework, we perform statistical significance analysis of the CI by implementing a surrogate framework. Briefly, we use the circular shift methodology for analyzing the p values of the hypothesis testing, aiming to detect significant values in CI for each region for each single participant (for further details in the surrogate creation and rational, see Deco, Vidaurre, et al., 2021). After computing these individual p values, we aggregate the p values for each brain area across the whole group of participants. The combination of different p values across subjects is addressed by implementing the Fisher’s method (Fisher, 1992). Here, we use a more sensitive methodology, namely, the Stouffer’s method (Stouffer et al., 1949), which sums the inverse normal transformed p values. After the aggregation of the pairs of p values across participants, we correct for multiple comparisons by using the traditional false discovery rate method of Benjamini and Hochberg (Hochberg & Benjamini, 1990). The result of the significance test across participants determines a binary vector M (with the dimension being the number of brain regions in a given parcellation) that indicates with ones or zeros whether the corresponding brain region is significant.

Our overall framework is schematically displayed in Figure 1 and is explained in detail in the Methods section. Briefly, we used whole-brain models that combine different sources of empirical information in order to replicate static and dynamic empirical properties observed in fMRI recordings. In this work, we focused on models where local dynamics are described by a Hopf bifurcation, that is, the capability to change the dynamic behavior from a fixed point governed by noise toward sustained oscillations. Specifically, we investigated two models with different levels of abstraction near the Hopf bifurcation point: a widely investigated phenomenological model built by coupled nonlinear Stuart–Landau oscillators (Deco & Jirsa, 2012; Deco, Kringelbach, et al., 2017; Deco et al., 2014; Freyer et al., 2012; Ghosh et al., 2008; Honey et al., 2007; Figure 1B); and an exact mean-field model, consisting of a system of two coupled ordinary differential equations (Montbrió et al., 2015; Figure 1C). This system of equations—often referred to as firing rate equations (FRE)—exactly describes the dynamics of the mean firing rate (R) and the mean membrane potential (V) of a large, heterogeneous population of all-to-all coupled QIF neurons. The heterogeneity within each population is represented by the half-width Δ and center η of the Lorentzian distribution standing for the neurons’ input currents, which are necessary for the exact derivation of the mean-field model. Consequently, the FRE corresponding to one population depends on only these parameters, and thus the mean firing rate and mean membrane potential dynamics are determined by the level of heterogeneity within the population (Devalle et al., 2017; see the Methods section). We numerically explored the dynamical landscape of an interacting pair of populations, one inhibitory and one excitatory, following the FRE derived by Montbrió et al. (2015), after performing a nondimensionalization of the FRE (Devalle et al., 2017; see the Methods section). For this exploration, we initially considered model parameters based on previous studies (Devalle et al., 2017; Montbrió et al., 2015), and we then tuned these values to obtain the two-dimensional parameter space determined by d_e − d_i (the rescaled half-width of the excitatory and inhibitory population, whereas the other parameters were fixed: η = $4030$⁠; j_ee = $5040$⁠; j_ie = $2040$⁠; j_ie = j_ii = −$2040$⁠) a Hopf bifurcation (representing excitation-inhibition-based oscillations, purple line in the phase diagram in Figure 1C). We included regional heterogeneity to modulate the bifurcation parameter, a, in the Hopf whole-brain model and the parameters d_e and d_i in the exact mean-field whole-brain model. The modulation was the same for both models, consisting of an additive term (referred to as the “bias”), and a term that is a constant multiplied by the regional heterogeneity (referred to as “scale”; Figure 1D).

First, we fitted the homogeneous model, that is, with zero bias and zero scale, to the empirical level of synchronization, quantified as the mean of the modulus of the Kuramoto order parameter (KoP). The quality of fit was defined as the absolute difference between the empirical and the modeled values (ErrKoP). Second, on top of the homogeneous model we exhaustively explored the values of the bias and scale, and determined a curve of minimum level of ErrKoP in the bias-scale parameter space (Figure 1E). Finally, we computed the model fit capacity to other static and dynamical properties along the iso-level curve of minimum ErrKoP. We evaluated the performance of the model with respect to some of the major outcome measures evaluated thus far in the literature; namely, the static FC and the GBC, but also with respect to a measure that quantifies the causal effect that each brain region has on the whole-brain of synchronization. To do so, we considered the dynamics of one brain region at a certain time i and then computed the causal effect that it has on the whole-brain level of synchronization in time i + 1, in terms of Granger causality (Figure 1F; see the Methods section for a more detailed explanation).

We assessed the homogeneous SL whole-brain model performance in reproducing empirical properties of resting-state data in a parcellation comprising 360 cortical areas (Glasser et al., 2016). We explored the global coupling parameter G with all regional bifurcation parameter a_i = −0.02 (i.e., heterogeneity is not considered). Figure 2A shows how well the model fits, as a function of G, by computing the correlation between model and empirical functional connectivity (FCfit) and absolute error of the model and empirical Kuramoto order parameter (ErrKoP) (see Methods for further details). We selected as the homogeneous model working point the minimum of ErrKoP, which shows a clear optimum at G = 3.23, while the FC fitting level asymptotically reaches a maximum.

We then studied how regional heterogeneity given by the T1w/T2w affects the fitting of the Kuramoto order parameter (ErrKoP) and the functional connectivity (FCfit). We also considered the fitting of causal ignition (CI) and global brain connectivity (GBC), both computed as the Pearson correlation between the empirical and simulated measures, referred to as CIfit and GBCfit. The CI metric is a measure of local causality reflecting how the activity of one brain region affects the whole-brain dynamics, while the GBC is a local metric that indicates how each brain region is functionally connected to the rest of the brain. We introduced this heterogeneity by modulating the regional bifurcation parameter a_i, at the optimal working point of the homogeneous model (G = 3.23), by a multiplicative factor, named scale, and additive term, named bias. Figure 2B shows the evolution of the fitting of the four measures as a function of the scale and bias. We identified an iso-level curve of ErrKoP, represented by red stars in the four matrices, computed as the value of bias that reaches the minimum value of ErrKoP for each value of scaling. Figure 2C presents the level of fitting of these four measures along the iso-level curve of ErrKoP starting from zero scale and bias, that is, the homogeneous case. We defined as the heterogeneous model’s optimum working point the scale and bias values where the maximal CIfit is reached within the iso-level ErrKoP curve (i.e., the maximum of the CIfit, purple, curve in Figure 2C). The distributions of fit statistics of CIfit and GBCfit across 50 runs of the homogenous model, heterogeneous model, and null hypothesis generated by spatially shuffling the regional heterogeneity are displayed in Figure 2D.

The spatial shuffling was carefully generated to preserve the spatial autocorrelation of the heterogeneous map, following the work of Burt et al. (2020). In the Supporting Information, Figure S1, we displayed the map and its corresponding surrogates. The heterogeneous model presents a significant improvement in the fitting capacity compared with the homogenous model and shuffled control (P < 0.001, Wilcoxon rank sum test with Bonferroni correction in all the comparisons). In summary, we can conclude that the inclusion of structural regional heterogeneities improves the model fitting beyond the inclusion of more model parameters (the shuffled control presents the same number of parameters as the heterogeneous case without improvements in the model-fitting capacity).

We evaluated the homogeneous exact mean-field model performance in reproducing empirical properties of the resting state. For a coarser scale parcellation suitable for mean-field whole-brain computational modeling, we used a Desikan–Killiany parcellation comprising 68 cortical areas (Desikan et al., 2006). We explored two parameters: the global coupling parameter, G, and the rescaled dispersion of the local level of heterogeneity of inhibitory population, d_i, with the other model parameters equal to all regions (η = $4030$⁠; d_e = 1; j_ee = $5040$⁠; j_ie = $2040$⁠; j_ie = j_ii = −$2040$⁠) to be near the Hopf bifurcation (see Figure 1 and the Methods section). It is remarkable that we have initially explored the parameter corresponding to the inhibitory populations following previous work that reveals the importance of the inhibitory network in facilitating stable macroscopic emergent dynamics in the brain (Hellyer et al., 2016), in providing a mechanism for brain oscillatory dynamics (Devalle et al., 2017), and in compensating the excess of long-range excitatory connections (Deco et al., 2014). To compare with fMRI empirical observables, the simulated firing activity of each excitatory pool of neurons was transformed into BOLD-like signals using the Balloon–Windkessel model (Deco et al., 2014; Stephan et al., 2007). We fitted the model to empirical data by computing the correlation between the model and empirical functional connectivity (FCfit, blue) and absolute error of the model and empirical Kuramoto order parameter (ErrKoP, orange). To do so, we exhaustively explored the parameter space defined by the global coupling parameter, G, and the dispersion of the local heterogeneity level of inhibitory population, d_i (Figure 3A). We selected the minimum value of ErrKoP as the optimal working point of the homogenous model found at G = 1.04 and d_i = 1.175, indicated with red stars in both matrices in Figure 3A. We replicated these results by exploring the level of heterogeneity in the excitatory population (d_e) and fixed the value of the inhibitory population (d_i = 1). The exploration shows a similar structure to fit the synchronicity and functional connectivity (compared with Figure 3A), presenting a similar optimum working point with similar fitting levels (see the Supporting Information, Figure S2). This result is aligned with the similar role these parameters play in the dynamical system, as seen in the diagram with the Hopf bifurcation displayed in Figure 1.

We then explored the impact of local T1w/T2w ratio heterogeneity in this model through the inclusion of two parameters, the scaling and bias, which equally modify d_i and d_e, and computed the GBCfit and CIfit for each bias-scale pair. Following the procedure that we performed with the phenomenological model, we found for each scale the bias that yields the minimum ErrKoP value, and we selected the best heterogeneous fit as the one that maximizes the fitting of CIfit (scale = 2.1 and bias = −0.7). We computed 50 times the CIfit and the GBCfit in the regional heterogeneity optimal working point and compared with the homogeneous model (scale = 0; bias = 0) and a spatial null model (generated by shuffling the regions of heterogeneity, preserving the spatial autocorrelation following Burt et al. [2020]; see the Supporting Information, Figure S1, for a rendering of the map and its surrogates). The boxplots show that the heterogeneous model yields the best fitting performance compared with the other two models in the CIfit and GBCfit (P < 0.001, Wilcoxon rank sum test with Bonferroni correction).

We investigated how the heterogeneous bifurcation parameters are organized in the parameter space with respect to the bifurcation point. We computed the bifurcation parameters for both models shaped by the scaling and bias obtained from the heterogenous fitting optimization procedure described in previous section. We found that the heterogenous bifurcation parameters, a_i, in the SL model are mostly located below the Hopf bifurcation point (Figure 4A, left panel). We also observed that the heterogenous bifurcation parameters, d_e and d_i, in the exact mean-field model are located below the bifurcation curve (Figure 4A, right panel). Interestingly, in both models the heterogenous bifurcation parameters are distributed around the homogeneous optimal working point (green line/point in SL model and exact mean-field model, Figure 4A).

We rendered the values of the heterogenous bifurcation parameters for both models onto a brain cortex. The bifurcation parameter, a, of the SL whole-brain model and the norm of the vector defined by each (d_e; d_i) pair are shown in the left and right panels of Figure 4B. Remarkably, despite the same brain areas in both models shifting toward high values and low values by construction, the distribution with respect to the bifurcation point shows different behavior.

We moved further and investigated the role of disease-specific regional heterogeneities to model Alzheimer’s disease (AD) patients. To do so, we assessed the homogeneous SL whole-brain model performance in reproducing empirical properties of resting-state data of Alzheimer’s in the fine parcellation with 360 brain regions (Glasser et al., 2016). We explored the global coupling parameter G with all regional bifurcation parameters a_i = −0.02. Figure 5A shows how well the model fits as a function of G in terms of FCfit and ErrKoP, and we defined the homogeneous model optimal working point as the minimum of ErrKoP, which shows a clear optimum at G = 2.8.

We then studied how disease-specific regional functional heterogeneity given by the node level of brain connectivity (GBC, see Methods) affects the fitting of the empirical properties. We introduced this heterogeneity as in the previous sections and computed the level of fitting of the FC, GBC, KoP, and a slightly different version of CI due to the sampling rate of patients’ data (see Methods). We exhaustively explored the bias and scale of the heterogeneous model, and in Figure 5B we show the evolution of the fitting of the four measures. As in previous sections, we identified an iso-level curve of ErrKoP, and we identified the optimal working point of the heterogeneous model as the maximum CIfit within the iso-level curve of ErrKoP.

Finally, we compared the performance of the heterogeneous model in fitting the empirical observables with the homogeneous model and the null hypothesis provided by the heterogeneity spatial shuffling (as in the previous sections). We also assessed the impact of including the functional heterogeneity given by the GBC of healthy controls (CNT) in the model-fitting performance of Alzheimer’s disease patients’ data. In Figure 5C, upper panel, we rendered onto a brain the AD-specific regional heterogeneity; in the lower panel, we showed the high correlation between regional functional heterogeneities of patients and controls. Figure 5D shows the statistical difference between the AD-specific heterogeneous model, the CNT-specific heterogeneous model, the homogeneous model, and the spatial shuffling (generated by shuffling the regions of heterogeneity, preserving the spatial autocorrelation following Burt et al. [2020]; see the Supporting Information, Figure S1, for a rendering of the map and its surrogates) for 50 runs of each model (P < 0.001, Wilcoxon rank sum test with Bonferroni correction). We extended these results to a different fitting metric, such as the Euclidian distance between FC matrices, to investigate to what extent these models can be useful to a different measure (see the Supporting Information, Figure S3). We noted that slight variations in the heterogeneity, such as those presented between the controls and disease functional maps (highly correlated, see Figure 5C), yield significantly different fitting levels. This high sensitivity could be related to the relationship between maps and the underlying network that couples the system. We computed the differences between the two sources of heterogeneity (GBC_CNT − GBC_AD) and compared them with the weighted degree of the underlying structural connectivity (Supporting Information, Figure S4). We found that both measures are correlated (R = 0.31, p < 0.001), showing that the regions where healthy controls present more GBC than AD patients correspond to highly connected regions. Thus, the differences in these regions could strongly impact the whole-brain dynamics, providing a possible explanation of the difference in the level of fitting between both heterogeneities. In summary, we demonstrated that the best model-fitting capability is reached when we include the corresponding disease-specific functional heterogeneity.

In this study, we extended previous work focused on investigating the dynamical consequences of local heterogeneities in asynchronous brain activity to the case of synchronous oscillatory behavior. We implemented phenomenological Stuart–Landau and exact mean-field whole-brain models of fMRI dynamics informed by T1w/T2w regional heterogeneity to evaluate the dynamical consequences of these variations in the oscillatory regime. The fit of the two models that are able to present oscillatory behavior with different levels of abstraction allows us to widely explore the static and dynamic implication of the inclusion of regional heterogeneities for modeling a large-scale dataset of resting-state fMRI recordings. We then identified that disease-specific regional functional heterogeneity also brings dynamical consequences within the oscillatory regime in fMRI recordings from patients with neurodegeneration (Alzheimer’s disease). In the following paragraphs, we discuss our findings in light of whole-brain modeling and regional heterogeneity within an oscillatory context, extending the findings of previous work in asynchronous models.

It has been demonstrated that cortical areas present multiple sources of regional variations, such as myeloarchitecture, cytoarchitecture, the density of neurotransmitter receptors, gene expressions, excitation/inhibition ratio, and interregional connectivity, among others (Deco et al., 2018; Deco, Kringelbach, et al., 2021; Deco et al., 2014; Huntenburg et al., 2018; Margulies et al., 2016; X. J. Wang, 2020). The relation between the spatial structure variation and the function has largely been studied, and the specific dynamical consequences of each of these heterogeneities have been investigated through whole-brain models. In previous work, Chaudhuri et al. (2015) included a hierarchical ordering of macaque cortical areas based on linearly scaling the excitatory input strength of each population into a mean-field model, changing the dynamical scenario from a single stable point toward a bistable state accordingly. In the same direction, Demirtaş et al. (2019) and Deco, Kringelbach, et al. (2021) included different sources of heterogeneity (comprising the T1w/T2w) to modulate the gain of the local neuronal response functions of the corresponding excitatory and inhibitory pools of each brain region. So far, these approaches have considered biophysical mean-field asynchronous models. Here, first, we introduced heterogeneity in the phenomenological SL whole-brain model by modifying the local bifurcation parameter of each region by the inclusion of a bias and a scaling that modulates the regional heterogeneity given by the T1w/T2w. In this sense, the heterogeneity can change the dynamical scenario for a single region from oscillatory toward noise dynamics, or vice versa. The relationship between the bifurcation parameter and the local excitation/inhibition ratio can be inferred by studying the bifurcation diagram of more biophysically realistic models, such as the Wilson–Cowan model (Cowan et al., 2016). In line with previous results on mean-field non-oscillatory models, we found that the inclusion of T1w/T2w regional heterogeneity has substantial implications on fitting static measures, such as the functional connectivity (FC) and the global brain connectivity (GBC) (Deco, Kringelbach, et al., 2021; Demirtaş et al., 2019). We also found that this heterogeneity has implications in the information transmission flow measured by the causal ignition within the NDTE framework (Deco, Vidaurre, et al., 2021). Specifically, this measure indicates the influence that each brain region presents over the whole level of brain synchronization.

We further investigated the implications that regional heterogeneities have on modeling brain activity by building a biophysically grounded mean-field model capable of reproducing the switch between oscillations and stable fixed points. In the last few years, several authors have modeled different observables from empirical fMRI data using Stuart–Landau whole-brain models (Deco & Kringelbach, 2020; Ghosh et al., 2008; Ipiña et al., 2020; Jobst et al., 2017). On the other hand, more sophisticated biophysical grounded models, such as mean-field models, were also successfully used to build whole-brain models capable of modeling empirical data from fMRI recordings (Breakspear et al., 2003; Deco & Jirsa, 2012; Deco, Kringelbach, et al., 2017; Deco et al., 2014; di Volo et al., 2019; Freyer et al., 2012; Ghosh et al., 2008; Herzog et al., 2020; Honey et al., 2007; Ponce-Alvarez et al., 2015). However, what level of abstraction is considered in the local dynamics to construct whole-brain models is an open question and focus of ongoing research. We hypothesize that one reason why the SL whole-brain model is one of the best models that fit to fMRI BOLD imaging despite its simplicity (Deco, Kringelbach, et al., 2017) is the capability of the Stuart–Landau model to capture the oscillatory nature of brain signals. In contrast, the spiking and standard mean field are asynchronous and unable to capture these oscillatory couplings. Furthermore, when these models are expressly set on the edge of the bifurcation point, they have demonstrated an excellent capability to reproduce metastability (Deco, Kringelbach, et al., 2017) and critical behavior of brain activity (Chialvo, 2010; Cocchi et al., 2017; Perl et al., 2021). In order to shed light on this question, we included oscillations in the mean-field model and compared them with the SL whole-brain model results.

We implemented an exact mean-field whole-brain model based on the mean-field description for populations of QIF neurons derived in Montbrió et al. (2015). We considered the local dynamics of each brain region to be described by a pool of inhibitory and excitatory neurons, coupled via structural connectivity with the other brain regions. We found that, as in the SL model, the best fit for this mean-field model is slightly below the Hopf bifurcation. We were able to fit the functional connectivity and the level of synchronization as well as other homogeneous whole-brain models based on mean-field local dynamics (Deco et al., 2018; Deco, Kringelbach, et al., 2021; Demirtaş et al., 2019; P. Wang et al., 2019). We investigated the implication of the regional heterogeneities, as in the SL model, when we model empirical observables by equally modulating the half-wide distribution of the variability within each population of inhibitory and excitatory neurons. We found that the T1w/T2w heterogeneity plays a crucial role in improving the level of fitting of GBC and causal ignition using this biophysical mean-field model capable of oscillatory dynamics.

We evaluated how the heterogeneous bifurcation parameters, that is, the ones that reach the best heterogenous level of fitting, were distributed with respect to the Hopf bifurcation. In both scenarios, the values are distributed around the bifurcation point, providing an indication that the optimal dynamical behavior is on the edge of criticality, between fluctuations and oscillations, as suggested by whole-brain models in previous research (Deco et al., 2013; Deco, Kringelbach, et al., 2017; Spiegler et al., 2016) and demonstrated in more detail in Freyer et al. (2012).

Finally, we investigated the dynamical consequences by including disease-specific regional heterogeneities based on functional fMRI recordings of Alzheimer’s disease patients. In general, neurodegenerative diseases such AD and bvFTD (behavioral variant frontotemporal dementia) are characterized by heterogeneity across different brain levels involving both structure and dynamics (Maestú et al., 2021; Mehta et al., 2013; Verdi et al., 2021), challenging the assumptions of spatially homogeneous models. Specifically, different measures on patients such as positron emission tomography (PET), which informs the accumulation of pathological proteins in the brain, such as amyloid-beta plaques (Nordberg, 2004) and neurofibrillary tau tangles (Schöll et al., 2016), or atrophy maps from fRMI, which may constitute signs of neurodegeneration (Schoonenboom et al., 2004), provide different sources of heterogeneity. However, insufficient mechanistic accounts of neurodegeneration prevent the development of integrative models of neurodegenerative diseases. We proposed investigating the heterogeneity implications in modeling AD data by following the approach proposed by Kong et al. (2021), who demonstrated that functional gradients improve the model-fitting capabilities of static and dynamic empirical observables from resting-state fMRI. We defined the functional heterogeneity as the GBC, which stands for how connected each brain region is with the rest of the brain in terms of functionality, and we computed for AD and healthy controls. Based on the results presented here and previous work, which successfully fit AD data with the SL model (Demirtaş et al., 2017), we used the SL model because it is less computationally intensive than the exact mean-field model and reaches a similar fitting level. We found that the AD functional heterogeneity in the model yielded a more accurate reproduction of the spatiotemporal structure of empirical FC, GBC, and CI. Note that this model largely overperformed the homogeneous model. Similar results were found by Kong et al. (2021), including functional gradient as regional heterogeneity in a mean-field model fitting the static and dynamic FC. We compared fitting improvements by including healthy control functional heterogeneity to fit the empirical observables for AD. Notably, despite that the two specific regional heterogeneities, healthy and AD, are highly correlated (r = 0.86), we found that the model with healthy control heterogeneity overperforms the homogeneous case but is not as good as the specific AD regional heterogeneity. We hypothesize that this fitting difference could be related to the fact that regions with the highest differences in heterogeneity between CNT and AD are the most connected by the underlying structure in the brain. Our results on AD fMRI data are proof of concept in fulfilling the gap of modeling neurodegeneration, providing mechanisms to include regional heterogeneity to build integrative models that allow the integration of different recording modalities.

In summary, we showed that the best fits are obtained by including the disease-specific functional heterogeneity from the AD patients, with poor results obtained when we used the GBC from healthy controls or the homogeneous model. This result prompts the need to discuss the relationship between this functional map and local dynamics. AD is linked to altered cellular energy metabolism (Gu et al., 2012), excitation/inhibition ratio (Lopatina et al., 2019; Maestú et al., 2021; Mehta et al., 2013), and neurotrophic factor release (Murer et al., 2001), impairing neural microcircuit function (Palop & Mucke, 2016). In particular, the excitation/inhibition ratio seems to play an important role in characterizing AD (Benussi et al., 2020; Cammisuli et al., 2022; Maestú et al., 2021). Considering that the Hopf bifurcation parameter and the local excitation/inhibition ratio can be related (which can be inferred by studying the bifurcation diagram of the Wilson–Cowan model [Cowan et al., 2016]), it is likely that this feature of cortical dynamics is being captured by the specific disease functional heterogeneity that directly modulates the local bifurcation parameter. Nevertheless, we acknowledge that with imaging technologies the number of high-resolution reference maps of brain structure and function are increasing (e.g., frequency gradients maps [Mahjoory et al., 2020] or the gene expression and protein density as used by Hansen et al. [2022]). In this direction, the comparison of our disease functional maps and these reference maps are not investigated here because it is out of the scope of this work. We focused on proposing a dynamic mechanism of the impact of including a map in oscillatory models, but extension to the evaluation of different maps will be considered in further studies. In particular, we propose leveraging the work of Markello et al. (2022), who developed an excellent tool to interpret structural and functional brain maps.

Overall, we demonstrated that two conceptually different whole-brain models presenting similar dynamic scenarios are equally good at fitting static and dynamic empirical data from resting-state fMRI recordings of healthy participants. Notably, the optimal working point of both models is at the edge of the Hopf bifurcation. We then showed that the dynamical consequences of the regional heterogeneities in the fitting capacity are similar in both models. To the extent that the SL whole-brain model performs similarly to an exact mean-field biophysically grounded whole-brain model, we conclude that the SL whole-brain model represents a suitable level of abstraction that captures one of the main dynamic behaviors of brain signals. Finally, we showed that functional specific-disease regional heterogeneities in the SL whole-brain model are relevant to improving the model fitting of empirical observables from fMRI recordings from participants with Alzheimer’s disease. We demonstrated that models with structural (Tw1/Tw2 ratio) or functional regional heterogeneities perform better than the homogeneous model. Furthermore, we showed that the specificity of the heterogeneity has a positive impact on the model’s performance; that is, the model performs better when the inclusion of the heterogeneity is related to the condition that has been modeled. A further direction is to include at the same time different classes of heterogeneities, structural and functional, specific and nonspecific, to improve model-fitting capacities but also to establish hierarchies between different heterogeneities. Steps in that direction were given by Kong et al. (2021) for healthy resting-state participants. We postulated that the combination of functional and structural sources of local heterogeneities could provide better model performance, and this framework can provide an analytical tool to determine the relevance of each class of heterogeneity.

Supporting information for this article is available at https://doi.org/10.1162/netn_a_00299.

Yonatan Sanz Perl: Conceptualization; Formal analysis; Investigation; Methodology; Software; Writing – original draft; Writing – review & editing. Gorka Zamora-Lopez: Investigation; Methodology; Writing – review & editing. Ernest Montbrió: Formal analysis; Investigation; Methodology; Writing – review & editing. Martí Monge-Asensio: Formal analysis; Visualization; Writing – review & editing. Jakub Vohryzek: Methodology; Writing – review & editing. Sol Fittipaldi: Data curation; Writing – review & editing. Cecilia González Campo: Data curation; Writing – review & editing. Sebastián Moguilner: Data curation; Writing – review & editing. Agustín Ibañez: Data curation; Writing – review & editing. Enzo Tagliazucchi: Formal analysis; Methodology; Writing – review & editing. B. T. Thomas Yeo: Methodology; Writing – review & editing. Morten Kringelbach: Conceptualization; Methodology; Software; Writing – review & editing. Gustavo Deco: Conceptualization; Formal analysis; Investigation; Methodology; Software; Supervision; Writing – original draft; Writing – review & editing.

Yonatan Sanz Perl, H2020 Marie Skłodowska-Curie Actions (https://dx.doi.org/10.13039/100010665), Award ID: 896354. Yonatan Sanz Perl, EU H2020 FET Flagship program and by European Union’s Horizon 2020, Award ID: 945539. Gorka Zamora-Lopez, EU H2020 FET Flagship program and by European Union’s Horizon 2020, Award ID: 945539. Ernest Montbrió, Ministerio de Educación y Formación Profesional (https://dx.doi.org/10.13039/501100020636), Award ID: PID2019-109918GB-I00/AEI/10.13039/501100011033. Agustín Ibañez, CONICET, Award ID: FONCYT- PICT (2017-1818, 2017-1820). Agustín Ibañez, Takeda, Award ID: CW2680521. Agustín Ibañez, Alzheimer’s Association, Award ID: SG-20-725707-ReDLat. Agustín Ibañez, National Institutes of Aging of the National Institutes of Health, Award ID: R01AG057234. Agustín Ibañez, ANID/FONDECYT Regular, Award ID: 1210195, 1210176, 1220995. Agustín Ibañez, ANID/FONDAP, Award ID: 15150012. Agustín Ibañez, ANID/PIA/ANILLOS, Award ID: ACT210096. Enzo Tagliazucchi, CONICET, Award ID: PICT-2019-02294. Enzo Tagliazucchi, ANID/FONDECYT, Award ID: 1220995. Morten Kringelbach, Danish National Research Foundation, Award ID: DNRF117. Gustavo Deco, EU H2020 FET Flagship program and by European Union’s Horizon 2020, Award ID: 945539. Gustavo Deco, Spanish Ministry of Science, Innovation, and Universities (MCIU), Award ID: ref. PID2019-105772GB-I00 MCIU AEI.

Hopf bifurcation:

The critical point where the system changes the solutions from fixed point to a self-sustained oscillatory behavior.

Global brain connectivity (GBC):

This measure stands for the mean level of functional connectivity that each brain region has with the rest of brain regions.

Functional connectivity (FC):

The pairwise Pearson’s correlation between fMRI time series.

Kuramoto order parameter (KoP):

This measure reflects the level of global synchronization of a system of N oscillators; when KoP = 1, the system is fully synchronized.

Causal ignition (CI):

This measure reflects the direct causality between the activity of one region and the global level of synchronization of the system.

Hopf whole-brain model:

Computational model describing the whole-brain dynamics, where the dynamics of each brain region is described by a Stuart-Landau oscillator.

Exact mean-field whole-brain model:

Model where the dynamics of each brain region are described by the mean-field equations exactly derived from a large population of all-to-all coupled quadratic integrate-and-fire (QIF) neurons.