# The impact of homeostatic inhibitory plasticity in a generative biophysical model

Iván Mindlin, Carlos Coronel-Oliveros, Jacobo D Sitt, Rodrigo Cofre, Andrea

Luppi, Thomas Andrillon, Yonatan Sanz Perl, Rubén Herzog

► To cite this version:

Iván Mindlin, Carlos Coronel-Oliveros, Jacobo D Sitt, Rodrigo Cofre, Andrea Luppi, et al. The impact of homeostatic inhibitory plasticity in a generative biophysical model. 2026. hal-05473610

HAL Id: hal-05473610

https://inria.hal.science/hal-05473610v1

Preprint submitted on 23 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
1

1 **Title: The impact of homeostatic inhibitory plasticity in a generative biophysical**

2 **model**

3

4 **Authors:** Iván Mindlin¹*, Carlos Coronel-Oliveros²,³,⁴, Jacobo D. Sitt¹, Rodrigo Cofré⁵,

5 Andrea Luppi⁶,⁷,⁸,⁹, Thomas Andrillon¹, Yonatan Sanz Perl⁹,¹⁰,¹¹,¹³, Rubén Herzog¹²,²,¹⁴*

6

7 **Affiliations**

8 1. Sorbonne Université, Institut du Cerveau - Paris Brain Institute - ICM, Inserm, CNRS, 75013, Paris, France

10 2. Latin American Brain Health Institute (BrainLat), Universidad Adolfo Ibañez, Santiago de Chile, Chile.

12 3. Trinity College Dublin, The University of Dublin, Dublin, Ireland.

13 4. Global Brain Health Institute (GBHI), Trinity College Dublin, Dublin, Ireland.

14 5. Université Côte d'Azur, INRIA CRONOS Team, Sophia Antipolis, France.

15 6. Department of Psychiatry, University of Oxford, Oxford, UK

16 7. St John's College, University of Cambridge, Cambridge, UK

17 8. Montreal Neurological Institute, McGill University, Montreal, CA

18 9. International Centre for Flourishing, Universities of Oxford (UK), Aarhus (Denmark) and Pompeu Fabra (Spain)

20 10. National Scientific and Technical Research Council (CONICET), CABA, Buenos Aires, Argentina

22 11. Department of Engineering, Universidad de San Andrés, Buenos Aires, Argentina

23 12. Instituto de Física Interdisciplinar y Sistemas Complejos (IFISC, UIB-CSIC), Campus UIB, Palma de Mallorca, Spain.

25 13. Center for Brain and Cognition, Computational Neuroscience Group, Universitat Pompeu Fabra, Spain

27 14. Department of Psychology, University of the Balearic Islands, Palma de Mallorca, Spain

29

30

31 * Corresponding authors: rherzoga(at)gmail(dot)com, mindlinster(at)gmail(dot)com

32

33 **Abstract**

34 A main characteristic of biological systems is their capacity to dynamically adapt to environmental
35 changes. In the brain, synaptic plasticity enables the strengthening or weakening of connections
36 between neurons, allowing neural circuits to adapt based on experience, learning, and environmental
37 changes. Yet, it is homeostatically regulated such that it avoids excessive proliferation of synaptic
38 contacts. These mechanisms can be studied with large-scale models of brain activity. Here, we embed
39 a biologically grounded inhibitory-homeostatic plasticity rule into the Dynamic Mean Field (DMF)
40 model, creating a Homeostatic Dynamic Mean Field (HDMF) model that dynamically tunes local
41 excitation-inhibition balance. Convergence of excitatory firing rates is reached by mapping a large
42 range of coupling strength to parameters of inhibitory synapses. The HDMF reproduces statistical
43 observables of brain activity as well as the original DMF, and can sustain neuromodulatory
44 perturbations without overhead computations. The HDMF can generate unprecedented sleep-like
45 slow-wave activity, which can also coexist with wake-like asynchronous dynamics, permitting to
46 model dissociated states of consciousness such as parasomnias. Together, these results show that a
47 single homeostatic rule broadens the stability and expressiveness of the DMF, providing a unified
48 platform for studying how local adaptive processes shape the diverse global dynamics of the human
49 brain.
50

51

Keywords

52

whole-brain model; dynamic mean field: slow-wave; homeostasis; plasticity;

53

neuromodulation

54

55

Introduction

56

57

Distinct brain dynamics underlie different states of arousal¹, consciousness²,³, cognition⁴, and health

58

conditions⁵. This varied landscape is grounded in the brain's remarkable capacity to self-regulate the

59

underlying physiological processes that induce them⁶. While some variables remain tightly controlled

60

leading to stable functioning, others can shift over time, giving rise to the wide array of observable

61

brain states⁷. At the microscopic level, plasticity serves as an adaptive mechanism, continuously fine-

62

tuning neural activity in response to a dynamic environment. However, synaptic changes are tightly

63

regulated by homeostasis, which constrains neuronal firing rates within functional ranges. This

64

regulation prevents pathological excitation-inhibition imbalances⁸, such as hyperactivity or

65

quiescence, thereby preserving the brain's capacity to sustain diverse and flexible dynamics

67

Biophysically detailed neuronal models are crucial for elucidating the mechanistic principles of neural

68

function, enabling researchers to simulate and analyze how specific dynamics and synaptic

69

interactions contribute to neuronal signaling and information processing. At a macroscopic level,

70

many models have been used to explain the wide variety of global brain dynamics observed

71

empirically⁹⁻¹³. Whole-brain models are pivotal in neuroscience as they enable the simulation of large-

72

scale neural dynamics by integrating anatomical connectivity data with computational frameworks,

73

providing insights into the mechanisms underlying brain functions. These models are instrumental in

74

understanding how distributed brain networks interact to produce complex behaviors and cognitive

75

processes. Models that incorporate plasticity mechanisms (in the form of local adaptation rules) can

76

support a remarkably broad spectrum of neural dynamics using a single unifying mechanism¹⁴⁻¹⁸. For

77

example, by strategically tuning these adaptive processes—such as homeostatic plasticity of

78

excitatory/inhibitory connections or intrinsic neuronal excitability—they can reproduce both wake-

79

like asynchronous firing and sleep-like slow-wave activity in a cohesive framework¹⁶. This way,

80

complex shifts in global brain state can emerge naturally from local plastic changes, allowing

81

networks to transition smoothly between different regimes while retaining stable yet flexible

82

functional organization¹⁵,¹⁷. This built-in adaptability may be crucial not only to describe healthy

83

conditions, like wakefulness and sleep, but also for investigating how pathological¹⁹,²⁰ or

84

pharmacologically altered states¹¹ arise when these mechanisms are disturbed.

85

86

While most such frameworks do not explicitly include homeostatic adaptation rules, they can still reproduce many different brain states by introducing targeted mechanisms for each phenomenon²¹.

87

88

The Dynamic Mean Field (DMF) model ²²⁻²⁴ is a well-established example: through careful tuning of

89

local excitability parameters and global coupling—often guided by empirical maps of

90

neurotransmitter expression or other empirical brain heterogeneities¹¹,²¹,²⁵⁻²⁷—it has replicated the

91

effects of pharmacological manipulations²⁸,¹¹ on brain-wide activity. Similarly, gradual changes in

92

structural connectivity weights have been used to simulate aging²⁹, by modeling the progressive

93

decline of synaptic or axonal integrity. Schizophrenia-like states have been reproduced via increasing

94

local self-excitation in specific subnetworks³⁰, while neurodegenerative disorders have been modeled

95

through excitation/inhibition imbalances in regions undergoing atrophy³¹ with specific

96

hypo/hyperexcitability spatial patterns. Collectively, these parameter-based approaches highlight how

97

a single model architecture can account for diverse observed states²¹, albeit typically without relying
on a unified adaptation mechanism to maintain stable firing across conditions. However, this requires exhaustive parameter searches²³ or restrictive heuristic solutions²⁴ to set the local inhibition, also known as Feedback Inhibitory Control (FIC). Additionally, the solutions obtained are static, which can constrain free parameters and limit the model's capacity to capture a broader range of observed brain dynamics.

As presented thus far, whole-brain models use individual mechanisms to control the dynamics that characterize a specific brain state. It is unclear how this approach can tackle conditions that are characterized by concurrent and distinct dynamical regimes. Dissociative states fall into this category, where asynchronous activity akin to wakefulness co-exists with slow waves typically present in sleep³²,³³. This paradoxical appearance of diverse dynamics bears resemblance to the notion of ‘chimeras’ from complexity science: when a subset in a system of identical parts, sustains disparate dynamics than the rest. In fact, chimeric states of brain dynamics are linked to cognitive functions³⁴. Although the brain is not composed of identical parts since it features spatial gradients³⁵ and use-dependent homeostatic regulation can alter dynamics in a regional way³⁶, we take inspiration from complexity science to call these dissociative states ‘chimeric-like’ and investigate how regional modulations of homeostatic processes underlie the emergence of concurrent distinct dynamical regimes.

Here, we embed a biologically grounded homeostatic inhibitory plasticity³⁷ into the DMF framework, creating the Homeostatic DMF (HDMF). We then systematically chart the model parameters' landscape to pinpoint an optimal solution that ensures homeostasis at the whole-brain level. Next, we test its adaptiveness and robustness when applying biologically plausible neuromodulation, revealing enhanced stability and predictability of the change on firing distributions. We then demonstrate that just by manipulating the global coupling the model exhibits sleep like slow-wave activity (SWA), an unprecedented feature for this model. Furthermore, we achieve chimeric-like states by introducing spatially heterogeneous homeostasis. Finally, we show that the HDMF enables the joint reproduction of the brain's functional connectivity and its dynamics in a subspace of parameters where both wake and sleep-like activity can be generated. This demonstrates how a unified adaptation mechanism can reconcile diverse brain states in a single, cohesive framework.

## Results

We extended the DMF model, a whole-brain approximation introduced by Deco et al²³ by incorporating a dynamic feedback inhibition control (FIC) mechanism with homeostatic plasticity, creating the HDMF model. The DMF model simulates neural activity by defining each region as consisting of interconnected excitatory (E) and inhibitory (I) neuron pools, with long-range interactions mediated solely between E pools via a structural connectivity matrix $C_{ij}$ (here we used AAL cortical atlas Figure 1A, see Methods). A global coupling parameter $G$ scales the influence of network interactions on each region. $G$ is a free parameter used to fit empirical data by simulating blood-oxygenation-level-dependent (BOLD)-like signals from the excitatory firing rates time series $r^E$ (Figure 1B). In its original version²³, the DMF used a static FIC that required time-consuming optimizations to achieve stable firing rates, but a heuristic structure-dependent linear solution was later proposed bypassing this computational bottleneck²⁴. This latter solution expresses the FIC as a function of the structural connectivity strength, $G$ and another free parameter $\alpha$, which represents the global E/I balance (see Methods). However, both approaches are static in nature (i.e. implemented as constant scalar values), imposing strong limitations in terms of functional repertoire and firing rates stability especially in the high $G$ (i.e. high coupled) regime. Our alternative dynamic solution for the
146 FIC values regulates the excitatory firing rates around a target firing rate ρ based on a learning rate (LR) and a decay ( decay) term (Figure 1A, B).

148

149

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

150

151 **Figure 1. Overview of the Homeostatic Dynamic Mean Field (HDMF).**

152 A) Each brain region (circles) is defined by a parcellation and their connections by a structural connectivity matrix (edges) reconstructed from in vivo diffusion tractography. The dynamics of each region consists of reciprocally connected excitatory (E, yellow) and inhibitory (I, blue) neural populations, while inter-regional connections are only through the E population.

156 B) The I->E connection, i.e the FIC variable $J_i$ is controlled by a plasticity mechanism which dynamically adjusts the E firing rate ($r^E$) to approach the target firing rate $\rho$ depending on LR and $\tau_{decay}$.

158 C) Example of time series of the simulated whole-brain average firing rates of an E (yellow) and I population (blue), and the respective average dynamic FIC (green) across all brain regions. Note that the FIC rapidly increases after sudden increases of the E firing rates, compensating for the excess firing rate, and thus reducing them. Note the delayed coupling between the E firing rates and the FIC.

162

163 **Homeostatic dynamic mean field model (HDMF)**

164

165 Here, we aimed to identify the regime in which an inhibitory plasticity rule enables effective homeostatic control within the DMF model, which we refer to as the HDMF. To this end, we computed the homeostatic fit, measured as the percentage mismatch between the average excitatory firing rate and the target firing rate ($\rho = 3.44 Hz^{23}$), as a function of the learning rate (LR = [10⁰, 10³] ms⁻¹) and the decay ( decay = [10², 10⁵] ms) (Figure 2A). To demonstrate the enhanced stability, the average mismatch was computed across a wide range of G values, from 0 to 8.5. The results, plotted in log-log scale, revealed a power-law relationship between the learning rate and decay that optimizes homeostasis, identifying the optimal parameters for the HDMF. Further, we found that the region of parameters for optimal homeostasis also corresponds to minimal variability in firing rates across the range of G, indicating enhanced stability regardless of the G value (Figure 2B) and target firing rates (Supplementary Figure 1). Note that without the decay term the average firing rate never converged to the target firing rate (Supplementary Figure 2). Then, in the optimal regime, the homeostatic plasticity rule can be expressed only in terms of the LR and $\rho$, as $\tau_{decay}$ can be derived
178 from the power-law fit. In biophysical terms, this relationship implies that the faster the synaptic 179 weight is learnt (high LR), the faster it needs to be forgotten (fast decay; small $\text{decay}$) in order to adapt 180 to the ongoing fluctuations of network inputs. From here on, whenever we use the HDMF, it means 181 that the $\text{decay}$ is derived for fix LR and $\rho$ using the power-law relation from the inset in Figure 182 2A.

183

184 To investigate the nature of the FIC values obtained through plasticity, we compared the average 185 firing rates as a function of G for three different models. First, the original DMF model that uses a 186 static FIC for each region, derived from a linear solution<sup>24</sup> or using pseudo gradient-descent 187 algorithms<sup>23</sup>. Second, the HDMF that dynamically changes FIC values to sustain a constant average 188 firing rate. Notably, one might expect that, when averaged over time, the time-varying FIC of the 189 HDMF could converge toward FIC values similar to those obtained through static methods. Thus we 190 also created a third model, the “mixed” DMF, that uses the FIC time-average values found from the 191 HDMF but does not dynamically adjust these values. For the same range of G values, we computed 192 the correlation and mean squared error (MSE) between the static FIC values obtained from the linear 193 solution and the time-averaged dynamic FIC values (Figure 2C). At very low G values— 194 corresponding to weak network input—the FIC dynamics were predominantly noise-driven, reducing 195 the correlation between the static and dynamic FIC values in this regime. For the remaining range we 196 found that while the correlation remained maximal, the MSE markedly increased near G=2.5, 197 indicating a difference between static and dynamic FIC values in this regime. Consistently with 198 previous results<sup>23,24</sup>, the DMF diverged into a high-firing rate regime near G=2.5 whereas the HDMF 199 maintained firing rates close to the target across the entire range of G values (Figure 2D). 200 Additionally, despite using FIC values derived from the HDMF, the mixed DMF model failed to 201 achieve stability near this critical point and remained well off the target rate -40% mismatch- for 202 higher G values. This highlights the necessity of dynamic FIC to compensate for rapid fluctuations in 203 firing rates.

204 Overall, these results demonstrated that the dynamics of the FIC, beyond the simple average, were 205 necessary for achieving stability in a wide range of G. In fact, in the low-to-mid coupling regime 206 (0.5<G<2.5) the heuristic linear solution for the static FIC is a particular case of the dynamic FIC, 207 suggesting that the dynamics of the latter could be represented by its steady-state value within this 208 regime. Beyond this regime, this relationship does not hold anymore.

209
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Figure 2. The local inhibitory plasticity controls whole-brain homeostasis.</figcaption></figure>

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

A) Average homeostatic fit (the percentage of mismatch between the average E firing rate and $\rho$) as a function of the learning rate and the decay for 4 different $\rho$ in log-log scale. Averages were taken across a range of G values from 0 to 8.5. Green dots show the minima for each row and the darker lines are the best linear fit, representing the optimal relationship between learning rate and decay to achieve homeostasis. Note that since the plot is in log-log scale, this relationship is a power-law. The coefficients of the linear equation inset in the panel are used to implement the HDMF at the desired target firing rate.

B) Same as in A, but showing the standard deviation in the same range of G values. This shows that the region of minimal mismatch is also the region with more stability across G values.

C) Correlation and the mean squared error (MSE) between the static FIC obtained by the linear solution (see Methods) and the one obtained after averaging the dynamic FIC over time, as a function of G. Note also that close to G=2.5 the MSE diverges, meaning that beyond this point the linear solution for the FIC is no longer the time-average of the dynamic FIC. Despite their high correlation, a static FIC is not able to achieve stability, demonstrating the relevant role for the dynamic FIC in achieving firing rate stability.

D) Percentage of firing rate mismatch as a function of G for the DMF (Static, green dashed line), for the HDMF (Dynamic, yellow line) and for a DMF where the FIC values corresponded to the average of the dynamic FIC (Mixed, black dots). Note that close to G=2.5 the DMF diverges into a high-firing rate regime, while the HDMF maintains the activity close to the target. The mixed model could not sustain stability near this divergence point. For higher G values, the mixed model stabilized but below the target rate, indicating over inhibition. For panels C and D, the DMF used $\alpha$=0.75 while the HDMF used LR=3.5. Both were simulated for 50000 time points using the same set of 30 different initial conditions. Lines and shadowed areas represent the mean and standard deviation across initial conditions, respectively.
235 **Enhanced stability under neuromodulation**

236 We investigated the stability of the HDMF model under neuromodulation, using the receptor map of
237 the 5HT-2A receptor obtained from positron emission tomography (PET) images<sup>38</sup>, which is widely
238 expressed in the cortex. Neuromodulation was simulated via response gain modulation<sup>27</sup>, where an
239 increase in the neuromodulatory gain elevates the slope of the F-I curve in both E and I pools<sup>11</sup>,
240 resulting in higher firing rates with the same input current (Figure 3A). As each region will have a
241 specific value of receptor expression, this mechanism will impact regions differently. A major
242 obstacle of this approach in the DMF is that it largely impacts the firing rates, thus losing biological
243 plausibility for high neuromodulatory gains.

244

245 Here, we compared the firing rate stability of the DMF model with the HDMF model in response to
246 increasing neuromodulatory gains. In the DMF model, firing rates increased with gain, but at high
247 gain values, the positive correlation with receptor density was highly variable, showing a wide spread
248 around the regression line (Figure 3B). In contrast, the HDMF model maintained firing rates closer to
249 the target (3.44 Hz), even with large neuromodulatory gains. Furthermore, the changes in firing rates
250 in the HDMF model exhibited a tight negative dependence on receptor density, demonstrating a more
251 predictable and controlled response (Figure 3C). Notably, this negative correlation is the result of the
252 homeostatic process: as more gain implies amplified responses for the same input currents, the
253 departure from the target firing rate will also amplify the dynamic FIC response, driving the
254 excitatory firing rate down on average in consequence. However, this offset was minor, in comparison
255 to the explosive effect observed with the DMF, and still provides an enhanced range for simulating
256 neuromodulation with biologically plausible firing rates. This stability was not limited to a particular
257 gain value or neurotransmitter distribution (Supplementary Figure 3), underscoring that the HDMF
258 supports a broader and more robust neuromodulatory regime than the DMF. By contrast, a mixed
259 model (as in the previous section) used as a control yielded inconsistent outcomes (Supplementary
260 Figure 3), supporting again that the dynamics of the FIC provide a stability that its steady-state value
261 only can't achieve.

262
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Figure 3. Homeostasis stabilizes the effect of increasing neuromodulation.</figcaption></figure>

A) Schematic representation of neuromodulation through response gain modulation in the simulation. As the neuromodulatory gain increases, the slope of the F-I curve increases, evoking higher firing rates with the same input current.

B) Average firing rate vs the receptor density for the DMF. Each dot represents a brain region. Colors denote 5 different neuromodulatory gains. Note that the firing rates are blown up with the gain and that the positive correlation between both variables is weak and noisy at higher gains.

C) Same as B, but for the HDMF. Note that despite the large increases in the neuromodulatory gain, the firing rates remain closer to the target firing rate.

Both models were simulated for 100 different initial conditions. DMF used α=0.75 and G=2.11; HDMF used LR=3.5 and G=2.55.

A potential mechanism for whole-brain slow wave generation

Recent advances in whole-brain modeling aim to reproduce different brain states through single parameter manipulations thus providing insightful and interpretable mechanisms about the transitions
281 of brain dynamics²¹. Slow waves, whose underlying bistable dynamics alternates periods of neural
282 activity and silence ³⁹, are central for diverse biological and cognitive functions ⁴⁰,⁴¹ and are present in
283 states such as deep sleep, general anesthesia and disorders of consciousness ³⁹. Previous attempts to
284 generate slow wave activity (SWA) in whole-brain models have employed different mechanisms ¹⁶,⁴²,⁴³
285 but it has never been done in the DMF, although this model has been widely used to study different
286 states of consciousness.

287

288 Here, we aimed to demonstrate that the HDMF model can generate SWA through modulation of the
289 learning rate LR parameter. Given a high global coupling G, as LR increased, the whole-brain average
290 of the excitatory firing rates transitioned from fully regular SWA to irregular activity. (Figure 4A).
291 This transition was further captured by the emergence of pronounced peaks in the power spectra at
292 frequencies consistent with empirical delta slow waves (1-2 Hz), confirming the presence of periodic
293 activity, opposed to the flat power across frequencies characteristic of irregular regimes (Figure 4B).

294

295 We then systematically explored the power spectra across a wide grid of LR and G values, identifying
296 for each parameter combination the average dominating frequency across all regions and their
297 corresponding power. This analysis revealed a broad landscape of delta band frequency (Figure 4C)
298 with high power (Figure 4D). As we increase the LR, frequency increases and power decreases, as a
299 consequence of faster learning dynamics. This reaches a critical point where variability peaks and the
300 system enters a general high frequency with lower power dominated state (Figure 4E-F). These
301 findings reveal a potential whole-brain mechanism of slow wave generation within the HDMF model.
302 The homeostatic plasticity is at the core of this mechanics: for the high-coupled regime, the excitatory
303 and synchronized inputs from the network abruptly increase a given region's firing rate, with the
304 consequence of an abrupt rise on the FIC. As a consequence, the firing rates decrease rapidly as well
305 and remain down until the FIC decays enough to trigger another active state. When LR is too high, the
306 compensation is much faster, generating more noise driven irregular dynamics with less dominance of
307 slow frequencies. The capacity to induce SWA through a biophysically grounded homeostatic rule
308 highlights the HDMF model's versatility and interpretability in capturing diverse brain states.
<figure><img src="image_11.png" /><figcaption>Figure 4. Slow waves emerge in the HDMF in the highly coupled regime.</figcaption></figure>

A) Whole-brain average excitatory firing rate for three different learning rate values (LR), with fixed high global coupling (G=3.5).

B) Power spectrum of the time series from A. Power around 0.5Hz increases with G, reaching its maximum value with the highest G. At high values of LR low frequency power is comparatively neglectable.

C) Average maximal frequency across regions as a function of LR and G. Colored squares show areas where points are taken for Fig 5A.

D) Average power of the maximal frequencies shown in C. Values in C and D are averages out of 8 different initial conditions.

E and F) Standard deviations for (C) and (D), respectively, across initial conditions. Note the increased variability near the transition between the slow-wave and irregular activity regimes.
322

323

Heterogeneous plasticity enables different and coexistent dynamics

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

Until now we explored the new dynamical repertoire of the DMF model when incorporating a homeostatic plasticity rule, which facilitates the emergence of SWA. However, this plasticity was uniformly applied across all regions, leading to fully synchronized dynamics. In contrast, biological plasticity is inherently heterogeneous, varying across the brain<sup>35,44</sup> –giving rise to distinct roles and computational capacities. Since the HDMF can self-regulate, we asked whether it could simultaneously sustain different dynamical regimes via heterogeneous plasticity. Motivated by the aforementioned chimeric-like states, we studied what happened to a model set in a system-wide asynchronous regime, when we change adaptiveness through the LR of some regions to induce in them a regional SWA regime, and, conversely, starting from system-wide SWA, shifting regions toward an asynchronous regime.

We created baseline simulations at fixed high coupling (G=3.55): a uniform low LR (LR=100; SWA regime) state and a uniform high LR (LR=1,500; asynchronous regime) (Figure 5A). To select which regions to change their LR value with respect to the starting global uniform LR value we used two criteria based on anatomical connectivity, namely node strength (STR) and s-core (S-CORE)<sup>45,46</sup>. Node strength is a local graph theoretical measure that corresponds to the sum of all the connection weights for a given region. S-core is a mesoscopic measure aiming to detect the maximal subgraph in which every node has strength -at least- *s*. When using STR we choose the *k* regions that have the highest average connectivity. For S-CORE we select the *k* regions that have the same or higher connectivity that are also connected between each other.

To characterize each region's local dynamics, we computed its peak frequency and corresponding power. Figure 5B shows these frequencies for both uniform and heterogeneous cases. For the baseline asynchronous regime, regions have high frequency (in the alpha band) and comparatively low power; in the SWA regime they exhibit low frequency (in the delta band) and comparatively high power. After heterogenizing a uniform asynchronous state (via STR and S-CORE) the modified regions shift towards lower frequency and higher power, approaching the SWA state but not reaching it. To quantify this shift we define the Regime Type Index (RTI) which measures proximity to each baseline frequency (see Methods for details). For an RTI of 1, the region's peak frequency matches the asynchronous regime, while for an RTI of -1, it matches the SWA regime. Inspired by chimera states studied in complex systems theory, we call the co-existence of regions with both high and low RTI values a chimeric-like state.

We repeated this procedure for a different number of regions going from a global high LR regime to a low LR regime and vice versa (Figure 5C-D). Triangles represent regions selected to have a different LR regime. When selecting nodes via the STR or RANDOM strategy and starting from an asynchronous regime, these could not reach a fully SWA regime. As the amount of modified regions increased, they reached a lower RTI, approaching a SWA regime. When the selection criteria was S-CORE, selected regions reached the lowest RTI, corresponding to a fully SWA regime. The difference in RTI distributions of the selected regions between both criteria held for group sizes 9 and 10 regions (*** means Cohen's d>0.8). In contrast, when starting from a uniform SW, modified regions reached the asynchronous RTI regardless of group size. Region selection criteria did not make a difference with respect to the frequency reached.
<figure><img src="image_12.png" /><figcaption>**Figure 5.** Regional heterogeneity in the learning rate (LR) gives rise to chimeric-like dynamics. **A)** Example of time series with high coupling strength (G=3.5) exhibiting SWA (LR=100) or asynchronous dynamics (LR=1500)</figcaption></figure>

**B)** Changes in region frequency and power when modifying specific region LR. Dots represent the maximum frequency and respective power of each region in a uniform SWA (yellow) or asynchronous (grey) regime. When we select a few regions (9, green triangles) to have the high value of LR, and the rest (81, green circles) at a low value of LR. Selected regions reach the asynchronous high LR regime, while SWA regions reduce their power. Regions were selected according to the STR criterion.

**C)** We expand the observation from **B** to bigger groups of selected regions. We quantify the difference in frequency from the baseline uniform models (SWA and asynchronous) to assess how much the dynamics of the model change for each region after heterogenization (see Methods for details on the Regime Type Index). For small groups of regions, the heterogeneous group frequency only reaches the frequency of a SWA uniform group if selected regions are a S-CORE group.. Trivially, this difference disappears as the selected group increases in size until it reaches SWA homogeneity.

**D)** Similar to **C**) but starting from a uniform SWA group and modifying regions to have asynchronous regimes. Nonetheless, regions reach immediately the same frequency as an uniform asynchronous regime. The S-CORE
387 and STR difference disappears. The number of regions was selected to match the sizes of the subnetworks found
388 by the S-CORE algorithm<sup>46</sup>.*** means Cohen's d>0.8

389

390 **HDMF fit to resting state fMRI data**

391 After showing the novel features of the HDMF in this section we confirm that these modifications do
392 not affect classic use of the DMF to reproduce slow fluctuating dynamics. We aimed to compare the
393 DMF model (with the heuristic linear solution for the FIC) with the HDMF model in terms of their fit
394 to empirical resting state fMRI data. To this end, we used resting state fMRI data from 13 healthy
395 subjects and computed two observables: i) the functional connectivity (FC, Figure 6C) matrices, and
396 ii) the functional connectivity dynamic (FCD) histograms (Figure 6D). These observables have been
397 typically used to fit whole-brain models to fMRI data<sup>11,24,21</sup>, as they can be reliably estimated from the
398 usual small sample sizes of fMRI recordings. We explored parameter space using a grid of 160 values
399 of the global coupling G (0.1 to 8.05 in steps of 0.05) and 150 values of either the learning rate (LR,
400 ranging from 1 to 10,000 in logarithmic space) or the excitation-inhibition balance parameter α
401 (ranging from 0.01 to 1, linearly spaced). For each parameter combination, we ran 13 simulations to
402 match the empirical sample size and computed both observables. Goodness of fit was quantified using
403 the correlation between the population-averaged empirical FC and the simulated FC, and the
404 Kolmogorov-Smirnov (KS) distance between the pooled empirical FCD histogram and the simulated
405 FCD histogram<sup>27</sup>. As shown in Supplementary Figure 4, the parameter regions yielding optimal fits
406 for FC and FCD were disjoint for both models. To identify parameter configurations that fit both
407 observables, we selected the top-performing parameter sets for FC and FCD separately and retained
408 only those belonging to the intersection of the two sets. Increasing the tolerance threshold increased
409 the size of this overlap, yielding parameter configurations that provided simultaneous, though not
410 individually optimal, fits to both observables (Figure 6A-B). We selected the top 13% (DMF) and top
411 14% (HDMF) overlap regions, which resulted in an equal number of parameter points for both models
412 (38 points each). For all parameter points within these overlap regions, we then performed 100
413 simulations and bootstrapped taking 13 samples each time to calculate the final goodness of fit. DMF
414 best fitted the FC with G=2.05 and α=0.72, FCD with G=1.95 α=0.70. HDMF best fitted the FC with
415 G=2.55 and LR=3.5, FCD with G=2.6 and LR=2.26.

416

417 Average FC correlation in the overlap area was higher for DMF than the HDMF (Pearson *r* = 0.36 ±
418 0.04, 95% CI [0.28, 0.44]; *r* = 0.33 ± 0.03, 95% CI [0.27, 0.39], respectively) with a large effect
419 (Cohen's d=-3.64) and KS distance was lower on average for the DMF compared with the HDMF
420 (KS=0.11±0.006, 95% CI [0.116,0.118] and KS=0.12±0.006, 95% CI [0.122, 0.124], respectively)
421 with a large effect (Cohen's d=0.91). Although the fit is better for the DMF we can see in both
422 observables that the difference does not come from a major change in how the model reproduces these
423 observables (Figure 6C), but rather a numerical difference. Different optimization strategies and
424 datasets may yield different performances. What we show here is that without distorting the purpose
425 of the DMF, the HDMF can additionally incorporate new behaviour, i.e. SWA, that was not possible
426 before. This is remarked in Figure 6D, where the overlap between FC and FCD fit (threshold >0.25)
427 colocalizes with a region in the parameter space where sleep-like SWA and wake-like irregular
428 activity emerge.

429

430
<figure><img src="image_13.png" /><figcaption>Figure 6. Slow wave activity falls within the biologically plausible regimes of the HDMF.</figcaption></figure>

A) Maximum FC correlation obtained within the overlap region as a function of the top X% tolerance threshold for FC and FCD fits, shown separately for the DMF (green) and the HDMF (orange). Stars indicate the tolerance thresholds selected for further analysis (13% for DMF, 14% for HDMF).

B) Minimum FCD KS distance obtained within the overlap region as a function of the top X% tolerance threshold for both models. Stars denote the selected tolerance thresholds.

C) Top: empirical FC matrix and simulated FC matrices from the DMF and HDMF at representative parameter points within the overlap region. Bottom: pooled histograms of FCD values for empirical data, DMF, and HDMF.

D) Parameter-space representation of the HDMF showing mean dominant frequency as a function of global coupling GGG and learning rate (LR). Contours indicate overlap regions corresponding to different tolerance thresholds (14%, 20%, 25%, and 30%). The star marks the parameter point selected within the overlap region.

## Discussion

### Overview of Main Contributions

One of the brain's astonishing features is the ability to switch between radically different cognitive states while maintaining baseline functions. This robustness is partly grounded in its capacity to adapt to external and internal fluctuations. In this study, we introduced a homeostatic plasticity rule into the DMF model to enhance its biological plausibility and prove how this can support a broader repertoire of dynamics. This self-tuned regulation led to a robust maintenance of firing rates in different scenarios such as high coupling or following changes after neuromodulatory drive. In addition, the
455 HDMF generated slow waves—a phenomenon not previously captured by the DMF—and supported
456 the coexistence of multiple dynamical states within the same system, thereby enabling the modeling
457 of dissociative conditions such as parasomnias³².

458 **Relevance of Homeostatic Plasticity for Whole-Brain Modeling**

459 We first tuned the plasticity rule to self-regulate (i.e., reach homeostasis with respect to the model's
460 target firing rates). We achieved this by finding a power-law relationship between the two control
461 parameters, such that simulated firing rates reached the target rate with a neglectable deviance. The
462 central premise of homeostatic plasticity is to regulate local or global excitability, preventing runaway
463 dynamics or quiescence⁴⁷,⁴⁸. Integrating such a mechanism into the DMF —a well established
464 biophysical whole-brain model— offers a biophysically grounded means to self-tune neuronal activity,
465 thereby reducing sensitivity to precise parameter configurations. Previous solutions for calculating
466 FIC include linear heuristics²⁴ and pseudo-gradient optimization²³. The former relies on a relationship
467 between structural connectivity and global coupling—providing mechanistic insight into E/I balance
468 that is useful for spontaneous activity but falls short under extrinsic perturbations such as
469 neuromodulation or external stimuli. In contrast, region-specific optimization overcomes this
470 limitation by using short simulations to adjust FIC values until stability is reached, thereby adapting
471 weights to any modeling condition. However, this approach incurs high computational cost and
472 compromises biological realism, as the brain can't compute its synaptic weights in advance. Within
473 the typical range of G values used in previous studies (0-2.5), the average FIC values produced by our
474 rule match those obtained with previous methods. This highly effective self-regulation may not be
475 present in pathological states. This opens an intriguing avenue for future exploration—one that first
476 requires us to lay the groundwork. Accordingly, we have integrated homeostatic regulation into the
477 fastDMF framework²⁴, in keeping with its mission to simplify model development, ensure
478 accessibility and move toward canonical models whose empirical validation can be more fully
479 explored.

480 By adding the homeostatic plasticity rule, the HDMF generated slow waves. In this regime the system
481 exhibits bistable dynamics, rapidly alternating between states of synchronous activity (UP) and
482 synchronous silence (DOWN). When global coupling was high, slow adaptiveness gave rise to slow
483 waves, whereas fast adaptiveness generated noisy activity. This coupling-adaptation mechanism is
484 compatible with the low-level mechanisms driven by activity-dependent potassium currents⁴⁹.
485 Additionally other modelling studies propose this generates as a bistable state between two
486 attractors⁵⁰. In vivo measurements observe the UP-DOWN states by building accumulative excitatory
487 recurrences that are controlled through a downscaling process of synaptic weight⁵¹. Transitions in the
488 sleep-cycle have also been modelled through the interactions of neurotransmitters that modulated
489 hyperpolarization and depolarization of each neuron⁵². The cumulative control of neuromodulation
490 between these mechanisms yields the overall effect of a self-sustaining UP-DOWN regime, or an
491 asynchronous noisy regime. In this work, the specific pathway by which this is achieved is abstracted
492 to the global adaptiveness that is tuned to the state of each region. What is crucial is that the capacity
493 to adapt to increased excitatory inputs allows the model to enter this regime that is prohibited to the
494 DMF model.

495

496

497 **Coexisting Dynamical States and Chimeric-like Activity**

498 The HDMF generated distinct dynamical regimes simultaneously, including UP-DOWN states and
499 asynchronous activity. This feature resonates with empirical observations of local sleep³⁶, partial
500 awakenings, and dissociative conditions³², where different brain areas can manifest oscillatory sleep-
501 like states that concur with asynchronous activity¹⁷. When this happens in the HDMF we refer to it as
502 a "chimeric-like" state. In complex systems, a "chimeric system" exhibits different dynamics across
503 identical regions⁵³. Although this is not the case for our model –regions have different connectivity
504 patterns and different learning rates–, it is neither the case for the brain. More broadly, several
505 modeling approaches highlight that embedding empirical heterogeneity is a fruitful direction for
506 capturing realistic dynamics²⁵,⁵⁴. In the HDMF heterogeneity in adaptation gave rise to different
507 dynamics that not only co-existed but influenced each other, and interacted effectively via structural
508 coupling. This capacity could help illuminate the role of SWA within asynchronous regimes, the use-
509 dependent regional changes in homeostatic regimes³⁶ and potentially lesioned areas of the brain.
510 Additionally it may explain how certain states (e.g., lucid dreams, false awakenings) could emerge
511 through specific pathways that shift local inhibition and excitatory drive⁵⁵,⁵⁶.

512 Performance of the HDMF Relative to Prior DMF Applications

513

514 We fitted both the DMF and HDMF to resting-state fMRI data. Different observables were optimally
515 fitted in distinct regions of the explored parameter space. For a fair comparison, we also included
516 DMF parameter points that are typically excluded due to excessively high firing rates. Although it has
517 been shown that not all observables can be optimally fitted at a single working point in the DMF²¹,
518 we chose to compare the models in regions where both FC and FCD were reproduced at an acceptable
519 level. While the HDMF did not surpass the DMF in overall performance for the current dataset, it
520 retained the ability to reproduce slow, fluctuating spontaneous activity. This modest reduction in
521 absolute performance is accompanied by the previously described dynamical behaviors and by
522 increased robustness.
523

523

524 The HDMF remained stable under neuromodulation across a range of intensities, whereas the static
525 FIC often failed to maintain firing rates within a biologically plausible regime. This was reflected in
526 the strong correlation between regional receptor density and resulting activity in the HDMF. In
527 contrast, the DMF and the mixed model –where the FIC is derived from the average of the dynamic
528 FIC– exhibited weaker correlations and pronounced hyperexcitability. The ability to sustain coherent
529 dynamics while simulating pharmacological interventions is particularly relevant for modeling
530 approaches aimed at developing digital twins to test pharmacological strategies²⁶,⁵⁷

531

532 Comparison to Other Model Families

533 The plasticity rule we implement is biologically inspired by how inhibitory synapses regulate the
534 balance between excitation and inhibition. Other studies have also implemented it for large-scale
535 networks. Stasinski et al.⁵⁸ show that by incorporating it there is a complexity increase in local
536 activity, as well as a more constrained and accurate space for optimal fitting of BOLD data. Hellyer et
537 al.⁵⁹ did the same using the Wilson-Cowan model, showing that FIC adaptation to specific target rates
538 can yield power-law distributed avalanche sizes. This suggested that the homeostatic system was at
539 criticality, explaining its enhanced flexibility. The AdEX model⁶⁰ implements a different form of
540 adaptation and can also generate UP-DOWN states. In this work we modified the aforementioned rule
541 by adding a *decay* parameter, which was necessary to reach homeostasis (Supplementary Figure 2).
542 Furthermore, the Frequency-Current curve of the DMF generates firing rates from currents through an
543 exponential-linear function that does not plateau, i.e. can increase indefinitely, contrary to typical
544 sigmoid functions used in other models, such as the Wilson-Cowan. The *decay* works as an adaptive
545 threshold that can promote potentiation or depression according to the current state, penalizing the
546 FIC in a magnitude-dependent way (the larger the FIC, the larger the penalization), which in turn
547 stabilizes the activity without imposing a saturation point. It is worth noticing that in contrast to other
548 models that externally impose saturation through a sigmoid, in the HDMF a comparable constraint
549 emerged endogenously through the adaptive threshold. In this sense, the regulation of activity was not
550 hardwired but self-organized within the system's dynamics. This control mechanism prevents
551 runaway excitation characteristic of pathological states such as epilepsy.
552 **Broader Applications and Future Directions**

553 Because the HDMF can generate slow-wave activity, it also provides a framework to introduce
554 explicit temporal structure into adaptive processes. This opens the possibility of modeling transient,
555 spatially localized slow-wave events resembling “local sleep” during otherwise wake-like activity⁶¹.
556 When combined with probe-based resting-state paradigms, such as experience sampling during mind
557 wandering⁶², this framework could be used to test whether specific temporal rules of adaptation align
558 with moments of mind blanking or attentional disengagement.

559 In addition, we showed that the HDMF can accommodate neuromodulatory influences as external
560 perturbations acting on internal network dynamics. This approach can be extended to other forms of
561 extrinsic input, including auditory or visual stimulation⁶³,⁶⁴. In this context, the HDMF provides a
562 setting in which slow resting-state fluctuations coexist with adaptive mechanisms that regulate
563 responses to external information, allowing resting-state dynamics and task-related processing to be
564 studied within the same model.

565 Beyond reproducing a diverse set of brain states, homeostatic inhibitory control in the HDMF could
566 prove valuable in personalized connectome-based modeling, where individual structural or PET
567 receptor maps may push local network dynamics into regimes of pathological hyper- or hypo-
568 activity⁶⁵,⁶⁶, or cortical development⁶⁷. Adding region-specific target firing rates or multi-timescale
569 plasticity might further enhance the model’s ability to simulate neurological disorders –such as
570 epilepsy, brain lesions or neurodegenerative conditions– and to explore therapeutic strategies.
571 Combining homeostatic and Hebbian-like rules could help capture learning-dependent processes
572 alongside baseline stability.

573
574 **Limitations**

575 Although the HDMF addresses important hurdles in large-scale brain modeling, several limitations
576 remain. First, our approach focuses on a single target firing rate, which may not reflect the multiple
577 concurrent homeostatic targets in the brain –such as metabolic constraints or different firing rates
578 across regions. Also, the LR- decay relationship need to be obtained for every combination of
579 structural connectivity and target firing rate. However, this can be done in a comparatively shorter
580 time than the original pseudo-gradient algorithm (which also lacks the FIC dynamics of the HDMF),
581 and a potential analytical solution could provide a shortcut for deriving this relationship. Second, we
582 tested primarily cortical connectivity, leaving open questions about subcortical and thalamocortical
583 contributions. Third, we did not explore interactions between homeostatic plasticity and other
584 plasticity mechanisms (e.g., Hebbian or STDP), which could be critical for understanding learning
585 and memory formation. Nevertheless, learning and memory may occur at much finer spatial scales
586 than the ones modelled with whole-brain models, making these limitations transversal to the field of
587 whole-brain modelling. Finally, further empirical validation, including clinical populations or
588 perturbation experiments, will be necessary to confirm the predictive utility of the HDMF across
589 diverse real-world contexts.

590 **Conclusions**

591 By adding a homeostatic plasticity rule to the DMF model we showed how this mechanism is
592 beneficial for ensuring stability and enhancing realism. This also allowed the modelisation of slow
593 waves, which is unprecedented in the DMF model. We finally propose that by adding heterogeneity in
594 homeostatic properties we can recreate chimeric-like states, suggesting a potential mechanism for
595 sleep dissociative states. Altogether, the HDMF illustrates how embedding homeostatic self-
596 regulation into whole-brain models can overcome fragility, extend dynamical repertoire, and provide
597 mechanistic links to complex brain states.
599 **Materials and methods**

600 **Participants**

601 **fMRI data**

602 The fMRI data used here was obtained from previous studies⁶⁸. Data were acquired from 13 healthy
603 volunteers (mean age ± s.d., 42.54 ± 13.64 years; 7 females; no history of neurological disorders) on a
604 3 T General Electric Signa System scanner. T2*-weighted whole-brain resting-state images were
605 obtained through a gradient-echo EPI sequence with axial orientation (200 volumes, 48 slices, slice
606 thickness: 3 mm, TR/TE: 2400ms/30 ms, voxel size: 3.4375 × 3.4375 × 3.4375 mm, flip angle:
607 90°, FOV: 220 mm²).

608 Resting-state fMRI data were preprocessed using FSL, following previously established procedures.
609 Briefly, data were processed with MELODIC (Multivariate Exploratory Linear Optimised
610 Decomposition into Independent Components)⁶⁹. Preprocessing steps included removal of the first
611 five volumes, motion correction using MCFLIRT⁷⁰, brain extraction with the Brain Extraction Tool
612 (BET)⁷¹, spatial smoothing with a 5-mm full-width at half-maximum Gaussian kernel, rigid-body
613 registration, and high-pass temporal filtering with a 0.01-Hz cutoff. Single-session independent
614 component analysis (ICA) with automatic dimensionality estimation was then performed. Noise
615 components, including motion-related and physiological artefacts, were identified and regressed out
616 on a subject-by-subject basis using FIX (FMRIB's ICA-based X-noiseifier⁷²). Finally, images were
617 co-registered to MNI space, and regional time series were extracted using the AAL atlas⁷³, yielding
618 signals from 90 cortical regions. All resulting regional time series were then converted to z-scores by
619 subtracting the mean and dividing by the standard deviation of each time series. These z-scored time
620 series formed the basis for downstream connectivity analyses.

621 **DWI data**

623 Diffusion-weighted MRI data were obtained from a previously collected cohort of 16 healthy right-
624 handed participants (5 women; mean age ± s.d., 24.8 ± 2.5 years) scanned at Aarhus University,
625 Denmark⁷⁴. The study protocol was approved by the local ethics committees, and all participants
626 provided written informed consent. Imaging was performed in a single session on a 3 T Siemens
627 Skyra scanner. Diffusion data were acquired with 62 nonlinear diffusion-encoding directions (b =
628 1500 s/mm²), interleaved with non-diffusion-weighted images (b = 0), using a repetition time of 9000
629 ms and an echo time of 84 ms. Images were collected with a voxel size of 1.98 × 1.98 × 2 mm³ and
630 reconstructed on a 106 × 106 matrix, with both anterior-to-posterior and posterior-to-anterior phase-
631 encoding directions. Motion and eddy-current correction was performed using MRtrix⁷⁵. Diffusion
632 images were co-registered to each participant's T1-weighted anatomical image and then to MNI
633 space. The AAL atlas was then transformed from MNI space to individual native diffusion space
634 using the inverse transformations, preserving discrete labels via nearest-neighbour interpolation. Fibre
635 orientation distributions were estimated using constrained spherical deconvolution, and whole-brain
636 probabilistic tractography was performed to reconstruct inter-regional pathways. Streamlines were
637 aggregated between atlas-defined regions to generate individual structural connectivity matrices,
638 which were subsequently averaged across participants to obtain a group-representative connectome
639 Further details are provided in the original study⁷⁴.

640

641 **PET-derived receptor density data**

642 We incorporated publicly available positron emission tomography (PET) receptor density maps from
643 the normative neurotransmitter atlas published by Hansen et al. ³⁸ These maps include whole-brain
644 estimates of multiple neurotransmitter receptor systems aggregated across healthy individuals,
allowing in-vivo quantification of receptor distributions. For the serotonin 2A receptor (5-HT₂A) specifically, the atlas is based on PET imaging using the agonist radioligand Cimbi-36, which binds with high affinity to 5-HT₂A receptors and has been validated as a marker of their density in the living human brain. PET tracer maps were registered to MNI space using the Neuromaps⁷⁶ library and parcellated using AAL atlas to extract mean receptor density values for each region denoted as $D_i$. Density values were normalised by dividing by the maximum value, ensuring that $\max(D_i) = 1$.

651

### Parcellation

For all neuroimaging data, we employed the AAL atlas with 90 parcels⁷³.

654

### Dynamic mean field (DMF) model

We used the DMF to simulate the brain activity of each brain region, as proposed in ⁷⁷. The DMF simulates whole-brain activity by defining the local dynamics of each region as two interconnected pools of excitatory (E) and inhibitory (I) neurons while long range interactions are only between the E pools of each region, adjusted by the structural connectivity $C_{ij}$ (for details, refer to Methods - Anatomical Connectivity). The model is based on the earlier work of ²², where excitatory synaptic currents $I_i^{(E)}$ are mediated by NMDA receptors, while inhibitory currents $I_i^{(I)}$ are mediated by GABA receptors.

663 In detail, the DMF is defined by the following system of coupled differential equations:

664

$$I_i^{(E)} = I_i^{\text{ext}} + W_E I_0 + w + \zeta J_{\text{NMDA}} S_i^{(E)} + G J_{\text{NMDA}} \sum_j C_{ij} S_j^{(E)} - J_i S_i^{(I)} \zeta$$

665

$$I_i^{(I)} = W_I I_0 + J_{\text{NMDA}} S_i^{(E)} - S_i^{(I)}$$

666

$$r_i^{(E)} = H^{(E)}(I \zeta i^{(E)}) \zeta = \zeta \zeta \frac{M_i(a_E I_i^{(E)} - b_E)}{1 - \exp(-d_E M_i(a_E I_i^{(E)} - b_E))}$$

667

$$r_i^{(I)} = H^{(I)}(I \zeta i^{(I)}) \zeta = \zeta \zeta \frac{M_i(a_I I_i^{(I)} - b_I)}{1 - \exp(-d_I M_i(a_I I_i^{(I)} - b_I))}$$

668

$$M_i = s D_i + 1$$

669

$$\frac{d S_i^{(E)}(t)}{dt} \zeta \frac{d S_i^{(E)}}{\tau_E} + (1 - S_i^{(E)}) \gamma r_i^{(E)} + \sigma v_i(t)$$

670

$$\frac{d S_i^{(I)}(t)}{dt} \zeta \frac{d S_i^{(I)}}{\tau_I} + r_i^{(I)} + \sigma v_i(t)$$

671

672

673

Here, $I_i^{(E,I)}$ (in nA) represents the total input current, and $r_i^{(E,I)}$ (in Hz) represents the respective firing rates. $S_i^{(E,I)}$ represent the synaptic gating variables of each pool. The neuronal nonlinear response functions, $H^{(E,I)}$, convert the total input currents received by the E and I pools into firing rates, $r_i^{(E,I)}$, based on the input-output function proposed by ⁷⁸. $H^{(E,I)}$ is also known as the F-I curve. The model employed specific gain factors, threshold currents, and constants to determine the shape of the curvature of $H^{(E,I)}$ around the threshold. s is called the neuromodulatory gain, which scales the impact of the receptor density $D_i$ in each region, which is mapped into $M_i^{27}$. G is a free parameter called the global coupling, which scales the influence of the network on each region. See Table 1 for the name and value of each model parameter.
### 684 Feedback inhibitory control (FIC)

The parameter $J_i$, referred to as FIC, defines the synaptic weight from the inhibitory (I) to the excitatory (E) population. As demonstrated by Deco et al²³ precise tuning of this parameter allows the excitatory firing rates to be clamped around 3.4 Hz, maintaining them within a biologically plausible regime. Because $J_i$ must be optimized for each value of the global coupling parameter G, this step substantially increases the computational cost of the simulations. However, Herzog et al.²⁴ proposed a heuristic linear solution for the FIC which derives the optimal value from a linear relationship between G and the anatomical connectivity strength $\beta$ as follows:

$$J_i = \alpha G \beta_i + 1 \qquad 692$$

where $\alpha$ is a free parameter –representing the global E/I balance– that depends on the parcellation and in the target firing rate, and $\beta_i = \sum_j^\square C_{ij}$, i.e. the anatomical connectivity strength. This solution achieves stability of firing rates up to critical value of G, where this solution cannot control the firing rates anymore and the system enters into an hyperactivity regime. See Herzog et al²⁴ for a detailed description of this approximation.

698

### 699 Inhibitory plasticity rule with decay

Following Hellyer et al.¹⁴, we implemented a local inhibitory plasticity rule to dynamically adapt $J_i$ such that the excitatory average firing rates converge to a target firing rate $\rho$. We extended this plasticity rule to include a decay term $\tau_{decay}$ on the synaptic weight. This rule follows:

$$\frac{d J_i}{dt} = \frac{-J_i}{\tau_{decay}} + LR r_i^{(I)} (r_i^{(E)} - \rho) \qquad 703$$

LR is the learning rate and represents the speed of the update of the synaptic weight by an active mechanisms while $\tau_{decay}$ represents an adaptive threshold mechanism: when excitatory firing rates exceed the target above a threshold, inhibition will be potentiated; and vice versa, firing rates below the threshold will induce a depression in inhibition.

708

### 709 Homeostatic dynamic mean field model (HDMF)

The inhibitory plasticity rule enables a dynamic control of the local inhibition based on the relationship between the E and I pool and the excess or lack of excitatory activity. However, for some combination of parameters the average excitatory firing rate will not converge to the target, but will converge to some other value, failing to achieve homeostasis. To fully achieve homeostasis, i.e. that the average firing rates converge to the expected target, we used the following heuristic relationship between the LR and $\tau_{decay}$:

$$\tau_{decay} = \chi LR^{\epsilon} \qquad 716$$

Then, $\tau_{decay}$ has a power-law relationship with LR controlled by the coefficient $\chi$ and the exponent $\epsilon$, which are obtained by numerical simulations (see Figure 2A). This implies that $\tau_{decay}$ depends on LR and then the plasticity rule can be fully expressed in terms of LR. As shown in Figure 2A, these parameters depend on the target firing rate but, regardless of this target, this relationship achieves maximal homeostasis in a wide range of G values, i.e. minimizes the difference between the target and the average firing rates. To distinguish the DMF from its version with homeostasis, we call the former homeostatic DMF (HDMF).

724
725 **Simulating BOLD signals**

726 To convert the simulated mean field activity into a BOLD signal, we used the generalized 727 hemodynamic model by Stephan et al.⁷⁹. The BOLD signal for each brain region i was calculated 728 based on the firing rate of the excitatory pools, $r_i^{(E)}$. In this framework, an increase in firing rate 729 caused a rise in a vasodilatory signal, $s_i$, which was also subject to auto-regulatory feedback. 730 Consequently, changes in blood inflow, $f_i$, were influenced by this signal, leading to variations in 731 blood volume, $v_i^\square$ and deoxyhemoglobin content, $q_i$. Also, for sustained firing increase, bold response 732 saturates ⁸⁰,⁸¹. We modeled this neurovascular coupling with baseline invariant sigmoid $g$. The 733 interactions between these biophysical variables were described by the following equations:

734 $z_i(t) = g(r_i^{(E)}(t)) - g(\rho)$
735 $g(u) = \frac{1}{1 + \exp[-0.25(u - 12)]}$.

736 $\frac{d s_i(t)}{dt} \approx 0.5 z_i(t) + 3 - s_i k_i - \gamma_i (f_i - 1)$

738 $\frac{d f_i(t)}{dt} \approx s_i$

739 $\tau_i \frac{d v_i(t)}{dt} \approx f_i - v_i^{1/\alpha}$

740 $\tau_i \frac{d q_i(t)}{dt} \approx \frac{f_i(1 - (1 - \zeta_i)^{1/f_i})}{\zeta_i} - \frac{q_i v_i^{1/\alpha}}{v_i}$

741

742 where $\zeta$ is the resting oxygen extraction fraction, $\tau$ is a time constant and $\alpha$ that represents the 743 resistance of the veins. Finally, the BOLD signal in each area i, $B_i$, is a static nonlinear function of 744 volume, $v_i$, and deoxyhemoglobin content, $q_i$, that comprises a volume-weighted sum of extra- and 745 intravascular signals:

746 $B_n = V_0 \dot{\zeta}$

748

749

750 All biophysical parameters were adopted from reference⁷⁹. To isolate the frequency range most 751 relevant for resting-state activity, both empirical and simulated BOLD signals were band-pass filtered 752 between 0.01 and 0.1 Hz.

753

754

755 **Fit to empirical data**

756 Model parameters were optimized selecting optimal points in an exploratory grid. The objective 757 function was defined either as: (i) the Pearson correlation between the simulated and population- 758 averaged empirical FC matrices, or (ii) the Kolmogorov-Smirnov (KS) distance between the 759 histogram of the pooled empirical Functional Connectivity Dynamics (FCD) and that of the simulated 760 FCD. Each simulation generated 460 seconds of BOLD signals sampled at TR = 2.4 s, yielding a 761 number of time points comparable to the empirical data. The optimization assumed a stochastic 762 objective function, allowing the algorithm to randomly initialize each simulation and estimate 763 variability due to different initial conditions.

764
765 Functional connectivity and functional connectivity dynamics

766 Functional connectivity (FC) was estimated by computing the Pearson correlation between the 767 empirical (simulated) fMRI time series. For the functional connectivity dynamics (FCD) we used the 768 same approach, but we computed the FC in a sliding window of 30 TR with 29 points of overlap. 769 Then, the upper triangular part of the FC of each window was transformed into a vector and the 770 Pearson correlation between all the time windows was computed, obtaining the FCD matrix. This 771 matrix captures the correlation of FC patterns across time. 772

773

774

775 Detection of oscillatory dynamics

776 We computed the power spectral density (PSD) of neuronal firing-rate traces using Welch's method 777 with a 1000 Hz sampling frequency, segment lengths of 4000 samples, and a 50% overlap. Each PSD 778 was evaluated only up to 100 Hz, and the maximum point of this spectral region was identified to 779 determine the dominant oscillation frequency and its corresponding power. This approach produced 780 four principal outputs: the distribution of frequencies (freqs), the PSD values for each frequency 781 (psd), the peak frequency (max freqs), and the maximal PSD within 0-100 Hz (max power), 782 providing a concise characterization of each signal's strongest oscillatory component. 783

784

785 Computing regime type index

786 To quantify how strongly a network endowed with region-specific plasticity deviated from either of 787 the two baseline uniform regimes, we defined the metric “Regime type index” (RTI) that yields a 788 value for each region. For each analysis we first identified a subset of N cortical regions that were 789 assigned a learning-rate parameter LR<sub>HET</sub> while all remaining areas retained the reference uniform 790 value LR<sub>UNI</sub>. One hundred independent simulations were then run with this heterogeneous vector. For 791 every run we extracted the dominant oscillation frequency of each region (see Detection of oscillatory 792 dynamics) and averaged these frequencies across simulations, yielding f<sub>HET</sub>. Each region uses two 793 reference values obtained from fully uniform networks f<sub>HIGH</sub> from the “asynchronous” 794 high-learning-rate regime and f<sub>LOW</sub> from the “slow-wave-like” low-learning-rate regime. The RTI 795 rescales the position of f<sub>HET</sub> between these anchors according to 796

797

$$
RTI = \zeta \frac{f_{HET} - f_{LOW}}{f_{HIGH} - f_{LOW}} \zeta - 1
$$

798

799 With this normalization the index equals -1 when the heterogeneous network is indistinguishable 800 from the low-learning-rate regime and +1 when it fully reproduces the high-learning-rate regime, 801 while intermediate values indicate partial convergence toward either regime. 802

803

804 Statistics

801 The distributions of goodness-of-fit values for the FC and FCD histograms were compared using Kolmogorov-Smirnov test.

803

804 References

805

806 1. Brown, E. N., Purdon, P. L. & Dort, C. J. V. General Anesthesia and Altered States of 807 Arousal: A Systems Neuroscience Analysis. Annu. Rev. Neurosci. 34, 601–628 (2011).
808 2. Luppi, A. I. et al. Consciousness-specific dynamic interactions of brain integration and functional diversity. *Nat. Commun.* **10**, 4616 (2019).

810 3. Sanz Perl, Y. et al. Nonequilibrium brain dynamics as a signature of consciousness. *Phys. Rev. E* **104**, 014411 (2021).

812 4. Varela, F., Lachaux, J.-P., Rodriguez, E. & Martinerie, J. The brainweb: Phase synchronization and large-scale integration. *Nat. Rev. Neurosci.* **2**, 229–239 (2001).

814 5. Dehaene, S., Changeux, J.-P. & Naccache, L. The Global Neuronal Workspace Model of Conscious Access: From Neuronal Architectures to Clinical Applications. in *Characterizing Consciousness: From Cognition to the Clinic?* (eds. Dehaene, S. & Christen, Y.) 55–84 (Springer, Berlin, Heidelberg, 2011). doi:10.1007/978-3-642-18015-6_4.

818 6. Wen, W. & Turrigiano, G. G. Keeping Your Brain in Balance: Homeostatic Regulation of Network Function. *Annu. Rev. Neurosci.* **47**, 41–61 (2024).

820 7. Foster, M. & Scheinost, D. Brain states as wave-like motifs. *Trends Cogn. Sci.* **28**, 492–503 (2024).

822 8. Chen, L., Li, X., Tjia, M. & Thapliyal, S. Homeostatic plasticity and excitation-inhibition balance: The good, the bad, and the ugly. *Curr. Opin. Neurobiol.* **75**, 102553 (2022).

824 9. Cofré, R. et al. Whole-Brain Models to Explore Altered States of Consciousness from the Bottom Up. *Brain Sci.* **10**, 626 (2020).

826 10. Coronel-Oliveros, C., Cofré, R. & Orio, P. Cholinergic neuromodulation of inhibitory interneurons facilitates functional integration in whole-brain models. *PLOS Comput. Biol.* **17**, e1008737 (2021).

829 11. Herzog, R. et al. A whole-brain model of the neural entropy increase elicited by psychedelic drugs. *Sci. Rep.* **13**, (2023).

831 12. Luppi, A. I. et al. Computational modelling in disorders of consciousness: Closing the gap towards personalised models for restoring consciousness. *NeuroImage* **275**, 120162 (2023).

833 13. Patow, G., Martin, I., Sanz Perl, Y., Kringelbach, M. L. & Deco, G. Whole-brain modelling: an essential tool for understanding brain dynamics. *Nat. Rev. Methods Primer* **4**, 1–2 (2024).
14. Hellyer, P. J., Jachs, B., Clopath, C. & Leech, R. Local inhibitory plasticity tunes macroscopic brain dynamics and allows the emergence of functional brain networks. *NeuroImage* **124**, 85–95 (2016).

15. Rocha, R. P., Koçillari, L., Suweis, S., Corbetta, M. & Maritan, A. Homeostatic plasticity and emergence of functional networks in a whole-brain model at criticality. *Sci. Rep.* **8**, 15682 (2018).

16. Goldman, J. S. et al. A comprehensive neural simulation of slow-wave sleep and highly responsive wakefulness dynamics. *Front. Comput. Neurosci.* **16**, (2023).

17. Stasinski, J. et al. Homeodynamic feedback inhibition control in whole-brain simulations. *PLOS Comput. Biol.* **20**, e1012595 (2024).

18. Rocha, R. P., Zorzi, M. & Corbetta, M. Role of homeostatic plasticity in critical brain dynamics following focal stroke lesions. *Sci. Rep.* **14**, 31631 (2024).

19. Coronel-Oliveros, C. et al. Viscous dynamics associated with hypoexcitation and structural disintegration in neurodegeneration via generative whole-brain modeling. *Alzheimers Dement.* **20**, 3228–3250 (2024).

20. Stefanovski, L. et al. Linking Molecular Pathways and Large-Scale Computational Modeling to Assess Candidate Disease Mechanisms and Pharmacodynamics in Alzheimer’s Disease. *Front. Comput. Neurosci.* **13**, (2019).

21. Luppi, A. I. et al. An engine for systematic discovery of cause-effect relationships between brain structure and function. 2025.04.05.647237 Preprint at https://doi.org/10.1101/2025.04.05.647237 (2025).

22. Wong, K.-F. & Wang, X.-J. A Recurrent Network Mechanism of Time Integration in Perceptual Decisions. *J. Neurosci.* **26**, 1314–1328 (2006).

23. Deco, G. et al. How Local Excitation–Inhibition Ratio Impacts the Whole Brain Dynamics. *J. Neurosci.* **34**, 7886–7898 (2014).

24. Herzog, R. et al. Neural mass modelling for the masses: Democratising access to whole-brain biophysical modelling with FastDMF. *Netw. Neurosci.* 1–38 (2024) doi:10.1162/netn_a_00410.

25. Perl, Y. S. et al. The impact of regional heterogeneity in whole-brain dynamics in the presence of oscillations. *Netw. Neurosci.* **7**, 632–660 (2023).
26. Mindlin, I. et al. Whole brain modelling for simulating pharmacological interventions on patients with disorders of consciousness. *Commun. Biol.* **7**, (2024).

27. Deco, G. et al. Whole-Brain Multimodal Neuroimaging Model Using Serotonin Receptor Maps Explains Non-linear Functional Effects of LSD. *Curr. Biol.* **28**, 3065-3074.e6 (2018).

28. Luppi, A. I. et al. Whole-brain modelling identifies distinct but convergent paths to unconsciousness in anaesthesia and disorders of consciousness. *Commun. Biol.* **5**, 1–15 (2022).

29. Gatica, M. et al. High-order functional redundancy in ageing explained via alterations in the connectome in a whole-brain model. *PLOS Comput. Biol.* **18**, e1010431 (2022).

30. Yang, G. J. et al. Functional hierarchy underlies preferential connectivity disturbances in schizophrenia. *Proc. Natl. Acad. Sci.* **113**, E219–E228 (2016).

31. Moguilner, S. et al. Biophysical models applied to dementia patients reveal links between geographical origin, gender, disease duration, and loss of neural inhibition. *Alzheimer's Res. Amp Ther.* **16**, (2024).

32. Antelmi, E. et al. From state dissociation to status dissociatus. *Sleep Med. Rev.* **28**, 5–17 (2016).

33. Sodré, M. E., Wießner, I., Irfan, M., Schenck, C. H. & Mota-Rolim, S. A. Awake or Sleeping? Maybe Both... A Review of Sleep-Related Dissociative States. *J. Clin. Med.* **12**, 3876 (2023).

34. Bansal, K. et al. Cognitive chimera states in human brain networks. *Sci. Adv.* **5**, eaau8535 (2019).

35. Huntenburg, J. M., Bazin, P.-L. & Margulies, D. S. Large-Scale Gradients in Human Cortical Organization. *Trends Cogn. Sci.* **22**, 21–31 (2018).

36. Andrillon, T. & Oudiette, D. What is sleep exactly? Global and local modulations of sleep oscillations all around the clock. *Neurosci. Biobehav. Rev.* **155**, 105465 (2023).

37. Vogels, T. P., Sprekeler, H., Zenke, F., Clopath, C. & Gerstner, W. Inhibitory Plasticity Balances Excitation and Inhibition in Sensory Pathways and Memory Networks. *Science* **334**, 1569–1573 (2011).
890 38. Hansen, J. Y. et al. Mapping neurotransmitter systems to the structural and functional organization of the human neocortex. *Nat. Neurosci.* **25**, 1569–1581 (2022).

892 39. Sanchez-Vives, M. V. Origin and dynamics of cortical slow oscillations. *Curr. Opin. Physiol.* **15**, 217–223 (2020).

894 40. Wilckens, K. A., Ferrarelli, F., Walker, M. P. & Buysse, D. J. Slow-Wave Activity Enhancement to Improve Cognition. *Trends Neurosci.* **41**, 470–482 (2018).

896 41. Sheybani, L. et al. Wake slow waves in focal human epilepsy impact network activity and cognition. *Nat. Commun.* **14**, 7397 (2023).

898 42. Marsh, B. et al. Emergent effects of synaptic connectivity on the dynamics of global and local slow waves in a large-scale thalamocortical network model of the human brain. *PLOS Comput. Biol.* **20**, e1012245 (2024).

901 43. Sacha, M., Tesler, F., Cofre, R. & Destexhe, A. A computational approach to evaluate how molecular mechanisms impact large-scale brain activity. *Nat. Comput. Sci.* **5**, 405–417 (2025).

903 44. Appelbaum, L. G., Shenasa, M. A., Stolz, L. & Daskalakis, Z. Synaptic plasticity and mental health: methods, challenges and opportunities. *Neuropsychopharmacology* **48**, 113–120 (2023).

905 45. Castro, S., El-Deredy, W., Battaglia, D. & Orio, P. Cortical ignition dynamics is tightly linked to the core organisation of the human connectome. *PLOS Comput. Biol.* **16**, e1007686 (2020).

907 46. Rubinov, M., Kötter, R., Hagmann, P. & Sporns, O. Brain connectivity toolbox: a collection of complex network measurements and brain connectivity datasets. *NeuroImage* **47**, S169 (2009).

909 47. Abeysuriya, R. G. et al. A biophysical model of dynamic balancing of excitation and inhibition in fast oscillatory large-scale networks. *PLOS Comput. Biol.* **14**, e1006007 (2018).

911 48. Soldado-Magraner, S., Seay, M. J., Laje, R. & Buonomano, D. V. Paradoxical self-sustained dynamics emerge from orchestrated excitatory and inhibitory homeostatic plasticity rules. *Proc. Natl. Acad. Sci. U. S. A.* **119**, e2200621119 (2022).

914 49. Compte, A., Sanchez-Vives, M. V., McCormick, D. A. & Wang, X.-J. Cellular and Network Mechanisms of Slow Oscillatory Activity (<1 Hz) and Wave Propagations in a Cortical Network Model. *J. Neurophysiol.* **89**, 2707–2725 (2003).
917 50. Holcman, D. & Tsodyks, M. The Emergence of Up and Down States in Cortical Networks.
*PLOS Comput. Biol.* **2**, e23 (2006).

918
919
51. Bartram, J. et al. Cortical Up states induce the selective weakening of subthreshold synaptic
inputs. *Nat. Commun.* **8**, 665 (2017).

921
922
52. Krishnan, G. P. et al. Cellular and neurochemical basis of sleep stages in the thalamocortical
network. *eLife* **5**, e18607 (2016).

923
924
53. Haugland, S. W. The changing notion of chimera states, a critical review. *J. Phys. Complex*.
**2**, 032001 (2021).

925
926
54. Patow, G., Martin, I., Sanz Perl, Y., Kringelbach, M. L. & Deco, G. Whole-brain modelling:
an essential tool for understanding brain dynamics. *Nat. Rev. Methods Primer* **4**, 53 (2024).

927
928
55. Aru, J., Siclari, F., Phillips, W. A. & Storm, J. F. Apical drive—A cellular mechanism of
dreaming? *Neurosci. Biobehav. Rev.* **119**, 440–455 (2020).

929
930
56. LaBerge, S., LaMarca, K. & Baird, B. Pre-sleep treatment with galantamine stimulates lucid
dreaming: A double-blind, placebo-controlled, crossover study. *PLOS ONE* **13**, e0201246 (2018).

931
932
57. Alnagger, N. L. N. et al. A Virtual Clinical Trial of Psychedelics to Treat Patients With
Disorders of Consciousness. *Adv. Sci. n/a*, e11780.

933
934
58. Stasinski, J. et al. Homeodynamic feedback inhibition control in whole-brain simulations.
*PLOS Comput. Biol.* **20**, e1012595 (2024).

935
936
937
59. Hellyer, P. J., Jachs, B., Clopath, C. & Leech, R. Local inhibitory plasticity tunes
macroscopic brain dynamics and allows the emergence of functional brain networks. *NeuroImage
124*, 85–95 (2016).

938
939
60. Goldman, J. S. et al. A comprehensive neural simulation of slow-wave sleep and highly
responsive wakefulness dynamics. *Front. Comput. Neurosci.* **16**, (2023).

940
941
61. Andrillon, T., Burns, A., Mackay, T., Windt, J. & Tsuchiya, N. Predicting lapses of attention
with sleep-like slow waves. *Nat. Commun.* **12**, 3657 (2021).

942
943
944
62. Mortaheb, S. et al. Mind blanking is a distinct mental state linked to a recurrent brain profile
of globally positive connectivity during ongoing mentation. *Proc. Natl. Acad. Sci.* **119**,
e2200511119 (2022).
945 63. Bekinschtein, T. A. et al. Neural signature of the conscious processing of auditory regularities. *Proc. Natl. Acad. Sci. 106*, 1672–1677 (2009).

947 64. Türker, B., Belloli, L., Owen, A. M., Naci, L. & Sitt, J. D. Processing of the same narrative stimuli elicits common functional connectivity dynamics between individuals. *Sci. Rep.* **13**, 21260 (2023).

950 65. Maestú, F., de Haan, W., Busche, M. A. & DeFelipe, J. Neuronal excitation/inhibition imbalance: core element of a translational perspective on Alzheimer pathophysiology. *Ageing Res. Rev.* **69**, 101372 (2021).

953 66. Ranasinghe, K. G. et al. Altered excitatory and inhibitory neuronal subpopulation parameters are distinctly associated with tau and amyloid in Alzheimer's disease. *eLife* **11**, e77850 (2022).

955 67. Zhang, S. et al. In vivo whole-cortex marker of excitation-inhibition ratio indexes cortical maturation and cognitive ability in youth. *Proc. Natl. Acad. Sci. U. S. A.* **121**, e2318641121 (2024).

958 68. Escrichs, A. et al. Unifying turbulent dynamics framework distinguishes different brain states. *Commun. Biol.* **5**, 1–13 (2022).

960 69. Beckmann, C. F. & Smith, S. M. Probabilistic independent component analysis for functional magnetic resonance imaging. *IEEE Trans. Med. Imaging* **23**, 137–152 (2004).

962 70. Jenkinson, M., Bannister, P., Brady, M. & Smith, S. Improved Optimization for the Robust and Accurate Linear Registration and Motion Correction of Brain Images. *NeuroImage* **17**, 825–841 (2002).

965 71. Smith, S. M. Fast robust automated brain extraction. *Hum. Brain Mapp.* **17**, 143–155 (2002).

966 72. Griffanti, L. et al. ICA-based artefact removal and accelerated fMRI acquisition for improved resting state network imaging. *NeuroImage* **95**, 232–247 (2014).

968 73. Tzourio-Mazoyer, N. et al. Automated Anatomical Labeling of Activations in SPM Using a Macroscopic Anatomical Parcellation of the MNI MRI Single-Subject Brain. *NeuroImage* **15**, 273–289 (2002).

971 74. Deco, G., Kringelbach, M. L., Jirsa, V. K. & Ritter, P. The dynamics of resting fluctuations in the brain: metastability and its dynamical cortical core. *Sci. Rep.* **7**, 3095 (2017).
973 75. Tournier, J.-D. et al. MRtrix3: A fast, flexible and open software framework for medical
974 image processing and visualisation. *NeuroImage* **202**, 116137 (2019).

975 76. Markello, R. D. et al. neuromaps: structural and functional interpretation of brain maps. *Nat. Methods* **19**, 1472–1479 (2022).
976

977 77. Deco, G. et al. How Local Excitation–Inhibition Ratio Impacts the Whole Brain Dynamics. *J. Neurosci.* **34**, 7886–7898 (2014).
978

979 78. Abbott, L. F. & Chance, F. S. Drivers and modulators from push-pull and balanced synaptic
980 input. in *Progress in Brain Research* vol. 149 147–155 (Elsevier, 2005).

981 79. Stephan, K. E., Weiskopf, N., Drysdale, P. M., Robinson, P. A. & Friston, K. J. Comparing
982 hemodynamic models with DCM. *NeuroImage* **38**, 387–401 (2007).

983 80. Tesler, F., Linne, M.-L. & Destexhe, A. Modeling the relationship between neuronal activity
984 and the BOLD signal: contributions from astrocyte calcium dynamics. *Sci. Rep.* **13**, 6451 (2023).

985 81. Rosa, M. J., Kilner, J. M. & Penny, W. D. Bayesian Comparison of Neurovascular Coupling
986 Models Using EEG-fMRI. *PLOS Comput. Biol.* **7**, e1002070 (2011).

987

988 Acknowledgments

989 We thank the financial support of the FLAG-ERA research funding organisation (project
990 ModelDXConsciousness).

991

992

993 Author contributions:

994 Conceptualization: Iván Mindlin, Rubén Herzog

995 Methodology: Iván Mindlin, Rubén Herzog, Andrea Luppi, Carlos Coronel

996 Investigation: Iván Mindlin, Rubén Herzog

997 Visualization: Iván Mindlin

998 Supervision: Rubén Herzog

999 Writing-original draft: Iván Mindlin, Rubén Herzog

1000 Writing-review & editing: Iván Mindlin, Carlos Coronel, Jacobo D. Sitt, Rodrigo Cofré,

1001 Andrea Luppi, Thomas Andrillon, Yonatan Sanz Perl and Rubén Herzog

1002

1003

1004 Competing interests

1005 Authors declare that they have no competing interests.

1006

1007 Data and materials availability

1008 Computational codes to run the HDMF and reproduce the results are available at

1009 https://github.com/Picardian14/fastHDMF-code.git. We’ve also made a pure Python

1010 implementation of the same model. https://github.com/carlosmig/Homo_DMF.
<table><thead><tr><th>Parameter name</th><th>Value</th><th>Description</th></tr></thead><tbody><tr><td>dt</td><td>0.1</td><td>DMF integration step (ms)</td></tr><tr><td>E</td><td>100</td><td>NMDA characteristic time (ms)</td></tr><tr><td>I</td><td>10</td><td>GABA characteristic time (ms)</td></tr><tr><td></td><td>0.641</td><td>Kinetic parameter of excitation</td></tr><tr><td></td><td>0.01</td><td>Noise standard deviation (nA)</td></tr><tr><td>J<sub>NMDA</sub></td><td>0.15</td><td>Excitatory synaptic coupling (nA)</td></tr><tr><td>I<sub>0</sub></td><td>0.382</td><td>Effective external input (nA)</td></tr><tr><td>W<sub>E</sub></td><td>1</td><td>External-to-excitatory coupling</td></tr><tr><td>W<sub>I</sub></td><td>0.7</td><td>External-to-inhibitory coupling</td></tr><tr><td>w<sub>+</sub></td><td>1.4</td><td>Local excitatory recurrence</td></tr><tr><td>a<sub>E</sub></td><td>0.16</td><td>Excitatory conductance</td></tr><tr><td>b<sub>E</sub></td><td>0.403</td><td>Excitatory threshold for nonlinearity</td></tr><tr><td>d<sub>E</sub></td><td>310</td><td>Excitatory nonlinear shape parameter</td></tr><tr><td>a<sub>I</sub></td><td>0.087</td><td>Inhibitory conductance</td></tr><tr><td>b<sub>I</sub></td><td>0.288</td><td>Inhibitory threshold for nonlinearity</td></tr><tr><td>d<sub>I</sub></td><td>615</td><td>Inhibitory nonlinear shape parameter</td></tr><tr><td>TR</td><td>2</td><td>BOLD signal sampling frequency (s)</td></tr><tr><td>dtt</td><td>0.001</td><td>BW integration step (s)</td></tr></tbody></table>