Emergence of metastability in frustrated oscillatory networks: the key role of hierarchical modularity

Oscillatory complex networks in the metastable regime have been used to study the emergence of integrated and segregated activity in the brain, which are hypothesised to be fundamental for cognition. Yet, the parameters and the underlying mechanisms necessary to achieve the metastable regime are hard to identify, often relying on maximising the correlation with empirical functional connectivity dynamics. Here, we propose and show that the brain’s hierarchically modular mesoscale structure alone can give rise to robust metastable dynamics and (metastable) chimera states in the presence of phase frustration. We construct unweighted 3-layer hierarchical networks of identical Kuramoto-Sakaguchi oscillators, parameterized by the average degree of the network and a structural parameter determining the ratio of connections between and within blocks in the upper two layers. Together, these parameters affect the characteristic timescales of the system. Away from the critical synchronization point, we detect the emergence of metastable states in the lowest hierarchical layer coexisting with chimera and metastable states in the upper layers. Using the Laplacian renormalization group flow approach, we uncover two distinct pathways towards achieving the metastable regimes detected in these distinct layers. In the upper layers, we show how the symmetry-breaking states depend on the slow eigenmodes of the system. In the lowest layer instead, metastable dynamics can be achieved as the separation of timescales between layers reaches a critical threshold. Our results show an explicit relationship between metastability, chimera states, and the eigenmodes of the system, bridging the gap between harmonic based studies of empirical data and oscillatory models.

1 Introduction

The macroscopic dynamics of the human brain are characterized by a balance and flexible switching between integrated and segregated activity (Fox and Raichle, 2007). The coexistence of both segregative and integrative tendencies is thought to be essential for healthy brain functioning and cognitive processing (Tognoli and Kelso, 2014; López-González et al., 2021; Wang et al., 2021; Capouskova et al., 2023). In coordination dynamics (Kelso, 2013), this coexistence is the defining feature of complex dynamical systems in the metastable regime (Tognoli and Kelso, 2014; Hancock et al., 2023a). In the last 2 decades, various indexes of metastability have been proposed to quantify these dynamical behaviours from empirical data (Hancock et al., 2023a). Recently, measures of metastability have been successful in characterizing changes in the brain’s metastable-like dynamical features in the presence of pathological or pharmacological alterations (Lee et al., 2018; Lord et al., 2019; Hancock et al., 2023b), brain injury (Hellyer et al., 2015), the aging brain (Cabral et al., 2017) [as well as the aging rat brain (Alteriis et al., 2024)], and after brain stimulation (Bapat et al., 2024).

Motivated by these empirical observations, numerous whole-brain dynamical models have been employed to reproduce the metastable features of macroscopic brain dynamics (Cabral et al., 2014; Hansen et al., 2015). A whole-brain model consists of a complex dynamical system in which each of the constituent elements interacts with other elements according to the structural connectivity of real brains. The most frequently used indexes of metastability are the variation of the Kuramoto Order Parameter (KOP) (Shanahan, 2010) and the variation of the leading eigenvector of the instantaneous phase locking matrix and its eigenvalue (Cabral et al., 2017; Alteriis et al., 2024) (see review (Hancock et al., 2023a) for an in-depth discussion). Whole-brain oscillatory based models, in particular, allow for an extensive exploration of the structural and dynamical parameter space, enabling researchers to identify key parameters that dictate the models’ behaviours in the metastable regime (Cabral et al., 2011; Deco et al., 2017; Torres et al., 2024). A large number of key control parameters have been suggested to play a vital role in the emergence of metastable-like features. For instance, heterogeneity in couplings, delays, and connectivity has been shown to be crucial for the emergence of metastable states that synchronize at frequencies lower than the average oscillator frequency (Cabral et al., 2022). Turning to the role of structure, graph theoretical properties of connectomes have been shown to alter indexes of metastability (Váša et al., 2015) and similar measures of functional complexity (Zamora-López et al., 2016).

In these models, indexes of metastability have been hypothesised to peak when the correlation with empirical dynamical functional connectivity (dFC) is maximised (Deco and Kringelbach, 2016) [however, see also (Pope et al., 2023)]. Thus, most studies primarily focus on the properties of a system in the metastable regime and how the system’s parameters change these observations. However, the criteria that should be satisfied by an oscillatory system to achieve the metastable regime, and the mechanisms underpinning metastability, are not fully understood. As suggested in (Cabral et al., 2022), understanding the mechanisms that give rise to the emergence of metastable states is essential to be able to make predictions about the appearance and duration of specific metastable modes, and eventually to design possible therapeutic interventions.

In dynamical systems theory, metastable-like dynamics have been observed in various systems of oscillators. Here, metastability is usually associated with instability of chimera states: symmetry breaking states characterised by the coexistence of coherent and incoherent patterns in systems of phase-lagged identical oscillators (Abrams and Strogatz, 2004; Panaggio and Abrams, 2015; Haugland, 2021). For instance, in (Shanahan, 2010), metastable chimera states arise as a result of winnerless competition between modules of identical oscillators to join the coalition of synchronized modules. Similar metastable behaviours can be found in (Bick, 2018) due to the addition of higher-order interactions. Chaos, turbulence, and other dynamics with metastable-like features have also been shown to arise in the case of heterogeneous couplings and phase lags in the two-population model (Bick et al., 2018). Thus, these studies suggest that the interplay between modular structures, coupling heterogeneities, and phase frustration is necessary for metastable chimera states to arise. The effects of non-local hierarchical topologies on chimera states have also been investigated. In these settings, the structural symmetries and the clustering coefficient have been found to be important parameters for the emergence of chimera states in networks of Van der Pol oscillators (Ulonska et al., 2016). Closer to the approach taken in this study, authors in (Makarov et al., 2019) analyzed a system of Kuramoto-Sakaguchi oscillators composed of subnetworks connected via hub nodes placed on a ring. The authors showed how the size and connectivity of the subnetworks promote a competition mechanism between scales which results in an expansion of the range of parameters within which chimera states arise when compared to the case without subnetworks. Hence, hierarchical network properties seem to promote the emergence of chimera states, but their role in achieving the metastable regime and the specific mechanisms enabled by hierarchical structures remain unclear.

Some insights have come from computational neuroscience studies that explored in more detail effects of the hierarchically modular structure of the brain (Meunier et al., 2010), such as Griffith phases (Rubinov et al., 2011; Moretti and Muñoz, 2013; Ódor et al., 2015). In particular, recent research suggested that including the brain’s mesoscale connectivity in whole-brain oscillatory models can help us identify and understand the mechanisms underlying metastable dynamics. In (Mackay et al., 2023), the authors generated spatially constrained random networks to approximate the mesoscale neocortex, allowing the implementation of heterogeneous phase delays in a model of non-identical Kuramoto oscillators. The authors observed how the topology of these constrained networks widens the range of critical couplings within which metastable behaviours arise, analogous to Griffith phases. In (Villegas et al., 2014) instead, the authors studied in detail real and synthetic hierarchically modular brain networks of oscillators in the absence of phase frustration. Whilst it is well known that the hierarchical network structure affects the dynamics of oscillatory systems in the path towards global synchronization (Arenas et al., 2006; Arenas and Díaz-Guilera, 2007; Villegas et al., 2022), in (Villegas et al., 2014) the authors observed that modules of identical units at the bottom of the hierarchy synchronize first (fast timescale dynamics). However, synchronization can be lost as interactions with other modules become more significant at longer timescales, giving rise to transient metastable states on the path towards synchronization. Thus, these observations highlight how the slow dynamics associated with the lowest eigenvalues of the graph Laplacian affect the faster timescale dynamics that dictate the oscillators’ behaviour within lower-layer modules. This study, however, did not consider the effect of phase frustrations, such as phase lags or phase delays, and whether robust metastable and chimera states can persist, preventing global synchronization even in the case of identical oscillators.

Taken together, these studies suggest that modular or hierarchically modular structures, combined with either higher-order (non-additive) interactions, or the heterogeneity of both couplings and delays, are crucial to achieve the metastable regime [but also see the case for nonlinear couplings (Haugland et al., 2015; Haugland, 2021)]. In this work, we investigate the more parsimonious hypothesis that the brain’s hierarchically modular mesoscale structure alone can give rise to robust chimera states and modulate metastable dynamics in the presence of phase-frustration. Under the assumption that dFC patterns in real brains are supported by the brain’s mesoscale structural connectivity, our hypothesis is motivated by the fact that dFC patterns can be constructed using the spectral information of structural brain networks (Atasoy et al., 2016; Atasoy et al., 2017). Inspired by the network structures used in (Villegas et al., 2014; Zamora-López et al., 2016; Kang et al., 2019; Fuscà et al., 2023), we construct 3-layer hierarchical networks using a variation of the nested Stochastic Block Model (nSBM) (Peixoto, 2014). Our variation of the nSBM is parametrized by two quantities: the average degree of the network, and a structural parameter controlling the ratio of connections between and within blocks at high hierarchical layers while keeping the average degree of the network conserved and homogeneous. In the absence of coupling and oscillator frequency heterogeneities, we show how the Kuramoto-Sakaguchi dynamics on these networks display robust metastable and chimera states at different hierarchical layers depending on the mesoscale structure of the network. Through explaining the mechanisms underpinning the emergence of metastable and chimera states we elucidate an explicit relationship between these states and the eigenmodes of the graph Laplacian. In particular, we show how the system we investigated may be reduced to the classic two-population model for a specific range of the structural parameter via the Laplacian Renormalization Group (LRG) flow. Our results suggest that, while the slow eigenmodes determine the functional organization in higher coarse-grained layers, the spectral gap between fast and slow eigenmodes affects the stability of cluster synchronization in the lower fine-grained layers. We conclude by pointing towards possible future extensions of this work to further bridge the gap between harmonic based studies of empirical dFC and oscillatory based whole-brain models.

2 Methods

2.1 Hierarchical network model choice

To test our hypothesis, we seek to construct hierarchically modular networks while avoiding the presence of structural heterogeneities such as network motifs, hubs and the rich club, as well as non-homogeneous degree distributions, all of which have already been studied in (Ulonska et al., 2016; Zamora-López et al., 2016; Krishnagopal et al., 2017). Hence, we seek to construct hierarchical networks with a single control parameter that smoothly varies the mesoscale organization of the network while maintaining the average degree of the network fixed and homogeneous across all vertices. Critical mesoscale structures, which, in the context of information-diffusion and synchronization dynamics, are identified as structures that emerge at a characteristic timescale in the dynamical process, are an important feature of hierarchically modular networks such as the brain (Villegas et al., 2022). Such structures naturally arise in community-structured networks such as the Stochastic Block Model (SBM) and its nested variations (nSBM) used in various computational neuroscience studies such as (Villegas et al., 2014; Kang et al., 2019; Villegas et al., 2022; Fuscà et al., 2023). Therefore, nSBMs are a good starting point for our study. In the next subsection, we introduce a single-parameter nSBM variation, allowing us to smoothly vary the system from a SBM to a nSBM for a given desired average degree.

2.2 Variation of the nSBM

Given an arbitrary partition $P$ of $N=|\mathcal{V}|$ nodes, where $\mathcal{V}$ is the set of nodes, a SBM is defined by assigning a probability of connection between all pairs of nodes in terms of $P$. Starting from an initial node partition into $B_1$ blocks (layer $l=1$), we obtain a $L$-layer nSBM (Peixoto, 2014) by defining a new SBM in each layer $l^{\prime }>l$ in terms of the partition at layer $l^{\prime }-1$ (Figure 1A). The last layer $L$ is defined as a single block for convenience. This process results in a partition matrix $P\in {\mathbb{N}}^{L\times N}$ (Figure 1B). Each column of the partition matrix $P$ assigns a node to a specific block in each layer $l$, where each row of $P$ corresponds to the partition in layer $l$. Thus, the system is parametrized by assigning a probability of connection to all $N(N-1)$ possible connections in terms of the partition matrix of dimension $L\times N$. We may reduce the dimension of the $[L\times N]$-dimensional parameter space by choosing specific probability of connection rules that depend on the nodes’ block membership in each layer. In our variation of the nSBM, see Figure 1C, we introduce the following constraints:

$•$ each block in layer $l$ contains the same number $n_l$ of layer $l-1$ blocks;

$•$ the probability of connection between two nodes in a block in layer $l$ is the same for all $B_l$ blocks in the same layer;

$•$ in layer $L$ there is only one block containing $n_L=2$ blocks from layer $L-1$.

The first two constraints ensure structural homogeneity at layer level for all layers, while the third constraint was chosen to model the brain’s two-hemisphere configuration, similar to previous studies (Villegas et al., 2014; Kang et al., 2019). We assign non-zero probabilities only for edges between nodes in the same block of any layer. Then, since layer $L$ is a single block, only $L$ probabilities are required given these constraints. For the case $L=3$, all possible edges are assigned a probability $p_1,p_2$, or $p_3$ as depicted in Figure 1C. However, we wish these probabilities to depend on a single control parameter $H\in [\mathrm{0,1}]$ such that $p_i\equiv p_i(H),i=\mathrm{1,2,3}$. To do this, we add three more requirements:

$•$ for communities in layer 1 to be almost fully connected, i.e., $p_1$ must be as close to unity as possible;

$•$ for layers 2 and 3 to be identical when $H=0$. Then, the nSBM reduces to the classical SBM with only intra- and inter-community connection probabilities, i.e., $p_2(0)=p_3(0)$, with $p_1(H)$ reducing to the intra-community connection probability and $p_2=p_3$ reducing to the inter-community connection probability;

$•$ for communities in layer 2 to be completely decoupled for $H=1$, such that $p_3(H=1)=0$.

These constraints allow to smoothly vary the system configuration from a 2-layer community structured network to a 3-layer hierarchically modular network while keeping blocks in layer 1 always almost fully connected. Finally, we obtain expressions for $p_1(H),p_2(H)$ and $p_3(H)$ by requiring the average degree of the network to remain constant with respect to varying $H$ (see derivation in Supplementary Appendix SA):

$p_1\left(H\right)=1-\left(\frac{1-H}{2}\right)\left(\frac{n_1\gamma }{n_1-1}\right),(1)$

$p_2\left(H\right)=\left(\frac{1+H}{2}\right)\gamma ,(2)$

$p_3\left(H\right)=\left(\frac{1-H}{2}\right)\gamma ,(3)$

where the factor $\gamma \in [\mathrm{0,1}]$ in the expressions for $p_2$ and $p_3$ can be set to control the average degree of the network. Given these constraints, the theoretical average degree of the network can be written as

$⟨k⟩=\left(n_1-1\right)p_1\left(H\right)+n_1\left(n_2-1\right)p_2\left(H\right)+n_2{n}_1\left(n_3-1\right)p_3\left(H\right),(4)$

allowing us to find $\gamma$ given a desired average degree $k$. Intuitively, we may see Equation 4 as the sum of the connection probabilities of any single row in Figure 1B. Additionally, we may also note how this equation shows the structural homogeneity in each layer; each node (block) has the same average degree as any other node (block) in the same layer. After plugging in the expressions for $p_i$, Equation 1–3, we obtain the following expression for the average degree of the network (see Supplementary Appendix SA):

$⟨k⟩=n_1+\gamma \left(n_1{n}_2-n_1\right)-1,\gamma \in \left[\mathrm{0,1}\right].$

Hence our model allows to choose an average degree in the range between $n_1-1$ and $n_1{n}_2-1$, i.e., the maximum number of connections in a layer 1 module and the maximum number of connections in a layer 2 module. This limitation is a result of the imposed model requirements. Note that the total number of nodes is determined by $n_l$, with $N=n_1{n}_2{n}_3$. In what follows, we used $n_3=2$, $B_3=1$, and $B_2=2$, to model the two-hemisphere division of the brain and to allow comparison of our results with the 2-population model analytical results (Abrams et al., 2008). Finally, following the model’s constraints, we also have that $B_1=2n_2$. Therefore, given the size of blocks in layers 1 and 2, specified by $n_1$ and $n_2$, our variation of the nSBM is fully parametrized by $k$ and $H$, which control the average degree and the density of connections between and within blocks in higher hierarchical layers of the resulting network, respectively. In Figure 2, we show three examples of adjacency matrices constructed using our variation of the nested SBM. The parameter $H$ changes the density of connections within and between the two populations of nodes, while preserving the average degree and the node’s degree distribution of the resulting networks. For easier understanding, in the remainder of the manuscript, we will refer to blocks of nodes in layer 2 as populations, and blocks of nodes in layer 1 as modules. Then, the subset of nodes belonging to a population ${\rho }_i,i=\mathrm{1,2}$ are formally defined as ${\rho }_i=\{v|P_{2v}=i,v\in \mathcal{V}\}$ where $P$ is the partition matrix and $\mathcal{V}$ is the set of all nodes. Instead, the subsets of nodes belonging to a module are defined as ${\mu }_i=\{v|P_{1v}=i,v\in \mathcal{V}\}$ and a module $i$ is said to belong to population ${\rho }_j$ if ${\mu }_i\subset {\rho }_j$. To ease comparisons with other studies and real brain structural connectivities, we also computed the modularity index

$Q=\frac{1}{2|\mathcal{E}|}{\displaystyle \sum _c}\left[e_c-{\left(\frac{{\sum }_{i\in c}{k}_i}{2|\mathcal{E}|}\right)}^2\right],$

where each $c$ in the summation denotes a community given a partition $P^{\prime }$, $e_c$ is the number of edges in $c$, $k_i$ is the degree of a node $i\in c$, and $|\mathcal{E}|$ is the number of edges of a network $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with $\mathcal{V}$ and $\mathcal{E}\in \mathcal{V}\times \mathcal{V}$ being the set of nodes and edges, respectively. In Figures 2D, E, we report the modularity index $Q_l$ of the networks generated with our variation of the nSBM using the partition of nodes in layers $l=\mathrm{1,2}$.

Figure 1

Figure 1. Example of a 3-layer hierarchical network. (A) At each layer $l$, we define a SBM with $B_l$ blocks containing $n_l$ blocks from the layer below (or nodes if $l=1$). The last layer is fixed to be a single block, $B_3=1$, containing $n_3=2$ blocks from layer 2. In this example in layer 1 (blue blocks ${\mu }_i$) we have $B_1=8,n_1=2$, in layer 2 (green blocks ${\rho }_i$) we have $B_2=2,n_2=4$, in layer 3 (red block) we have $B_3=1,n_3=2$. (B) Corresponding partition matrix. (C) Using our constraints, the 3-layer hierarchical network is fully parametrized by the connection probabilities $p_i,i=\mathrm{1,2,3}$ (blue, green, red).

Figure 2

Figure 2. (top) Example adjacency matrices for $H=0.0$ (A) $H=0.5$ (B) and $H=1$ (C). Each network was generated using $n_1=16,n_2=8$ and exemplifying average degree $k=51.2$. Note how for $H=0.0$ the generated network is similar to a SBM while for $H=1.0$ the two populations are decoupled. (bottom) Modularity index $Q_l$ calculated over partitions $\{P_{l,i}|i\in \mathcal{V}\}$ for $l=1$ (D) and $l=2$ (E) as a function of $H$ and $k$.

2.3 Dynamical model and measures

We consider the underlying network to be static, and associate each node to a Kuramoto-Sakaguchi oscillator. The phase ${\theta }_i$ of each oscillator $i$ is governed by the following equation of motion:

$\frac{\mathrm{d}{\theta }_i}{\mathrm{d}t}={\omega }_i-K{\displaystyle \sum _j^N}{A}_{ij}\sin \left({\theta }_j-{\theta }_i-{\alpha }_{ij}\right)(5)$

where $A={\{A_{ij}\}}_{i,j\in \mathcal{V}}$ is the adjacency matrix, ${\omega }_i$ is the natural frequency of oscillator $i$, ${\alpha }_{ij}$ is the phase lag between oscillators $i$ and $j$, and $K$ is the coupling constant. As in the classical models for chimera states and metastability in oscillatory networks (Abrams et al., 2008; Shanahan, 2010; Panaggio and Abrams, 2015), we used identical oscillators by choosing ${\omega }_i=\omega =1$. We also fixed the phase lag to be equal to zero if $i,j$ are in the same module at layer $l=1$, and ${\alpha }_{ij}=\alpha$ otherwise. The zero phase lag for within-module connections at the lowest hierarchical layer was chosen as an approximation of shorter phase delays between oscillators in the same module of real networks embedded in a metric space. As standard practice in the dynamical systems’ literature (Panaggio and Abrams, 2015), we parametrized the phase lag $\alpha$ using the lag parameter $\beta =\pi /2-\alpha$. In this manuscript, we report results for $\beta =0.1$ which exemplifies our results. Since the oscillators are identical (i.e., they have the same natural frequency), the coupling constant effectively acts as a time re-scaling factor, leaving results qualitatively invariant. Here, numerical simulations were performed using a coupling constant $K=50/k$ normalized by the average degree of the network, allowing us to compare networks of the same size for varying degree. At time $t=0$, oscillator phases ${\theta }_i$ were randomized in $U(-\pi ,\pi )$. All simulations were performed using the Euler method with steps $\mathrm{\Delta }t$ of size 0.001 for a total of 55,000 steps with 5,000 steps of relaxation. We confirmed that this step size was appropriate by verifying that results were unchanged when using smaller Euler step sizes (results not shown).

2.3.1 Kuramoto Order Parameter

To quantify the degree of synchrony of the whole system we used the Kuramoto Order Parameter (KOP)

$Z={\mathrm{R}\mathrm{e}}^{i\psi }=\frac{1}{N}{\displaystyle \sum _j^N}\exp \left(i{\theta }_j\right),(6)$

where the real part $R\in [\mathrm{0,1}]$ quantifies the degree of synchronization, with $R\sim 1/\sqrt{N}$ suggesting incoherent motion and $R\sim 1$ full synchronization.

2.3.2 Local Kuramoto Order Parameter

In a similar fashion, we defined the local KOP for a subset $\mu \subset N$ of oscillators as

$Z_{\mu }=R_{\mu }{e}^{i{\psi }_{\mu }}=\frac{1}{|\mu |}{\displaystyle \sum _{j\in \mu }}\exp \left(i{\theta }_j\right).(7)$

Similarly to Equation 6, the real part of this measure, $R_{\mu }\in [\mathrm{0,1}]$, quantifies the degree of coherence of the subset of oscillators $\mu$ in isolation.

2.3.3 Measures for metastability

Oscillatory systems are said to be in the metastable regime when they exhibit spontaneous flexible switching between segregative (incoherence) and integrative (synchronization) tendencies¹ (Kelso, 2013). To quantify this switching behaviour, here we used the metastability index ${\sigma }_{\text{met}}$, introduced in (Shanahan, 2010). In its original form, ${\sigma }_{\text{met}}$ is obtained by computing the average standard deviation of the local KOPs (Equation 7). Since the 3-layer hierarchical model is equipped with a partition matrix defining distinct blocks in each layer, the index of metastability can be calculated in each layer $l$. Consequently, here we define a distinct index ${\sigma }_{\text{met}}^l$ for each layer $l$. For instance, consider the set of blocks in layer 1, then the degree of metastability of the layer 1 blocks is ${\sigma }_{\text{met}}^1={⟨\sigma (R_{{\mu }_i}(t))⟩}_{i\in [1,B_1]}$, where $B_1$ is the number of blocks in the first layer, ${\mu }_i$ with $i\in [1,B_1]$ identifies the set of oscillators that make up a block $i$ in the first layer, and $\sigma$ is the standard deviation. Similarly, we define ${\sigma }_{\text{met}}^2$ for layer 2 blocks (i.e., the populations ${\rho }_1$ and ${\rho }_2$). For layer 3 instead, the metastability index is just the variation of the whole system KOP since all oscillators belong to the same block, i.e., ${\sigma }_{\text{met}}^3=\sigma (R)$.

2.3.4 Measures for chimera states

The term chimera state is generally invoked to denote the coexistence of coherent and incoherent patterns in systems of identical oscillators. Over the past decade, various measures have been defined to characterize different types of chimera states (Kemeth et al., 2016; Haugland, 2021). In the case of identical sinusoidally-coupled oscillators with fixed amplitude, the local KOPs can be used to measure the degree of phase-coherence (Panaggio and Abrams, 2015) between the system’s parts. However, since an analytical characterization of chimera states is only possible in the limit of an infinite number of oscillators, an arbitrary threshold needs to be introduced, ${\delta }_1$, to numerically assess the presence of chimera states (Kemeth et al., 2016). For instance, if a system self-organizes into two subsets ${\rho }_1,{\rho }_2$ of oscillators with different degrees of synchronization at some time $t$, i.e., $d(t)=|R_{{\rho }_1}(t)-R_{{\rho }_2}(t)|>{\delta }_1$, we conclude that system is in a chimera state. Here, due to the similarity with the classic two-population model (Abrams et al., 2008), we focus on chimera states in layer 2, which can be fully characterized using the degree of synchronization of the oscillator subsets ${\rho }_1$ and ${\rho }_2$ (or populations) defined in Section 2.2. Specifically, we numerically identified chimera states by computing the mean and standard deviation of the difference between the populations’ local KOPs in time:

$\bar{d}=\frac{1}{T-r}{\displaystyle \sum _{t=r+1}^T}d\left(t\right)(8)$

$\sigma \left(d\right)=\sqrt{\frac{1}{T-r}{\displaystyle \sum _{t=r+1}^T}\left(d\left(t\right)-\bar{d}\right)},(9)$

where $T$ is the number of steps in the simulation and $r$ the number of relaxation steps. These measures were then averaged over 100 distinct initial conditions. Using these measures and an additional threshold ${\delta }_2$, we can distinguish three types of chimera states:

$•$ stable chimera state (Panaggio and Abrams, 2015): a chimera state in which the local KOPs are different and remain constant, i.e., $\bar{d}>{\delta }_1$ and $\sigma (d)<{\delta }_2$ (Equations 8, 9);

$•$ breathing chimera state (Panaggio and Abrams, 2015): a chimera state in which one of the local KOPs oscillates, i.e., $\bar{d}>{\delta }_1$ and $\sigma (d)>{\delta }_2$;

$•$ metastable and alternating chimera state (Shanahan, 2010; Buscarino et al., 2015; Feng et al., 2023): the state of a system in which the degree of synchronization of both subsets of oscillators varies irregularly or in an alternate manner, i.e., $\bar{d}<{\delta }_1$ and $\sigma (d)>{\delta }_2$;

where thresholds ${\delta }_1$ and ${\delta }_2$ are chosen to be 3 standard deviations higher than the baseline mean $\bar{d}$ and standard deviation $\sigma (d)$ calculated in the case of no mesoscale structure, i.e., $H=0$ (see also sensitivity analysis in Supplementary Figure SA10). These measures allow us to identify the regions of the parameter space in which symmetry breaking states occur in the second layer. This classification scheme is similar to that described in (Kemeth et al., 2016), with the main difference being the correlation measure. While the local curvature measure used in (Kemeth et al., 2016) is a more general correlation measure, its application to our case is not practical as it would require approximating the phase of a modules or population using the local order parameters to detect symmetry breaking states in a specific layer.

2.4 Data availability

All code and notebooks related to this work can be downloaded from (Caprioglio, 2024).

3 Results

3.1 Metastable dynamics are detected at the whole-system level

To detect metastable dynamics and chimera states, we numerically analyzed the dynamics of Kuramoto-Sakaguchi oscillators (governed by Equation 5) on 3-layer hierarchical networks at different levels of observation. We started by examining the dynamics at the whole-system level, presenting the KOP and ${\sigma }_{\text{met}}^3$ for networks of size $N=512$ with $n_1=16,n_2=8$ in Figure 3. Our analysis of the whole parameter space (Figures 3B, C) suggests that, while there is some degree of variability depending on the average degree of the network, $k$, the degree of synchronization of the system and the index of metastability mostly depend on the structural parameter $H$. Thus, we analysed in more detail our results for fixed average degree $k=N/5=51.2$, which exemplifies the results obtained for both smaller and larger average degrees (Figure 3A). Metastable dynamics are detected only when mesoscale structures are well defined ($H>0.3$, see also Figure 2D) in good agreement with our hypothesis. In contrast, in the absence of well-defined populations (i.e., low values of $H$) the system displays high coherence and a low index of metastability. As the connectivity within populations increases further, metastability reaches a local maximum around $H=0.5$ while the whole-system degree of synchronization decreases sharply. Beyond $H>0.7$, the density of connections between populations is very low (see Figure 2D). Consequently, the metastable index increases sharply while the KOP approaches 0.5, indicating that the system has decoupled into two independent populations. While metastable dynamics can already be detected at the whole-system level in the presence of well-defined populations, global measures do not provide insights into why metastability starts to increase only for values of $H>0.3$.

Figure 3

Figure 3. (A) Degree of synchronization, $R$ (blue), and index of metastability, ${\sigma }_{\text{met}}^3$ (purple), in the $l=3$ (whole-system) layer as a function of the structural parameter $H$. (B) Degree of synchronization $R$ as a function of $H$ and the average degree of the network $k$. (C) Index of metastability ${\sigma }_{\text{met}}^3$ as a function of $H$ and $k$. Results were obtained using $n_1=16$, $n_2=8$, and $k=51.2$, averaged over 100 seeds.

3.2 Self-organized chimera states are detected at the population layer

To understand why metastable dynamics are detected as the connectivity within population increases, we analysed the populations’ local KOPs, $R_{{\rho }_i}$, and the populations’ metastability index, ${\sigma }_{\text{met}}^2$ for varying values of $H$. We computed the values of $\bar{d}$ and $\sigma (d)$ to identify chimera states (Figure 4C) as illustrated in the Methods section. Independently of the value of $k$, we identified 5 regions of the parameter space, $i-v$, in which distinct states are detected in the population layer. Symmetry breaking states are detected in regions $ii-iv$, with stable $(ii)$, breathing $(iii)$, metastable and alternating $(iv)$ chimera states. When stable and breathing chimera states are detected, we define as stable population the subset of oscillators displaying the higher mean local KOP, and as unstable population the subset of oscillators with a lower degree of synchronization. As an example, in Figure 4$(iii)$, we identified ${\rho }_2$ as the stable population and ${\rho }_1$ as the unstable (breathing in this case) population. In Figure 4, we present our results for $k=51.2$.

Figure 4

Figure 4. (A) Mean degree of synchronization $R_{{\rho }_i},i=\mathrm{1,2}$ for varying values of $H$. (B) Populations’ metastability index, ${\sigma }_{\text{met}}^2$, as a function of $H$. (C) (Blue line) mean difference between the populations’ local KOP $\bar{d}$ and (orange line) standard deviation of the difference between local KOPs, $\sigma (d)$ as a function of $H$. Red dashed lines separate the regions of the parameter space in which symmetry-breaking states are detected. $(i-iv)$ Examples of the evolution in time of the populations’ local KOPs and the whole-system KOP (black dashed lines) in distinct dynamical regimes. $(i)$ For low values of $H$ the system displays coherent dynamics and the two populations display the same degree of synchronization. $(ii,iii)$ For intermediate values of $H$ stable and breathing chimera states emerge, respectively. $(iv)$ Metastable chimera states emerge for $H=0.5$. Results obtained using $n_1=16$, $n_2=8$, and $k=51.2$, averaged over 100 seeds.

In good agreement with our results of Section 3.1, metastability was found to reach its minimum in region $i$, where the absolute difference between local KOPs ($\bar{d}$ from Equation 8) was below the threshold ${\delta }_1$. This confirms that the system reaches a symmetric coherent global state (for instance, see Figure 3$(i)$) in the absence of well-defined mesoscale structures. In this region, both populations display a high degree of synchronization, equal to that of the whole system, $R=0.84\pm 0.04$, with some fluctuations due to phase lag frustration.

In the same region of the parameter space in which metastable dynamics are initially detected at the whole-system layer, we encounter stable and breathing chimera states (regions $ii$ and $iii$) in the second layer. Within these regions, the system self-organizes into a symmetry breaking state in which one of the two populations is more synchronized than the other, such that $\bar{d}>{\delta }_1$. While in region $ii$ the unstable population displays a stable degree of synchronization $(\sigma (d)<{\delta }_2)$, in region $iii$ the unstable population’s local KOP presents oscillatory dynamics $(\sigma (d)>{\delta }_2)$. This result suggests that the initial increase in metastability ${\sigma }_{\text{met}}^3$ detected at the whole-system level reflects the presence of stable and breathing chimeras in the second layer.

Metastable and alternating chimeras are detected in the range $0.49<H<0.66$ in which $\bar{d}<{\delta }_1$ and $\sigma (d)>{\delta }_2$. In the vicinity of $H=0.5$, both local KOPs display a similar degree of synchronization and the populations’ metastability, ${\sigma }_{\text{met}}^2$, is maximised. When the metastability index is at its maximum, populations (de) synchronize in a non-predictable manner. For instance, in Figure 4$(iii)$, when a population desynchronizes the other populations’ response can either be an increase or decrease in synchronization. In the same region, as the parameter $H$ is increased further, the number of connections between the two populations decreases and alternating chimera states can also found (see Supplementary Figure SA12).

In region $v$ both populations display a high local KOP but they are not synchronized between each other, indicating once again that beyond $H\sim 0.7$ the two populations become effectively uncoupled. In Figure 4$(i-iv)$, we show examples of the evolution of the local KOPs over time for four exemplifying values of the structural parameter $H$.

3.3 Chimera states can be supported by metastable modules

The remaining question is whether modules in the first layer also exhibit metastable dynamics, and whether modules in different populations display distinct dynamics depending on the symmetry-breaking states detected in layer 2.

To answer the first question, we computed the degree of metastability ${\sigma }_{\text{met}}^1$ for varying values of $H$ and $k$ (Figure 5A). While in layers 2 and 3, metastable dynamics were dependent on the structural parameter $H$, our analysis shows that in layer 1 the dominant parameter is the average degree of the network $k$. For systems of size $n_1=16$, $n_2=8$, the modules’ degree of metastability was zero for values of $k$ below $\sim 30$, the same range of values in which modularity in layer 1 is the highest (Figure 2E). In this regime, all modules remain fully synchronized for all values of $H$ (see Figure 5A$(a)$ for an example). As $k$ increases and modularity in layer 1 decreases, modules start to display metastable dynamics modulated by the value of $H$ (see Figure 5A$(b)$). Likewise, as in the previous subsection, for fixed $k>30$, the metastability index calculated for all modules peaks at $H=0.5$. However, as $k$ is increased further, the metastability index slightly increases only around $H=0.5$ while slightly decreasing for all other values of $H$ (see Supplementary Figure SA13). Thus, our analysis indicates that while the average connectivity of the network and modularity in layer 1 determines whether spontaneous (de) synchronization manifests or not in modules in layer 1, the presence of a mesoscale structure modulates the degree of metastability.

Figure 5

Figure 5. (A) Left: Metastability index ${\sigma }_{\text{met}}^1$ calculated over all modules ${\mu }_i$ as a function of $H$ and $k$. Right: example of $R_{{\mu }_i}$ dynamics in time for a single simulation with average degree $k=51.2$ (top) and $k=21$ (bottom). (B) Average local KOP within each population $j=\mathrm{1,2}$, i.e., $⟨R_{{\mu }_i}⟩$ for ${\mu }_i\subset {\rho }_j$, for fixed $k=51.2$ and varying values of $H$. (C) Metastability index within each population $j=\mathrm{1,2}$, i.e., ${\sigma }_{\text{met}}(R_{{\mu }_i})$, for ${\mu }_i\subset {\rho }_j$, for fixed $k=51.2$ and varying values of $H$, and overall metastability of the first layer blocks ${\sigma }_{\text{met}}^1$ (black dashed line). Red dashed lines separate the regions of the parameter space in which symmetry-breaking states are detected in layer 2.

To address our second question, we computed the average local KOP and the degree of metastability of modules within each population separately, with fixed $k=51.2$ and varying values of $H$ (Figures 5B, C). Here, the metastability index of modules in population ${\rho }_j$ was defined as ${\sigma }_{\text{met}}(R_{{\mu }_i})={⟨\sigma (R_{{\mu }_i})⟩}_{{\mu }_i\subset {\rho }_j}$. We divided the parameter range into the same regions $i-v$ as presented in the previous section (see Figure 4C). In the absence of a well-defined mesoscale structure (region $i$), modules in each population display the same index of metastability and the same average degree of synchronization. In contrast, in regions $ii$ and $iii$, modules in the stable and unstable populations display different degrees of synchronization and metastability. The stable population displays a higher average KOP and lower metastability while the modules in the unstable population exhibit the opposite. As expected, in regions $iv$ and $v$, modules in both populations display the same index of metastability, as well as the same degree of synchronization. Therefore, it is only in the presence of stable and breathing chimeras in layer 2 that modules in distinct populations display different degrees of metastability, reflecting the stability or instability of the populations to which the modules belong.

3.4 Relationship between symmetry-breaking states and characteristic timescales of the system

In the previous sections, we provided numerical evidence for the emergence of chimera and metastable states in hierarchically modular networks at different levels of observation. While in layer 2 different types of chimera states were identified depending on the mesoscale structure of the network (determined by $H$), in layer 1 we showed how metastability arises as the average connectivity of the network reaches a certain threshold. In this section, we aim to explain these observed phenomena by analyzing the spectral information of the graph Laplacian.

The systems we investigated are structurally homogeneous in each layer (uniform degree distribution). Nonetheless, a peculiar property of diffusive-like dynamics on hierarchically modular networks, such as synchronization dynamics, is the richness of distinct information-diffusion pathways across different timescales (Villegas et al., 2022). In the absence of phase frustration, blocks of oscillators synchronize hierarchically, from layer 1,until the whole system is synchronized (Arenas et al., 2006). As shown in (Arenas et al., 2006), the linearized Kuramoto model (without phase lags) reduces to a linear diffusion model governed by the graph Laplacian. The graph Laplacian is defined as $L=D-A$ with eigenvalues ${\lambda }_i$ ranked as $0={\lambda }_1<{\lambda }_2<\dots <{\lambda }_N$, where $A$ is the adjacency matrix of the network and $D$ is a diagonal matrix with $D_{ii}={\sum }_j^N{A}_{ij}$. Crucially, authors in (Arenas et al., 2006) showed how the gaps between eigenvalues, which separate blocks of nodes at distinct hierarchical layers, are associated with the relative difference in characteristic timescales between hierarchical layers. When the spectral gaps are small, the relative difference between hierarchical timescales is also small, and blocks at distinct hierarchical layers are not well defined. Conversely, when the spectral gaps are large, the opposite is true. Owing to the similarities with the networks studied in (Arenas et al., 2006), our Laplacian spectrum analysis (depicted in Figure 6) also reveals similar spectral gaps depending on the value of $k$ and $H$. In particular, we identify two spectral gaps:

$•$ the first spectral gap is found to be between ${\lambda }_{B_2}^{-1}$ and ${\lambda }_{B_2+1}^{-1}$, separating the slowest mode associated with the characteristic timescale of layer 3 and the slow modes associated with layer 2,

$•$ the second spectral gap is found to be between ${\lambda }_{B_1}^{-1}$ and ${\lambda }_{B_1+1}^{-1}$, separating the slow modes associated with the characteristic timescales of layer 2 and the fast modes associated with layer 1,

where $B_1$ and $B_2$ are the number of blocks in layers 1 and 2, respectively, and the inverse of each eigenvalue ${\lambda }_i$ is indicative of the characteristic timescale of mode $i$, i.e., ${\lambda }_i^{-1}\sim {\tau }_i$. Informed by these observations, we hypothesised that (a) the eigenmodes associated with the slower timescales of the system, ${\lambda }_j$ for $j\in [1,B_1]$, determine the state of the system in layers 2 and 3, and that (b) the second spectral gap between fast and slow eigenmodes affects the synchronizability and metastable dynamics of modules in layer 1.

Figure 6

Figure 6. Gaps between consecutive eigenvalues ${\lambda }_2,{\lambda }_3,\dots {\lambda }_N$ are affected by the structural parameter $H$ and the average degree $k$. In this model, increasing $k$ reduces the size of the second spectral gap, between ${\lambda }_{B_1}$ and ${\lambda }_{B_1+1}$ (where $B_1=16$ is the number of layer 1 blocks). On the other hand, increasing $H$ increases the size of the first spectral gap (inset figure), between ${\lambda }_{B_2}$ and ${\lambda }_{B_2+1}$ (where $B_2=2$ is the number of layer 2 blocks). Each eigenvalue was obtained by averaging over 10 randomly generated 3-layer networks with $n_1=16,n_2=8$.

To test hypothesis (a), we adopted the Laplacian Renormalization Group (LRG) flow approach recently introduced in (Villegas et al., 2022; Villegas et al., 2023). In summary, this approach allows us to rescale the system up to a timescale $1/{\lambda }^{*}\sim {\tau }^{*}$ by discarding the fast modes associated with eigenvalues ${\lambda }_i>{\lambda }^{*}$. Then, the resulting coarse grained structure is obtained by using the slow modes associated with ${\lambda }_i<{\lambda }^{*}$ only. We adopt the bra-ket notation used in (Villegas et al., 2023), such that $|{\lambda }_i\rangle$ denotes the column eigenvector associated with eigenvalue ${\lambda }_i$ while $\langle {\lambda }_i|$ is its corresponding row vector. Then, the dot product can be written as $⟨{\lambda }_j|{\lambda }_i⟩$ while the outer product is indicated by $|{\lambda }_i\rangle \langle {\lambda }_j|$. Using this notation, we rewrote the graph Laplacian in terms of its eigenmodes $|{\lambda }_i\rangle$, i.e.,

$L={\displaystyle \sum _i^N}{\lambda }_i|{\lambda }_i\rangle \langle {\lambda }_i|,$

and constructed the time-rescaled Laplacian at ${\tau }^{*}=1/{\lambda }^{*}=1/{\lambda }_{B_1+1}$, as

$L^{\prime }={\displaystyle \sum _{i<B_1+1}^N}{\lambda }_i|{\lambda }_i\rangle \langle {\lambda }_i|,(10)$

discarding the $N-B_1$ fast modes of the system. Next, we constructed the adjacency matrix $A^{\prime }$ from the resulting time-rescaled Laplacian $L^{\prime }$. This was done by aggregating nodes within each of the $B_1$ modules. Specifically, we defined $B_1$ orthogonal vectors $|\alpha \rangle ={({\alpha }_1,{\alpha }_2,\dots {\alpha }_N)}^{\prime }$, $\alpha =\mathrm{1,2},\dots B_1$, using the partition of the 3-layer hierarchical networks at layer 1, i.e., ${\alpha }_i=1$ if $P_{1\alpha }=\alpha$ and ${\alpha }_i=0$ otherwise. Each element of the $B_1\times B_1$ adjacency matrix $A^{\prime }$ was obtained by projecting the rescaled Laplacian (Equation 10) onto the new $B_1$-dimensional basis defined by the vectors $|\alpha \rangle$, i.e., $A_{\alpha \beta }^{\prime }=-⟨\alpha |L^{\prime }|\beta ⟩$, with $A_{\alpha \alpha }^{\prime }=0$ to avoid self interactions. The resulting weighted fully-connected adjacency matrix $A^{\prime }$ (an example is shown in Figure 7A) shares similar properties with the weighted adjacency matrix of the analytically solvable two-population model. Specifically, pairs of nodes in the same population are more strongly connected than pairs of nodes in distinct populations.

Figure 7

Figure 7. (A) Left: example of populations’ dynamics in a stable chimera state in layer 2, where $A$ was obtained using the 3-layer variation of the nSBM with $H=0.2$, $n_1=16,n_2=8$, and $k=21$. Right: the weighted coarse-grained adjacency matrix $A^{\prime }$ is obtained by time-rescaling the Laplacian associated with the original adjacency matrix $A$. The symmetry breaking parameter $a^{\prime }=0.23$, obtained by computing the normalized difference between intra-population $K_1$ (red squares) and inter-population $K_2$ couplings of the coarse-grained adjacency matrix $A^{\prime }$, is also associated with stable chimera states in the two-population model (see (Abrams et al., 2008) and results therein). (B) The symmetry breaking parameter $a^{\prime }$ for varying $H$ and six exemplifying average degrees $k$. (C) Log-log plot displaying the cutoff point when metastability in layer 1, ${\sigma }_{\text{met}}^1$, starts to sharply decrease as the size of the second spectral gap, ${\lambda }_{B_1+1}-{\lambda }_{B_1}$, increases.

While $A^{\prime }$ loses the fine-grained structural information in layer 1, it preserves the structural information in layer 2 which depends on the structural parameter $H$. This allows a numerical comparison of $a$, the symmetry breaking parameter of the two-population model, with the structural parameter $H$, where $a$ is defined as the difference between the interaction strength between pairs of nodes in the same population and that of between pairs of nodes in different populations. From the reconstructed adjacency matrix $A^{\prime }$, we obtained a numerical approximation of the symmetry breaking parameter $a^{\prime }$ by computing the normalized difference between the average strength of interaction between nodes in the same population, $K_1$, and nodes in distinct populations, $K_2$, as illustrated in Figure 7A. Next, we investigated the relationship between $a^{\prime }$ and $H$ for different average degrees of the network (see Figure 7B). Our results suggest that there is a linear relationship between $a^{\prime }$ and $H$ when the average degree is low, and a nonlinear relationship as $k$ increases. Notably, we found that when the relationship is linear or approximately linear, the symmetry-breaking parameter $a^{\prime }$ predicts well the range of $H$ in which stable and breathing chimeras are detected (for instance, see the case for $k=21$ in Figure 7A and Supplementary Figure SA11). However, we found this not to be the case for larger average degrees, such as $k=51.2$ used in the previous sections, when the relationship between $a^{\prime }$ and $H$ was found to be nonlinear. For instance, using $k=51.2$ we found stable and chimera states in the range $0.31<H<0.49$ of the structural parameter, which corresponds to the range $0.8<a^{\prime }<0.9$ of the numerical approximation of the symmetry breaking term, which is well outside the range of values in which stable and breathing chimeras appear in the two-population model (Abrams et al., 2008). This result suggests that, as the second spectral gap decreases (equivalently, as the average degree of the network increases), the slow modes of the graph Laplacian alone do not describe well the behaviours we observe in layer 2 and that the fast modes contributions should be taken into account.

To test hypothesis (b), we computed the size of the second spectral gap, ${\lambda }_{B_1+1}-{\lambda }_{B_1}$, and the degree of metastability ${\sigma }_{\text{met}}^1$ of layer 1 modules for varying $k$. As depicted in the log-log plot in Figure 7C, the degree of metastability ${\sigma }_{\text{met}}^1$ decreases sharply for values of the spectral gap ${\lambda }_{B_1+1}-{\lambda }_{B_1}>2$, indicating a clear cutoff point in the vicinity of ${\lambda }_{B_1+1}-{\lambda }_{B_1}\approx 2$ for systems of size $n_1=16,n_2=8$ with fixed value of $H=0.5$. The same behaviour was found for any other value of $H$. This result suggests that when the timescale separation between hierarchical layers 1 and 2 is large enough and modules are well defined, modules of oscillators remain completely synchronized, such that each subset effectively behaves as a single macro-oscillator. In contrast, when there is no clear timescale separation between hierarchical layers 1 and 2, the slow modes of the system (associated with ${\lambda }_i$, $i<B_1$) affect the synchronization dynamics within modules, resulting in random (de)synchronization patterns that persist in the system.

3.5 Robustness against structural perturbations and heterogeneous natural frequencies

Our numerical analyses have shown the relationship between different symmetry-breaking states and the system’s hierarchically modular structure. In layers 2 and 3, chimera and metastable states are detected as the populations become more defined (with increasing $H$). In layer 1 instead, we found the opposite relationship: as modules become less defined for increasing $k$ (as shown by decreasing spectral gap and modularity) metastable dynamics are detected. To measure the robustness of these observed dynamics, we analysed systems generated with fixed values of $H$ and $k$ when perturbed in two separate ways: by randomly rewiring $r$ edges in the network (Figure 8), and by including heterogeneous frequencies ${\omega }_i$ drawn from a normal distribution $\mathcal{N}(1,\delta \omega )$ with mean 1 and increasing standard deviation $\delta \omega$ (Figure 9). We report results for systems constructed using $H=0.5$, since most peaks in metastability in all layers were found for this value, and using two exemplifying average degrees $k=51.2$ and $k=21$, since metastability in layer 1 was found to depend on the size of the second spectral gap, itself controlled by the average degree of the network.

Figure 8

Figure 8. Structural perturbation analysis for fixed $H=0.5$ and for $k=51.2$ (top row) and $k=21$ (bottom row). In each layer $l$, we compute the average degree of synchronization $R_l$ of all blocks in the layer (left column) and the metastability index ${\sigma }_{\text{met}}^l$ (middle column). The modularity index $Q_l$ is computed with respect to the partitions in layer 1 and 2 (right column). Results were obtained for systems of size $n_1=16,n_2=8$ and averaged over 100 seeds.

Figure 9

Figure 9. Impact of heterogeneous frequencies with increasing $\delta \omega$ for $H=0.5$ and for $k=51.2$ (top row) and $k=21$ (bottom row). In each layer $l$, we compute the average degree of synchronization $R_l$ of all blocks in layer $l$ (left column) and the metastability index ${\sigma }_{\text{met}}^l$ (right column). Results were obtained for systems of size $n_1=16,n_2=8$, using normalized coupling $K/k$ with $K=50$, and averaged over 100 seeds.

By randomly rewiring $r$ edges in the network, the average degree of the network remains fixed while modularity in layers 2 and 1 is progressively lost (Figure 4). We found that the system starts to lose its modular features when $r=10^3$ random rewirings are performed, i.e., $10\%$ of the edges at most since the expected number of edges is $\mathcal{E}=kN$, and completely loses modularity (as well as the spectral gaps) when $10^4$ rewirings are performed. As expected, measures of metastability and synchrony behave differently in layer 1 compared to layers 2 and 3. While the system loses its metastable features in layers 2 and 3 when $r=10^3$ random rewirings are performed, in layer 1 metastability slightly decreases in the case of $k=51.2$ and sharply increases in the case of $k=21$. This analysis reflect the fact that the presence of a mesoscale structure in layer 2 and 3 is necessary for metastable dynamics to emerge in the upper layers. Modules of oscillators in layer 1, instead, are characterized by the absence of phase lags between pairs of oscillators in the same module, and therefore display metastable dynamics only when their modular structure is less defined, or, equivalently, when the second spectral gap disappears (which is the case for $k=51.2$). While not the focus of this study, the presence of a high degree of metastability in a random network ($10^4$ rewirings) with non-overlapping subsets of oscillators without phase-lagged interaction prompted further investigation. In Supplementary Figure SA14, we report an example of a system of oscillators constructed using the configuration model with a uniform degree distribution $k=51$. The random network system is only characterized by the presence of a partition of oscillators, similar to the partition in layer 1 we studied here, in which only pairs of oscillators that do not belong to the same subset $c$ in the partition have phase-lagged interaction. Interestingly, results confirm the presence of metastable chimera states very similar to those displayed in the Shanahan model (Shanahan, 2010) despite the absence of heterogeneous couplings and modular structures.

With the addition of natural frequencies, the system’s coupling becomes an important factor in the dynamics. It is well known that oscillatory systems display various degrees of metastability within the range of critical couplings (Villegas et al., 2014). However, we seek to investigate the behaviour of the system away from the critical synchronization point to assess in what range of values of $\delta \omega$ the results we reported in the previous section with $K=50/k$ can hold with the addition of heterogeneous frequencies. In the case of $k=51.2$ $(k=21)$, the analysis reported in Figure 9 shows that for values of $\delta \omega <2$ $(\delta \omega <10)$ metastable dynamics at all layers remain unchanged. As the heterogeneity of frequencies is increased further, all blocks of oscillators in each layer display incoherent dynamics as each average local KOP $R_l$ reaches the inverse of the square root of the number of oscillator in a block, i.e., ${\lim }_{\delta \omega \gg ⟨\omega ⟩}{R}_l={\left({\prod }_{k=1}^i{n}_k\right)}^{-1/2}$ (for instance, $R_1$ reaches $1/\sqrt{n_1}=1/\sqrt{16}$). Consequently, the indexes of metastability for high frequency heterogeneities reflect this incoherence. Peaks in metastability are also detected at values of $\delta \omega$ for which the slope of decreasing $R_l$ is maximal (dashed red lines in Figure 9), in good agreement with previous studies (Villegas et al., 2014).

4 Discussion

One of the hallmarks of complex systems is the presence of hierarchically modular structures: systems composed of interrelated substructures whose interaction across scales is thought to be fundamental for robust and efficient information processing (Simon, 1991). The identification of interrelated subsystems and their role in shaping the dynamics and in the emergence of macroscopic collective patterns is an important ongoing research endeavour in complexity science (Jensen, 2023) and particularly in network neuroscience (Bassett and Gazzaniga, 2011). In this work, we hypothesised that the interplay between nested structures alone can give rise to chimera and metastable states in minimal brain-inspired models of phase-frustrated oscillatory systems. First, we provided evidence for the emergence of symmetry-breaking states in the absence of heterogeneous couplings, delays, and structural heterogeneities. Second, we revealed the relationships between the symmetry-breaking states we observed and the spectral properties typical of hierarchically modular networks.

Our analysis highlighted two distinct pathways towards achieving the metastable regime in different layers. First, metastable dynamics were detected in layers 2 and 3 at a point in which stable and breathing chimeras cease to exist. While it was possible to explain the emergence of stable and breathing chimera states by integrating out the fast modes of the system and reducing the network to a structure similar to that of the two-population model, the precise mechanism that leads to metastable chimera states remains unclear. Due to the small size of the rescaled system, we cannot exclude the possibility that the behaviours we observed are the result of finite size effects. Additionally, a good agreement with the results in (Abrams et al., 2008) was only possible for low values of $k$ and the presence of a large second spectral gap. This result suggests that, when the gap between fast and slow modes is small, the contribution from the fast modes should not be discarded. To reach a conclusive answer, an analytical approach is likely necessary. In (Abrams et al., 2008), the two-population model was solved via the Ott-Antonsen dimensionality reduction approach (Ott and Antonsen, 2008), by showing the existence of inhomogeneous solutions to the dynamical system of equations in the limit of an infinite number of nodes in each population. The inhomogenous solutions correspond to the stable and breathing chimera states which are very similar to the chimeras observed in our model. However, the two-population model does not predict the emergence of the metastable chimera states we observed for a wide range of the parameter space. Due to the nested nature of the order parameters in our model and the absence of heterogeneous couplings, a similar analytical approach is not immediate. Informed by our observations, we suggest that a synergetic-based study of hierarchical systems of oscillators could be an interesting future research approach (Zheng et al., 2024). As an intuitive example, in Section 3.3 we found that chimera states detected in layer 2 coexist with highly metastable modules in layer 1. Crucially, these modules display a high average degree of synchronization (Figure 5B) despite the high degree of metastability. This result suggests a synergetic interpretation of the behaviours we observed, in which the spontaneous (de) synchronization patterns of the modules may be interpreted as fast fluctuations of the “enslaved” fine-grained order parameters of the system (Haken, 2004).

Second, the emergence of spontaneous (de) synchronization patterns in layer 1 was found to be related to the size of the second spectral gap. When there is enough separation of timescales and modules are structurally well defined (high modularity), each oscillator within a module has time to integrate with other units in the same module the phase-frustrated inter-module interactions. When the separation of timescales is not large enough, however, modules do not have time to integrate those interactions, resulting in the spontaneous desynchronization of oscillators within the module. This behaviour has been observed in Sections 3.3, 3.4, in which we systematically increased the average degree of the network, leading to the loss of modularity in layer 1, as well as in Section 3.5, in which modularity was lost by randomly rewiring edges in the network. We also considered the more parsimonious hypothesis that spontaneous (de) synchronization is caused by increasing the average degree of the network since this increases the number of phase-frustrated interactions in the system. However, as the average degree of the network is increased beyond a threshold value, layer 1 modules do not display increasing metastability; instead metastability reaches a plateau for all degrees $k>30$ (Supplementary Figure SA13). Additionally, our structural perturbation analysis in Section 3.5 revealed the presence of highly metastable dynamics in layer 1 even in the absence of any modular structure (see also Supplementary Figure SA14). This begs the question of the impact of the presence of phase lags only in the upper layers on the mechanism underlying the metastable dynamics we observed in layer 1. In the absence of phase lags and any other source of frustration or symmetry-breaking term, the system possesses a single solution, the synchronization manifold. Therefore, no metastable dynamics can be detected in this case (not shown). In contrast, the addition of phase lagged interactions in layer 1, i.e., ${\alpha }_{ij}=\alpha$ for all pairs $i,j$, has some non-trivial effects. Using our variation of the nested stochastic block model, our results do not hold if the number of oscillators in the first layer is small (for instance using $n_1=16$). However, if the number of oscillators is increased, for instance to $n_1=64$, preliminary results (Supplementary Figure SA15) show that similar metastable behaviours in the first layer can be observed when the second spectral gap decreases to a minimum threshold while the average connectivity of the network increases, i.e., the same mechanism underlying metastability in the first layer we observed here. This suggests that our results are not exclusive to systems where phase lags are present only in the upper layers.

An important distinction evident in our model is that between metastability at the critical synchronization point and metastability away from the critical synchronization point. In general, at the phase transition between coherence and incoherence indexes of metastability detect the presence of large fluctuations in the global and local order parameters (Arenas et al., 2008; Villegas et al., 2014; Mackay et al., 2023) indicative of the metastable regime. In fact, these behaviours can be observed in our robustness analysis in Section 3.5 with the addition of heterogeneous frequencies. More specifically, in Figure 9 we comment on the observation that at the transition between coherence and incoherence, at around $\delta \omega \approx 10$, the metastability indexes detect large fluctuations in the order parameters at all layers. Crucially, these metastable behaviours are different in nature from the fluctuations we observed in Sections 3.1–3.3. The systems we studied here are characterized by the presence of identical oscillators and a homogeneous coupling. In this case, one may expect a system to always reach the synchronized state given enough time. However, as we have discussed, fluctuations in the order parameters can still arise due to the instability of chimera states in layer 2 (Section 3.2) or due to the spontaneous (de) synchronization of modules in layer 1 when the timescale separation between slow and fast modes is below a critical threshold (Section 3.3). This distinction reflects the important observation highlighted in (Hancock et al., 2023a) that metastability is a dynamical regime rather than a mechanism. Therefore, as we also point out here, multiple distinct mechanisms can lead to the metastable regime.

Our systematic numerical analysis across different levels of observation allowed us to detect and compare metastable dynamics and chimera states in different hierarchical layers. However, it also highlighted the limitations of the variation of the whole-system KOP as a measure of metastability. For example, while we detected metastability at the whole-system level when the structural parameter ranged between 0.3 and 0.49, layer 2 analysis revealed the presence of stable and breathing chimera states instead. At the same time, metastable dynamics in layer 1 remained entirely undetected in our analysis at the whole-system level. In computational neuroscience, the variation of the whole-system KOP is often used to tune whole-brain dynamical models, and has been hypothesized to peak when the correlation with dFC is maximised (Deco and Kringelbach, 2016). In the model we studied, and which allows a more systematic analysis of the indexes of metastability due to the model’s inherent modular construction, the peaks in metastability were all found to be in the same region of the parameter space, except for the value of metastability of the modules in the unstable population (Figure 5C). Therefore, our findings suggest that a thorough analysis across scales should be performed when possible, and that the index of metastability should be properly defined in terms of the relevant order parameters of the system. This approach may help resolve inconsistencies in the current literature, where maximising the variation of the whole-system order parameter has not necessarily aligned with the best empirical fit (Pope et al., 2023).

When dealing with empirical structural brain networks, such analysis is often not possible due to the lack of well defined communities and consequently the absence of clearly identified relevant local order parameters of the system. Both in our model and the Shanahan model (Shanahan, 2010), where the metastability index was originally formulated, the systems exhibit a clear modular structure, inherent to the system’s construction. Consequently, the issue of identifying the relevant local KOPs is absent. Novel approaches based on the harmonic decomposition of structural connectomes (Atasoy et al., 2016; Atasoy et al., 2017) have been proposed [see also the eigen-microstate approach (Zheng et al., 2024)], which we suggest could overcome this problem. Within this framework, dynamical functional connectivities can be characterized by computing the contribution of each eigenmode of the system at each point in time (Luppi et al., 2024; Zheng et al., 2024). Inspired by the harmonic brain modes framework and the work of Villegas and colleagues (Villegas et al., 2014; Villegas et al., 2022; Villegas et al., 2023), we analysed our systems using the spectral information of the graph Laplacian. Our analysis in layer 1 showed that when the relative separation of timescales between modes is large enough the fast eigenmodes associated with dynamics in the first layer do not contribute to the dynamics of the system, since modules remain completely synchronized at all times. In layer 2, instead, we found that the slow modes of the system encode the relevant information about the system’s behaviours at higher coarse-grainings. In particular, the time-rescaled system reduced to the two-population model for a specific range of parameters. These results highlight the effectiveness of the LRG approach to identify the critical mesoscale structures with a characteristic timescale, and suggest that a similar approach may be used to identify the relevant order parameters of the system at different levels of observation, as well as the spectral profiles that may facilitate hierarchical integration (Wang et al., 2021) in oscillatory systems.

The limitations imposed on the system have allowed us to study the effects of the presence of a mesoscale structure in isolation, however, it is unclear if the behaviours we observed would persist with the addition of structural heterogeneities such as core-periphery structures, and, in particular, the presence of hubs. The rich club has been shown to change the path towards synchronization (Gómez-Gardeñes et al., 2007) when compared to random networks: rather than small clusters synchronizing first, the subset of hubs and the nodes more strongly connected to them synchronize first. The combination of rich club and modular structure has also been investigated in (Arenas and Díaz-Guilera, 2007; Zamora-López et al., 2016), highlighting the role hubs play in the formation of isolated synchronized communities and their contribution to the overall functional complexity. It remains unclear how these structural heterogeneities affect systems of phase-lagged or phase-delayed oscillators and, in particular, away from the critical regime. As noticed in (Arenas and Díaz-Guilera, 2007), the dynamical equation corresponding to a hub node is a topological average of the phases of its neighbours. In the presence of frustration, it is unclear if the frustrations get cancelled out, promoting integration and synchronization, or have the opposite effect. The answer to this will likely depend on the choice between uniform or heterogeneous lags, or, with the addition of an underlying metric space, heterogeneous delays.

Finally, we believe that the study of metastability in hyperbolic networks (Krioukov et al., 2010) of oscillators might help address these limitations. Hyperbolic networks have been shown to display many of the universal properties of real-world networks (Krioukov et al., 2010), including clustering, hubs, rich club, and modularity (Kovács and Palla, 2021; Balogh et al., 2023). Crucially, hyperbolic networks come equipped with an underlying metric space, allowing a geometric renormalization procedure (García-Pérez et al., 2018) which has been shown to detect the possible organizing principles of the brain (Zheng et al., 2020) across scales.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author contributions

EC: Conceptualization, Formal Analysis, Methodology, Software, Writing–original draft, Writing–review and editing. LB: Conceptualization, Funding acquisition, Methodology, Supervision, Writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. The PhD studentship of the first author was fully funded by Moogsoft Ltd. The funder had no role in the research.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnetp.2024.1436046/full#supplementary-material

Footnotes

¹Note the distinction between metastability and multistability. While in multistable systems a trajectory can escape a basin of attraction only due to external perturbations or noise, trajectories in metastable systems do this spontaneously (Hancock et al., 2023a).

References

Abrams, D. M., Mirollo, R., Strogatz, S. H., and Wiley, D. A. (2008). Solvable model for chimera states of coupled oscillators. Phys. Rev. Lett. 101 (8), 084103. doi:10.1103/physrevlett.101.084103

PubMed Abstract | CrossRef Full Text | Google Scholar

Abrams, D. M., and Strogatz, S. H. (2004). Chimera states for coupled oscillators. Phys. Rev. Lett. 93 (17), 174102. doi:10.1103/physrevlett.93.174102

PubMed Abstract | CrossRef Full Text | Google Scholar

Alteriis, G. de, MacNicol, E., Hancock, F., Ciaramella, A., Cash, D., Expert, P., et al. (2024). EiDA: a lossless approach for dynamic functional connectivity; application to fMRI data of a model of ageing. Imaging Neurosci. 2, 1–22. doi:10.1162/imag_a_00113

CrossRef Full Text | Google Scholar

Arenas, A., and Díaz-Guilera, A. (2007). Synchronization and modularity in complex networks. Eur. Phys. J. Special Top. 143 (1), 19–25. doi:10.1140/epjst/e2007-00066-2

CrossRef Full Text | Google Scholar

Arenas, A., Díaz-Guilera, A., Kurths, J., Moreno, Y., and Zhou, C. (2008). Synchronization in complex networks. Phys. Rep. 469 (3), 93–153. doi:10.1016/j.physrep.2008.09.002

CrossRef Full Text | Google Scholar

Arenas, A., Díaz-Guilera, A., and Pérez-Vicente, C. J. (2006). Synchronization reveals topological scales in complex networks. Phys. Rev. Lett. 96 (11), 114102. doi:10.1103/physrevlett.96.114102

PubMed Abstract | CrossRef Full Text | Google Scholar

Atasoy, S., Deco, G., Kringelbach, M. L., and Pearson, J. (2017). Harmonic brain modes: a unifying framework for linking space and time in brain dynamics. Neurosci. 24 (3), 277–293. doi:10.1177/1073858417728032

CrossRef Full Text | Google Scholar

Atasoy, S., Donnelly, I., and Pearson, J. (2016). Human brain networks function in connectome-specific harmonic waves. Nat. Commun. 7 (1), 10340. doi:10.1038/ncomms10340

PubMed Abstract | CrossRef Full Text | Google Scholar

Balogh, S. G., Kovács, B., and Palla, G. (2023). Maximally modular structure of growing hyperbolic networks. Commun. Phys. 6 (1), 76. doi:10.1038/s42005-023-01182-4

CrossRef Full Text | Google Scholar

Bapat, R., Pathak, A., and Banerjee, A. (2024). Metastability indexes global changes in the dynamic working point of the brain following brain stimulation. Front. Neurorobotics 18, 1336438. doi:10.3389/fnbot.2024.1336438

PubMed Abstract | CrossRef Full Text | Google Scholar

Bassett, D. S., and Gazzaniga, M. S. (2011). Understanding complexity in the human brain. Trends Cognitive Sci. 15 (5), 200–209. doi:10.1016/j.tics.2011.03.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Bick, C. (2018). Heteroclinic switching between chimeras. Phys. Rev. E 97 (5), 050201. doi:10.1103/physreve.97.050201

PubMed Abstract | CrossRef Full Text | Google Scholar

Bick, C., Panaggio, M. J., and Martens, E. A. (2018). Chaos in Kuramoto oscillator networks. Chaos Interdiscip. J. Nonlinear Sci. 28 (7), 071102. doi:10.1063/1.5041444

PubMed Abstract | CrossRef Full Text | Google Scholar

Buscarino, A., Frasca, M., Gambuzza, L. V., and Hövel, P. (2015). Chimera states in time-varying complex networks. Phys. Rev. E 91 (2), 022817. doi:10.1103/physreve.91.022817

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabral, J., Castaldo, F., Vohryzek, J., Litvak, V., Bick, C., Lambiotte, R., et al. (2022). Metastable oscillatory modes emerge from synchronization in the brain spacetime connectome. Commun. Phys. 5 (1), 184. doi:10.1038/s42005-022-00950-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabral, J., Hugues, E., Sporns, O., and Deco, G. (2011). Role of local network oscillations in resting-state functional connectivity. NeuroImage 57 (1), 130–139. doi:10.1016/j.neuroimage.2011.04.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabral, J., Kringelbach, M. L., and Deco, G. (2014). Exploring the network dynamics underlying brain activity during rest. Prog. Neurobiol. 114, 102–131. doi:10.1016/j.pneurobio.2013.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabral, J., Vidaurre, D., Marques, P., Magalhães, R., Silva Moreira, P., Miguel Soares, J., et al. (2017). Cognitive performance in healthy older adults relates to spontaneous switching between states of functional connectivity during rest. Sci. Rep. 7 (1), 5135. doi:10.1038/s41598-017-05425-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Capouskova, K., Zamora-López, G., Kringelbach, M. L., and Deco, G. (2023). Integration and segregation manifolds in the brain ensure cognitive flexibility during tasks and rest. Hum. Brain Mapp. 44 (18), 6349–6363. doi:10.1002/hbm.26511

PubMed Abstract | CrossRef Full Text | Google Scholar

Caprioglio, E. (2024). Emergence of metastability in frustrated oscillatory networks: the key role of hierarchical modularity. Available at: https://github.com/EnricoCaprioglio/Emergence-of-metastability-in-frustrated-oscillatory-networks.

Google Scholar

Deco, G., and Kringelbach, M. L. (2016). Metastability and coherence: extending the communication through coherence hypothesis using A whole-brain computational perspective. Trends Neurosci. 39 (3), 125–135. doi:10.1016/j.tins.2016.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Deco, G., Kringelbach, M. L., Jirsa, V. K., and Ritter, P. (2017). The dynamics of resting fluctuations in the brain: metastability and its dynamical cortical core. Sci. Rep. 7 (1), 3095. doi:10.1038/s41598-017-03073-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, P., Yang, J., Wu, Y., and Liu, Z. (2023). Alternating chimera states in complex networks with modular structures. Chaos Interdiscip. J. Nonlinear Sci. 33 (3), 033136. doi:10.1063/5.0132072

CrossRef Full Text | Google Scholar

Fox, M. D., and Raichle, M. E. (2007). Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging. Nat. Rev. Neurosci. 8 (9), 700–711. doi:10.1038/nrn2201

PubMed Abstract | CrossRef Full Text | Google Scholar

Fuscà, M., Siebenhühner, F., Wang, S. H., Myrov, V., Arnulfo, G., Nobili, L., et al. (2023). Brain criticality predicts individual levels of inter-areal synchronization in human electrophysiological data. Nat. Commun. 14 (1), 4736. doi:10.1038/s41467-023-40056-9

PubMed Abstract | CrossRef Full Text | Google Scholar

García-Pérez, G., Boguñá, M., and Serrano, M. Á. (2018). Multiscale unfolding of real networks by geometric renormalization. Nat. Phys. 14 (6), 583–589. doi:10.1038/s41567-018-0072-5

CrossRef Full Text | Google Scholar

Gómez-Gardeñes, J., Moreno, Y., and Arenas, A. (2007). Paths to synchronization on complex networks. Phys. Rev. Lett. 98 (3), 034101. doi:10.1103/physrevlett.98.034101

PubMed Abstract | CrossRef Full Text | Google Scholar

Haken, H. (2004). Synergetics: introduction and advanced topics. Berlin Heidelberg: Springer.

Google Scholar

Hancock, F., Rosas, F. E., McCutcheon, R. A., Cabral, J., Dipasquale, O., and Turkheimer, F. E. (2023b). Metastability as a candidate neuromechanistic biomarker of schizophrenia pathology. PLOS ONE 18 (3), e0282707. doi:10.1371/journal.pone.0282707

PubMed Abstract | CrossRef Full Text | Google Scholar

Hancock, F., Rosas, F. E., Zhang, M., Mediano, P. A. M., Luppi, A., Cabral, J., et al. (2023a). Metastability demystified — the foundational past, the pragmatic present, and the potential future. Preprints. doi:10.20944/preprints202307.1445.v1

CrossRef Full Text | Google Scholar

Hansen, E. C., Battaglia, D., Spiegler, A., Deco, G., and Jirsa, V. K. (2015). Functional connectivity dynamics: modeling the switching behavior of the resting state. NeuroImage 105, 525–535. doi:10.1016/j.neuroimage.2014.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Haugland, S. W. (2021). The changing notion of chimera states, a critical review. J. Phys. Complex. 2 (3), 032001. doi:10.1088/2632-072x/ac0810

CrossRef Full Text | Google Scholar

Haugland, S. W., Schmidt, L., and Krischer, K. (2015). Self-organized alternating chimera states in oscillatory media. Sci. Rep. 5 (1), 9883. doi:10.1038/srep09883

PubMed Abstract | CrossRef Full Text | Google Scholar

Hellyer, P. J., Scott, G., Shanahan, M., Sharp, D. J., and Leech, R. (2015). Cognitive flexibility through metastable neural dynamics is disrupted by damage to the structural connectome. J. Neurosci. 35 (24), 9050–9063. doi:10.1523/jneurosci.4648-14.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Jensen, H. J. (2023). Complexity science: the study of emergence. New York, NY: Cambridge University Press.

Google Scholar

Kang, L., Tian, C., Huo, S., and Liu, Z. (2019). A two-layered brain network model and its chimera state. Sci. Rep. 9 (1), 14389. doi:10.1038/s41598-019-50969-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Kelso, J. A. S. (2013). “Coordination dynamics,” in Encyclopedia of complexity and systems science. New York: Springer, 1–41.

CrossRef Full Text | Google Scholar

Kemeth, F. P., Haugland, S. W., Schmidt, L., Kevrekidis, I. G., and Krischer, K. (2016). A classification scheme for chimera states. Chaos Interdiscip. J. Nonlinear Sci. 26 (9), 094815. doi:10.1063/1.4959804

PubMed Abstract | CrossRef Full Text | Google Scholar

Kovács, B., and Palla, G. (2021). The inherent community structure of hyperbolic networks. Sci. Rep. 11 (1), 16050. doi:10.1038/s41598-021-93921-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Krioukov, D., Papadopoulos, F., Kitsak, M., Vahdat, A., and Boguñá, M. (2010). Hyperbolic geometry of complex networks. Phys. Rev. E 82 (3), 036106. doi:10.1103/physreve.82.036106

PubMed Abstract | CrossRef Full Text | Google Scholar

Krishnagopal, S., Lehnert, J., Poel, W., Zakharova, A., and Schöll, E. (2017). Synchronization patterns: from network motifs to hierarchical networks. Philosophical Trans. R. Soc. A Math. Phys. Eng. Sci. 375 (2088), 20160216. doi:10.1098/rsta.2016.0216

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, W. H., Doucet, G. E., Leibu, E., and Frangou, S. (2018). Resting-state network connectivity and metastability predict clinical symptoms in schizophrenia. Schizophrenia Res. 201, 208–216. doi:10.1016/j.schres.2018.04.029

PubMed Abstract | CrossRef Full Text | Google Scholar

López-González, A., Panda, R., Ponce-Alvarez, A., Zamora-López, G., Escrichs, A., Martial, C., et al. (2021). Loss of consciousness reduces the stability of brain hubs and the heterogeneity of brain dynamics. Commun. Biol. 4 (1), 1037. doi:10.1038/s42003-021-02537-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Lord, L. D., Expert, P., Atasoy, S., Roseman, L., Rapuano, K., Lambiotte, R., et al. (2019). Dynamical exploration of the repertoire of brain networks at rest is modulated by psilocybin. NeuroImage 199, 127–142. doi:10.1016/j.neuroimage.2019.05.060

PubMed Abstract | CrossRef Full Text | Google Scholar

Luppi, A. I., Uhrig, L., Tasserie, J., Signorelli, C. M., Stamatakis, E. A., Destexhe, A., et al. (2024). Local orchestration of distributed functional patterns supporting loss and restoration of consciousness in the primate brain. Nat. Commun. 15 (1), 2171. doi:10.1038/s41467-024-46382-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Mackay, M., Huo, S., and Kaiser, M. (2023). Spatial organisation of the mesoscale connectome: a feature influencing synchrony and metastability of network dynamics. PLOS Comput. Biol. 19 (8), e1011349. doi:10.1371/journal.pcbi.1011349

PubMed Abstract | CrossRef Full Text | Google Scholar

Makarov, V. V., Kundu, S., Kirsanov, D. V., Frolov, N. S., Maksimenko, V. A., Ghosh, D., et al. (2019). Multiscale interaction promotes chimera states in complex networks. Commun. Nonlinear Sci. Numer. Simul. 71, 118–129. doi:10.1016/j.cnsns.2018.11.015

CrossRef Full Text | Google Scholar

Meunier, D., Lambiotte, R., and Bullmore, E. T. (2010). Modular and hierarchically modular organization of brain networks. Front. Neurosci. 4, 200. doi:10.3389/fnins.2010.00200

PubMed Abstract | CrossRef Full Text | Google Scholar

Moretti, P., and Muñoz, M. A. (2013). Griffiths phases and the stretching of criticality in brain networks. Nat. Commun. 4 (1), 2521. doi:10.1038/ncomms3521

PubMed Abstract | CrossRef Full Text | Google Scholar

Ódor, G., Dickman, R., and Ódor, G. (2015). Griffiths phases and localization in hierarchical modular networks. Sci. Rep. 5 (1), 14451. doi:10.1038/srep14451

PubMed Abstract | CrossRef Full Text | Google Scholar

Ott, E., and Antonsen, T. M. (2008). Low dimensional behavior of large systems of globally coupled oscillators. Chaos Interdiscip. J. Nonlinear Sci. 18 (3), 037113. doi:10.1063/1.2930766

PubMed Abstract | CrossRef Full Text | Google Scholar

Panaggio, M. J., and Abrams, D. M. (2015). Chimera states: coexistence of coherence and incoherence in networks of coupled oscillators. Nonlinearity 28 (3), R67–R87. doi:10.1088/0951-7715/28/3/r67

CrossRef Full Text | Google Scholar

Peixoto, T. P. (2014). Hierarchical block structures and high-resolution model selection in large networks. Phys. Rev. X 4 (1), 011047. doi:10.1103/physrevx.4.011047

CrossRef Full Text | Google Scholar

Pope, M., Seguin, C., Varley, T. F., Faskowitz, J., and Sporns, O. (2023). Co-evolving dynamics and topology in a coupled oscillator model of resting brain function. NeuroImage 277, 120266. doi:10.1016/j.neuroimage.2023.120266

PubMed Abstract | CrossRef Full Text | Google Scholar

Rubinov, M., Sporns, O., Thivierge, J. P., and Breakspear, M. (2011). Neurobiologically realistic determinants of self-organized criticality in networks of spiking neurons. PLoS Comput. Biol. 7 (6), e1002038. doi:10.1371/journal.pcbi.1002038

PubMed Abstract | CrossRef Full Text | Google Scholar

Shanahan, M. (2010). Metastable chimera states in community-structured oscillator networks. Chaos Interdiscip. J. Nonlinear Sci. 20 (1), 013108. doi:10.1063/1.3305451

PubMed Abstract | CrossRef Full Text | Google Scholar

Simon, H. A. (1991). “The architecture of complexity,” in Facets of systems science. Springer US, 457–476.

CrossRef Full Text | Google Scholar

Tognoli, E., and Kelso, J. A. S. (2014). The metastable brain. Neuron 81 (1), 35–48. doi:10.1016/j.neuron.2013.12.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Torres, F. A., Otero, M., Lea-Carnall, C. A., Cabral, J., Weinstein, A., and El-Deredy, W. (2024). Emergence of multiple spontaneous coherent subnetworks from a single configuration of human connectome-coupled oscillators model. bioRxiv. doi:10.1101/2024.01.09.574845

CrossRef Full Text | Google Scholar

Ulonska, S., Omelchenko, I., Zakharova, A., and Schöll, E. (2016). Chimera states in networks of Van der Pol oscillators with hierarchical connectivities. Chaos Interdiscip. J. Nonlinear Sci. 26 (9), 094825. doi:10.1063/1.4962913

PubMed Abstract | CrossRef Full Text | Google Scholar

Váša, F., Shanahan, M., Hellyer, P. J., Scott, G., Cabral, J., and Leech, R. (2015). Effects of lesions on synchrony and metastability in cortical networks. NeuroImage 118, 456–467. doi:10.1016/j.neuroimage.2015.05.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Villegas, P., Gabrielli, A., Santucci, F., Caldarelli, G., and Gili, T. (2022). Laplacian paths in complex networks: information core emerges from entropic transitions. Phys. Rev. Res. 4 (3), 033196. doi:10.1103/physrevresearch.4.033196

CrossRef Full Text | Google Scholar

Villegas, P., Gili, T., Caldarelli, G., and Gabrielli, A. (2023). Laplacian renormalization group for heterogeneous networks. Nat. Phys. 19 (3), 445–450. doi:10.1038/s41567-022-01866-8

CrossRef Full Text | Google Scholar

Villegas, P., Moretti, P., and Muñoz, M. A. (2014). Frustrated hierarchical synchronization and emergent complexity in the human connectome network. Sci. Rep. 4 (1), 5990. doi:10.1038/srep05990

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, R., Liu, M., Cheng, X., Wu, Y., Hildebrandt, A., and Zhou, C. (2021). Segregation, integration, and balance of large-scale resting brain networks configure different cognitive abilities. Proc. Natl. Acad. Sci. 118 (23), e2022288118. doi:10.1073/pnas.2022288118

PubMed Abstract | CrossRef Full Text | Google Scholar

Zamora-López, G., Chen, Y., Deco, G., Kringelbach, M. L., and Zhou, C. (2016). Functional complexity emerging from anatomical constraints in the brain: the significance of network modularity and rich-clubs. Sci. Rep. 6 (1), 38424. doi:10.1038/srep38424

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, M., Allard, A., Hagmann, P., Alemán-Gómez, Y., and Serrano, M. Á. (2020). Geometric renormalization unravels self-similarity of the multiscale human connectome. Proc. Natl. Acad. Sci. 117 (33), 20244–20253. doi:10.1073/pnas.1922248117

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, Z., Xu, C., Fan, J., Liu, M., and Chen, X. (2024). Order parameter dynamics in complex systems: from models to data. Chaos Interdiscip. J. Nonlinear Sci. 34 (2), 022101. doi:10.1063/5.0180340

CrossRef Full Text | Google Scholar