# Metastable oscillatory modes emerge from synchronization in the brain spacetime connectome

Joana Cabral, Francesca Castaldo, Jakub Vohryzek, Vladimir Litvak, Christian Bick, Renaud Lambiotte, Karl Friston, Morten L Kringelbach, Gustavo Deco

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

## To cite this version:

Joana Cabral, Francesca Castaldo, Jakub Vohryzek, Vladimir Litvak, Christian Bick, et al. Metastable oscillatory modes emerge from synchronization in the brain spacetime connectome. Communications Physics, 2022, 5 (1), pp.184. 10.1038/s42005-022-00950-y. hal-04106033

HAL Id: hal-04106033

https://amu.hal.science/hal-04106033

Submitted on 25 May 2023

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Metastable oscillatory modes emerge from synchronization in the brain spacetime connectome

Joana Cabral ID 1,2,3,4,14✉, Francesca Castaldo ID 5,14, Jakub Vohryzek²,⁶, Vladimir Litvak ID ⁵, Christian Bick ID 7,8,9,10, Renaud Lambiotte ID ⁹, Karl Friston ID ⁵, Morten L. Kringelbach ID 1,2,3,11 & Gustavo Deco 6,12,13

A rich repertoire of oscillatory signals is detected from human brains with electro- and mag-
netoencephalography (EEG/MEG). However, the principles underwriting coherent oscillations
and their link with neural activity remain under debate. Here, we revisit the mechanistic
hypothesis that transient brain rhythms are a signature of metastable synchronization,
occurring at reduced collective frequencies due to delays between brain areas. We consider a
system of damped oscillators in the presence of background noise – approximating the short-
lived gamma-frequency oscillations generated within neuronal circuits – coupled according to
the diffusion weighted tractography between brain areas. Varying the global coupling strength
and conduction speed, we identify a critical regime where spatially and spectrally resolved
metastable oscillatory modes (MOMs) emerge at sub-gamma frequencies, approximating the
MEG power spectra from 89 healthy individuals at rest. Further, we demonstrate that the
frequency, duration, and scale of MOMs – as well as the frequency-specific envelope functional
connectivity – can be controlled by global parameters, while the connectome structure
remains unchanged. Grounded in the physics of delay-coupled oscillators, these numerical
analyses demonstrate how interactions between locally generated fast oscillations in the
connectome spacetime structure can lead to the emergence of collective brain rhythms
organized in space and time.

<footnote n="1">¹Life and Health Sciences Research Institute (ICVS), School of Medicine, University of Minho, Braga, Portugal. ²Centre for Eudaimonia and Human Flourishing, Linacre College, University of Oxford, Oxford, UK. ³Center for Music in the Brain, Department of Clinical Medicine, Aarhus University, Aarhus, Denmark. ⁴ICVS/3B's - Portuguese Government Associate Laboratory, Braga/Guimarães, Portugal. ⁵Wellcome Centre for Human Neuroimaging, University College London, Queen Square Institute of Neurology, London, UK. ⁶Center for Brain and Cognition, Computational Neuroscience Group, Universitat Pompeu Fabra, Barcelona, Spain. ⁷Department of Mathematics, Vrije Universiteit Amsterdam, Amsterdam, The Netherlands. ⁸Amsterdam Neuroscience - Systems & Network Neuroscience, Amsterdam, The Netherlands. ⁹Mathematical Institute, University of Oxford, Oxford, UK. ¹⁰Department of Mathematics, University of Exeter, Exeter, UK. ¹¹Department of Psychiatry, University of Oxford, Oxford, UK. ¹²Institució Catalana de la Recerca i Estudis Avançats (ICREA), Barcelona, Spain. ¹³Department of Neuropsychology, Max Planck Institute for Human Cognitive and Brain Sciences, Leipzig, Germany. ¹⁴These authors contributed equally: Joana Cabral, Francesca Castaldo. ✉email: joanacabral@med.uminho.pt</footnote>
The human brain is one of the most complex networks in nature, exhibiting a rich repertoire of activity patterns organized not only in space and time but also in the frequency domain. Indeed, rhythmicity is a central property of brain function—and perhaps of all biotic self-organization: from fast gamma activity in neurons to the life-cycle itself¹⁻⁴. Within the broad range of oscillations emerging at frequencies between 0.05 Hz and 500 Hz, the oscillations detected extracranially with electro- and magnetoencephalography (EEG/MEG) in resting humans typically peak between 0.5 and 30 Hz, being categorized as delta (~0.5–4 Hz), theta (~4–8 Hz), alpha (~8–13 Hz), and beta (~13–30 Hz)⁵. Notably, these oscillations lock in phase over long distances, generating metastable spatial topographies lasting up to a few hundred milliseconds⁶⁻⁸.

Falling significantly below the range of frequencies generated in local neuronal networks by feedback inhibition (>35 Hz, in the gamma-frequency range), it is generally agreed that sub-gamma oscillatory activity does not have a purely local origin and is associated with synchronization between distant neural assemblies⁹⁻¹⁴. Notably, there is a relation between the distance over which synchronization is observed and the frequency of the synchronized oscillations¹⁵⁻¹⁷. Specific brain circuitries, including among others the thalamocortical loop, have been proposed to play a role in the generation of rhythmic activity¹⁸⁻²⁰, which appears disrupted in neurological/neuropsychiatric disorders¹,²¹. Still, the fundamental mechanisms driving the spontaneous emergence of short-lived spatially and spectrally resolved oscillatory patterns remain unclear⁹,²²⁻²⁴.

Given the spatial distance and the finite propagation speed, interactions between brain areas are intrinsically time-delayed, which can manifest in network activity in the frequency domain. Indeed, delay-coupled limit-cycle oscillators have been demonstrated to synchronize at frequencies slower than the natural frequency of the oscillators, leading to a form of *collective* frequency emerging from synchronization mechanisms²⁵,²⁶. Briefly, when *N* phase oscillators—with natural frequency *ω*—are coupled together with a time delay *τ*, they synchronize at a delay- and interaction-dependent collective frequency *Ω* given by *Ω* = *ω* / (1 + *K*\**N*\**τ*), where *K* is the global coupling strength²⁵. However, this phenomenon has so far only been demonstrated for networks of limit-cycle oscillators²⁵, and it is unclear how it generalizes to systems where oscillations are not self-sustained, but instead emerge only transiently.

Computational models have proved helpful for demonstrating how the brain's complex network structure of long axonal projections connecting brain areas—the so-called structural

connectome²⁷—can shape brain activity in space and time²⁸⁻³⁶. Particularly, simulations of oscillatory units interacting in the connectome reveal a critical regime where different subsets of units temporarily synchronize and desynchronize, leading to transiently correlated activity across spatially segregated units²⁰,³⁰,³¹,³⁷. This reinforces the hypothesis that long-range functional connectivity between brain areas is driven by synchronization mechanisms²⁴,³⁸⁻⁴². Importantly, when considering realistic time delays in the Kuramoto model of coupled phase oscillators, periods of increased synchrony are accompanied by increased power at slower frequencies, generating spatially-organized band-limited power fluctuations similar to the ones captured with MEG²⁰. While these numerical results revealed the critical role of time delays to generate collective oscillations at reduced frequencies, it remains to be verified whether this phenomenon holds in the more realistic setting, wherein local oscillations have fluctuating amplitude—which is neglected in the Kuramoto model—as observed empirically in electrophysiological recordings of neural activity⁴³,⁴⁴. Furthermore, understanding the parameters that control the duration, size and occupancy of collective oscillations is crucial to inform the prediction of therapeutic strategies aimed at modulating dysfunctional oscillatory brain activity.

To address these fundamental questions, we build a phenomenological brain network model with realistic connectivity and time delays, where each node is described by a Stuart-Landau oscillator operating in the subcritical regime, i.e., responding to a stimulus with an oscillation with decaying amplitude³⁴,³⁵,⁴⁵. As the amplitude dynamics introduces an additional degree of complexity, it needs to be verified if the analytic predictions made for coupled limit-cycle oscillators⁴⁶ (valid for phase oscillators or supercritical Stuart-Landau oscillators) still hold⁴⁷. Selecting 40 Hz as a typical frequency of gamma oscillations, we set all units with identical natural frequency to exclude additional effects of frequency dispersion⁴⁸,⁴⁹, and perturb all units with uncorrelated white noise, considering that units resonate at their natural frequency in the presence of background noisy activity⁵⁰. Assuming the generalizability of collective synchronization frequencies to delay-coupled damped oscillators, we hypothesize to identify a critical range of global model parameters (global coupling and conduction speed) where metastable synchronization generates the transient emergence of sub-gamma collective oscillations, approximating features of human MEG recordings.

## Results

**Dynamical regimes of the brain network model.** The reduced brain network model comprises *N* = 90 dynamical units

<figure><img src="image_4.png" /><figcaption>**Fig. 1 Reducing the brain to a system of oscillators coupled in Connectome spacetime structure.** a The phenomenological brain network model consists in *N* = 90 nodes representing brain regions with links representing diffusion tracts between them. b A Stuart-Landau oscillator in the subcritical regime responding to perturbation (vertical arrow) with an oscillation with decaying amplitude. c In the supercritical regime, the Stuart-Landau oscillator enters a limit-cycle (with constant amplitude), approximating a phase oscillator.</figcaption></figure>
representing anatomically defined brain areas coupled according to a normative structural connectome of the human brain (see Methods—Structural Connectome) with reciprocal (i.e., bidirectional/symmetric) coupling $C_{NxN}$ and distance $D_{NxN}$ matrices (Fig. 1a). Each unit is described by a Stuart–Landau oscillator operating in the subcritical (underdamped) regime, such that when perturbed it decays to a fixed-point equilibrium with a damped oscillation at a natural frequency $\omega$ (Fig. 1b), in contrast with the supercritical regime, where the oscillations are in a limit cycle (Fig. 1c, see Methods and Supplementary Note 1).

To verify that novel frequencies emerge purely from delayed interactions, the natural frequency of all units is set at $\omega = 40$ Hz (representing the resonant frequency of isolated neural masses driven by feedback inhibition) and each unit is perturbed with uncorrelated white noise. The model—represented mathematically by a system of stochastic delay coupled differential equations—is solved numerically for two parameter ranges: the global coupling strength, $K$, which scales all pairwise connections, and the mean conduction delay, $\langle\tau\rangle$, which scales the time delays between units in proportion to the diffusion tract lengths (Fig. 2, see Methods for details).

The synchrony degree of the system, evaluated using the Kuramoto Order Parameter (KOP), is modulated by the global coupling strength $K$: for weak coupling, the synchrony is low, and all units exhibit oscillations close to the natural frequency $\omega$ (Fig. 2a). In the critical range between incoherence and full synchrony, periods of weakly stable synchronization drive slow fluctuations in the KOP (Fig. 2b). For sufficiently strong coupling, all units tend to synchronize at a global collective frequency $\Omega$, which, in the presence of time-delays, is distinct from the natural frequency $\omega$ (Fig. 2c).

Observing the levels of synchrony and metastability across the range of parameters explored (Fig. 2d, e), we find that the critical value of $K$ above which the system can synchronize increases logarithmically with the mean delay, in line with analytic predictions for coupled oscillators with heterogeneous delays⁵¹ (see Supplementary Note 2 and Supplementary Fig. 5). When synchronization occurs in the presence of delays, we observe a sharp decrease in the global peak frequency (Fig. 2f), closely approximating the analytic prediction given by $\Omega = \omega/(1+K*N^*(\tau))$ (Fig. 2g, see also Supplementary Note 2 and Supplementary Fig. 6).

These findings serve to verify that the phenomenon of synchronization at reduced collective frequencies is not restricted to coupled phase oscillators and generalizes to units in the subcritical regime, where damped oscillations emerge in response to perturbation (Supplementary Fig. 7). Further, it demonstrates that the peak frequency of synchronization can be predicted analytically from global variables such as the mean natural frequency $\omega$, the number of units $N$, the coupling strength $K$, and the mean delay $\langle\tau\rangle$. The robustness of this prediction to distributed natural frequencies is reported in Supplementary Fig. 8.

**Simulations reveal spectral features of human brain activity.** One characteristic feature of MEG (and EEG) signals from healthy humans at rest is the transient emergence of oscillations in the alpha frequency range (~8–13 Hz), resulting in a peak in the power spectrum whose prominence varies strongly across people (see Fig. 3a for the normalized power spectrum of MEG signals from 89 healthy young adults resting with eyes open from the Human connectome Project (HCP) open-source database; details in Methods section, individual power spectra reported in Supplementary Fig. 9).

We find that the brain network model approximates the average MEG power spectrum of awake resting subjects within

the critical region of high metastability where synchronization occurs at reduced collective frequencies (comparing Fig. 3c with Fig. 2e, f). In detail, for each pair of model parameters we calculate the squared Euclidean distance between the power spectrum of the simulated signals (Fig. 3b) and the MEG power spectrum averaged across all sensors and all subjects (Fig. 3a), revealing the greatest disparity when no delays are considered or if the global coupling is too weak (see Methods for details).

Given the observed (and well-established) variability between MEG power spectra across individuals (Fig. 3a), we investigate the extent to which this variability can be associated with changes in global model parameters, while keeping the structural connectivity unchanged. To do so, we identify the pair of model parameters that approximates the individual MEG power spectra of each of the 89 participants, falling in 29 pairs of parameters (white asterisks in Fig. 3c, see also Supplementary Fig. 10). Notably, this reveals a confined region in parameter space for a range of average delays $\langle\tau\rangle$ of 2–11 milliseconds, with slight changes in the coupling strength and conduction speed maximizing the fit to individual MEG power spectra, while the structural connectivity remains unchanged. These results do not exclude the role of individual variability in structural connectivity across subjects but reveal additional parameters that modulate a network’s frequency spectrum. This serves to demonstrate that the same connectome structure can support distinct activity patterns depending on global model parameters, with longer/ shorter time delays and stronger/weaker coupling inducing shifts in the peak frequency and modulating the distribution of power across the spectrum (Fig. 3b).

**Metastable oscillatory modes emerge from weakly stable cluster synchronization.** In the range of parameters where the model optimally approximates the power spectrum of MEG signals, fluctuations in the magnitude of the order parameter are driven by metastable cluster synchronization. In other words, when the coupling is strong, but not sufficiently strong to stabilize full synchronization, some subsets of units that are more strongly connected together (i.e., clusters/communities) can engage in partially-synchronized modes that remain stable for a short period in time. Given the presence of time delays, these clusters do not synchronize at the natural frequency of the individual units ($\omega = 40$ Hz), but instead synchronize at slower cluster-specific collective frequencies, leading to the emergence of metastable oscillatory modes (MOMs) at sub-gamma frequencies.

To detect the occurrence of MOMs and characterize them in space and time, we band-pass filter the simulated signals in four frequency bands (delta 0.5–4 Hz, theta 4–8 Hz, alpha 8–13 Hz, and beta 13–30 Hz). In Fig. 4, a colored shade is added when the amplitude in each frequency band exceeds 5 standard deviations of the amplitude in that same frequency range detected in simulations without delays (see Supplementary Note 3, Supplementary Fig. 11).

As shown in Fig. 4, we find that MOMs are structured both in space and in time. Specifically, the units synchronizing together exhibit the simultaneous emergence of an oscillation at the same collective frequency, leading to the vertical alignment of shaded areas, particularly visible for the alpha frequency range in Fig. 4a. Notably, for different sets of parameters, the configuration of Fig. 4a changes strongly. Indeed, while for very weak coupling, almost no supra-threshold oscillations are detected (Supplementary Fig. 12), for stronger coupling, globally synchronized supra-threshold oscillations emerge transiently in the delta band (Fig. 5a). For longer delays, oscillations are detected with a less definitive temporal alignment between brain areas (Supplementary Fig. 13).
<figure><img src="image_5.png" /><figcaption>d. Synchrony mean(KOP)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>b. Intermediate Global Coupling (K = 10)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>e. Metastability std(KOP)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>c. Strong Global Coupling (K = 50)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>f. Global Peak Frequency (Hz)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>g. Predicted Global Frequency (Hz)</figcaption></figure>

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

**Fig. 2** Collective oscillations emerge at reduced frequencies from time-delayed synchronization. The system of N = 90 coupled oscillators, Z, was simulated for 50 s in the presence of white noise, varying only two global parameters: the Global Coupling K (increasing exponentially to better capture the effect of delays) and the conduction speed, which scales the Mean Conduction Delay. **a–c** To illustrate the effect of the coupling strength in the frequency of synchronization, the collective signal given by $\sum_{n=1}^{N}Z_{n}$, with N = 90 is reported for three levels of global coupling, keeping the same mean conduction delay of 5 milliseconds. The corresponding power spectra are reported on the right of each plot, and the Kuramoto Order Parameter (KOP) is reported below. For weak coupling (**a**) the simulated signal exhibits oscillations peaking close to the node's natural frequency. For intermediate coupling (**b**), weakly stable synchronization generates transient oscillations at reduced frequencies. For strong coupling (**c**), global synchronization becomes more stable, and all units are entrained in a collective oscillation at a reduced frequency. For intermediate coupling, fluctuations in the order parameter are indicative of metastability. **d–g** For each simulation across the parameters explored, we report: (**d**) the mean of the KOP (referred to as Synchrony); (**e**) the standard deviation of the KOP (referred to as Metastability$^{87}$); (**f**) the peak frequency of the simulated collective signal; (**g**) the synchronization frequency predicted analytically, showing agreement with simulation results for sufficient synchrony.
<figure><img src="image_12.png" /><figcaption>**Fig. 3** Approximation of human magnetoencephalography (MEG) power spectra (PS) in a critical range of parameters. **a** MEG power spectra from 89 healthy young adults resting with eyes open from the open-source database of the Human connectome Project. The average power spectrum across individuals is reported in blue. **b** For each pair of parameters, the power spectra of the simulated signals (averaged across units and normalized between 0 and 80 Hz) is reported. **c** Squared Euclidean distance between the MEG power spectrum averaged across all sensors and subjects and the power spectrum of the simulated signals. Asterisks indicate the sets of parameters that optimally approximate the MEG power spectra of each of the 89 individuals (size scaled according to the number of subjects in each point).</figcaption></figure>

Furthermore, the power at sub-gamma frequencies is found to correlate strongly with the instantaneous phase synchronization evaluated by the KOP over time ($r = 0.7595$ and $r = 0.8247$ for Figs. 4b and 5b correspondingly). This demonstrates that the emergence of oscillations at sub-gamma frequencies in the simulations is modulated by fluctuations in the synchrony degree.

We further define quantitative metrics to characterize the MOMs emerging at different frequency bands for different sets of model parameters in terms of their duration (i.e., consecutive time that the power remains above threshold), their size (i.e., the number of units simultaneous displaying power above threshold) and occupancy (i.e., the proportion of time that the power is detected above threshold). As can be seen in Fig. 6, in the range of parameters where optimal fits to MEG data are obtained (Optimal Range), the alpha MOMs last longer, recruit more units and occur more often. Importantly, we demonstrate that global parameters, such as the coupling strength and the conduction speed, modulate the spatiotemporospectral properties of the whole system in a non-trivial way, while the dynamics at the local level and the underlying structural network remain unchanged.

The implicit sensitivity to global model parameters is illustrated in Fig. 7, where the emergence of supra-threshold oscillations in different frequency bands is represented in the brain at a single time point for five distinct sets of parameters. The evolution over time is shown in Supplementary Movie 1.

**Frequency-specific functional connectivity.** To link with studies of functional connectivity in MEG, we further investigate how the model parameters modulate the correlation between the amplitude envelopes across frequency bands. To do so, we band-pass filter the signals in each frequency band, extract the amplitude of the Hilbert transform and report the envelope correlation matrices in Fig. 8 for each frequency band and for four representative sets of model parameters. For weak coupling, the envelope correlations are close to zero (Pearson's correlation coefficient $cc < 0.1$ for all pairs of brain areas), indicating that the

coupling is insufficient to drive functional connections between brain areas. For global parameters in the optimal range (here $K=10$ and $\langle\tau\rangle = 3$ ms), different brain areas exhibit correlated envelopes, with stronger correlations (up to $cc = 0.78$) being detected in the alpha frequency range. In contrast, for strong coupling the functional connectivity in the alpha band is reduced (maximum pairwise correlation of $cc = 0.25$), while the envelopes of delta and theta oscillations are strongly correlated across the brain (up to $cc = 0.89$). Keeping the optimal range of global coupling, $K=10$, but increasing the delays to an average of $\langle\tau\rangle = 20$ ms, envelope functional connectivity is detected mostly in the delta frequency range. This illustrates how, given the same underlying spacetime network structure (i.e., the matrices of coupling weights $C$ and distances $D$), changes in global parameters strongly affect the envelope functional connectivity patterns at different frequency bands.

To illustrate the level of functional connectivity across the brain, next to each correlation matrix in Fig. 8, we represent each area as a sphere placed at its centre of gravity and colored according to the strongest correlation with any other brain area. This shows that, for the optimal range of parameters, the areas exhibiting the strongest functional connectivity in the alpha band are distributed mostly in posterior and dorsal cortical areas, aligning with empirical observations of stronger functional connectivity in the alpha band in the visual and somatomotor systems. However, it is important to consider that the specific spatial configuration of functional connections is inherently dependent on the resolution and topology of the structural connectome, which is known to depend on the parcellation scheme and on the brain parts (i.e., cortical, subcortical) considered. In Supplementary Methods 1, we perform the same analysis on data simulated using a structural connectome including 200 cortical-only brain areas⁵². These results show that, while the phenomenology of MOMs is robust to changes in the parcellation scheme, the spatial specificity across frequency bands is sensitive to the parcellation scheme considered (Supplementary Figs. 16 and 17). Most importantly, this analysis illustrates how
<figure><img src="image_13.png" /><figcaption>**Fig. 4** Sub-gamma oscillations emerge from weakly stable cluster synchronization. **a** An example of the simulated signals in all 90 units plotted over 25 s, each representing a brain area from a brain parcellation template, filtered below 30 Hz to highlight the sub-gamma oscillatory activity typically detected with magnetoencephalography (MEG). Shades indicate the time points of increased power in the delta (yellow), theta (red), alpha (blue) and beta (green) frequency bands. For each frequency band, the threshold was defined as five standard deviations (STD) of the amplitude—in the same frequency bands—when no delays were considered. For the simulations, the resonant frequency, $\omega_0$, of all units was set to 40 Hz, the conduction speed was tuned such that the average delay between units, $\langle\tau\rangle$, was 3 milliseconds (ms) and the global coupling strength was set to $K=10$. **b** The mean amplitude envelope (blue) of the filtered signals shown in (a) correlates with a Pearson's correlation coefficient $r = 0.7595$ with the phase synchronization evaluated by the Kuramoto Order Parameter (orange, right y-axis).</figcaption></figure>

frequency-specific functional connectivity patterns depend sensitively on global variables modulating the distributed dynamics, while the structural connectivity remains unchanged.

Discussion

This work addresses the physical mechanisms underlying brain rhythms detected empirically, employing a reductionist perspective to ground the inner complexity of encephalographic signals to universal theoretical principles$^{5,54}$. Approaching the problem from a macroscopic perspective, we focus on the emergent properties of interacting dynamical units, where the collective ensemble engages in functionally relevant activity patterns that cannot be inferred from the isolated units alone$^{55-58}$.

Specifically, we first demonstrate the generalizability of a synchronization mechanism described for networks of delay-coupled *limit-cycle* oscillators to networks of delay-coupled *damped* oscillators (i.e., in the subcritical range of a Hopf bifurcation). This is important for the neuroscience field, since empirical electrophysiological recordings show that local field oscillations in the gamma-frequency band are *not* limit-cycle oscillations (as considered in previous models using the Kuramoto of coupled

oscillators$^{20}$), but instead emerge only transiently. Therefore, the substantial reduction of brain areas to phase oscillators in Cabral et al. (2014) has raised concerns on the generalizability of the proposed mechanism to more realistic settings, given the demonstrated importance of considering the amplitude dynamics on the connectivity between phases$^{46,47,59,60}$.

Subsequently, we extend on previous brain network modelling works by demonstrating that the synchronization frequency can be approximated analytically from global model parameters, namely the number of units, the mean coupling strength, the average time delay between units, and the mean natural frequency of the units. Regarding the latter, we show that, in the presence of delays, the system is less sensitive to the spread of frequencies across units, in line with theoretical predictions$^{25}$ (Supplementary Fig. 8).

These insights are crucial to explain the macroscopic spatiotemporally organized oscillatory signals detected with EEG/MEG at sub-gamma frequencies, without explicitly introducing these oscillations in the model$^{61}$. Here, we consider that only gamma-frequency oscillations can be generated at the local neuronal level, with power at other frequencies resulting purely from synchronization with time delays. Furthermore, we demonstrate
<figure><img src="image_14.png" /><figcaption>**Fig. 5** Global delta waves emerge for strong coupling. **a** The simulated signals in all 90 units plotted over 25 s, each representing a brain area from a brain parcellation template, filtered below 30 Hz to focus on the sub-gamma oscillatory activity typically detected with magnetoencephalography (MEG). Shades highlight the time points of increased power in the delta (yellow), theta (red), alpha (blue) and beta (green) frequency bands. For each frequency band, the threshold was defined as five standard deviations (STD) from the amplitude in the same frequency bands when no delays were considered. These simulations were performed setting the resonant frequency of all units ω₀ = 40 Hz, the average delay between units, ⟨τ⟩ = 3 milliseconds (ms) and the global coupling strength was increased to K = 50 with respect to the simulations shown in Fig. 4. **b** The mean amplitude envelope (blue) of the filtered signals shown in (a) correlates with *r* = 0.8247 with the phase synchronization evaluated by the Kuramoto Order parameter (orange, right y-axis).</figcaption></figure>

the impact of global model parameters in the modulation of frequency-specific collective oscillations emerging across space and time. The detailed characterization of metastable oscillatory modes in terms of number of units synchronizing together, duration and occupancy provides a new framework to analyze collective brain oscillations complementary to frequency-specific envelope functional connectivity analysis.

Our hypothesis is endorsed using a phenomenological brain network model, reduced to its key essential ingredients to allow efficient numerical approximations to analytic predictions, but at the same time sufficiently complex to allow a fair approximation of MEG spectral features. The deliberate reductionist perspective inherent in this brain network model is intended to link with theoretical works on delay-coupled oscillatory systems<sup>25,51,62,63</sup>. Towards this end, we consider identical units with same natural frequency, same damping coefficient and same noise level, coupled in the structural connectome. Therefore, we focus solely on the effects of global variables, namely the global coupling *K* and the mean conduction delay ⟨τ⟩ in the emerging synchronization phenomena. To establish the construct validity of our numerical simulations, we show that the peak synchronization frequency can be approximated by the analytic prediction derived for

synthetic networks of coupled Kuramoto oscillators with time delays<sup>25</sup>. Further, in line with theoretical predictions<sup>51,64</sup>, we find that the complex spacetime topology of the structural connectome widely expands the critical border between incoherence and global synchrony where fluctuations in the order parameter are indicative of metastability<sup>25</sup>. Despite its simplicity, this model provides a robust framework to test a theoretically grounded mechanistic scenario for the spontaneous formation of frequency-specific long-range coherence in complex networks.

While the investigation of mechanistic principles and control parameters benefits from reduced complexity, adding heterogeneity is certainly needed to improve the fitting to real brain activity from individuals in different conditions. Building up on these fundamental aspects, additional degrees of complexity can be added to the model, namely by considering more fine-grained connectome structures, considering non-homogeneous intrinsic frequencies and damping properties, or even replacing the noisy input by dynamic concentration patterns to mimic local neuro-modulatory effects. Further, given the potential generalizability of this synchronization mechanism, we expect our analysis may provide valuable insight to interpret some of the complex self-organizing phenomena emerging in more realistic biophysical
<figure><img src="image_15.png" /><figcaption>**Fig. 6 Characterization of metastable oscillatory modes (MOMs) emerging from the system.** For different Global Coupling strength (K) and Conduction Delays $\langle\tau\rangle$, MOMs are characterized in terms of duration (i.e., consecutive time that the power remains above threshold), size (i.e., the number of units simultaneous displaying power above threshold) and occupancy (i.e., the proportion of time that the power is detected above threshold over the entire simulation), for each frequency band. This demonstrates that the same network structure, i.e., the connectome, can exhibit different oscillatory modes organized in space and in time, depending on global parameters of the system. In the critical range of parameters (Optimal Range), oscillations in the alpha frequency band emerge more frequently and involve more units. Globally synchronized delta oscillations—as typically observed in states of reduced consciousness—are associated to an increase in the global coupling strength (Strong Coupling). Error bars represent 1 standard deviation. See also Supplementary Movie 1.</figcaption></figure>

models of neural networks<sup>65,66</sup> for which a precise analytic prediction cannot be solved.

Our findings reinforce the idea that conduction delays—often neglected in network models of whole brain activity due to the added complexity—play a crucial role in shaping the frequency spectrum of coupled oscillatory systems. Although the frequency of the oscillations considered herein is relatively fast with respect to the ultra-slow fluctuations (<0.1 Hz) detected with functional Magnetic Resonance Imaging (fMRI), it is important to highlight that metastable synchronization drives power fluctuations on ultra-slow timescales, and therefore, even relatively short time delays can significantly modulate spontaneous activity at ultra-slow time-scales. We note that for the numerical integration of stochastic delay differential equations to be stable and align with analytic predictions, the time step for numerical integration needs to be sufficiently small and a running history needs to be saved for the length of the maximum delay between units, which significantly increases the computation times when compared to simulations where delays are neglected (here the numerical results were found to stabilize for $dt \le 10^{-4}$ s, see Supplementary Figs. 18–20).

The discovery of multistability in systems of delay-coupled oscillators, initially described in 1999 by Young and Strogatz<sup>62</sup> and extended to heterogeneous delays in 2009 by Lee, Ott and Antonsen<sup>51</sup>, was crucial to develop the theoretical hypothesis behind this work, opening grounds to speculate that this phenomenon may be related to the maintenance of the right balance

between integration and segregation in living brains<sup>67,68</sup>. Beyond the range where the model best approximates healthy awake brain activity, we find that higher coupling enhances global order, where the whole brain displays slow coherent oscillations in the delta-range (0.5–2 Hz), which nicely approximate the most powerful brain rhythms detected during unconscious states such as slow-wave sleep, coma or anaesthesia. On the other hand, operating at weaker coupling hinders the formation of MOMs at sub-gamma frequencies, altering the spectral profile similarly to what is observed in M/EEG recordings of patients with neuropsychiatric disorders associated to disconnection, such as schizophrenia, where the power in alpha appears to be significantly reduced<sup>69–72</sup>. Such abnormal interactions within cortico-subcortical oscillatory networks may emerge from specific local deregulation or neural circuit disruption<sup>73</sup>. However, how a local change may alter the communication between brain-areas and brain network dynamics remains an open question. Overall, these results are aligned with recent works proposing that spontaneous transitions between multiple space-time patterns on complex networks provide a solid theoretical framework for the interpretation of the non-stationary but recurrent macroscopic patterns emerging spontaneously in brain activity, and ultimately supporting brain function<sup>74,75</sup>. From a technical perspective, it may be surprising that this kind of itinerant dynamics emerges under symmetrical coupling between nodes; in the sense that asymmetric coupling is normally required for breaking detailed balance—and engendering
<figure><img src="image_16.png" /><figcaption>**Fig. 7** *Metastable Oscillatory Modes (MOMs) emerge transiently from interactions in the Connectome spacetime structure only for sufficient coupling and conduction times.* Each brain area is represented as a sphere located at its centre of gravity. A color code is used to highlight the brain areas with power exceeding five standard deviations from the baseline power at a given time point. This image is a still frame from Supplementary Movie 1. While the structural connectome is the same for all simulations, MOMs only emerge at reduced frequencies in the presence of Conduction Delays ($\langle\tau\rangle$) and for sufficient Coupling strength (K).</figcaption></figure>

stochastic chaos of the sort described above. However, the dynamics of each node are generated with asymmetric Jacobians, suggesting that symmetry breaking of intrinsic connectivity is a sufficient condition for the nonequilibrium dynamics that characterize real brains.

While metastability appears to be crucial for brain function, the specific role of MOMs to support cognitive functions remains unclear$^{4,24,68,76,77}$. One possibility is that the areas engaged in a MOM are directly involved in long-range functional integration, but another is that these areas are inhibited by entering in a collective low-energy mode$^{13,39}$. Shedding some light on this open question, we find that synchronization with delays induces not only a shift to slower frequencies but also a decrease in amplitude, in line with theoretical studies reporting amplitude death in systems with distributed delays$^{78}$ (see the vertical axes in Fig. 2a-c, top panels). From a ‘metabolic’ perspective, this shows that MOMs can be approached as ‘low-energy modes’ with respect to high power gamma oscillations, providing a physical explanation for the emergence of the so-called ‘idle rhythms’$^{79}$. Although the functional implications of this mechanism are beyond the scope of this work, we expect it will provide fertile grounds for the formulation of novel falsifiable predictions to be further tested. Moreover, these findings give room to further investigations of how local perturbations can affect the spatiotemporal dynamics on the macroscopic scale, to gain insight on the mechanisms of action of perturbative strategies such as transcranial magnetic stimulation or deep brain stimulation.

## Methods

**Ethics statement.** All human data used in this study is from the public repository of the Human Connectome Project (HCP)$^{80}$ (https://www.humanconnectome.org), which is distributed in compliance with international ethical guidelines.

**Structural connectome.** The NxN matrices of structural connectivity, C, and distances, D, used for the network model were derived from a probabilistic tractography-based normative connectome provided as part of the leadDBS toolbox (https://www.lead-dbs.org/)⁸¹. This normative connectome was generated from diffusion-weighted and T2-weighted Magnetic Resonance Imaging (MRI) from 32 healthy participants (mean age 31.5 years old ± 8.6, 14 females) from the HCP. The diffusion-weighted MRI data was recorded for 89 min on a specially-designed MRI scanner with more powerful gradients then conventional scanners. The dataset and the acquisition protocol details are available in the Image & Data Archive under the HCP project (https://ida.loni.usc.edu/). DSI Studio (http://dsi-studio.labsolver.org) was used to implement a generalized q-sampling imaging algorithm to the diffusion data. A white-matter mask, derived from the segmentation of the T2-weighted anatomical images, was used to co-register the images to the b0 image of the diffusion data using the SPM12 toolbox (https://www.fil.ion.ucl.ac.uk/spm/software/spm12/). Within the white-matter mask, 200,000 most probable fibres were sampled for each participant. Then, fibres were transformed to the standard Montreal Neurological Institute (MNI) space applying a nonlinear deformation field derived from the T2-weighted images via a diffeomorphic registration algorithm$^{82}$. The individual tractograms were then aggregated into a joint dataset in MNI standard space resulting in a normative tractogram representative of a healthy young adult population and made available in the leadDBS toolbox$^{81}$.

The NxN matrices were computed from the normative tractogram using the Automated Anatomical Labelling (AAL) parcellation scheme$^{83}$ with N = 90 cortical and subcortical areas, by calculating the number of tracts, C(n,p), and mean tract length, D(n,p), between the voxels belonging to each pair of brain areas n and p. Further details on the structural matrices in the AAL and other parcellation schemes are reported in Supplementary Methods 1 and Supplementary Fig. 14.
<figure><img src="image_17.png" /><figcaption>**Fig. 8** Influence of global model parameters in frequency-specific envelope functional connectivity patterns. For 4 simulations obtained with different Global Coupling strength (K) and Conduction Delays (τ), we report the frequency-specific functional connectivity (FC) estimated as the correlation matrices of the envelopes of signals band-pass filtered in the delta (δ), theta (θ), alpha (α) and beta (β) frequency bands. The colormap limits of the matrices are scaled by the maximum absolute correlation and centred at zero. Next to each matrix, each of the N = 90 brain areas is represented as a sphere placed in its centre of gravity and colored according to the maximum envelope FC to any other brain area (same colorbar applied to all spheres, scaled between -1 and 1).</figcaption></figure>

**MEG power spectra from healthy participants.** The power spectra from human resting-state MEG signals were also downloaded from the HCP database as a FieldTrip structure in a MATLAB file. The MEG power spectra are provided for 89 healthy participants at rest (mean age 28.7 years old, range 22–35, 41 female) distinct from the 32 participants from which the structural connectomes were derived, but with similar age range and gender ratio. Resting-state MEG signals were recorded on a Magnes 3600 MEG (4D NeuroImaging) with 248 magnetometers for 6 min and the “powavg” pipeline was used to obtain the power spectrum of the resting-state MEG data in each MEG sensor. Briefly, the signals were segmented, Hanning-tapered, Fourier-transformed and the power spectrum was averaged over all segments. Notch filters were applied to remove the power line noise (cut-off frequencies 59–61 Hz and 119–121 Hz). Additional details are explained in the HCP reference manual (https://humanconnectome.org/storage/app/media/documentation/s1200/HCP_S1200_Release_Reference_Manual.pdf). The MEG power spectra were averaged across the 248 sensors to obtain a power spectrum representative of each subject.

**Brain network model.** The Stuart–Landau (SL) equation (first term in Eq. 1) is the canonical form to describe the behavior of a system near an Andronov–Hopf bifurcation, i.e. exhibiting the birth of an oscillation from a fixed point$^{40,84}$. In other words, it is used to describe systems that have a static fixed point (like a resting spring), but respond to perturbation with an oscillation, which may damped or self-sustained depending on the operating point of the system with respect to the bifurcation (see Supplementary Note 1 and Supplementary Figs. 1–4). This model allows to describe complex-systems behavior among several applications, bridging the gap between the simplicity of the Kuramoto model and the extensiveness of the phase-amplitude frameworks$^{85,86}$. It describes how the oscillator behaves both when it is weakly attracted to a limit cycle (displaying only damped oscillations in response to perturbation) and, on the other hand, when it is purely restricted to a limit cycle (oscillations remain self-sustained).

Our analysis is based on a system of N = 90 SL oscillators coupled in the connectome, considering both the connectivity strength, $C_{np}$, and the conduction delays, $\tau_{np}$, between each pair of brain areas n and p. The conduction delays are defined in proportion to the fiber lengths between brain areas, assuming an homogenous conduction speed v, such that $\tau_{np} = D_{np}/v$, where $D_{np}$ is the length of the fibres detected between brain areas n and p. To simulate how the activity in

node n is affected by the behavior of all other nodes p ($p \in N \land p \neq n$), we describe the interaction between nodes in the form:

$$ \frac{dZ_n}{dt} = Z_n[a + i\omega - |Z_n^2|] + K \sum_{p \neq n}^{N} C_{np} [Z_p(t - \tau_{np}) - Z_n(t)] + \beta\eta_1 + i\beta\eta_2, \forall n \in N, \quad (1) $$

where the complex variable $Z_n(t)$ describes the state of the n-th oscillator at time t. The first term in Eq. (1) describes the intrinsic dynamics of each unit, the second term describes the input received from coupled units and the last terms represent uncorrelated white noise (see Supplementary Note 2 for detailed analysis of the model).

With this approach, we consider that the first term of Eq. 1 represents the natural excitability of neuronal assemblies, where $\omega = 2\pi * f_f$ is the angular frequency, with $f_f$ as the fundamental frequency in Hertz. For our proof of concept, we set all nodes with identical natural frequency $\omega_0 = 2\pi * 40Hz$, representing the undifferentiated ability of a neural mass to engage in gamma-frequency oscillations.

The parameter a determines the position of the each unit with respect to the limit cycle. For $a > 0$, a stable limit cycle appears via a supercritical Hopf bifurcation, while when $a < 0$ there is only a stable fixed point at the origin $Z_n = 0$, so the bifurcation point is at $a = 0$. Importantly, if a is sufficiently close to the bifurcation, the system is still weakly attracted to the limit cycle and damped oscillations emerge in response to external input, with a decay time scaled by a. In this work, we pick a value of $a = -5$ for all nodes, such that a single input drives a damped oscillation decaying after ~1 s, approximating the slowest decay time-constants of inhibitory receptors (Supplementary Fig. 4) (rdecay(GABA_B) ≈ 500–1000 ms). In Supplementary Note 2 and Supplementary Fig. 7 we show that our results are qualitatively similar for a broad range of a values, both positive and negative, thus demonstrating the generalizability of synchronization at collective frequencies to coupled oscillatory systems with fluctuating amplitude, be they damped or self-sustained. We note that this mechanism only fails when the units have an overdamped response (exponential decay without oscillation), which, in this case, only occurred for $a = -500$. Thus it is of great interest in future research to investigate whether the local bifurcation parameters can be tuned based on sensitive observables to fit the MEG data of different individuals in different conditions.

The second term represents the total input received from other brain areas, scaled by parameter K, which sets the strength of all network interactions with respect to the intrinsic node dynamics. Because we wish to focus on the nonlinear
phenomena introduced by time delays, we model the node-to-node interactions using a particular linear diffusive coupling, as the simplest approximation of the general coupling function, considering delayed interactions. Here, the signal of node n at time t is calculated with respect to the activity of all other nodes p at time t − τ$_{np}$ (where τ$_{np}$ is the time delay between n and p), scaled by the relative coupling strength given by C$_{np}$.

The third term of Eq. 1 represents the added uncorrelated noise to each unit (with real and imaginary components η₁ and η₂). In this analysis, the system is perturbed with uncorrelated white noise, where η₁ and η₂ are independently drawn from a Gaussian distribution with mean zero and standard deviation β = 0.001 (integrated as β*√dt).

In this framework, our whole-brain network model is purely bottom-up (i.e., not inferred from the MEG data we aim at explaining). For a qualitative comparison with the literature in delay-coupled oscillatory systems<sup>25,51,62</sup>, we explore the network dynamics as a function of the coupling strength K and the mean delay 〈τ〉 = 〈D〉/ν, where 〈D〉 is the mean length of all fibres detected between each pair of brain areas. For each set of parameters, the model is solved numerically for 50 s with an integration step dt = 10⁻⁴ s.

**Kuramoto order parameter.** To evaluate the global synchrony of the simulated network activity over time, we use the Kuramoto order parameter [Eq. 2]:

$$R(t)e^{i\theta(t)} = \frac{1}{N} \sum_{n=1}^{N} e^{i\theta_n(t)} \quad (2)$$

where θ$_{n}$ (t) is the phase of each node, given by the argument of Z$_{n}$. The temporal evolution of the Order Parameter R(t) provides an instantaneous measure of the degree of synchrony of the ensemble. Since we add noise in the simulations, we first band-pass filter the signals Z$_{n}$ around the peak frequency of the ensemble. A steady order parameter indicates a stable solution (be it asynchronous, 〈R(t)〉~0 or synchronous 〈R(t)〉~1), whereas fluctuations in the order parameter are indicative of Metastability, driven by constant transitions between different weakly stable solutions<sup>65,87</sup>. For the analysis in parameter space, we take the mean 〈R(t)〉 as a measure of the global synchronization while the standard deviation STD(R(t)) indicates how much R(t) fluctuates in time<sup>87</sup>.

**Analytic prediction of the collective frequency of synchronization.** Previous theoretical studies have shown analytically that coupled oscillatory networks with homogeneous delays can find stable solutions at multiple collective frequencies Ω. Let us consider the Kuramoto transition in a population of phase oscillators defined by Eq. 3:

$$\dot{\theta}_n(t) = \omega_0 + K \sum_{p=n}^{N} C_{np} \sin[\theta_p(t - \tau_{np}) - \theta_n(t)], \quad (3)$$

and the fully synchronized, uniformly rotating one-cluster state θ$_{n}$ = ... = θ$_{N}$ = Ωt. Substituting this expression into Eq. 3 we obtain<sup>25,26,62,78,88</sup>:

$$\Omega = \omega_0 - KN\sin(\Omega\tau) \quad (4)$$

where ω₀ corresponds to the nodes' intrinsic frequency and τ is the homogeneous time delay between nodes. As K is increased and full synchrony is approached, the system finds an equilibrium point at the lowest stable solution for Eq. 4, which is given by:

$$\Omega_{\min} = \omega_0/(1 + KN\tau) \quad (5)$$

Note that, for collective oscillations to emerge, the global coupling K needs to be sufficiently strong such that the synchronized solutions are at least weakly stable. To approximate the analytic prediction from Eq. (5), the coupling matrix was normalized by its mean, such that <C> = 1.

**Model performance.** We perform a parameter space exploration by tuning the two free parameters K and 〈τ〉. We choose to increase K exponentially as a power of 10 from 10⁻¹ to 10¹.⁷ in steps of 10⁰.¹, to ensure a range that covers from weak to strong coupling. 〈τ〉 is explored in the range from 0 ms to 30 ms in steps of 1 ms.

We measure the fitting between the empirical sensor MEG PS for each of the 89 subjects and the simulated PS for each pair of parameters as the squared Euclidean distance, resulting in one fitting value for each subject. This can be regarded as a maximum likelihood procedure under the assumption of Gaussian observation noise.

**Metastable oscillatory modes.** To detect MOMs and characterize them in space and time, we band-pass filter the simulated signals in each frequency band and obtain the corresponding amplitude envelopes using the Hilbert transform for each band. Subsequently, we consider that a node (or brain area) engages in a MOM if the amplitude increases five standard deviations above the baseline amplitude in that frequency range. We define the baseline threshold considering the simulations with the optimal K but with zero delays. Since some areas are more coupled together than others, even with “zero delays” these areas may exhibit more power across frequencies that is purely due to noisy interactions. Therefore, we define a different threshold for each node and each band.

**Envelope functional connectivity patterns.** Following standard procedures to estimate frequency-specific functional connectivity in empirical source-projected MEG data<sup>42</sup>, we first band-pass filter the simulated signals in each frequency band of interest, compute the analytic signal using the Hilbert transform and then calculate the correlation matrices between the amplitude (i.e. the absolute value) of the analytic signal. This is done in one optimal point (K = 10, 〈τ〉 = 3 ms), for weak coupling (K = 0.1, 〈τ〉 = 3 ms), strong coupling (K = 50, 〈τ〉 = 3 ms), no delays (K = 10, 〈τ〉 = 0 ms) and long delays (K = 10, 〈τ〉 = 20 ms). The same analysis performed using N = 200 units is shown in Supplementary Methods 1 and Supplementary Fig. 17.

## Data availability

Human neuroimaging data used in this study were provided by the Human Connectome Project (HCP)<sup>80</sup> (https://www.humanconnectome.org), WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. The normative connectomes were computed from Human Connectome Project data and included as part of the leadDBS toolbox<sup>81</sup> (https://www.lead-dbs.org/). The matrices computed from the normative connectomes used for simulations, together with the MEG power spectra from 89 individuals, are publicly available in .mat format at: https://github.com/fcast7/Hopf_Delay_Toolbox. Simulated data are available from the corresponding author on reasonable request. Supplementary Notes, Supplementary Methods and Supplementary Video 1 are available with this paper.

## Code availability

All simulations were performed in MATLAB2021b. The codes used in this study are publicly available at: https://github.com/fcast7/Hopf_Delay_Toolbox.

Received: 19 November 2021; Accepted: 20 June 2022; Published online: 15 July 2022

## References

1. Uhlhaas, P. J. & Singer, W. Neural synchrony in brain disorders: relevance for cognitive dysfunctions and pathophysiology. *Neuron* **52**, 155–168 (2006).

2. Engel, A. K., Fries, P. & Singer, W. Dynamic predictions: oscillations and synchrony in top-down processing. *Nat. Rev. Neurosci.* **2**, 704–716 (2001).

3. Buzsáki, G. & Draguhn, A. Neuronal oscillations in cortical networks. *Science* **304**, 1926–1929 (2004).

4. Singer, W. Synchronization of cortical activity and its putative role in information processing and learning. *Annu. Rev. Physiol.* **55**, 349–374 (1993).

5. Buzsáki, G. *Rhythms of the Brain*. (Oxford University Press, 2006).

6. Baker, A. P. et al. Fast transient networks in spontaneous human brain activity. *eLife* **3**, e01867 (2014).

7. Vidaurre, D. et al. Spectrally resolved fast transient brain states in electrophysiological data. *NeuroImage* **126**, 81–95 (2016).

8. Michel, C. M. & Koenig, T. EEG microstates as a tool for studying the temporal dynamics of whole-brain neuronal networks: a review. *NeuroImage* **180**, 577–593 (2018).

9. Jensen, O., Spaak, E. & Zumer, J. M. In *Magnetoencephalography* 359–403 (Springer, 2014).

10. Sherman, M. A. et al. Neural mechanisms of transient neocortical beta rhythms: Converging evidence from humans, computational modeling, monkeys, and mice. *Proc. Natl Acad. Sci.* **113**, E4885–E4894 (2016).

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

12. Hipp, J. F., Hawellek, D. J., Corbetta, M., Siegel, M. & Engel, A. K. Large-scale cortical correlation structure of spontaneous oscillatory activity. *Nat. Neurosci.* **15**, 884–890 (2012).

13. Palva, S. & Palva, J. M. Discovering oscillatory interaction networks with M/EEG: challenges and breakthroughs. *Trends Cogn. Sci.* **16**, 219–230 (2012).

14. Traub, R. D., Whittington, M. A., Stanford, I. M. & Jefferys, J. G. A mechanism for generation of long-range synchronous fast oscillations in the cortex. *Nature* **383**, 621–624 (1996).

15. Schnitzler, A. & Gross, J. Normal and pathological oscillatory communication in the brain. *Nat. Rev. Neurosci.* **6**, 285–296 (2005).

16. Von Stein, A. & Sarnthein, J. Different frequencies for different scales of cortical integration: from local gamma to long range alpha/theta synchronization. *Int. J. Psychophysiol.* **38**, 301–313 (2000).

17. Nunez, P. L. & Srinivasan, R. *Electric fields of the brain: the neurophysics of EEG*. (Oxford University Press, 2006).
18. Bhattacharya, S., Cauchois, M. B., Iglesias, P. A. & Chen, Z. S. The impact of a closed-loop thalamocortical model on the spatiotemporal dynamics of cortical and thalamic traveling waves. *Sci. Rep.* **11**, 1–19 (2021).

19. Freyer, F. et al. Biophysical mechanisms of multistability in resting-state cortical rhythms. *J. Neurosci.* **31**, 6353–6361 (2011).

20. Cabral, J. et al. Exploring mechanisms of spontaneous functional connectivity in MEG: How delayed network interactions lead to structured amplitude envelopes of band-pass filtered oscillations. *NeuroImage* **90**, 423–435 (2014).

21. Llinás, R. R., Ribary, U., Jeanmonod, D., Kronberg, E. & Mitra, P. P. Thalamocortical dysrhythmia: a neurological and neuropsychiatric syndrome characterized by magnetoencephalography. *Proc. Natl Acad. Sci.* **96**, 15222–15227 (1999).

22. Vidaurre, D. et al. Spontaneous cortical activity transiently organizes into frequency specific phase-coupling networks. *Nat. Commun.* **9**, 1–13 (2018).

23. O'Neill, G. C. et al. Dynamics of large-scale electrophysiological networks: a technical review. *NeuroImage* **180**, 559–576 (2018).

24. Friston, K. J. Transients, metastability, and neuronal dynamics. *NeuroImage* **5**, 164–171 (1997).

25. Niebur, E., Schuster, H. G. & Kammen, D. M. Collective frequencies and metastability in networks of limit-cycle oscillators with time delay. *Phys. Rev. Lett.* **67**, 2753–2756 (1991).

26. Atay, F. M., Jost, J. & Wende, A. Delays, connection topology, and synchronization of coupled chaotic maps. *Phys. Rev. Lett.* **92**, 144101 (2004).

27. Sporns, O., Tononi, G. & Kotter, R. The human connectome: a structural description of the human brain. *PLoS Comput. Biol.* **1**, e42 (2005).

28. Honey, C. J., Kotter, R., Breakspear, M. & Sporns, O. Network structure of cerebral cortex shapes functional connectivity on multiple time scales. *Proc. Natl Acad. Sci. USA* **104**, 10240–10245 (2007).

29. Ghosh, A., Rho, Y., McIntosh, A. R., Kotter, R. & Jirsa, V. K. Cortical network dynamics with time delays reveals functional connectivity in the resting brain. *Cogn. Neurodyn.* **2**, 115–120 (2008).

30. Deco, G., Jirsa, V., McIntosh, A. R., Sporns, O. & Kotter, R. Key role of coupling, delay, and noise in resting brain fluctuations. *Proc. Natl Acad. Sci. USA* **106**, 10302–10307 (2009).

31. Cabral, J., Hugues, E., Sporns, O. & Deco, G. Role of local network oscillations in resting-state functional connectivity. *NeuroImage* **57**, 130–139 (2011).

32. Cabral, J., Kringelbach, M. & Deco, G. Functional Connectivity dynamically evolves on multiple time-scales over a static Structural Connectome: Models and Mechanisms. *NeuroImage* **160**, 84–96 (2017).

33. Deco, G. & Kringelbach, M. L. Turbulent-like dynamics in the human brain. *Cell Rep.* **33**, 108471 (2020).

34. Deco, G. et al. Single or Multi-Frequency Generators in on-going brain activity: a mechanistic whole-brain model of empirical MEG data. *NeuroImage* **152**, 538–550 (2017).

35. 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).

36. Tewarie, P. et al. How do spatially distinct frequency specific MEG networks emerge from one underlying structural connectome? The role of the structural eigenmodes. *NeuroImage* **186**, 211–220 (2019).

37. Roberts, J. A. et al. Metastable brain waves. *Nat. Commun.* **10**, 1–17 (2019).

38. Ponce-Alvarez, A. et al. Resting-state temporal synchronization networks emerge from connectivity topology and heterogeneity. *PLoS Comput. Biol.* **11**, e1004100 (2015).

39. Engel, A. K., Gerloff, C., Hilgetag, C. C. & Nolte, G. Intrinsic coupling modes: multiscale interactions in ongoing brain activity. *Neuron* **80**, 867–886 (2013).

40. Cocchi, L., Gollo, L. L., Zalesky, A. & Breakspear, M. Criticality in the brain: a synthesis of neurobiology, models and cognition. *Prog. Neurobiol.* **158**, 132–152 (2017).

41. Deco, G. et al. Rare long-range cortical connections enhance human information processing. *Curr. Biol.* **31**, 4436–4448. e4435 (2021).

42. Brookes, M. J. et al. Measuring functional connectivity using MEG: methodology and comparison with fcMRI. *NeuroImage* **56**, 1082–1104 (2011).

43. Buhl, E. H., Tamas, G. & Fisahn, A. Cholinergic activation and tonic excitation induce persistent gamma oscillations in mouse somatosensory cortex in vitro. *J. Physiol.* **513**, 117–126 (1998).

44. Sanchez-Vives, M. V. & McCormick, D. A. Cellular and network mechanisms of rhythmic recurrent activity in neocortex. *Nat. Neurosci.* **3**, 1027 (2000).

45. Selivanov, A. A. et al. Adaptive synchronization in delay-coupled networks of Stuart-Landau oscillators. *Phys. Rev. E* **85**, 016201 (2012).

46. Daffertshofer, A. & van Wijk, B. C. On the Influence of Amplitude on the Connectivity between Phases. *Front. Neuroinform.* **5**, 6 (2011).

47. Ashwin, P., Coombes, S. & Nicks, R. Mathematical frameworks for oscillatory network dynamics in neuroscience. *J. Math. Neurosci.* **6**, 1–92 (2016).

48. Strogatz, S. H. & Mirollo, R. E. Stability of incoherence in a population of coupled oscillators. *J. Stat. Phys.* **63**, 613–635 (1991).

49. Petkoski, S., Iatsenko, D., Basnarkov, L. & Stefanovska, A. Mean-field and mean-ensemble frequencies of a system of coupled oscillators. *Phys. Rev. E* **87**, 032908 (2013).

50. Samanta, H. S., Bhattacharjee, J. K., Bhattacharyay, A. & Chakraborty, S. On noise induced poincaré-andronov-Hopf bifurcation. *Chaos: Interdiscip. J. Nonlinear Sci.* **24**, 043122 (2014).

51. Lee, W. S., Ott, E. & Antonsen, T. M. Large coupled oscillator systems with heterogeneous interaction delays. *Phys. Rev. Lett.* **103**, 044101 (2009).

52. Schaefer, A. et al. Local-Global Parcellation of the Human Cerebral Cortex from Intrinsic Functional Connectivity MRI. *Cereb. Cortex* **28**, 3095–3114 (2018).

53. Pikovsky, A., Kurths, J., Rosenblum, M. & Kurths, J. *Synchronization: a universal concept in nonlinear sciences*. (Cambridge university press, 2003).

54. Strogatz, S. *Sync: The emerging science of spontaneous order*. (Penguin UK, 2004).

55. Winfree, A. T. Biological rhythms and the behavior of populations of coupled oscillators. *J. Theor. Biol.* **16**, 15–42 (1967).

56. Haken, H. *Information and Self-Organization - A Macroscopic approach to Complex Systems*. (Springer, 1988).

57. Mirollo, R. E. & Strogatz, S. H. Synchronization of pulse-coupled biological oscillators. *SIAM J. Appl. Math.* **50**, 1645–1662 (1990).

58. Strogatz, S. H. & Stewart, I. Coupled oscillators and biological synchronization. *Sci. Am.* **269**, 102–109 (1993).

59. Daffertshofer, A., Ton, R., Kringelbach, M. L., Woolrich, M. & Deco, G. Distinct criticality of phase and amplitude dynamics in the resting brain. *NeuroImage* **180**, 442–447 (2018).

60. Siems, M. & Siegel, M. Dissociated neuronal phase-and amplitude-coupling patterns in the human brain. *NeuroImage* **209**, 116538 (2020).

61. Deco, G. et al. Single or multiple frequency generators in on-going brain activity: a mechanistic whole-brain model of empirical MEG data. *NeuroImage* **152**, 538–550 (2017).

62. Yeung, M. K. S. & Strogatz, S. H. Time Delay in the Kuramoto Model of Coupled Oscillators. *Phys. Rev. Lett.* **82**, 648–651 (1999).

63. Earl, M. G. & Strogatz, S. H. Synchronization in oscillator networks with delayed coupling: a stability criterion. *Phys. Rev. E* **67**, 036204 (2003).

64. Wildie, M. & Shanahan, M. Metastability and chimeras states in modular delay and pulse-coupled oscillator networks. *Chaos* **22**, 043131 (2012).

65. Bick, C., Goodfellow, M., Laing, C. R. & Martens, E. A. Understanding the dynamics of biological and neural oscillator networks through exact mean-field reductions: a review. *J. Math. Neurosci.* **10**, 1–43 (2020).

66. Izhikevich, E. M. & Edelman, G. M. Large-scale model of mammalian thalamocortical systems. *Proc. Natl Acad. Sci. USA* **105**, 3593–3598 (2008).

67. Tononi, G., Sporns, O. & Edelman, G. M. A measure for brain complexity: relating functional segregation and integration in the nervous system. *Proc. Natl Acad. Sci.* **91**, 5033–5037 (1994).

68. Tognoli, E. & Kelso, J. A. The metastable brain. *Neuron* **81**, 35–48 (2014).

69. Cabral, J., Kringelbach, M. L. & Deco, G. Functional graph alterations in schizophrenia: a result from a global anatomic decoupling. *Pharmacopsychiatry* **45**(Suppl 1), S57–S64 (2012).

70. Cabral, J., Hugues, E., Kringelbach, M. L. & Deco, G. Modeling the outcome of structural disconnection on resting-state functional connectivity. *NeuroImage* **62**, 1342–1353 (2012).

71. Schartner, M. M., Carhart-Harris, R. L., Barrett, A. B., Seth, A. K. & Muthukumaraswamy, S. D. Increased spontaneous MEG signal diversity for psychoactive doses of ketamine, LSD and psilocybin. *Sci. Rep.* **7**, 46421 (2017).

72. Carhart-Harris, R. L. et al. The entropic brain: a theory of conscious states informed by neuroimaging research with psychedelic drugs. *Front. Hum. Neurosci.* **8**, 20 (2014).

73. Goriely, A., Kuhl, E. & Bick, C. Neuronal oscillations on evolving networks: dynamics, damage, degradation, decline, dementia, and death. *Phys. Rev. Lett.* **125**, 128102 (2020).

74. Bick, C. Heteroclinic switching between chimeras. *Phys. Rev. E* **97**, 050201 (2018).

75. Ansmann, G., Lehnertz, K. & Feudel, U. Self-induced switchings between multiple space-time patterns on complex networks of excitable units. *Phys. Rev. X* **6**, 011030 (2016).

76. Rabinovich, M. I., Simmons, A. N. & Varona, P. Dynamical bridge between brain and mind. *Trends Cogn. Sci.* **19**, 453–461 (2015).

77. Roberts, J. A. et al. Metastable brain waves. *Nat. Commun.* **10**, 1056 (2019).

78. Atay, F. M. Distributed delays facilitate amplitude death of coupled oscillators. *Phys. Rev. Lett.* **91**, 094101 (2003).

79. Ritter, P., Moosmann, M. & Villringer, A. Rolandic alpha and beta EEG rhythms' strengths are inversely related to fMRI-BOLD signal in primary somatosensory and motor cortex. *Hum. Brain Mapp.* **30**, 1168–1187 (2009).

80. Van Essen, D. C. et al. The WU-Minn Human Connectome Project: an overview. *NeuroImage* **80**, 62–79 (2013).

81. Horn, A. et al. Lead-DBS v2: Towards a comprehensive pipeline for deep brain stimulation imaging. *NeuroImage* **184**, 293–316 (2019).
82. Ashburner, J. A fast diffeomorphic image registration algorithm. *NeuroImage* **38**, 95–113 (2007).

83. 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).

84. Andronov, A. A., Vitt, A. A. & Khakin, S. E. *Theory of oscillators*. (Dover Mathematics, 1987).

85. Aranson, I. S. & Kramer, L. The world of the complex Ginzburg-Landau equation. *Rev. Mod. Phys.* **74**, 99–143 (2002).

86. Olinger, D. J. A low-dimensional model for chaos in open fluid flows. *Phys. Fluids A: Fluid Dyn.* **5**, 1947–1951 (1993).

87. Shanahan, M. Metastable chimera states in community-structured oscillator networks. *Chaos* **20**, 013108 (2010).

88. Sethia, G. C., Sen, A. & Atay, F. M. Clustered chimera states in delay-coupled oscillator systems. *Phys. Rev. Lett.* **100**, 144102 (2008).

## Acknowledgements

J.C. is funded by the Portuguese Foundation for Science and Technology grants UIDB/50026/2020, UIDP/50026/2020 and CEECIND/03325/2017, Portugal. F.C. is funded by the EU-project euSNN European School of Network Neuroscience (MSCA-ITN-ETN H2020-860563). The Wellcome Centre for Human Neuroimaging is supported by core funding from Wellcome [203147/Z/16/Z]. J.V. is supported by the EU H2020 FET Proactive project Neurotwin grant agreement no. 101017716. R.L. acknowledges support from EPRSC Grants No. EP/V013068/1 and EP/V03474X/1. C.B. acknowledges support from the Engineering and Physical Sciences Research Council (EPSRC) through the grant EP/T013613/1. MLK is supported by the Center for Music in the Brain, funded by the Danish National Research Foundation (DNRF117), and the Centre for Eudaimonia and Human Flourishing, funded by the Pettit Foundation and Carlsberg Foundation. G.D. is supported by the Spanish national research project (PID2019-105772GB-I00 MCIU AEI), funded by the Spanish Ministry of Science, Innovation and Universities (MCIU), State Research Agency (AEI); HBP SGA3 Human Brain Project Specific grant agreement 3 (945539), funded by the EU H2020 FET Flagship program; SGR Research Support Group (reference 2017 SGR 1545), funded by the Catalan Agency for Management of University and Research Grants (AGAUR); Neurotwin Digital twins for model-driven non-invasive electrical brain stimulation (grant agreement 101017716), funded by the EU H2020 FET Proactive program; euSNN (grant agreement 860563), funded by the EU H2020 MSCA-ITN Innovative Training Networks; The Emerging Human Brain Cluster (CECH) (001-P-001682) within the framework of the European Research Development Fund Operational Program of Catalonia 2014–2020; Brain-Connects: Brain Connectivity during Stroke Recovery and Rehabilitation (201725.33), funded by the Fundacio La Marato TV3; and Corticity, FLAG-ERA JTC 2017 (reference PCI2018-092891), funded by the MCIU, AEI.

## Author contributions

J.C., G.D., and M.L.K. conceived and designed the analysis. J.C., F.C., J.V., G.D. contributed with data and analysis tools. J.C. and F.C. performed the simulations and analysis. V.L., K.F., R.L., C.B., M.L.K., and G.D. supervised the analysis. J.C. wrote the original draft. F.C. wrote the Supplementary Information. All authors reviewed and edited the final paper and Supplementary Information.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s42005-022-00950-y.

**Correspondence** and requests for materials should be addressed to Joana Cabral.

**Peer review information** *Communications Physics* thanks Fatihcan Atay and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022