Stable but flexible brain activity patterns arise from local balance of excitation and inhibition in the resting cortex

To approximate the anatomical connections between brain areas, we use a normative connectome from 32 healthy participants (mean age = 31.5 years old ± 8.6, 14 females) generated as part of the Human Connectome Project (HCP; https://www.humanconnectome.org↗). In this work, we focus on the cortical regions of the AAL atlas, yielding a 78 × 78 matrix. Tract lengths between brain regions were also extracted to be used in the computation of conduction delays between network nodes. For more detailed information, refer to Castaldo et al. (2023).

To analyze network dynamics and FC properties in empirical data, we use resting-state BOLD time series from 99 healthy unrelated subjects from the HCP dataset. Data were processed using the pipeline described in Castaldo et al. (2023), yielding 99 time series of size 78 areas × 1200 TR, roughly corresponding to 14.5 min of data for each subject (TR = 0.72 s).

To model the activity of individual cortical regions, comprising excitatory (r^E) and inhibitory (r^I) neural masses, we make use of the Wilson–Cowan model (Cowan, Neuman, & Van Drongelen, 2016; Wilson & Cowan, 1972):τEriEtdt=−riEt+FEGEcEEriEt−cEIriIt+GEIiext+ξtτIriItdt=−riIt+FIcIEriEt+ξt(1)

Where c^EE represents the recurrent excitatory coupling, c^EI represents the coupling from inhibitory to excitatory populations, and c^IE represents the excitatory coupling unto inhibitory populations. G^E is a parameter introduced in Páscoa dos Santos and Verschure (2025), allowing for the scaling of all excitatory inputs to the excitatory neural mass (i.e., recurrent excitation and I^ext). The default values of each parameter can be consulted in Supporting Information Table S1. Unless stated otherwise, ξ(t) represents Gaussian noise with mean 0 and variance 0.01. Here, Iiext describes the incoming input to each node i from the rest of the network and can be written as follows:Iiext=C∑j=1NWijrjEt−τij+P(2)where C is a scaling factor referred to as global coupling, W_ij is the SC matrix, and τ_ij represents the conduction delay between nodes i and j. Here, similarly to Abeysuriya et al. (2018), Castaldo et al. (2023), and Páscoa dos Santos et al. (2023), P represents a background input imparted on the nodes and is set at 0.31 so that the default uncoupled node is close to the Hopf bifurcation. The time constants τ^E and τ^I were tuned so that the uncoupled node oscillates at ∼40 Hz. Furthermore, we define F^E/I(x) as the input–output function of the excitatory and inhibitory populations, respectively, using the following equation:FE/Ix=11+e−x−μE/IσE/I(3)

As opposed to the common approach to E–I homeostasis (Abeysuriya et al., 2018; Castaldo et al., 2023; Deco et al., 2014; Hellyer et al., 2016; Naskar et al., 2021; Páscoa dos Santos et al., 2023), we do not implement dynamical homeostasis. Instead, we follow the framework developed in Páscoa dos Santos and Verschure (2025) to mathematically estimate the parameter values that ensure the maintenance of mean excitatory activity at a target firing rate (ρ) for different levels of input (I^ext; Hengen, Lambo, Van Hooser, Katz, & Turrigiano, 2013; Ma et al., 2019), in line with the governing rules of the E–I homeostasis (Turrigiano, 2011). That said, to compute the homeostatic parameters of nodes embedded in a network, it is necessary to provide an estimate of I^ext for each node. To do this, we rely on the fact that E–I homeostasis works in timescales much slower than neural dynamics (Gainey & Feldman, 2017; Turrigiano, 2011; Wen & Turrigiano, 2024), ensuring timescale separation. Therefore, to compute the homeostatic parameters, we consider the average external input received by each node at the homeostatic set point, when the average activity across the network corresponds to ρ, leading to the following expression:Iiextt=Pt+C∑jWijrjEtIiextt=P+C∑jWijrjEtIiext=P+Cρ∑jWij(4)

Here, we explore the behavior of models with the following distinct mechanisms of the E–I homeostasis, based on experimental results (Desai et al., 1999; Maffei et al., 2004; Maffei & Turrigiano, 2008; Turrigiano, 2011; Turrigiano et al., 1998; Wen & Turrigiano, 2021, 2024) and described in detail in Páscoa dos Santos and Verschure (2025):G^E homeostasis: synaptic scaling of excitatory synapses in PY neurons.c^EI homeostasis: synaptic scaling of inhibitory synapses in PY neurons.μ^E homeostasis: plasticity of the intrinsic excitability of PY neurons, modulating the "firing threshold" (i.e., μ^E) of F^E(x).μ^E + σ^E homeostasis: plasticity of the intrinsic excitability of PY neurons, modulating both the "firing threshold" (μ^E) and "slope" (σ^E) of F^E(x) in a synergistic manner.G^E + c^EI homeostasis: synaptic scaling of excitatory and inhibitory synapses in PY neurons.G^E + c^EI + μ^E homeostasis: synaptic scaling of excitatory and inhibitory synapses and intrinsic excitability of PY neurons, implemented at the level of μ^E.G^E + c^EI + μ^E + σ^E homeostasis: synaptic scaling of excitatory and inhibitory synapses and intrinsic excitability of PY neurons, implemented at the level of μ^E and σ^E.

The equations to compute the set-point parameters corresponding to each form of homeostasis can be consulted in SI Methods.

From the raw activity of the excitatory populations, we extract simulated BOLD signals by using a forward hemodynamic model (Friston, Mechelli, Turner, & Price, 2000; Friston, Harrison, & Penny, 2003), described in detail in Páscoa dos Santos et al. (2023). This model represents the coupling between model activity and blood vessel diameter, which affects blood flow, blood volume, and deoxyhemoglobin content. After passing model activity through the hemodynamic model, the output is downsampled to a sampling period of 0.72 s to align simulated signals with empirical data.

For each type of homeostasis and combination of hyper-parameters (C, ρ, and mean delay), we evaluate if the homeostatic set point—that is, a mean firing rate of ρ for all nodes—can be maintained across the network. To do this, we run 15-s simulations and compute the mean firing rate of each node in the last 10 s of the simulation. We consider the given network fixed point to be invalid if there is any node in the network for which the mean firing rate differs from ρ by more than 1% (i.e., ∣rE−ρ∣ρ>0.01). The dynamical behavior of invalid models is examined in detail in SI Text. While it is possible to evaluate the stability of a network analytically (see Ghosh, Rho, McIntosh, Kötter, & Jirsa, 2008), the characteristic polynomial of our model (from which the eigenvalues can be extracted) does not have a closed-form solution due to the highly nonlinear nature of the Wilson–Cowan equations. For this reason, we base our exploration of model validity on the analysis of the numerical solution of the system.

We perform model optimization by exploring all combinations of the global coupling (C), mean delay, and target firing rate (ρ) within the following ranges of variation: C: [0, 10], mean delay: [0, 40] ms, and ρ: [0.05, 0.2]. Within these ranges, we selected 19 logarithmically spaced values for C, 16 values for ρ in steps of 0.01, and 41 mean delays in steps of 1 ms. For each simulation, we initialize the model with the homeostatic parameters corresponding to the combination of C and ρ and record model activity for 30 min. To evaluate model performance against empirical data, we use the following properties of FC.▪ Static FC: 78 × 78 matrix of correlations between BOLD time series across all network nodes. Modeled FC matrices were compared with group-averaged empirical FC by computing the correlation coefficient (r_FC) and mean squared error (MSE_FC) between their upper triangular elements.FCD: matrix of correlations between the upper triangular part of FC matrices computed in windows of 80 samples with 80% overlap (Deco et al., 2017, 2021). Model results are compared with empirical data by computing the Kolmogorov–Smirnov (KS) distance between the distributions of values in the respective FCD matrices (KS_FCD).

BOLD signals are filtered between 0.008 and 0.08 Hz before computing FC and FCD. In addition, we compute a cross-feature fitting score as rFC−MSEFC−KSFCD−SminSmax−Smin to evaluate the performance of models in simultaneously reproducing static and dynamic properties of empirical FC (Castaldo et al., 2023; Páscoa dos Santos et al., 2023). S_max = 1 represents the maximum value of r_FC − MSE_FC − KS_FCD (1–0–0) while S_min = −4.0024 represents its minimum (−1–2.0024–1). The value of 2.0024 represents the maximum possible MSE between an FC matrix and our empirical FC matrix, considering that FC can only vary between −1 and 1. With this expression, we obtain a cross-feature fitting score normalized between 0 and 1, indicating how well models represent FC and FCD.

All simulations were performed by solving the system's equations numerically using the Euler method with an integration time step of 0.2 ms.

To evaluate the magnitude of synchrony and the degree of dynamic switching between network states (i.e., metastability), we compute the Kuramoto Order Parameter (KOP; Kuramoto, 1975; Shanahan, 2010) of BOLD timeseries, filtered between 0.008 and 0.08 Hz. More specifically, the KOP (R(t) in Equation 5) represents the degree of phase alignment within a set of coupled oscillators at a given point in time and can be calculated as follows:RteiΦt=1N∑n=1Neiθnt(5)where θ_n(t) represents the instantaneous Hilbert phase of node n at time t. Synchrony and metastability are defined, respectively, as the mean and standard deviation of R(t) over time.

Previous studies have demonstrated that the architecture of structural networks has a significant impact on the complexity of the functional interactions between nodes (Arsiwalla & Verschure, 2016; Sporns et al., 2002; Zamora-López et al., 2016). To explore this, we measure the complexity of simulated FC matrices using the method introduced in Zamora-López et al. (2016). In short, we quantify complexity by measuring how strongly the distribution of FC weights deviates from an equivalent uniform distribution using the following formula:C=1−1CM∑m=1M∣pmFCij−1M∣(6)where M represents the number of bins used to evaluate the FC distribution. The term CM=2m−1m is added for normalization, corresponding to the deviation of a Dirac-delta function from the uniform distribution. Therefore, FC complexity is 0 when all entries of the FC matrix are the same and 1 when they are uniformly distributed. While there are other methods to evaluate the complexity of a matrix, Zamora-López et al. (2016) demonstrate that their formulation is more robust to changes in bin sizes, compared with other methods quantifying the entropy or variance of FC matrices. Here, we use a bin size of 0.05.

To investigate the effects of the E–I homeostasis on slow local dynamics, we employ the reduced Wong–Wang model of cortical populations (Deco et al., 2014). In this model, instead of simulating fast α-Amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA) synapses, excitatory synapses have a considerably slower decay time (100 ms), consistent with the timescale of N-methyl-D-aspartate (NMDA) synapses. The model dynamics are described as follows:IiE=GEWEI0+GEw+JNMDASiE+CGEJNMDA∑jWijSjE−JGABASiIIiI=WII0+JNMDASiI−JGABASiIriE=FEIiEriI=FIIiIdSiEtdt=−SiEτE+1−SiEγriE+ξtdSiItdt=−SiIτI+riI+ξt(7)where IiE/I represents the total input to excitatory or inhibitory neuronal population i, riE/I represents its firing rate, and SiE/It represents an average synaptic gating variable. In this case, we have also added the parameter G^E, which allows for the simultaneous modulation of all excitatory synapses unto the excitatory neuronal population. The function ξ(t) represents the additive Gaussian noise with mean 0 and variance 0.01 (nA). H^E/I(t) represents the input–output function of neuronal populations and is written as follows:FE/Ix=aE/Ix−bE/I1−exp−dE/IaE/Ix−bE/I(8)

The parameter values, adapted from Deco et al. (2014) and their physiological meaning are described in detail in Supporting Information Table S3. Importantly, this model does not consider the dynamics of AMPA synapses and τ^E = 100 ms (the time constant of excitatory synapses) represents the slow decay time of NMDA (Deco et al., 2014). Therefore, in line with Deco et al. (2014) and given the focus of the Wong–Wang model on slower dynamics, interareal conduction delays are neglected in this model.

Similarly to the Wilson–Cowan model, we derived the equations allowing for the mathematical estimation of the steady-state value of each model parameter corresponding to the different empirically based mechanisms of homeostasis. More specifically, we investigate the following forms of the E–I homeostasis:G^E homeostasis: synaptic scaling of excitatory synapses in PY neurons.J_GABA homeostasis: synaptic scaling of inhibitory synapses in PY neurons.b^E homeostasis: plasticity of the intrinsic excitability of PY neurons, modulating the "firing threshold" (i.e., b^E) of F^E(x).a^E homeostasis: plasticity of the intrinsic excitability of PY neurons, the "slope" of (σ^E) of F^E(x) in a synergistic manner. In this case, since changing a^E has a minimal impact on the value of F^E(x) at 0, we employ this form of homeostasis to simulate the plasticity of intrinsic excitability described in Desai et al. (1999); Nataraj, Le Roux, Nahmani, Lefort, and Turrigiano (2010); and Wen and Turrigiano (2021). For more details, refer to Supporting Information Figure S9.G^E + J_GABA homeostasis: synaptic scaling of excitatory and inhibitory synapses in PY neurons.G^E + J_GABA + b^E homeostasis: synaptic scaling of excitatory and inhibitory synapses and intrinsic excitability of PY neurons, implemented at the level of b^E.G^E + J_GABA + a^E homeostasis: synaptic scaling of excitatory and inhibitory synapses and intrinsic excitability of PY neurons, implemented at the level of a^E.

For more details on the derivation of the mathematical expressions and the effects of the E–I homeostasis at the circuit level, refer to the SI Methods.

Similarly to the procedure employed for the Wilson–Cowan model, we perform model optimization by running one simulation for each combination of the global coupling (C) and target firing rate (ρ) within the following ranges of variation: C: [0, 1] and ρ: [1, 10] Hz. Within these ranges, we selected 25 logarithmically spaced values for C and 10 values for ρ. For each simulation, we initialize the model with the homeostatic parameters corresponding to the combination of C and ρ and record model activity for 30 min. In Deco et al. (2014), the target firing rate is fixed at 3.0631 Hz; the model is considered to have reached the homeostatic fixed point when the mean firing rate is between 2.63 and 3.55 Hz, corresponding to an interval of ∼15% around the target. Therefore, here, we consider a simulation valid when the error between mean firing rates and the target is smaller than 15% for all nodes in the network.

Finally, we detect the optimal combination of C and ρ maximizing the fitting score for each form of homeostasis. To be able to statistically compare models, we run 10 simulations with each combination of optimal parameters and compute the respective cross-feature fitting score and metastability.

To explore the impact of noise in our model, we run simulations of the network with homeostasis of G^E, c^EI, and μ^E +σ^E with optimal parameters (C = 3.59, ρ = 0.12) and across all mean delays yielding valid simulations. For each simulation, we vary the variance of Gaussian noise (ξ(t)), corresponding to 29 logarithmically spaced values between 10⁻⁵ and 10−13. We then analyze model performance, synchrony, metastability, and functional complexity across all models with varying levels of noise.

Following the framework applied in Páscoa dos Santos et al. (2023), we investigate the effects of focal lesions in model dynamics and functional interactions by following the following protocol. Using the model with all mechanisms of homeostasis at the optimal point (C = 3.59; ρ = 0.12), we compute the homeostatic parameters by using the methodology from Páscoa dos Santos and Verschure (2025; SI Methods). From the balanced model, we extract 30 min of the BOLD signal corresponding to the activity in the prelesion (i.e., healthy) period. Then, we apply a structural lesion by removing all connections to and from a chosen node and extract 30 min of activity corresponding to the early acute period of stroke recovery. We then compute the new homeostatic parameters, allowing for the restoration of E–I balance in the lesioned connectome, using the prelesion homeostatic values as G0E, c0EI, and μ0E (SI Methods for the equations). From the rebalanced model, we extract 30 min of model activity corresponding to the chronic period. We follow this protocol for all possible single-node lesions and extract the following metrics from the healthy, acute, and chronic simulations: difference of FC from baseline, FC-SC correlation, FC modularity, FC complexity, synchrony, metastability, and FCD distributions. For a more detailed description of how each feature is computed, refer to SI methods. While, in the main text, we present the effects of lesions on models with the combination of all mechanisms of homeostasis, we also report the effects of lesions on models with G^E and G^E +c^EI homeostasis (Supporting Information Figures S27–S28).

All simulations and analyses were run using in-house scripts in Python. The code can be consulted in https://gitlab.com/francpsantos/wc_network_homeostasis↗.

We start our analysis by studying how distinct forms of homeostasis, based on the plasticity of excitation, inhibition, and intrinsic excitability (see the Methods section and Figure 1), shape the performance of our model in reproducing empirical FC and FCD obtained from resting-state BOLD signals, which are strongly modulated by E–I balance (Castaldo et al., 2023; Páscoa dos Santos et al., 2023). We simulate the activity of a large-scale model of the neocortex with different combinations of free parameters (C, ρ, and mean delay) and mechanisms of homeostasis (Figure 1). For each simulation, we compare empirical and simulated FC matrices and FCD distributions and extract a cross-feature fitting score to evaluate model performance. The score is calculated as rFC−MSEFC−KSFCD−SminSmax−Smin, where r_FC is the Pearson's correlation between empirical and simulated FC matrices, MSE_FC is their mean squared error, KS_FCD is the Kolmogorov–Smirnov distance between FCD distributions, and S_max/S_min is the maximum/minimum of r_FC − MSE_FC − KS_FCD. With this, we obtain a cross-feature fitting score between 0 and 1, reflecting how well models represent the features of interest (see the Methods section for more details). Furthermore, our results remain similar when employing a different method to compute the cross-feature fitting score (Supporting Information Figure S12).

First, it is relevant to point out that the main parameters impacting model performance are C and ρ, while the mean delay of interareal conduction does not have a strong impact, provided it is long enough to avoid widespread pathological synchronization (Supporting Information Figures S13–S16 for the full parameter spaces). Therefore, we present the performance scores for each combination of C and ρ, averaged across mean delays (Figure 2A). For each mechanism of homeostasis, we observed a narrow region in the parameter space where empirical resting-state FC and FCD emerge. To compare the optimal performance of each model, we select the combination of C and ρ yielding the highest fitting score while ensuring that simulations are valid for at least 10 different values of mean delay. Although conduction delays do not strongly impact the fitting score (Supporting Information Figures S13–S16), we enforce this constraint to ensure that models are robust to changes in conduction velocity without losing the ability to maintain stable firing rates. That said, we present the distribution of cross-feature fitting scores across mean delays for the optimal point of each model (Figure 2B, Table 1) and examples of activity, FC, and FCD distributions from empirical and simulated data at the optimal point of each model (Figure 2C). The best scores are obtained for models with G^E homeostasis (0.855 ± 0.009) and G^E +c^EI homeostasis (0.851 ± 0.019), followed by the homeostasis of G^E, c^EI, and μ^E +σ^E (i.e., all; 0.849 ± 0.009). The values of r_FC, MSE_FC, and KS_FCD corresponding to each optimal score can be consulted in Table 1. While the performance of the model with G^E homeostasis was significantly better than the one with all mechanisms of homeostasis (p = 0.020; Mann–Whitney U test, false discovery rate [FDR] correction), the effect size is not substantial (Cohen's d = 0.281). Furthermore, the optimal target firing rates for the three best models (G^E: 0.12; G^E +c^EI: 0.1; all: 0.12) are close to the bifurcation between damped and sustained oscillations (Supporting Information Figure S1), in line with theoretical (Cabral et al., 2014; Deco et al., 2017; Freeman, 2005) and experimental (Ma et al., 2019) results suggesting that the neocortex operates at the edge-of-bifurcation.

Our results indicate that models with E–I homeostasis can reproduce the emergence of empirical networks and their dynamics at the ultra-slow timescales characteristic of resting-state BOLD, depending on the interaction between global coupling (C) and target firing rate (ρ). In addition, empirical data are best approximated when accounting for mechanisms of homeostasis that modulate the bifurcation point (i.e., the ones that modulate G^E or σ^E; Páscoa dos Santos & Verschure, 2025), so that areas with stronger inputs are poised farther away from the bifurcation (Figure 1B and Supporting Information Figure S1). This effect compensates for the higher input fluctuations in strongly connected nodes (Supporting Information Figure S17), ensuring that such nodes are not constantly entering the regime of sustained oscillations. Therefore, even though this effect leads to heterogeneous distances to the bifurcation across cortical areas, we argue that our results still align with the edge-of-bifurcation theory of brain dynamics (Arsiwalla & Verschure, 2016, 2017; Cabral et al., 2014; Deco & Jirsa, 2012; Deco et al., 2017; Freeman, 2005; Freeman & Holmes, 2005; Schirner et al., 2022) when input fluctuations are accounted for. More importantly, this bifurcation modulation, observed in models with the homeostasis of excitation (G^E) and excitability slopes (σ^E; Figure 1B and Supporting Information Figure S1), is essential for the emergence of FC and FCD (Figure 2).

Given that the metastability of large-scale cortical networks is shaped by circuit-level edge-of-bifurcation dynamics (Deco et al., 2017; Freeman, 2005; Freeman & Holmes, 2005; Naskar et al., 2021), controlled by E–I homeostasis (Páscoa dos Santos & Verschure, 2025; Wen & Turrigiano, 2024), it is relevant to explore how distinct mechanisms of homeostasis shape the collective dynamics of the resting-state cortex. To do this, we measure the synchrony and metastability (see the Methods section) of the simulated BOLD signals of models under different homeostatic mechanisms and combinations of C and ρ. Again, the conduction delays do not significantly impact either synchrony or metastability (Supporting Information Figures S18–S19), and, for this reason, we present results averaged across mean delays.

We observe that synchrony generally follows increases in C and ρ, reaching a maximum before the networks become unstable (Supporting Information Figure S20). To explore how synchrony is shaped by E–I homeostasis, we compare the optimal models of each mechanism of homeostasis (i.e., best cross-feature fit) with empirical data (Figure 3A, Table 1). While some of the models that best reproduce empirical FC and FCD generally display moderate levels of synchrony, models based on the homeostasis of G^E +c^EI (0.707 ± 0.041) or G^E (0.661 ± 0.027) show significantly higher synchrony than empirical resting-state BOLD dynamics (0.537 ± 0.112; p < 0.001, Mann–Whitney U test, FDR corrected). In contrast, networks based on the combination of all mechanisms of homeostasis, which reproduce empirical FC and FCD to a similar degree (Figure 2B), show levels of synchrony similar to empirical data (0.580 ± 0.021; p = 0.108; Mann–Whitney U test, FDR corrected).

Regarding metastability, our results reveal a region of maximal metastability determined by the interaction between C and ρ (Supporting Information Figure S21), coinciding with the region of optimal fit to empirical data (Figure 2A). Furthermore, the mechanisms of homeostasis that best fit empirical FC and FCD also show levels of metastability close to empirical values (0.194 ± 0.025; Figure 3B, Table 1). More specifically, we found enhanced metastability in models with the homeostasis of G^E +c^EI (0.203 ± 0.011), G^E (0.201 ± 0.009), and the combination of all mechanisms of the E–I homeostasis (0.191 ± 0.008), with no significant difference from empirical values (G^E +c^EI: p = 0.254; G^E: p = 0.254; All: p = 0.521; Mann–Whitney U test, FDR corrected). Importantly, empirical levels of both synchrony and metastability can only be reproduced when including all mechanisms of homeostasis (Figure 3 and Supporting Information Figure S22). In addition, we found a pronounced positive correlation between the metastability of models and their fitting scores, not only in the model with all mechanisms of homeostasis (Figure 3C) but as a general principle (Supporting Information Figure S23). Therefore, our results confirm previous studies proposing that the resting-state dynamics of cortical networks might correspond to a regime of maximal metastability (Deco et al., 2017; Naskar et al., 2021; Yang, Shew, Roy, & Plenz, 2012).

In summary, our analysis of network dynamics suggests that empirical FC networks and dynamics emerge more prominently in networks with moderate synchrony and high metastability, which can only be replicated when combining three mechanisms of the E–I homeostasis (G^E, c^EI, and μ^E +σ^E). More importantly, our results show that, through the homeostatic regulation of firing rates, cortical networks self-organize toward a regime of maximal metastability where the empirical resting-state networks emerge.

The complexity of spontaneous activity in the human brain is strongly related to the topology of its structural networks (Arsiwalla & Verschure, 2016; Sporns et al., 2002; Zamora-López et al., 2016). However, the architecture of the connectome is likely not the only source of complexity in cortical networks. Given the relevance of the E–I homeostasis at both the local and global scales, determining how local circuits engage in collective dynamics, E–I balance might be key in the emergence of complex functional interactions between cortical areas. Therefore, we study how different mechanisms of the E–I homeostasis impact FC complexity, using the methodology presented in Zamora-López et al.(2016; see the Methods section).

Our results on FC complexity suggest that, for all mechanisms of homeostasis, it is maximized in the region of the parameter space corresponding to heightened metastability and optimal fit to empirical FC and FCD (Supporting Information Figures S24 and S25). Comparing complexity in the models that best fit empirical FC and FCD data (Figure 4A, Table 1), our results indicate that the combined homeostasis of G^E, c^EI, and μ^E (0.797 ± 0.020) maximizes FC complexity, although the resulting values are significantly higher than in empirical data (p = 0.014; Mann–Whitney U test, FDR corrected). Conversely, some of the homeostatic mechanisms with the best fit to empirical FC and FCD are associated with levels of complexity significantly lower than empirical values (G^E: 0.640 ± 0.017, p < 0.001; G^E +c^EI: 0.565 ± 0.052, p < 0.001; Mann–Whitney U test, FDR-corrected). However, models relying on the combination of all forms of homeostasis reach levels of FC complexity (0.750 ± 0.020; p = 0.026) consistent with empirical data. Therefore, accounting for multiple mechanisms of homeostasis is necessary to reproduce not only synchrony and metastability but also functional complexity. Importantly, we found a strong correlation between metastability and complexity, not only for the model with all mechanisms of homeostasis (Figure 4B) but also across models (Supporting Information Figure S23). This suggests that heightened metastability of resting-state brain dynamics (Deco et al., 2017; Naskar et al., 2021) results in highly complex activity patterns that reflect on the variability of FC. More importantly, this requires the combined action of multiple mechanisms of the E–I homeostasis at the circuit level.

Having investigated the relationship between E–I homeostasis and complexity we compare complexity in models with the empirical and shuffled structural connectomes to assert the relevance of network topology (Arsiwalla & Verschure, 2016; Sporns et al., 2002; Zamora-López et al., 2016). Since shuffling preserves the distribution of edge weights, the complexity of the original and shuffled structural matrices is the same (Supporting Information Figure S26). That said, our results show that the complexity of FC in models with the empirical connectome is significantly higher than those with its shuffled counterpart (original: 0.750 ± 0.020; shuffled: 0.439 ± 0.042; p < 0.001; Mann–Whitney U test; Figure 4C). Therefore, in line with Arsiwalla and Verschure (2016); Sporns et al. (2002); and Zamora-López et al. (2016), the complexity of functional interactions also depends on the topology of SC. However, we point out that, given the relationship between metastability and FC complexity, this result could be an effect of decreased metastability in models with a shuffled connectome (Figure 4D). Therefore, our results suggest that while E–I homeostasis maximizes metastability and complexity by regulating edge-of-bifurcation dynamics in local cortical circuits, the role of the structural connectome in constraining their interactions should not be overlooked.

In summary, our results indicate that the processes underlying the emergence of metastable dynamics also maximize the complexity of spontaneous activity in the cortex. More importantly, our results suggest that multiple mechanisms of homeostasis (modulating excitatory and inhibitory synapses and intrinsic excitability through μ^E and σ^E) are required to reproduce the various dynamical signatures of cortical networks (Table 1). Therefore, we propose that the role of multiple mechanisms of the E–I homeostasis does extend beyond the local circuit level (Gainey & Feldman, 2017; Wen & Turrigiano, 2024), shaping the collective dynamics of large-scale cortical networks.

The version of the Wilson–Cowan implemented to simulate the local dynamics of cortical areas is based on the fast dynamics of gamma oscillations (∼40Hz). These rhythms are generated by the recurrent loop between E–I populations, with fast inhibition provided by parvalbumin-positive (PV) interneurons and fast excitation through α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic (AMPA) receptors (Buzsáki & Wang, 2012; Wang, 2010; Wang & Buzsáki, 1996). While the regulation of AMPA receptors has an important participation in the homeostasis of excitatory synapses (Turrigiano, 2011), experimental results show that the slower NMDA receptors are also modulated to maintain the homeostasis of firing rates (Watt, et al., 2000). Therefore, it is relevant to investigate if the effects of the local E–I homeostasis on the emergence of FC, FCD, and collective metastable dynamics are limited to its impact on fast gamma rhythms, related to AMPA receptors, or can also be observed when accounting for the slower NMDA dynamics.

To investigate this, we implemented a large-scale model with local dynamics based on the reduced Wong–Wang model (Deco et al., 2014), where excitatory synapses have slower time constants based on the timescales of NMDA receptors (Figure 5A). First, following the framework of Páscoa dos Santos and Verschure (2025), previously implemented for the Wilson–Cowan model, we derive the analytical expressions to calculate the parameter values at the homeostatic fixed point for each empirically derived mechanism of homeostasis. For more details on the derivation, refer to the SI Methods and Supporting Information Figure S27. Importantly, since the timescale of excitatory synapses (100 ms) is substantially slower than the timescale of inhibition (10 ms), the Wong–Wang model does not have a Hopf bifurcation. Therefore, the local dynamics dependent on NMDA receptors remain in the stable regime for all tested combinations of external input and target firing rate, regardless of the combination of homeostatic mechanisms (Supporting Information Figure S28), as opposed to the strong bifurcation modulation observed for the Wilson–Cowan model (Figures 1B and Supporting Information Figure S1).

With that in mind, we explore the performance of large-scale models based on Wong–Wang local dynamics by studying the dynamics of networks with different combinations of target firing rate (ρ) and global coupling (C). In line with the implementation of Deco et al. (2014), we do not account for interareal conduction delays, due to their negligible impact on the slow NMDA dynamics. Furthermore, following Deco et al. (2014), we consider a simulation valid if the firing rate is within 15% of the target, as opposed to the more stringent condition of 1% imposed on the Wilson–Cowan model. For each combination of parameters, we analyze how the model reproduces empirical FC and FCD by computing the cross-fitting scores. Our results indicate that the score depends on the interaction between C and ρ, similarly to the Wilson–Cowan model (Figure 5B). The optimal values of C and ρ and the respective scores can be consulted in Supporting Information Table S1. When comparing results at the optimal point for each model, our results show that there is no significant difference between the fitting scores of models with distinct mechanisms of homeostasis (Figure 5C). This is due to the reduced impact of the E–I homeostasis on the local NDMA dynamics (Supporting Information Figure S28), when compared with its effect on the dynamics of faster gamma rhythms, dependent on AMPA receptors (Figures 1B and Supporting Information S1; Páscoa dos Santos & Verschure, 2025). More importantly, the performance of any of the networks based on the Wong–Wang model is significantly lower than the performance of the Wilson–Cowan model with the homeostasis of excitation (G^E), inhibition (c^EI) and excitability slope (σ^E), and threshold (μ^E; Figure 5C). In line with this result, even though we also found a strong relationship between metastability and the emergence of FC and FCD when using the Wong–Wang model (Supporting Information Figure S29), networks based on slow NMDA dynamics were not able to reach the levels of metastability observed in empirical data (Figure 5D).

That considered, our results suggest that the resting-state collective dynamics of the cortex, measured with BOLD fMRI, more strongly reflect the interaction between gamma rhythms generated by fast excitation and inhibition than the slower dynamics of NDMA receptors. Therefore, given the strong impact of the E–I homeostasis on local gamma dynamics, the collective dynamics of the cortex at rest are strongly dependent on the mechanisms employed for the homeostasis of firing rates.

Synaptic transmission is known that one to be noisy in brain networks (Rusakov, Savtchenko, & Latham, 2020). This variability can have a relevant role in spontaneous brain dynamics (Pinneo, 1966; Rusakov et al., 2020), underlying phenomena such as stochastic resonance (Deco, Jirsa, & McIntosh, 2011; Deco, Jirsa, McIntosh, Sporns, & Kötter, 2009; Schirner et al., 2022). In addition, systems with multiple weakly stable attractor states display noise-driven transitions between attractors, corresponding to multistable dynamics (Deco & Jirsa, 2012; Schirner et al., 2022), which can be confounded with metastability (Kelso, 2012). In metastable systems, collective dynamics emerge as a consequence of the interaction between network topology and node dynamics, largely independently of noise. However, it is not yet clear if resting-state brain activity is more reflective of a metastable or multistable system (Kelso, 2012; Schirner et al., 2022). For that reason, we explore the behavior of our network under different levels of noise (see the Methods section). Given that the model with all mechanisms of homeostasis better reproduces empirical dynamics, connectivity, and complexity, we take the optimal model with all forms of homeostasis (C = 3.59, ρ = 0.12) and simulate BOLD signals under varying levels of noise (see the Methods section). We then analyze model performance, collective dynamics, and FC complexity to assess how they are influenced by the variability of neuronal noise. Since we did not observe an interaction between noise and conduction delays (Supporting Information Figure S30), we present the average results across mean delays (Figure 6). Our results demonstrate that there is no significant effect of noise on model performance, collective dynamics, and FC complexity for lower noise levels. Conversely, as the variance of noise increases beyond these levels, model performance deteriorates. This effect is to be expected, given that the magnitude of noise becomes comparable with the intrinsic fluctuations in node dynamics (node activity is constrained between 0 and 1 by the sigmoid activation function). Furthermore, we observe a small peak in synchrony at levels of noise close to 0.1, which could suggest stochastic resonance (Deco et al., 2009), in line with studies showing noise-induced synchronization due to the structure of the human connectome (Pang, Gollo, & Roberts, 2021). However, given that this peak does not correspond to the emergence of empirically relevant dynamics and there are no significant effects on the remaining network features, we do not explore this further. That said, given that model performance and dynamics are largely independent of noise, our results favor the emergence of metastable dynamics, as opposed to multistability (Kelso, 2012; Schirner et al., 2022). For this reason, we consider that the collective dynamics of our models are most likely a product of the emergence of metastability, shaped by edge-of-bifurcation dynamics at the circuit level (regulated by E–I homeostasis) and the topology of interareal interactions (constrained by the structural connectome).

Our results establish the key role of local E–I balance and its homeostasis for the emergence of empirically valid large-scale networks and dynamics in cortical networks. E–I homeostasis can also be a significant driver of the reemergence of large-scale FC properties in the lesioned cortex (Chakraborty et al., 2023; Páscoa dos Santos & Verschure, 2022; Páscoa dos Santos et al., 2023; Vattikonda et al., 2016). However, our earlier results (Páscoa dos Santos et al., 2023) suggest that scaling of inhibition, the homeostatic mechanism most commonly implemented in large-scale models (Abeysuriya et al., 2018; Castaldo et al., 2023; Deco et al., 2014; Hellyer et al., 2016; Naskar et al., 2021), is insufficient by itself for the recovery of metastable dynamics. Therefore, we investigate if the conjugation of multiple forms of homeostasis identified here can potentiate the reemergence of not only FC but also network dynamics following structural lesions. To do this, we apply structural lesions in the optimal model: including homeostasis of excitation (G^E), inhibition (c^EI), and intrinsic excitability (μ^E and σ^E) with C = 3.59, ρ = 0.12, and mean delay = 40 ms, corresponding to a conduction velocity (3.9 m/s) close to the values reported in previous studies (Lemaréchal et al., 2022; Petkoski & Jirsa, 2022). We extract activity before the lesion (healthy), after the lesion (acute), and after node parameters are adapted to restore local E–I balance (chronic). For more details on the simulation protocol and analysis (Páscoa dos Santos et al., 2023), refer to the Methods section and SI Methods. All p values presented in this section relate to the result of the Wilcoxon rank-sum test and are FDR corrected.

Starting with FC, we measure the Euclidean distance between healthy FC and FC in the acute and chronic periods (Figure 7A). We observed that, even though FC patterns are disrupted in the acute period (10.751 ± 4.900), significant recovery toward prelesion FC configurations (5.800 ± 1.573; p < 0.001) occurs in the chronic period. Consistent with stroke literature (Zhang et al., 2017), the model displays an acute decrease in their correlation between structural and functional networks (−10.97 ± 13.34%; p < 0.001; Figure 7B). However, this disruption was recovered in the chronic period through the action of the E–I homeostasis (1.33 ± 6.10%; p = 0.059). In addition, we observe an acute reduction of modularity (−6.00 ± 10.39%; p < 0.001; Figure 7C), consistent with that observed in stroke patients (Gratton, Nomura, Pérez, & D'Esposito, 2012; Siegel et al., 2018). In line with previous results (Páscoa dos Santos et al., 2023; Siegel et al., 2018), the recovery of the E–I balance drives the reemergence of modularity toward prelesion levels in the chronic period (2.41 ± 12.44%; p = 0.257). Therefore, our model can replicate the effects of lesions on FC and its subsequent recovery through the focal E–I homeostasis. We further explore the effects of lesions on functional complexity (Figure 7D), revealing a minimal decrease in the acute period (−1.86 ± 6.21%; p = 0.046) and a return to baseline after recovery of the E–I balance (−0.62 ± 3.81%; p = 0.240). Therefore, our results suggest that, as opposed to the FC-SC relation and FC modularity, FC complexity may be minimally disrupted by lesioning the connectome.

Beyond reproducing the effects of stroke, and subsequent recovery, on FC (Gratton et al., 2012; Páscoa dos Santos et al., 2023; Siegel et al., 2016; Vattikonda et al., 2016; Zhang et al., 2017), we explore how lesions disrupt network dynamics. Starting with synchrony, our results indicate that the system tends toward hypersynchrony in the acute period (9.31 ± 13.44%; p < 0.001; Figure 7E), in line with literature on the effects of brain injury (Hillary & Grafman, 2017). More importantly, synchrony is recovered toward prelesion levels in the chronic period (p < 0.001), although a small, but significant, difference from baseline remains (−2.19 ± 4.88%; p < 0.001). Regarding metastability, we observe a significant decrease postlesion (−5.30 ± 7.32%; p < 0.01; Figure 7F), followed by a significant recovery (p < 0.001) to prelesion levels in the chronic period (−1.08 ± 5.17%; p = 0.060). Finally, given the relationship between metastability and FCD (Deco et al., 2017; Schirner et al., 2022), we analyze FCD distributions in the simulated healthy, acute, and chronic periods. In line with Páscoa dos Santos et al. (2023), we observe an acute disruption, with the mean FCD correlation values shifting from 0.429 in the healthy period to 0.522 postlesion (KS distance = 0.356; Figure 7F). Importantly, the combination of multiple forms of homeostasis can drive the recovery of FCD distributions to prelesion levels (mean correlation = 0.434; KS distance = 0.021), as opposed to the plasticity of inhibition alone, as attested by our previous results (Páscoa dos Santos et al., 2023).

Finally, we have also performed the same analysis steps in models with the G^E or G^E +c^EI homeostasis, showing that such models are either unable to reproduce the studied poststroke disruptions in FC and network dynamics, or such disruptions were not recovered through the E–I homeostasis (Supporting Information Figures S31–S32).

In summary, the conjugation of multiple mechanisms of homeostasis not only supports metastable dynamics and empirical FC networks but can also drive their reemergence following lesions to the connectome, which is not the case in models with homeostasis of excitation or inhibition only. Therefore, the multifactor homeostasis of E–I balance can also be considered a robust means of self-organization in large-scale cortical networks.