Metastable Dynamics Emerge from Local Excitatory-Inhibitory Homeostasis in the Cortex at Rest

The human cortex displays highly metastable dynamics at rest, underlying the spontaneous exploration of large-scale network states. This metastability depends on edge-of-bifurcation dynamics at the circuit level, which emerge due to the local control of firing rates through multiple mechanisms of excitatory-inhibitory (E-I) homeostasis. However, it is unclear how the distinct forms of homeostasis contribute to the metastability of large-scale cortical networks. Here, we propose that individual mechanisms of E-I homeostasis contribute uniquely to the emergence of metastable dynamics and resting-state functional networks and test that hypothesis in a large-scale model of the human cortex. We show that empirical networks and dynamics can only be reproduced when accounting for multiple mechanisms of E-I homeostasis. More specifically, while the homeostasis of excitation and inhibition enhances metastability, the complementary regulation of intrinsic excitability ensures moderate levels of synchrony, maximizing the complexity of functional networks. Furthermore, the modulation of distance-to-bifurcation by the homeostasis of excitation and intrinsic excitability supports collective dynamics by compensating for strong input fluctuations in strongly connected areas. Altogether, our results show that cortical networks self-organize toward maximal metastability through the multi-factor homeostatic regulation of E-I balance, which controls local edge-of-bifurcation dynamics. Therefore, the functional benefits of combining multiple homeostatic mechanisms transcend the circuit level, supporting the rich spontaneous dynamics of large-scale cortical networks.

C omplex systems, such as the human brain, show rich collective dynamics beyond the behavior of their units (1)(2)(3).One of such collective behaviors in cortical networks is metastability, underlying the spontaneous exploration of the state-space of network configurations (4-10), which has a pivotal role in cognitive function (11)(12)(13)(14).That said, the key ingredients underlying the emergence of metastable dynamics are not well understood.In this context, the connectome can have a pivotal role, either due to symmetry-breaking caused by its intrinsic inhomogeneity (9) or due to the conduction delays between cortical areas (10,15).In addition, studies show that topologies typical of the mammalian connectome evolve in networks optimized for complexity (16) and how connectome topology maximizes the complexity of intrinsic activity patterns (17,18), establishing a link between the structure of cortical networks and their complex and metastable spontaneous dynamics.
On a different note, it has been suggested that a necessary condition for the emergence of metastability is that the units making up the system are poised near a phase-transition, or bifurcation point (6,19).Indeed, local edge-of-bifurcation dynamics are thought to be behind the rich patterns of activity observed in restingstate dynamics of the neocortex (6,17,(19)(20)(21)(22)(23)(24).More specifically, when multiple small-scale phase transitions co-occur across the cortex, they can self-reinforce and lead to the emergence of collective activity patterns (19).Therefore, while the emergence of metastable dynamics requires the interaction between cortical areas, which is constrained by the connectome, local circuit dynamics must also be key in its genesis.For this reason, we propose that excitatory-inhibitory (E-I) homeostasis, which regulates edge-of-bifurcation dynamics in local cortical circuits (19,(25)(26)(27)(28)(29), could be a fundamental principle for understanding the cortex and its collective dynamics (3,30).

Significance Statement
Experimental studies have consistently shown that cortical circuits maintain a precise homeostasis of excitatory-inhibitory (E-I) balance, thereby optimizing local dynamics.While it is well-established that multiple homeostatic mechanisms are involved in this local regulation, it remains unclear how each contributes to the large-scale dynamics of cortical networks.This study presents evidence that, through E-I homeostasis, the cortex can self-organize towards a regime of highly complex and metastable spontaneous dynamics.Crucially, we demonstrate that this results from the synergistic action of multiple homeostatic mechanisms.Our findings advance our understanding of E-I homeostasis as a process of self-organization, demonstrating its key role in the maintenance of metastable dynamics in large-scale cortical networks.

D R A F T
for a myriad of operations performed by cortical networks (27,(36)(37)(38).For this reason, cortical neurons are equipped with mechanisms of E-I homeostasis, which contribute to the maintenance of E-I balance in the face of perturbations in activity (39)(40)(41)(42).These mechanisms include synaptic scaling of excitatory (39,40,43) and inhibitory synapses (35,40,43,44) and the adaptation of intrinsic excitability (40,45,46).Furthermore, recent results suggest that the conjugation of multiple forms of E-I homeostasis not only ensures stable firing rates but can also regulate local activity toward criticality (28,47), providing a mechanism for the control of edge-of-bifurcation dynamics in cortical circuits (29).Therefore, the maintenance of balanced E-I interactions at the circuit level is likely a necessary condition for the emergence of the metastable patterns observed in large-scale cortical dynamics.
Several computational studies have explored how local E-I balance ensures flexible collective dynamics (48,49), which support the emergence of functional connectivity networks (49)(50)(51).Furthermore, while E-I balance shapes spatiotemporal dynamics across timescales, its influence is particularly strong in the ultra-slow dynamics characteristic of blood-oxygen-level-dependent (BOLD) signals (52).In addition, E-I homeostasis can be a mechanism of selforganization, driving the recovery of functional connectivity (FC) following lesions to the connectome (53)(54)(55), and the re-emergence of functional properties such as modularity (56).
There is a common caveat to large-scale modeling studies accounting for E-I balance: the exclusive focus on the homeostasis of inhibitory synapses.Importantly, this neglects the role of the combined action of multiple modes of homeostasis in the maintenance of circuit function, especially edge-ofbifurcation dynamics (28,47).Importantly, our results suggest that this might be relevant at larger scales, since the homeostasis of inhibition alone is not sufficient for the recovery of collective dynamics, such as metastability, after lesions to the connectome (56).Furthermore, while the combination of multiple modes of homeostasis ensures stable firing rates by modulating the distance to the bifurcation of local circuits (29), it is not clear how this effect impacts large-scale dynamics.
Here, we aim to investigate the hypothesis that relying on distinct mechanisms of homeostasis at the circuit level is essential to support the emergence of large-scale metastable dynamics in the human cortex.To do this, we study how largescale models with local E-I homeostasis reproduce empirical FC and FC dynamics (FCD), network dynamics, and the complexity of functional networks.More importantly, we explore the dependence of these collective patterns not only on E-I balance but also on the specific mechanisms employed to maintain it.In this context, our results demonstrate that the only models able to reproduce empirical FC, FCD, network dynamics, and complexity are the ones accounting for the combined homeostasis of excitation, inhibition, and intrinsic excitability.Finally, we also show that networks with these mechanisms of homeostasis can reorganize to recover not only functional networks but also metastable dynamics following focal lesions.Therefore, our results demonstrate that the maintenance of E-I balance not only underlies the emergence of metastable spontaneous dynamics in the cortex, but is also a robust mechanism of self-organization, especially when based on multiple complementary mechanisms of E-I homeostasis.

Results
To study how E-I homeostasis impacts the resting-state dynamics of the neocortex we constructed a large-scale model considering some of its key features together with local regulation of E-I balance (Figure 1).The dynamics of cortical areas are modeled with the Wilson-Cowan model of coupled excitatory-inhibitory populations with oscillations in the gamma-band (∼40 Hz).Long-range connectivity and conduction delays between areas follow empirical structural connectivity data derived from diffusion-tensor imaging.Critically, we implement multiple empirically identified modes of local E-I homeostasis (29,41,47), which ensure that activity remains close to a target firing rate (ρ) by regulating excitatory and inhibitory synapses and the intrinsic excitability of excitatory populations in each cortical area (29) (Figure 1).To optimize model performance, we explore different combinations of free parameters, namely the local target firing rate (ρ), the global coupling (C), which scales the structural connectivity, and the mean conduction delay, which depends on the length of white-matter tracts and the conduction velocity.For each simulation, we first evaluate if E-I homeostasis can maintain stable firing rates across all areas -i.e. the difference between mean local firing rates and the target (ρ) does not exceed 1% -keeping models with stable parameter configurations for further analysis.For a detailed exploration of model stability, refer to the Supporting Information (SI Text and Figs S1 -10).Finally, we extract simulated BOLD signals from stable networks by using the Balloon-Windkessel hemodynamic model (see Methods) and analyze network dynamics based on the BOLD time series.

Modulation of Distance to Bifurcation at Circuit Level Opti-
mizes Fit to Empirical Data.We start our analysis by studying how distinct forms of homeostasis, based on the plasticity of excitation, inhibition, and intrinsic excitability (see Methods and Figure 1) shape the performance of our model in reproducing empirical FC and FCD obtained from restingstate BOLD signals, which are strongly modulated by E-I balance (52,56).
We simulate the activity of a large-scale model of the neocortex with different combinations of free parameters (C, ρ, and mean delay) and modes of homeostasis (Figure 1).For each simulation, we compare empirical and simulated FC matrices and FCD distributions and extract a crossfeature fitting score to evaluate model performance.The score is calculated as rF C − M SEF C − KSF CD , where rF C is the Pearson's correlation between empirical and simulated FC matrices, M SEF C their mean squared error and KSF CD the Kolmogorov-Smirnov distance between FCD distributions (see Methods).First, it is relevant to point out that, similarly to model stability (Figs S1-10), the main parameters impacting model performance are C and ρ, while the mean delay of inter-areal conduction does not have a strong impact, provided it is long enough to avoid widespread pathological synchronization (Figs S11-14 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 mode of homeostasis, We extract simulated BOLD signals for further analysis and benchmarking against empirical resting-state recordings.

D R A F T
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 present the distribution of cross-feature fitting scores across mean delays for the optimal point of each model (Figure 2B, Table 1).The best scores are obtained for models with G E homeostasis (0.274±0.046) and G E + c EI homeostasis (0.253±0.094), followed by homeostasis of G E , c EI and µ E + σ E (i.e.All) (0.243±0.044).The values of rF C , M SEF C , and KSF CD corresponding to each optimal score can be consulted in SI Table 2.While the performance of the model with G E homeostasis was significantly better than the one with all modes of homeostasis (p=0.020;Mann-Whitney u-test, 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 (Fig S15 ), in line with theoretical (19,21,22) and experimental (28) 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 restingstate BOLD, depending on the interaction between global coupling (C) and target firing rate (ρ).In addition, empirical data is best approximated when accounting for modes of homeostasis that modulate the bifurcation point (i.e. the ones that modulate G E or σ E ) (29), so that areas with stronger inputs are poised farther away from the bifurcation (Fig S15).This effect compensates for the higher input fluctuations in strongly connected nodes (Fig S16 ), 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 (6,17,(19)(20)(21)(22)(23)(24) 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 ) (Fig S15 ), is essential not only for network stability (Fig S10 ) but also for the emergence of FC and FCD (Figure 2).

Resting-State Functional Connectivity Emerges in Networks
with Metastable Dynamics.Given that the metastability of large-scale cortical networks is shaped by circuit-level edgeof-bifurcation dynamics (6,19,22,49), controlled by E-I homeostasis (29,47), it is relevant to explore how distinct modes of homeostasis shape the collective dynamics of the resting-state cortex.To do this, we measure the synchrony and metastability (see Methods) of the simulated BOLD signals of models under different modes of homeostasis and combinations of C and ρ.Again, the conduction delays do not significantly impact either synchrony or metastability (Figs S17-18) 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 (Fig S19).To explore how synchrony is shaped by E-I homeostasis, we compare the optimal models of each mode 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.112) 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 modes 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 ρ (Fig S20 ), coinciding with the region of optimal fit to empirical data (Figure 2A).Furthermore, the modes 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 homeostasis of G E + c EI (0.  synchrony and metastability can only be reproduced when including all modes of homeostasis (Figure 3 and S21).In addition, we found a pronounced positive correlation between the metastability of models and their fitting scores, not only in the model with all modes of homeostasis (Figure 3C) but as a general principle (Fig S22).Therefore, our results confirm previous studies proposing that the resting-state dynamics of cortical networks might correspond to a regime of maximal metastability (22,49,57).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 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 selforganize toward a regime of maximal metastability where the empirical resting-state networks emerge.
Functional complexity is shaped by E-I balance, metastability, and connectome topology.The complexity of spontaneous activity in the human brain is strongly related to the topology of its structural networks (16)(17)(18).However, the architecture of the connectome is likely not the only source of complexity in cortical networks.Given the relevance of 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 modes of E-I homeostasis impact FC complexity, using the methodology presented in (18) (see Methods).
Our results on FC complexity suggest that, for all modes of homeostasis, it is maximized in the region of the parameter space corresponding to heightened metastability and optimal fit to empirical FC and FCD (Fig S24).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 modes of homeostasis 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 modes of homeostasis is necessary not only to reproduce synchrony and metastability but also functional complexity.Importantly, we found a strong correlation between metastability and complexity, not only for the model with all modes of homeostasis (Figure 4B) but also across models (Fig S22).This suggests that heightened metastability of resting-state brain dynamics (22,49) results in highly complex activity patterns that reflect on the variability of FC.More importantly, this requires the combined action of multiple mechanisms of E-I homeostasis at the circuit level.
Having investigated the relationship between E-I homeostasis and complexity compare complexity in models with the empirical and shuffled structural connectomes to assert the rel- evance of network topology (16)(17)(18).Since shuffling preserves the distribution of edge weights, the complexity of the original and shuffled structural matrices is the same (Fig S25).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 (16)(17)(18), the complexity of functional interactions also depends on the topology of structural connectivity.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 modes of E-I homeostasis does extend beyond the local circuit level (42,47), shaping the collective dynamics of large-scale cortical networks.

Flexible Network Dynamics Emerge Independently of Noise.
Synaptic transmission is known to be noisy in brain networks (58).This variability can have a relevant role in spontaneous brain dynamics (58,59), underlying phenomena such as stochastic resonance (24,60,61).In addition, systems with multiple weakly stable attractor states display noisedriven transitions between attractors, corresponding to multistable dynamics (20,24), which can be confounded with metastability (62).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 (24,62).For that reason, we explore the behavior of our network under different levels of noise (see Methods).Given that the model with all modes 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 Methods).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 (Fig s26 ), we present the average results across mean delays (Figure 5).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 to 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 (60).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 (24,62).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 inter-areal interactions (constrained by the structural connectome).

Multiple Modes of Homeostasis Drive the Re-emergence of Functional Networks and Collective Dynamics after Struc-
tural Lesions.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 (53)(54)(55)(56).However, our earlier results (56) suggest that scaling of inhibition, the homeostatic mechanism most commonly  implemented in large-scale models (48)(49)(50)(51)(52), 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 re-emergence 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 close to 3.9 m/s (63).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 (56) refer to the Methods and SI Methods.All p-values presented in this section relate to the result of the Wilcoxon-Ranked 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 6A).We observed that, even though FC patterns are disrupted in the acute period (10.751±4.900)significant recovery toward pre-lesion FC configurations (5.800±1.573;p<0.001) occurs in the chronic period.Consistent with stroke literature (64), the model displays an acute decrease in their correlation between structural and functional networks (-10.972±13.338%;p<0.001) (Figure 6B).However, this disruption was recovered in the chronic period through the action of E-I homeostasis (1.334±6.103%;p=0.059).In addition, we observe an acute reduction of modularity (-5.999±10.387%;p<0.001) (Figure 6C), consistent with that observed in stroke patients (65,66).In line with previous results (56,66), the recovery of E-I balance drives the reemergence of modularity toward pre-lesion levels in the chronic period (2.412±12.435%;p=0.257).Therefore, our model can replicate the effects of lesions on FC and its subsequent recovery through focal E-I homeostasis.We further explore the effects of lesions on functional complexity (Figure 6D), revealing a minimal decrease in the acute period (-1.857±6.213%;p=0.046) and a return to baseline after recovery of E-I balance (-0.620±3.808%;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 functional connectivity (53,56,64,65,67), we explore how lesions disrupt network dynamics.Starting with synchrony, our results indicate that the system tends toward hypersynchrony in the acute period (9.313±13.435%;p<0.001) (Figure 6E), in line with literature on the effects of brain injury (68).More importantly, synchrony is recovered toward pre-lesion levels in the chronic period (p<0.001),although a small, but significant, difference from baseline remains (-2.19±4.881%;p<0.001).Regarding metastability, we observe a significant decrease post-lesion (-5.304±7.32;p<0.01) (Figure 6F), followed by a significant recovery (p<0.001) to pre-lesion levels in the chronic period (-1.08±5.165%;p=0.060).Finally, given the relationship between metastability and FCD (22,24) (Fig S22 ), we analyze FCD distributions in the simulated Healthy, Acute, and Chronic periods.In line with (56), we observe an acute disruption, with the mean FCD correlation values shifting from 0.429 in the Healthy period to 0.522 post-lesion (KS distance=0.356)(Figure 6F).Importantly the combination of multiple forms of homeostasis can drive the recovery of FCD distributions to pre-lesion levels (mean correlation=0.434;KS distance=0.021), as opposed to the plasticity of inhibition alone, as attested by our previous results (56).
In summary, the conjugation of multiple modes of homeostasis not only supports metastable dynamics and empirical FC networks but can also drive their re-emergence following lesions to the connectome, which is not the case in models with homeostasis of inhibition only.Therefore, the multifactor homeostasis of E-I balance can also be considered a robust means of self-organization in large-scale cortical.

Discussion
Building on the observation that multiple modes of E-I homeostasis are essential for the regulation of local circuit dynamics (47), we extend this hypothesis to large-scale networks, exploring how distinct modes of E-I homeostasis impact the global dynamics and functional networks of the human cortex.Our results demonstrate that the homeostatic regulation of local E-I balance supports the emergence of resting-state FC and FC dynamics, particularly when involving forms of homeostasis that modulate the distance to the bifurcation between damped and sustained oscillations (29) (Fig S15).In addition, we show that the resulting heterogeneity in the distance to the bifurcation across cortical areas supports network stability and the emergence of empirical FC and FCD by compensating for heterogeneities in the magnitude of input fluctuations.We further demonstrate that the optimal working points of models generally correspond to moderate levels of synchrony and heightened metastability, confirming that cortical networks might operate in a regime of enhanced metastability.More importantly, empirical levels of synchrony and metastability can only be replicated by models that incorporate multiple modes of homeostasis.Furthermore, we study the complexity of large-scale functional interactions, stressing that it not only depends on the topology of structural networks (16)(17)(18) but also on local E-I balance.Importantly, empirical levels of complexity were more closely approximated by the model with all forms of homeostasis enabled.Altogether, we confirm the hypothesis that the benefits of the interaction between multiple forms of E-I homeostasis extend beyond the local circuit level (40)(41)(42)(43)47), contributing to the regulation of large-scale dynamics toward a regime of maximal metastability and complexity, which underlies the FC networks observed in the resting-state cortex.Finally, we also demonstrate that the combined action of multiple modes of homeostasis can drive the re-emergence of not only FC but also network dynamics following focal lesions to the connectome (56), evidencing the robustness of E-I homeostasis as a self-organization process in cortical networks.

D R A F T
To our knowledge, we are the first to explore the combined action of diverse modes of E-I homeostasis in large-scale models, given that most studies focus on the homeostasis of inhibition (48)(49)(50)(51)(52)56).At the circuit level, the synergy between the synaptic scaling of excitation and inhibition and the plasticity of intrinsic excitability are essential for the maintenance of not only stable firing rates but also edge-ofbifurcation dynamics (i.e.criticality) (42,47,69).Our results add to this perspective by suggesting that these benefits extend to larger scales, supporting network dynamics and FC networks.More specifically, our results show that networks with synaptic scaling of excitation and inhibition reach heightened levels of global metastability, characteristic of the cortex at rest (6,10,19,22,24).Regarding the plasticity of intrinsic excitability, it is not yet clear from experimental data if it acts through the modulation of firing thresholds only (µ E ) (40) or the concerted adaptation of both the threshold and slope of neuronal excitability (µ E +σ E ) (45,46,70).Here, we suggest that the second option is more plausible given its contribution to ensuring the moderate levels of synchrony observed in empirical data, while still allowing for highly metastable dynamics.In this line, it should be highlighted that recent results suggest that, while the plasticity of intrinsic excitability might not be mobilized in the adult cortex, being instead pivotal during development (46).However, it has also been suggested that homeostasis of intrinsic excitability can also be mobilized following substantial perturbations that cannot be compensated for by the homeostasis of excitation and inhibition (41).This aligns with evidence from the poststroke brain, where processes typical of the developmental phases, such as neurogenesis, are also recruited for recovery (71).Therefore, we argue that the diverse mechanisms of homeostasis have specific contributions to cortical dynamics at both the circuit and global scales which may play a role in both the healthy and diseased brain.
In light of our results, it is relevant to reflect on how exactly local E-I balance contributes to the large-scale dynamics of the cortex.Here, we focus on the edge-of-bifurcation theory for cortical dynamics (6,(19)(20)(21)(22)(23)(24).Under this perspective, local circuits at the edge-of-bifurcation can undergo spontaneous phase transitions, which, when coincident across cortical areas, can self-reinforce, leading to the emergence of the collective behaviors typical of the cortex (19).In our model, these collective dynamics are strongly regulated by the maintenance of stable mean firing rates (ρ) through E-I homeostasis (29).For example, increases in global coupling (C), which may push node dynamics beyond the bifurcation, can be compensated by lowering ρ.For this reason, model performance and dynamics are heavily dependent on the interaction between C and ρ.On the other hand, the heterogeneity in the in-degrees of the connectome not only leads to stronger inputs in highly connected nodes, but also higher input fluctuations (Fig S16).Therefore, if all cortical areas were to be equally distant from the bifurcation, the highly connected ones would be more likely to engage in sustained oscillations due to strong input fluctuations that push them beyond the bifurcation.However, this effect can be compensated for by forms of homeostasis (i.e.excitation and excitability slopes) that maintain stable firing rates by scaling the distance to the bifurcation (Fig S15) (29), poising highly connected nodes farther away from the bifurcation point.This adjustment is critical to ensure that cortical areas exhibit edge-of-bifurcation dynamics regardless of the variability of their inputs.More so, this aligns with previous studies showing that distance to criticality varies across cortical areas (72) and that global metastable dynamics are better reproduced in models with heterogeneous distances to the bifurcation (22).Importantly, we demonstrate that such heterogeneity can emerge naturally from the maintenance of stable population firing rates via multiple modes of homeostasis.In this line, we suggest that future studies investigate how the heterogeneity in distance-to-bifurcation across cortical areas (22,72) relates to their structural connectivity and in microcircuitry (73).
Our results indicate that empirical FC dynamics emerge in a state of maximal metastability (Fig S22 ), confirming that FCD reflects the exploration of the state space of network configurations potentiated by metastable brain dynamics (9,10,19,22,24).However, it has been argued that metastability is generally derived from abstract Hopf-bifurcation models, with no direct relationship to physiological processes (24).In contrast, we demonstrate that these theoretical concepts are directly applicable in a model of the neocortex where the generation of oscillations is based on E-I interactions and bifurcation control implemented by E-I homeostasis.This confirms that the concept of metastability is compatible with both the architecture of the human cortex and its processes of self-organization.Importantly, our results also show that the collective dynamics of the human cortex most likely reflect metastability and not multistability (62), making a further argument in favor of the hypothesis that metastability is one of the key dynamical features of large-scale cortical networks (6,9,10,24).
Our analysis further shows a relationship between metastable dynamics and the complexity of functional interactions.In complex systems, the interactions between units are dictated by both the structural scaffold and circuit dynamics (74).In the cortex, the architecture of the connectome plays an important role in the emergence of complex spontaneous activity (16,17) and its graph properties are required to support the complexity observed in empirical FC networks (18).By exploring the relationship between complexity and metastability, we demonstrate that FC complexity not only reflects the topology of the connectome, but is also tightly related to metastable dynamics (Figures 4B and  S22), suggesting it reflects a balanced exploration of network states (9, 10).When averaged over long time scales, this exploration leads to a distribution of FC weights that is closer to uniformity, revealing a dynamic balance between integration and segregation.For this reason, we propose that the complexity of functional interactions can be another signature of local E-I balance, which allows for the emergence of metastable collective dynamics (6,19,21,22) through the regulation of local edge-of-bifurcation dynamics (29).That said, networks with homeostasis of excitation (G E ) or excitation and inhibition (G E + c EI ), associated with high metastability (Figure 3B), display levels of complexity considerably lower than found in empirical data (Figure 4B).However, such models also show overly synchronized dynamics (Figure 3A).Conversely, when the plasticity of intrinsic excitability (µ E + σ E ) is introduced, synchrony is maintained at empirical levels and complexity enhanced.Therefore, we propose that the plasticity of intrinsic excitability has an important role in maintaining complex interactions by avoiding hyper-synchrony while still allowing for highly metastable dynamics.

D R A F T
To further validate our model, we also examine the role of E-I homeostasis in the self-organization of the cortex following focal lesions, following the framework of (56).We replicate disruptions in FC properties such as SC-FC coupling (64) and modularity (66), characteristic of the post-stroke brain, which are then recovered through cortical re-organization, driven by standard mechanisms of homeostatic plasticity, in line with previous results (56).However, the model in (56) was equipped with the plasticity of inhibition only, which proved insufficient for the recovery of metastable dynamics and FCD.Here, we demonstrate that networks with multiple modes of E-I homeostasis (excitatory and inhibitory synapses and intrinsic excitability) can not only replicate the disruptions and subsequent recovery of macroscale FC properties, but also return metastability and FCD to pre-lesion levels through the restoration of E-I balance (Figure 6).That said, not all mechanisms of homeostasis may be necessary for the recovery of network dynamics.Therefore, we apply the same lesion protocol in networks with either G E or G E + c EI homeostasis (Figs S27-28).Critically, those are either unable to recover network dynamics or do not reproduce empirical post-stroke deficits in FC (64)(65)(66).Therefore, while the homeostasis of inhibition might strongly contribute to the recovery of FC networks, as demonstrated in (56), the plasticity of excitation and intrinsic excitability is critical for ensuring that network dynamics can be recovered following damage to structural networks.
Our results show that, by imposing stronger constraints on the local circuitry of cortical networks, E-I homeostasis potentiates the flexible exploration of macroscale network states.While this effect appears paradoxical, it is in line with the concept of "constraints that deconstrain" (75).This idea stems from the work of Kirschner and Gerhart (76), who suggest that some constraining mechanisms are preserved through evolution because they deconstrain other processes that are advantageous for the organisms.In neuroscience, this concept has been applied to explain how the layered architecture of brain networks (77, 78) ensures both robust motor control and high adaptability in behavior (75).We propose that the tight control of E-I balance by multiple mechanisms of homeostasis is another example of a "constraint that deconstrains", by supporting the emergence of metastability which endows the cortex with the flexibility needed to explore its state-space during rest (6)(7)(8)(9)(10)19).
There are, however, some limitations to point out in our model.The first is that we do not explore homeostatic plasticity at the level of inhibitory neurons.Recent studies have shown that the plasticity of excitatory connections (28) or intrinsic excitability (42) in parvalbumin-positive (PV) interneurons can also contribute to E-I homeostasis.Therefore, future models should consider the implementation of E-I homeostasis in inhibitory neurons, delving into the interaction between the homeostatic set-points of excitatory and inhibitory populations.Another limitation relates to the fact that, while synaptic scaling of excitation is ubiquitous across cortical layers (40,43,44,46,79), other modes of plasticity, such as the scaling of inhibition, might be absent from deeper layers (40,43,46,79,80).Therefore, we suggest that E-I homeostasis should be explored in large-scale models accounting for the laminar structure of the neocortex (81,82).Furthermore, we implement a homogeneous target firing rate across the network, in line with previous modeling studies (48)(49)(50)(51).While in vivo studies of E-I homeostasis are mostly based on the activity of visual or somatosensory cortices (28,79), different cortical areas might be regulated toward area-specific set points, possibly following hierarchical gradients of cortical organization (73).For this reason, future models should explore the possibility of not only implementing heterogeneous target firing rates, but also the interaction between cortical hierarchy and the homeostasis of E-I balance.Finally, we tune the Wilson-Cowan model to generate gamma oscillations, which have been related to fluctuations in BOLD signals (21,52,60,83) and transient rhythms across frequency bands (15).While there is empirical evidence for gamma resonance through the recurrent interactions between PY and PV neurons (84,85), some studies suggest that cortical networks also generate alpha oscillations in deeper layers (86,87).Therefore, the extension of our framework to the laminar architecture of the cortex (81, 82) might also provide new insight into how E-I balance interacts with different cortical rhythms, fast and slow inhibition (29) and how it shapes cross-frequency interactions.
To conclude, our results demonstrate the critical role of firing-rate homeostasis, which regulates edge-of-bifurcation dynamics at the circuit level, as a key piece underlying the selforganization of resting-state cortical activity toward a regime of heightened metastability where resting-state networks emerge.More importantly, we show that this reflects the joint action of multiple homeostatic mechanisms employed to maintain E-I balance, suggesting that their influence extends beyond the neuronal and circuit level, underlying the rich spontaneous dynamics of cortical networks.

Materials and Methods
Structural Connectivity Data.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 Anatomic Automatic Labeling (AAL) atlas, yielding a 78x78 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 (52).
BOLD fMRI Data.To analyze network dynamics and FC properties in empirical data, we use resting-state blood-oxygenation-leveldependent (BOLD) time series from 99 healthy unrelated subjects from the HCP dataset.Data was processed using the pipeline described in (52), yielding 99 time series of size 78 areas x 1200 TR, roughly corresponding to 14.5 minutes of data for each subject (TR=0.72 s).

Large-Scale Model.
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 (88,89):

D R A F T
Where c EE represents the recurrent excitatory coupling, c EI the coupling from inhibitory to excitatory populations and c IE the excitatory coupling unto inhibitory populations.G E is a parameter introduced in (29), 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 SI Table 1.Unless stated otherwise, ξ(t) represents Gaussian noise with mean 0 and variance 0.01.Here, I ext i describes the incoming input to each node i from the rest of the network and can be written as: where C is a scaling factor referred to as global coupling, W ij is the structural connectivity matrix and τ ij represents the conduction delay between nodes i and j.Here, similarly to (51,52,56), 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: Excitatory-Inhibitory Homeostasis.As opposed to the common approach to E-I homeostasis (48-52, 56), we do not implement dynamical homeostasis.Instead, we follow the framework developed in (29) 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 ) (28, 79), in line with the governing rules of E-I homeostasis (41).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 (41,42,47), 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: Here, we explore the behavior of models with the following distinct modes of E-I homeostasis, based on experimental results (39-41, 43, 45-47) and described in detail in ( 29): • 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 .• All modes of 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 mode of homeostasis can be consulted in SI Methods.

Hemodynamic Model.
From the raw activity of the excitatory populations, we extract simulated BOLD signals by using a forward hemodynamic model (90,91), described in detail in (56).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.72s to align simulated signals with empirical data.
Network Stability.For each type of homeostasis and combination of hyper-parameters (C, ρ, and mean delay), we evaluate if the homeostatic set point -i.e. a mean firing rate of ρ for all nodes -can be maintained across the network.To do this, we run 15-second simulations and compute the mean firing rate of each node in the last 10 seconds of the simulation.We consider the given network fixed point to be unstable if there is any node in the network for which the mean firing rate differs from ρ by more than 1% (i.e. The nature of model instabilities is examined in detail in SI Text.While it is possible to evaluate the stability of a network analytically (see (92)), the characteristic polynomial of our model (from which the eigenvalues can be extracted) does not have a closed-form solution due to the highly non-linear nature of the Wilson-Cowan equations.For this reason, we base our exploration of model stability on the analysis of the numerical solution of the system.• Static FC: 78x78 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 F C ) and mean squared error (M SE F C ) between their upper triangular elements.

Model
• FC Dynamics (FCD): matrix of correlations between the upper-triangular part of FC matrices computed in windows of 80 samples with 80% overlap (22,93).Model results are compared to empirical data by computing the Kolmogorov-Smirnov distance (KS F CD ) between the distributions of values in the respective FCD matrices.
BOLD signals are filtered between 0.008 and 0.08 Hz before computing FC and FCD.In addition, we compute a crossfeature fitting score as r F C − M SE F C − KS F CD to evaluate the performance of models in simultaneously reproducing static and dynamic properties of empirical FC (52,56).All simulations were performed by solving the system's equations numerically using the Euler method with an integration time step of 0.2 ms.

Synchrony and Metastability.
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) (94,95) of BOLD timeseries, filtered between 0.008 and 0.08 Hz.More specifically, the KOP (R(t) in Eq 5) represents the degree of phase alignment within a set of coupled oscillators at a given point in time and can be calculated as: e iθn(t) [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.Functional Complexity.Previous studies have demonstrated that the architecture of structural networks has a significant impact on the complexity of the functional interactions between nodes (16)(17)(18).To explore this, we measure the complexity of simulated FC matrices using the method introduced in (18).In short, we quantify complexity by measuring how strongly the distribution of FC weights deviates from an equivalent uniform distribution using the following formula:

D R A F T
where M represents the number of bins used to evaluate the FC distribution.The term C M = 2 m−1 m 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, (18) demonstrate that their formulation is more robust to changes in bin sizes, compared to other methods quantifying the entropy or variance of FC matrices.Here, we use a bin size of 0.05.

Effects of Noise.
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 stable simulations.For each simulation, we vary the variance of Gaussian noise (ξ(t)), corresponding to 29 logarithmically spaced values between 10 −5 and 10 − 1 3 .We then analyze model performance, synchrony, metastability, and functional complexity across all models with varying levels of noise.
Lesion Simulation Protocol.Following the framework applied in (56), we investigate the effects of focal lesions in model dynamics and functional interactions by following the following protocol.Using the model with all modes of homeostasis at the optimal point (C = 3.59; ρ = 0.12), we compute the homeostatic parameters by using the methodology from (29) (SI Methods).From the balanced model, we extract 30 minutes of BOLD signal corresponding to activity in the pre-lesion (i.e.healthy) period.Then, we apply a structural lesion by removing all connections to and from a chosen node and extract 30 minutes 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 pre-lesion homeostatic values as G E 0 , c EI 0 and µ E 0 (SI Methods for the equations).From the re-balanced model, we extract 30 minutes 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 modes of homeostasis, we also report the effects of lesions on models with G E and G E + c EI homeostasis (Figs S27-28).
Code.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.

D R A F T
In the second case, E-I homeostasis modulates both the firing threshold µ E and the sensitivity of the activation function σ E in a coordinated manner so that the value of F E (0), or the baseline firing rate under no external input, is maintained.To do this, we consider that both parameters can be related to each other as σ E = Kµ E , where K is a constant.For more details on the derivation, refer to (1).Here, we use the default parameters to compute K, so that the values of both parameters can be obtained through the following expression: Synaptic Scaling of Excitatory and Inhibitory Synapses.Considering that plasticity of excitation and inhibition operate under the same timescale, it is possible to estimate the steady-state values of both parameters given the initial values G E 0 and c EI 0 .For a detailed derivation of the expressions and the extension of the framework to the case of different timescales, please consult (1).
Synaptic Scaling of Excitation and Inhibition and Plasticity of Intrinsic Excitability.Following the previous implementations, we implement homeostasis of excitation and inhibition in conjunction with the two forms of homeostasis of intrinsic excitability.First, for plasticity of µ E , we have: Similarly, for plasticity of µ E and σ E , we obtain: corresponds to a stable fixed point or spiral, the mean activity of the excitatory population is precisely equal to the fixed-point r E .However, in systems in a limit-cycle regime, the average firing rate differs from r E f ixed .For this reason, in (1), we derive a method to obtain an estimate of the steady-state model parameters as a function of ρ and I ext , based on the following iterative algorithm, with ϵ = 10 −10 .For more details, refer to (1).
1. Define upper and lower limits for r E f ixed (r E lower = 0 and r E upper = 1).

For r
, compute the homeostatic parameters of the system under I ext using the equations in the previous sections with r E = r E f ixed .
3. Solve the system numerically for 5 seconds and compute the average firing rate ⟨r E ⟩ from the last 2 seconds of activity.

4.
If |⟨r E ⟩ − ρ| < ϵ, take r E f ixed as the fixed point corresponding to ρ and save the homeostatic model parameters.

Analysis of FC Properties in Lesion Simulations.
To analyze the impact of structural lesions on FC and the subsequent recovery through E-I homeostasis, we follow the procedure introduced in (2).To analyze changes in FC patterns, we make use of FC distance, which roughly quantifies the magnitude of disruption, following (2,3).Furthermore, to measure disruptions in the macroscale architecture of FC, we measure the correlation between structural and functional connectivity (4) and modularity (5,6), both of which are impacted in stroke patients.Furthermore, the recovery of modularity is thought to be relevant for cognitive function (6).Below, we describe how each of these metrics is computed from simulated FC matrices.The simulation protocol is described in the Methods section of the main text.
FC Distance.To measure the dissimilarity between FC matrices at baseline, acute, and chronic periods, we follow (3), defining FC distance as the Frobenius norm of the difference between two given matrices.
FC-SC Correlation.Given the results of (4), showing a decoupling between functional and structural connectivity in stroke patients, correlating with motor function, we test this biomarker at baseline, acute, and chronic periods, by computing the Pearson's correlation coefficient between the upper triangles of FC and SC matrices.
Modularity.Modularity measures the degree to which a network follows a modular structure, with dense connections within functional clusters and sparser ones between them.Modularity (Q) was calculated using the formula defined in (7): where M is a set of non-overlapping modules (groups of nodes) in the network and euv is the proportion of edges in the network that connect nodes in module u with nodes in module v. Similarly to previous studies, (2,6), we define modules a priori to avoid biasing the modularity measure by using modularity maximization to detect community structure.That said, modules were from the empirical FC data, by using a clustering algorithm described in (2), resulting in 6 clusters.Following the original formulation of modularity (7), FC matrices were transformed into unweighted graphs by applying a density threshold, through which only a percentage of the strongest connections are kept and their weights set to 1. Lesioned nodes were removed from the network before computing modularity, similarly to (2,6).

Francisco Páscoa dos Santos and Paul FMJ Verschure
Model Stability for Different Modes of Homeostasis.For each type of homeostasis, we run short simulations (15 seconds) with all combinations of parameters and evaluate the deviation of the mean node activity from the target firing rate ρ.If any node deviates by more than 1%, we consider the fixed point unstable.
To start, we illustrate the impact of each parameter on model stability for models under G E homeostasis and discuss the principles underlying the stability or instability of networks in different regions of the parameter space represented by the combination of the free parameters (C, rho and mean delay).We start by using models with G E homeostasis to illustrate the nature of instability in our model (Figure S1).However, the general principles apply to other types of homeostasis, as will be illustrated later.The respective parameter spaces can be consulted in Figure S6.
The first main conclusion is that, as the target firing rate ρ increases, the models become increasingly more unstable and, for any value of ρ higher than 0.15, no combination of parameters resulted in a stable system.Given the results of Chapter 4, this result is easily interpreted.As more nodes enter the limit cycle regime (S15) the presence of sustained oscillations leads to a tendency of nodes to synchronize in a self-reinforcing manner, generating run-away activity and resulting in the instability of the fixed point solution.Therefore, to guarantee the stability of models, the target firing rate ρ should be close to the threshold after which some nodes enter the limit cycle regime, corresponding to ρ = 0.12.
To delve into this topic in more detail, we examine the behavior of the model when ρ = 0.14.While most of the parameter space corresponds to an unstable solution of the system, systems can be stable if the global coupling is high enough (Figure S1).To study this effect, we simulated models with the same ρ (0.12) and mean delay (10 ms), but with low (C = 0.75) and high (C = 7.5) global coupling (Figures S2 and S3, respectively).For each, we ran simulations with the predicted homeostatic G E , no noise, and different initial conditions for node activity r E (0, ρ or 2ρ).In addition, we ran simulations with dynamical and different homeostatic timescales, to ensure that the instabilities were inherent to the model and not a consequence of errors in the estimation of homeostatic parameters.Starting with the lower global coupling (C = 0.75) (Figure S2), we observe that, regardless of the initial conditions, the model with predicted G E tends toward global synchronization and, since long-range connections are excitatory, the oscillations self-reinforce and quickly reach activity levels close to saturation.How fast this happens depends on the initial conditions.Can this be fixed by implementing explicit homeostatic plasticity of G E , instead of using the predicted weights?Our results indicate that such models, having the same tendency to synchronize, will exhibit short periods of high activity departing from ρ, which are then quenched by homeostatic plasticity.In this case, the speed of this quenching depends on the timescale of homeostatic plasticity (Figure S2).For this reason, we consider these models to still be unstable, indicating that the instability is an inherent property of systems with this combination of parameters and not due to errors in the estimation of parameters.Indeed, when observing the behavior of individual nodes in the homeostatic fixed points (Figure S2), we see that 52s of them are in the limit cycle regime, which explains the tendency of the model for over-synchronizing.
But why are models with ρ = 0.14 when the global coupling is increased?In Figure S3, we show one such case where, regardless of the initial conditions, all nodes in the model converge toward the target firing rate.Following the principles introduced before, we show that the homeostatic solution to this system corresponds to only one node in the limit cycle regime.
This phenomenon occurs due to how model dynamics are shaped by G E homeostasis (Figure S3).Shortly, for a given value of ρ, as the external current is increased, the models transition from a limit cycle to a stable spiral regime.Therefore, given that I ext scales with the global coupling, increasing the global coupling will bring nodes out of the limit cycle regime, contradicting the tendency to hypersynchrony.This effect is only present for some types of homeostasis (Figure S15) which has implications for the stability of models.Conversely, in models with µ E or c EI homeostasis, when ρ ≥ 0.12 all nodes are necessarily in the limit cycle regime and, therefore, have a greater tendency toward instability.
That said, it should be pointed out that, for target firing rates lower than 0.12, the general effect of increasing the global coupling is the opposite, disrupting the stability of models.However, in this case, the mechanisms of instability are different.
In Figure S4 we present the behavior of the model in one such case, where the model settles in different fixed points depending on the initial conditions.If the initial r E is higher than the target firing rate, the model will settle in a stable configuration where mean activity is higher than ρ.Conversely, the opposite happens when the model is initialized with r E = 0.More importantly, the target fixed point of the system is only stable when the initial conditions correspond to r E = ρ across all nodes.These results indicate that, for this combination of parameters, the model is bistable, with the target corresponding to the saddle point of the system, thus being unstable.Accordingly, when implementing explicit homeostatic plasticity (Figure S4), the system will constantly switch between the two stable attractors.Since the target firing rate is different from ρ in both, homeostatic plasticity will be permanently frustrated, leading to the constant switching between attractors.For this reason, even though the model can be stable in this case, the stable solutions do not correspond to the target firing rate and, therefore, we consider this to represent an unstable solution of the system.While we performed this analysis in the model with G E homeostasis, the same principles apply, for example, to models with G E + c EI homeostasis (Figure S9).
Our results so far suggest that the main parameters shaping model stability are the global coupling and ρ.However, there are instances where the mean delay also plays a significant role.As an example, we take C = 4.5 and ρ = 0.11 (close to the bifurcation after which nodes with weaker inputs start entering the limit cycle regime) and examine model activity for 4 different values of mean delay (Figure S5).Starting with the model with no delays (mean delay=0), we observe that the system tends toward hyper-synchrony, similar to what was described for high values of ρ.As the mean delay is increased to higher values, this tendency disappears and the model can settle in the target fixed point.Here, our results demonstrate the role of delayed interactions in avoiding widespread synchronization, through the disruption of phase-relations between oscillating systems.However, we observe a range of mean delays centered around 26 milliseconds, for which the model can also be unstable (Figure S5).In this case, the target fixed point is only weakly stable, since, when approached from one direction, dynamics settle at ρ.However, if the initial conditions are higher than the target, some of the nodes engage in sustained oscillations potentiated by their recurrent interactions.While this case is not as severe as the model with no delays, it still represents a solution of the system that does not correspond to the target ρ.But why does this happen specifically for this range of delays?For this combination of C and ρ, most of the nodes have a natural frequency of oscillation between 45 and 60 Hz (Figure S5).More importantly, the nodes with the lower input, which are the ones closer to the limit cycle regime (S15), oscillate at around 45 Hz, corresponding to a period of ∼22.5 ms.Therefore, when the mean delay is close to this value, nodes in the system can still synchronize, albeit with a phase shift close to a full cycle, which is sufficient to throw the system outside of the desired equilibrium.Accordingly, as the mean delay is further increased, the model becomes stable again (Figure S5).
While we performed this analysis in the model with G E homeostasis, the same principles apply, for example, to models with c EI (Figure S7) or G E + c EI (Figure S8) homeostasis.
Finally, for each mode of homeostasis, we plot the percentage of stable simulations as a function of each free parameter (C, ρ, and mean delay) and the percentage across all simulations (Figure S10).S10D).We suggest that this results from the modulation of the bifurcation point toward higher ρ in response to increases in external inputs, which is characteristic of these modes of homeostasis ( 1 Colors represent the correlation coefficient between empirical and simulated FC matrices for each combination of C and ρ.

Fig. S12. FC fitting of models with different modes of homeostasis
Colors represent the mean squared error between empirical and simulated FC matrices for each combination of C and ρ.

Fig. 1 .
Fig. 1.Model Architecture and Excitatory-Inhibitory Homeostasis To simulate the restingstate dynamics of the cortex, we built a large-scale model based on the healthy human connectome.Cortical areas are modeled with the Wilson-Cowan model of coupled excitatory-inhibitory populations, tuned to oscillate at 40 Hz.Based on experimental studies, we implement distinct forms of local E-I homeostasis to maintain the mean excitatory activity at a target firing rate (ρ) based on the scaling of excitatory synapses (G E ), inhibitory synapses (c EI ), or the plasticity of intrinsic excitability, either through the adaptation of the firing threshold (µ E ) or the threshold and slope (µ E and σ E ) of the input-output function of neural populations.In addition, we explore models with different combinations of these homeostatic mechanisms and of the free parameters (global coupling C, target firing rate (ρ), and mean delay).

Fig. 3 .
Fig. 3. Networks with multiple modes of homeostasis reproduce empirical network dynamics (A) Comparison between synchrony in empirical data and models at the optimal point (best cross-feature fitting score) for each mode of homeostasis.Simulated distributions correspond to models with optimal C and ρ and all values of mean delay yielding a stable solution.Brackets represent a significant difference between models and empirical data from a Mann-Whitney U-test (* p<0.05; ** p<0.01; *** p<0.001).All pvalues were FDR-corrected.(B) Same as (A), for metastability.(C) Relationship between fitting score and metastability in the model with homeostasis of excitation, inhibition, and intrinsic excitability (µ E and σ E ).Gray dots show the results of all simulations while red crosses represent simulations at the optimal point (C = 3.59, ρ = 0.12).The vertical bar and shaded area display the mean and standard deviation of metastability in empirical BOLD signals.

Fig. 4 .
Fig. 4. Functional complexity reflects E-I homeostasis, metastability, and connectome topology (A) Comparison between functional complexity in empirical data and models at the optimal point (best cross-feature fit) for each mode of homeostasis.Distributions correspond to models with optimal C and ρ and all values of mean delay yielding a stable system.Brackets represent a significant difference between models and empirical data from a Mann-Whitney U-test (* p<0.05; ** p<0.01; *** p<0.001).All pvalues were FDR-corrected.(B) Relationship between complexity and metastability in the model with homeostasis of excitation, inhibition, and intrinsic excitability (µ E and σ E ).Gray dots show the results of all simulations while red crosses represent simulations at the optimal point (C = 3.59, ρ = 0.12).The horizontal and vertical bars and shaded areas display the mean and standard deviation of complexity and metastability, respectively, in empirical BOLD signals.(C) Comparison between FC complexity in models with the original and shuffled structural connectomes and homeostasis of G E , c EI , µ E , and σ E .Simulations were performed with C = 3.59 and ρ = 0.12 and all mean delays yielding stable solutions.Asterisks represent the significance of a Mann-Whitney U-test between samples (*** p<0.001).(D) Same as C, for metastability.

Fig. 5 .
Fig. 5. Changes in network features with the variance of noise.Values relate to models with C = 3.59 and ρ = 0.12 and all mean delays between 0 and 40 ms corresponding to stable network solutions.Solid lines represent the average and shaded areas the standard deviation across mean delays.

Fig. 6 .
Fig. 6.Disruption and recovery of network properties after focal lesion in models with homeostasis of excitation (G E ), inhibition (c EI ) and intrinsic excitability (µ E + σ E ) (A) Distance from baseline of FC matrices in acute (post-lesion) and chronic (post-lesion with homeostatic recovery of E-I balance) simulations (B) Change from baseline (in percentage) of FC-SC correlation in the acute and chronic simulations (C) Same as B, for FC Modularity (D) Same as B, for FC complexity (E) Same as B, for synchrony (D) Same as B, for metastability (F) Distributions of FCD correlations in all stages of our lesion simulations.In addition, we present the Kolmogorov-Smirnov distance between distributions at baseline and acute/chronic simulations.In all plots, p-values show the significance of the Wilcoxon Ranked-sum test (* p<0.05; ** p<0.01; *** p<0.001).All p-values were FDR-corrected.
Optimization.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 minutes.To evaluate model performance against empirical data, we use the following properties of FC.
Excitatory Populations.Homeostatic plasticity of intrinsic excitability, which modulates the parameters of the activation function F E (x), can be implemented in two manners, considering empirical data(1).The first is at the level of the firing threshold of excitatory populations, µ E :

2 of 35 Francisco
Páscoa dos Santos and Paul FMJ Verschure Computation of Homeostatic Parameters as a Function of the Target Firing Rate ρ.When the fixed point of the Wilson-Cowan model ) (Fig S15).More specifically, when C is increased, amplifying not only the magnitude of external inputs, but also their fluctuations (Fig S16), nodes can compensate for the stronger input fluctuations by being farther away from the bifurcation.Conversely, with the homeostasis of c EI or µ E , nodes with higher inputs are more likely to engage in sustained oscillations (Fig S15), enhancing network instabilities.In short, the bifurcation modulation effect afforded by the homeostasis of excitation and the slope of intrinsic excitability is an important feature of local cortical circuits, ensuring the stability of large-scale dynamics.Francisco Páscoa dos Santos and Paul FMJ Verschure

Fig. S1 .
Fig. S1.Model stability under different combinations of parameters for G E Homeostasis Red represents unstable models, where the mean activity of at least one of the nodes differed from the target firing rate ρ by more than 1%.Conversely, blue represents models where all nodes have mean activity close to the target.

Fig. S2 .
Fig. S2.NetworkBehavior for C = 0.75, ρ = 0.14 and mean delay 10 ms On the top left, we present the stability parameter space for ρ = 0.14, where the star represents the current parameters.On the bottom left, we present the dynamical portrait for G E homeostasis at node level, with red and black dots representing network nodes in the limit cycle and stable spiral regimes, respectively.On the right, we present the activity of networks with predicted G E values, for different initial conditions, and of networks with dynamical G E homeostasis with different time constants.Black lines represent r E of each node, red lines the average across nodes and blue lines the target firing rate ρ.

Fig. S5 .
Fig. S5.Network behavior for C = 4.5, ρ = 0.11 and different mean delays.On the top left, we present the stability parameter space of the model for ρ = 0.11.On the top right, we present the natural frequency of oscillation of nodes in the network as a function of the average external input they receive.On the bottom, we display network activity for models with mean delay 0, 12, 26 and 40 ms, and with different initial conditions.Black lines represent r E of each node, red lines the average across nodes and blue lines the target firing rate ρ.

Fig. S6 .
Fig. S6.Stability of models with different modes of homeostasis (A) c EI Homeostasis (B) µ E Homeostasis (C) µ E + σ E Homeostasis (D) G E + c EI Homeostasis (E) G E + c EI + µ E Homeostasis (F) G E + c EI + µ E + σ E Homeostasis.Blue and red colors represent stable and unstable simulations, respectively.

Fig. S7 .
Fig. S7.Behavior of the network with homeostasis of c EI .We present the activity of models with predicted c EI values and different initial conditions: r E = 0 (Top), r E = ρ (Middle) and r E = 2ρ (Bottom)

Fig. S11 .
Fig. S11.FC fitting of models with different modes of homeostasis(A) G E Homeostasis (B) c EI Homeostasis (C) µ E Homeostasis (D) µ E + σ E Homeostasis (E) G E + c EI Homeostasis (F) G E + c EI + µ E Homeostasis (G) G E + c EI + µ E + σ E Homeostasis.Colors represent the correlation coefficient between empirical and simulated FC matrices for each combination of C and ρ.
Fig. S13.FCD fitting of models with different modes of homeostasis(A) G E Homeostasis (B) c EI Homeostasis (C) µ E Homeostasis (D) µ E + σ E Homeostasis (E) G E + c EI Homeostasis (F) G E + c EI + µ E Homeostasis (G) G E + c EI + µ E + σ E Homeostasis.Colors represent the Kolmogorov-Smirnov distance between empirical and simulated FCD distributions for each combination of C and ρ.

Fig. S14 .
Fig. S14.Fitting performance of models with different modes of homeostasis(A) G E Homeostasis (B) c EI Homeostasis (C) µ E Homeostasis (D) µ E + σ E Homeostasis (E) G E + c EI Homeostasis (F) G E + c EI + µ E Homeostasis (G) G E + c EI + µ E + σ E Homeostasis.Colors represent the fitting score (r F C − M SE F C − KS F CD )for each combination of C and ρ.

Fig. S15 .
Fig. S15.Phase portraits of the Wilson-Cowan model under different modes of homeostasis Dashed lines represent the transition between the stable fixed point and stable spiral regimes and solid lines represent the Andronov-Hopf bifurcation between damped and sustained oscillations (i.e.limit cycle).Refer to (1) for more detail on the formulation of each mode of homeostasis and the analysis of the dynamics of Wilson-Cowan nodes under the distinct mechanisms of homeostasis.

Fig. S16 .
Fig. S16.Relationship between node degree, average inputs and input fluctuations.(A) Relationship between node degree (i.e sum of connectivity weights of a node) and the average input received by nodes.(B) Relationship between node degree (i.e sum of connectivity weights of a node) and fluctuations (i.e. standard deviation) of inputs received by a node.These results were obtained from a 60s simulation of the model with homeostasis of G E , c EI and µ E + σ E at the optimal point (C = 3.59, ρ = 0.12) and with a mean delay of 40 ms.

Fig. S17 .
Fig. S17.Synchrony of models with different modes of homeostasis(A) G E Homeostasis (B) c EI Homeostasis (C) µ E Homeostasis (D) µ E + σ E Homeostasis (E) G E + c EI Homeostasis (F) G E + c EI + µ E Homeostasis (G) G E + c EI + µ E + σ E Homeostasis.Colors represent the mean of the Kuramoto Order Parameter (i.e.synchrony) for each combination of C and ρ.

Fig. S18 .
Fig. S18.Metastability of models with different modes of homeostasis(A) G E Homeostasis (B) c EI Homeostasis (C) µ E Homeostasis (D) µ E + σ E Homeostasis (E) G E + c EI Homeostasis (F) G E + c EI + µ E Homeostasis (G) G E + c EI + µ E + σ EHomeostasis.Colors represent the standard deviation of the Kuramoto Order Parameter (i.e.metastability) for each combination of C and ρ.

Fig. S19 .
Fig. S19.Parameter spaces representing network synchrony (average Kuramoto Order Parameter across time) for each combination of C and ρ, averaged across mean delays, for models with distinct modes of homeostasis

Fig. S20 .
Fig. S20.Parameter spaces representing network metastability (standard deviation of Kuramoto Order Parameter across time) for each combination of C and ρ, averaged across mean delays, for models with distinct modes of homeostasis

Fig. S21 .
Fig. S21.Synchrony and metastability in models with distinct modes of homeostasis and empirical dataThe center of the crosses represents the mean synchrony and metastability in models at the optimal working point (C and ρ) of each mode of homeostasis.Conversely, horizontal and vertical bars represent the standard deviation of synchrony and metastability, respectively.Note that the model with all modes of homeostasis is the one that better approximates empirical levels of both synchrony and metastability.

Fig. S22 .
Fig. S22.Relationship between metastability, model performance and complexity for all modes of homeostasis (A) Relationship between metastability and the KS distance between empirical and simulated FCD distributions (B) Relationship between metastability and fitting scores (C) Relationship between metastability and functional complexity.Each data point corresponds to a single simulation.Blue lines and shaded areas represent the mean and standard deviation of metastability and complexity in empirical data.

Fig. S23 .
Fig. S23.Functional complexity of models with different modes of homeostasis(A) G E Homeostasis (B) c EI Homeostasis (C) µ E Homeostasis (D) µ E + σ E Homeostasis (E) G E + c EI Homeostasis (F) G E + c EI + µ E Homeostasis (G) G E + c EI + µ E + σ E Homeostasis.Colors represent the complexity of FC matrices for each combination of C and ρ.

Fig. S24 .
Fig. S24.Parameter spaces representing functional complexity for each combination of C and ρ, averaged across mean delays, for models with distinct modes of homeostasis

Fig. S25 .
Fig. S25.Functional complexity in models with the original and shuffled connectome.(Left) Structural connectivity matrices and respective weight distributions and complexity, corresponding to the original and shuffled connectomes.(Right) Functional connectivity matrices generated by models with the homeostasis of G E , c EI , µ E , and σ E with the original and shuffled connectomes and respective weight distributons.We used the optimal models (C = 3.59, ρ = 0.12) with a mean delay of 40 ms.

Fig. S26 .
Fig. S26.Parameter space of model with all modes of homeostasis with different noise levels.We present the results for all combinations of mean delay and noise variance for models with C = 3.59 and ρ = 0.12 and across all mean delays yielding a stable solution of the model.Note that the empty areas in the parameter space relate to mean delays for which |⟨r E ⟩ − ρ|/ρ > 0.01 in at least one node.

Fig. S27 .
Fig. S27.Disruption and recovery of network properties after focal lesion in models with homeostasis of excitation (G E ) (A) Distance from baseline of FC matrices in acute and chronic simulations (B) Change from baseline (in percentage) of FC-SC correlation in the acute and chronic simulations (C) Same as B, for FC Modularity (D) Same as B, for FC complexity (E) Same as B, for synchrony (D) Same as B, for metastability (F) Distributions of FCD correlations in all stages of our lesion simulations.In addition, we present the Kolmogorov-Smirnov distance between distributions at baseline and acute/chronic simulations.P-values represent the result of a Wilcoxon Ranked-sum test.For 1-sample tests, asterisks represent the significance of a Wilcoxon Ranked sum test (* p<0.05; ** p<0.01; *** p<0.001).All p-values were FDR-corrected.

Fig. S28 .
Fig. S28.Disruption and recovery of network properties after focal lesion in models with homeostasis of excitation (G E ) and inhibition (c EI ) (A) Distance from baseline of FC matrices in acute and chronic simulations (B) Change from baseline (in percentage) of FC-SC correlation in the acute and chronic simulations (C) Same as B, for FC Modularity (D) Same as B, for FC complexity (E) Same as B, for synchrony (D) Same as B, for metastability (F) Distributions of FCD correlations in all stages of our lesion simulations.In addition, we present the Kolmogorov-Smirnov distance between distributions at baseline and acute/chronic simulations.P-values represent the result of a Wilcoxon Ranked-sum test.For 1-sample tests, asterisks represent the significance of a Wilcoxon Ranked sum test (* p<0.05; ** p<0.01; *** p<0.001).All p-values were FDR-corrected.

Table 1 . Optimal parameters, fitting scores, dynamics and complexity for each mode of E-I homeostasis in comparison to empirical data Fitting
scores correspond to r F C − M SE F C − KS F CD , where r F C is the Pearson's correlation between matrices, M SE F C their mean squared error and KS F CD the Kolmogorov-Smirnov distance between FCD distributions.Values are presented as mean±sd.Values in bold represent a

Fig. 2. Cross-feature model performance for different modes of homeostasis
(1) results suggest that, across the various mechanisms of homeostasis, the main parameters shaping network stability are the global coupling C and target firing rate ρ.However, stability is reduced in models with homeostasis of inhibition (c EI ) (19.26% stable simulations), firing threshold (µ E ) (14.80%), or the combined homeostasis of G E , c EI and µ E (20.47%), all of which strongly rely on the modulation of "additive/subtractive" model parameters(1).Conversely, homeostasis mechanisms involving excitatory synapses (G E ) or the slope of neural excitability (σ E ) (i.e."multiplicative parameters") enhance network stability (Figure