Collective Activity Bursting in a Population of Excitable Units Adaptively Coupled to a Pool of Resources

We study the collective dynamics in a population of excitable units (neurons) adaptively interacting with a pool of resources. The resource pool is influenced by the average activity of the population, whereas the feedback from the resources to the population is comprised of components acting homogeneously or inhomogeneously on individual units of the population. Moreover, the resource pool dynamics is assumed to be slow and has an oscillatory degree of freedom. We show that the feedback loop between the population and the resources can give rise to collective activity bursting in the population. To explain the mechanisms behind this emergent phenomenon, we combine the Ott-Antonsen reduction for the collective dynamics of the population and singular perturbation theory to obtain a reduced system describing the interaction between the population mean field and the resources.

1 Introduction

Complex dynamical networks are indispensable for modeling many processes in nature, technology, and social sciences (Strogatz, 2001; Boccaletti et al., 2006; Arenas et al., 2008; Yanchuk et al., 2021). In realistic situations, collective dynamics in such networks is affected by the constraints on available resources from the environment (Roberts et al., 2014; Kroma-Wiley et al., 2021), resulting in complex dynamical phenomena, especially if the systems are self-organized to operate close to criticality (Levina et al., 2007). Often, additional resource dynamics gives rise to adaptive mechanisms such as frequency adaptation (Fuhrmann et al., 2002; Taylor et al., 2010; Ha and Cheong, 2017; Kroma-Wiley et al., 2021), delay adaptation (Fields, 2015; Park and Lefebvre, 2020), or various forms of homeostatic plasticity in neuronal systems (Zierenberg et al., 2018).

Compared with other somatic cells, neurons have a very high energy consumption (Attwell and Laughlin, 2001) and are highly sensitive to energy limitations affecting their cellular metabolic processes. Hence, the availability of metabolic resources, their dynamics and their interplay with the neuronal activity are important factors for the overall performance of neural networks and their homeostasis (Vergara et al., 2019). Dynamical networks with resource constraints have been in the focus of recent studies (Taylor et al., 2010; Roberts et al., 2014; Virkar et al., 2016; Nicosia et al., 2017; Song et al., 2020; Kroma-Wiley et al., 2021). In particular, in (Song et al., 2020) it has been investigated how phase synchronization between the mutually uncoupled system elements depends on the interaction with the environment. A mini-review (Roberts et al., 2014) has highlighted the importance of reciprocal coupling between neuronal activity and metabolic resources in self-organizing and maintaining neuronal operation near criticality, and has also presented a general slow-fast formulation for the case where resources change slowly relative to neural activity. In (Virkar et al., 2016), a discrete two-layer model has been proposed to describe a mechanism by which metabolic resources are distributed to neurons via glial cells. An example of frequency adaptation in Kuramoto model was provided in (Taylor et al., 2010), reproducing certain phenomena that are not qualitatively accounted for the classical Kuramoto model, such as long waiting times before reaching synchronization. In (Nicosia et al., 2017), neuronal dynamics and nutrient transport were assumed to be bidirectionally coupled, such that the allocation of the transport process at one layer depends on the degree of synchronization in the other and vice versa. In (Kroma-Wiley et al., 2021), a system of coupled Kuramoto oscillators that consume or produce resources depending on their oscillation frequency was considered.

Inspired by the mechanisms for the interaction of a neuronal network with a population of glial cells, the studies (Fields, 2015; Lücken et al., 2017; Park and Lefebvre, 2020) introduced models of networks with adaptive time-delays.

Of particular interest are adaptive networks in which connectivity changes are related to intrinsic nodal dynamics (Gross and Blasius, 2008; Berner, 2021). For example, these types of networks can model synaptic neuronal plasticity (Meisel and Gross, 2009; Markram et al., 2011), chemical (Jain and Krishna, 2001; Kuehn, 2019), epidemic (Gross et al., 2006), biological, and social systems (Horstmeyer and Kuehn, 2020). A paradigmatic example of adaptively coupled phase oscillators gained considerable interest recently (Gutiérrez et al., 2011; Kasatkin et al., 2017; Berner et al., 2019a; Berner et al., 2019b; Berner et al., 2020; Feketa et al., 2020; Berner et al., 2021). This type of phase oscillator models seems to be useful for predicting and describing phenomena in more realistic and detailed models (Popovych et al., 2015; Lücken et al., 2016; Röhr et al., 2019) as well as for the understanding of collective phenomena such as multicluster states (Berner et al., 2019a; Berner et al., 2019b) or recurrent synchronization (Thiele et al., 2022).

In the present paper, we consider coupled excitable units (Lindner et al., 2004; Izhikevich, 2007), characterized by a linearly stable rest state susceptible to finite-amplitude perturbations. Excitable systems act as nonlinear threshold-like elements, such that applying a sufficiently small perturbation gives rise to a small-amplitude linear response, while a perturbation exceeding a certain threshold may trigger a large-amplitude nonlinear response. A classical example for the excitability feature are neurons (Ermentrout and Kopell, 1986; Izhikevich, 2007) which respond to a supra-threshold stimulation by emitting a spike. Apart from neuronal systems, excitability is important for other living cells (Scialla et al., 2021), lasers (Yanchuk et al., 2019; Terrien et al., 2021), chemical reactions (Chigwada et al., 2006), machine learning (Ceni et al., 2019), and many other fields. A variety of phenomena, including resonances, oscillations, patterns and waves, are caused by the interplay of excitability and noise (Pikovsky and Kurths, 1997; Neiman et al., 1999; Pototsky and Janson, 2008; Franović et al., 2015; Bačić et al., 2018a; Bačić et al., 2018b; Franović et al., 2018; Zheng and Pikovsky, 2018; Bačić and Franović, 2020; Franović et al., 2020) or time-delay (Brandstetter et al., 2010; Klinshov et al., 2016).

As a prototype of excitable local dynamics, we consider active rotators, paradigmatic for type I excitability (Shinomoto and Kuramoto, 1986; Park and Kim, 1996; Lindner et al., 2004; Osipov et al., 2007; Dolmatova et al., 2017; Franović et al., 2020; Klinshov et al., 2021). Active rotators have been used to study interacting excitable systems with noise (Lindner et al., 2004), synchronization in the presence of noise (Shinomoto and Kuramoto, 1986; Park and Kim, 1996; Dolmatova et al., 2017; Klinshov et al., 2021), the interplay of noise and an adaptive feedback (Franović et al., 2020), effects of an adaptive network structure (Thamizharasan et al., 2021), co-effects of noise, coupling, and adaptive feedback (Bačić et al., 2018b; Song et al., 2020) or delayed feedback (Yanchuk et al., 2019) and the impact of higher-order Fourier modes (Ronge and Zaks, 2021), to name but a few.

An important ingredient of our model is the multiscale structure of the dynamics, whereby the processes at the pool of resources are assumed to occur much slower than the dynamics of excitable units at the nodes. Utilizing this feature, we apply the methods of singular perturbation theory (Desroches et al., 2012; Kuehn, 2015) to first study the fast dynamics (layer dynamics) for fixed resource levels with the Ott-Antonsen approach, and then reduce the problem to the slow dynamics of resources.

Our main result consists in demonstrating how the adaptive interaction between a population of excitable units with a pool of resources gives rise to collective activity bursting. Such emergent dynamics is characterized by alternating episodes of stationary and oscillating behavior of the macroscopic order parameter. We describe the mechanisms behind the activity bursting and indicate parameter regions where this phenomenon can be reliably observed. So far, collective bursting phenomena have been considered to emerge due to time-varying neuronal inputs (Stoop et al., 2002), the interplay of external input and homeostatic plasticity (Zierenberg et al., 2018), or synaptic short-term plasticity (Gast et al., 2020). In these studies, possible implications for healthy and diseased brain states have been drawn. Moreover, the important role of bursting phenomena for the understanding of brain-organ interactions have been highlighted in the perspectives article (Ivanov, 2021). Our study complements recent research on emergent bursting dynamics in brain and organ systems by providing a simple and analytically tractable model generating collective activity bursting.

Our paper is organized as follows. In Section 2 we lay out the model of a heterogeneous population of excitable units adaptively coupled to a pool of resources, while in Section 3 we introduce the main phenomenon of collective activity bursting. Section 4 and Section 5 concern the analysis of the system’s multiscale dynamics within the framework of singular perturbation theory, first elaborating on the layer problem and then using the reduced problem to explain the mechanism of collective bursting and the origin of multistability in the full system. Section 6 proposes two different approaches to induce switches between the coexisting collective regimes, whereas Section 7 provides our concluding remarks and outlook.

2 Model

We consider a system of N coupled active rotators (Strogatz, 1994) with a Kuramoto-type coupling given by,

${\dot{\varphi }}_k=I_k\left(\mathit{r}\right)-\sin {\varphi }_k+\frac{\sigma }{N}\sum _{j=1}^N\sin \left({\varphi }_j-{\varphi }_k\right),(1)$

where ϕ_k ∈ [0, 2π), k = 1, … , N are the local phase variables, and σ is the coupling strength. While providing a simplified description of local dynamics, active rotators manifest the excitability feature crucial to neuronal activity (Strogatz, 1994; Izhikevich, 2007), and are similar to the model of theta neurons (Luke et al., 2013; Laing, 2014) paradigmatic for type I neural excitability. Note that more detailed models of neuronal dynamics, such as those of Morris-Lecar (Morris and Lecar, 1981) and Wang-Buzsáki (Wang and Buzsáki, 1996), also belong to this excitability class. External inputs I_k(r(t)) = r₁ (t) + r₂ (t) ν_k received by each unit comprise of a homogeneous component r₁ (t), acting identically at all the units, and a heterogeneous component, where the variability is due to parameters ν_k drawn from a normalized Gaussian distribution ${\nu }_k\in \mathcal{N}(\mathrm{0,1})$. Recall that in models of coupled active rotators, terms I_k are classically interpreted as local bifurcation parameters describing individual oscillation frequencies. Nevertheless, here I_k (t) at each moment follow a Gaussian distribution $g(I)=\mathcal{N}(r_1,r_2^2)$, such that the local velocities of the units are modulated by coupling to r₁ and r₂. The latter modulation can be seen as describing an interaction with the resources from the environment (Song et al., 2020; Kroma-Wiley et al., 2021) summarized within the two-component resource variable r = (r₁, r₂). In the context of neuroscience such modulation of local velocities is reminiscent of frequency adaptation of neuronal spiking (Fuhrmann et al., 2002; Ha and Cheong, 2017) due to a limited amount of metabolic resources affecting e.g. neurotransmitters.

Adaptation of spiking activity is a slow process compared to spike emission (Ha and Cheong, 2017), which should be reflected in the dynamics of metabolic resources r (t). In fact, a model involving such a separation of time scales has recently been proposed to describe the interplay of energy consumption and activity in neuronal populations (Roberts et al., 2014). Here we introduce a simple model of dynamical resources based on the Hopf normal form. We consider r as a complex variable, i. e, r = r₁ + ir₂, which satisfies the dynamical equation

$\dot{\mathit{r}}=ϵf\left(\mathit{r}-\mathit{s},\mathrm{\lambda }\right),(2)$

$\dot{\mathrm{\lambda }}=-{ϵ}^{\prime }\left(\mathrm{\lambda }-{\mathrm{\lambda }}_0-\gamma A\left(t\right)\right),(3)$

with activity

$A\left(t\right)≔\frac{1}{N}\sum _{j=1}^N{\dot{\varphi }}_j,(4)$

The metabolism describing function given by f (r, λ) = r (λ + iω − |r|²), the frequency ω and the resource base level given by s = s₁ + is₂. Small parameters ϵ ≪ 1 and ϵ′ ≪ 1 are introduced to account for the scale separation between the fast spiking dynamics of units and the slowly adapting dynamics of the resources. Note that we consider the case ϵ′ = ϵ throughout the paper.

System Eq. 2, 3 that describes the interaction between the two resources undergoes a supercritical Hopf bifurcation at λ = 0. This allows for the interpretation of the resource dynamics as being inactive if λ < 0, when it possesses a stable focus at s, or as active if λ > 0, when it displays a stable limit cycle. In other words, in the inactive states, the resource dynamics lies stationary at the resource base level s = s₁ + is₂, while for active states, the resource dynamics is attracted to a periodic orbit that encircles the resource base level. We further assume that the dynamics of metabolic resources adapts to the activity A (t) of the population, see Eq. 4. In particular, the adaptation dynamics Eq. 3 can be regarded as a feedback mechanism whereby due to a feedback loop, an activated neuronal population may activate the pool of resources which in turn may further activate or even deactivate the neuronal population. The adaptation strength is described by parameter γ which controls the impact of the population’s dynamics on the dynamics of resources. Throughout the paper, we keep γ = 0.5. In case of no spiking activity, i.e., if A (t) = 0 or γ = 0, the resource dynamics is inactive and the corresponding resource activity variable λ settles to the rest level λ₀. In the remainder, the level λ₀ = −0.05 is assumed to correspond to a stable steady state at s. Due to the dynamical interplay between the metabolic resources and the neuronal population, the activity variable λ (t) may change in time. Accordingly, the state of the resources may change between active (periodic attractor) and inactive (stationary state). To describe the coherence of the population dynamics, we use the complex order parameter Z defined by

$Z\left(\mathit{\varphi }\left(t\right)\right)=\frac{1}{N}\sum _{j=1}^N{e}^{\text{i}{\varphi }_j\left(t\right)}=R\left(\mathit{\varphi }\left(t\right)\right)e^{\text{i}\mathrm{\Theta }\left(\mathit{\varphi }\left(t\right)\right)},(5)$

where R is the Kuramoto order parameter, and Θ is the mean phase (Bick et al., 2020).

Summarizing, we have proposed a multiscale model of a heterogeneous population of active rotators, featuring local excitability and spike frequency adaptation as two important ingredients of typical neuronal activity, coupled to a pool of resources that slowly adjusts its dynamics to the activity of the population. Figure 1 provides an illustration of our model.

FIGURE 1

FIGURE 1. Schematic for the two-layer model consisting of a heterogeneous population of excitable units (green) and interacting pool of resources described by an adaptive Stuart-Landau oscillator (purple). The heterogeneity ν_i of the excitable units are randomly drawn from a distribution n (ν).

3 Collective Activity Bursting

In this section, we briefly introduce the phenomenon of collective activity bursting induced by an adaptive coupling to resources. A more detailed analysis of the phenomenon will be performed in the subsequent sections.

In Figure 2, we show a simulation result of a system consisting of N = 5,000 active rotators adaptively coupled to a pool of resources as described by Eqs 1–3. The emergent collective dynamics within the population is represented by the macroscopic variables A (t) and R (t). The dynamics within the resource pool is characterized by the activity variable λ (t). We observe that the population of active rotators displays a recurrent temporal formation of bursts in the macroscopic activity A (t) followed by periods of inactivity. Such episodes of macroscopic activity and inactivity correspond to episodes of a rapidly and slowly varying order parameter R (t), respectively. Switching between the different regimes is equally well visible in the evolution of the resource variable λ (t) showing the pattern of recurrent activation (λ > 0) and deactivation (λ < 0).

FIGURE 2

FIGURE 2. Collective activity bursting in system Eqs 1–3. Three panels show the time traces of the population activity A (t) (green), the order parameter R (t) (red) and the resource activity variable λ (t) (blue) from left to right, respectively. The trajectory is obtained from a random initial condition for a system of N = 5,000 active rotators and parameters: σ = 5, ϵ = 0.05, s₁ = 0.97, s₂ = 1.2, ω = 0.2, λ₀ = −0.05, γ = 0.5.

We note that this recurrent switching between macroscopic activity and inactivity is due to the adaptive feedback provided by the dynamical resources and can not be observed in a system of active rotators alone. In fact, active rotators are a paradigmatic model for excitable systems, supporting regimes of either activity ${\dot{\varphi }}_i>0$ or inactivity ${\dot{\varphi }}_i=0$ depending on parameters such as the input currents I_i, see e.g. (Franović et al., 2020) for more details. The slow adaptation of the input currents caused by the resource dynamics, however, provides a mechanism to switch between the two regimes. In the following sections, we systematically describe the emergence of collective activity bursting by making use of the separation of timescales between the dynamics of the population and the resources. The slow-fast analysis within singular perturbation theory, see e.g. (De Maesschalck and Wechselberger, 2015; Kuehn, 2015), allows for a splitting of multiscale dynamics into a so-called layer dynamics of the fast variables and an averaged dynamics for the slow variables.

The layer dynamics of system Eqs 1–3 consists of a population of actively coupled rotators with input currents drawn for a Gaussian distribution $\mathcal{N}(r_1,r_2^2)$. The subsequent analysis of the layer equation in Section 4 provides us with a clear mapping for the regimes of population activity and inactivity. Building on this, we analyse the full system Eqs 1–3 and show that the collective activity bursting emerge close to criticality, i.e., the boundary between activity and inactivity of the layer dynamics. We also describe regimes of multistability between activity bursting and inactivity, and provide insights into perturbations that give rise to transitions between different states.

4 Layer Dynamics: Heterogeneous Population of Active Rotators

The fast subsystem, describing the evolution of the original slow-fast problem Eqs 1–3 on the fast timescale, comprises of a heterogeneous assembly of N globally coupled active rotators

$\dot{{\varphi }_k}=r_1+r_2{\nu }_k-\sin {\varphi }_k+\frac{\sigma }{N}\sum _{j=1}^N\sin \left({\varphi }_j-{\varphi }_k\right).(6)$

In the absence of adaptation of the resource variables r₁ and r₂, the local dynamics ${\dot{\varphi }}_k=I_k-\sin {\varphi }_k$ depends on external input I_k = r₁ + r₂ν_k which may be seen as an effective bifurcation parameter mediating the transition between an excitable (I_k ≲ 1) and oscillatory regime (I_k > 1) via a SNIPER (saddle-node infinite period) bifurcation at |I_k| = 1. In the singular limit ϵ → 0, system (Eq. 6) defines the layer problem, where r₁ and r₂ are treated as additional system parameters.

According to classical singular perturbation theory (De Maesschalck and Wechselberger, 2015; Kuehn, 2015), the layer problem describes solutions of the multiscale system Eqs 1–3 on a timescale much shorter than 1/ϵ, where the variables r₁ and r₂ do not change significantly. In particular, it can describe fast (rapidly changing) segments of the solutions.

4.1 Ott-Antonsen Approach for the Layer Dynamics

We analyze the layer problem by determining the stability of stationary solutions of the layer dynamics and their bifurcations within the framework of Ott-Antonsen theory (Ott and Antonsen, 2008; Ott and Antonsen, 2009). We start by rewriting the layer dynamics in terms of complex order parameter (Eq. 5), which leads to

${\dot{\varphi }}_k=I_k-\sin {\varphi }_k+\sigma \mathrm{I}\mathrm{m}\left(Z\left(t\right)e^{-i{\varphi }_k}\right).(7)$

In the thermodynamic limit N → ∞, the state of the population can be described by the probability density h (ϕ, I, t), which satisfies the normalization condition ${\int }_0^{2\pi }h(\varphi ,I,t)\mathrm{d}\varphi =g(I)$, see e.g. (Omel’chenko and Wolfrum, 2012; Omel’chenko and Wolfrum, 2013). The continuity equation for h (ϕ, I, t) then reads

$\frac{\partial h}{\partial t}+\frac{\partial }{\partial \varphi }\left(hv\right)=0,(8)$

where the velocity is given by v = I − sin ϕ + σIm (Z(t)e^−iϕ). According to Ott-Antonsen ansatz (Ott and Antonsen, 2008; Ott and Antonsen, 2009), the long-term dynamics of Eq. 8 settles onto an invariant manifold of the form

$h\left(\varphi ,I,t\right)=\frac{g\left(I\right)}{2\pi }\left\{1+\sum _{n=1}^{\infty }\left[{\bar{z}}^n\left(I,t\right)e^{\text{i}n\varphi }+z^n\left(I,t\right)e^{-\mathrm{i}n\varphi }\right]\right\},(9)$

where z (I, t) is the local order parameter, connected with the global complex order parameter (Eq. 5) via

$Z\left(t\right)={\int }_{-\infty }^{\infty }g\left(I\right)z\left(I,t\right)dI.(10)$

Inserting Eq. 9 into Eq. 8, one obtains the Ott-Antonsen equation for the layer dynamics

$\dot{z}=\frac{1}{2}\left(1-z^2\right)+iIz+\frac{\sigma }{2}Z-\frac{\sigma }{2}\bar{Z}{z}^2,(11)$

where bar denotes the complex conjugate.

4.2 Stationary Solutions of the Layer Dynamics

To find stationary solutions of Eqs 10, 11, we first write the local order parameter in polar form z (I, t) = ρ(I, t)e^iϑ(I,t). Separating for the real and imaginary parts, Eq. 11 becomes

$\begin{array}{cc}\hfill \dot{\rho }& =\frac{1}{2}\left(1-{\rho }^2\right)B\cos \mathrm{\Phi },\hfill \\ \hfill \rho \dot{\mathrm{\Phi }}& =I\rho -\frac{1}{2}\left(1+{\rho }^2\right)B\sin \mathrm{\Phi },\hfill \end{array}(12)$

where the new variables B, β and Φ are given by

$\begin{array}{cc}\hfill B\left(t\right)e^{i\beta \left(t\right)}& =1+\sigma R\left(t\right)e^{i\mathrm{\Theta }\left(t\right)},\hfill \\ \hfill \mathrm{\Phi }& =\vartheta -\beta .\hfill \end{array}(13)$

From Eq. 13, it follows that B and β are related with the macroscopic order parameter Eq. 5 via

$\begin{array}{cc}\hfill B& =\sqrt{1+{\sigma }^2{R}^2+2\sigma R\cos \mathrm{\Theta }},\hfill \\ \hfill \tan \beta & =\frac{\sigma R\sin \mathrm{\Theta }}{1+\sigma R\cos \mathrm{\Theta }}.\hfill \end{array}(14)$

Note that the local dynamics can be rewritten in terms of B as ${\dot{\varphi }}_k=I_k-B\sin ({\varphi }_k-\beta )$, suggesting that B may be understood as an effective excitability parameter that describes how local excitability is changed by the impact of interactions. As a consequence, the structure of stationary solutions of the Ott-Antonsen system Eq. 12 depends on the relation between |I_k| and B, such that a population splits into two groups comprised of excitable (|I| < B) or oscillating units (|I| > B). In particular, the stationary solutions (ρ^∗, Φ^∗) are given by

$\begin{array}{cc}\hfill \left({\rho }^{\ast },{\mathrm{\Phi }}^{\ast }\right)& =\left(1,\arcsin \frac{I}{B}\right),\hfill \\ \hfill \left({\rho }^{\ast },{\mathrm{\Phi }}^{\ast }\right)& =\left(1,\pi -\arcsin \frac{I}{B}\right),\hfill \end{array}(15)$

for the excitable (inactive) group, and

$\left({\rho }^{\ast },{\mathrm{\Phi }}^{\ast }\right)=\left(\frac{|I|-\sqrt{I^2-B^2}}{B},\frac{\pi }{2}\mathrm{sign}\left(I\right)\right)(16)$

for the oscillating (active) group. An explicit expression for B can be obtained by invoking the self-consistency relation between the global and local order parameter (Eq. 10). Inserting the results for the stationary local and global order parameter [using Eq. 10, Eq. 15, Eq. 16, and the first equation from Eq. 13] and separating for the real and imaginary parts, one ultimately arrives at a self-consistency equation for B (Lafuerza et al., 2010; Klinshov and Franović, 2019)

$p\left(B\right)=B^2-2\sigma p_2\left(B\right)+\frac{{\sigma }^2}{B^2}\left(p_1^2\left(B\right)+p_2^2\left(B\right)\right)-1=0,(17)$

where p₁(B) and p₂(B) are given by

$\begin{array}{ll}\hfill p_1\left(B\right)& =r_1-\underset{|I|>B}{\int }Ig\left(I-r_1\right)\sqrt{1-{\left(\frac{B}{I}\right)}^2}\mathrm{d}I,\hfill \\ \hfill p_2\left(B\right)& =\underset{|I|<B}{\int }g\left(I-r_1\right)\sqrt{B^2-I^2}\mathrm{d}I.\hfill \end{array}(18)$

Having determined B, the stationary local and global order parameters can be obtained using the relations $R=\sqrt{p_1^2+p_2^2}/B$ and Θ = arctan (p₁/(p₂ − σR²)), which follow from Eqs 10–14.

For a fixed coupling strength σ, the function p(B) may have from one to three roots, depending on the mean value r₁ and the standard deviation r₂ of the distribution of intrinsic parameters I_k. The examples in Figure 3 illustrate how the number of solutions of Eq. 17 changes between one and three for fixed σ = 5, r₁ = 0.97 under increasing r₂. We refer to the stationary solutions by the corresponding B values, which we arrange in decreasing order B₁ > B₂ > B₃. Recalling the arguments above, one sees that the larger B value implies a prevalence of excitable over oscillating units within the local structure of the stationary state. This is evinced by the left column of Figure 4 which shows the dependence of the local order parameter z(I). Typically, the state B₁ comprises of a clear majority of excitable units, corresponding to a coherent domain z = 1, and may thus be referred to as a homogeneous stationary state. The two remaining stationary states B₂ and B₃ are heterogeneous in the sense that they involve a mixture of excitable and asynchronously oscillating units.

FIGURE 3

FIGURE 3. Changes in form and the number of roots of the function p (B) given by Eq. 17 under variation of r₁ and r₂ for fixed σ = 5. The function p (B) has three roots for (r₁, r₂) = (0.9, 2) (blue line; roots indicated by letters) and a single root for (r₁, r₂) = (1.1, 2) (red) and (r₁, r₂) = (0.9, 2.25) (green).

FIGURE 4

FIGURE 4. Local structure and spectra of stationary solutions B₁ (a), B₂ (b) and B₃ (c) of Ott-Antonsen equation Eq. 11 for σ = 5 and (r₁, r₂) = (0.9, 2). (A) shows the dependencies of the local order parameter on the input z (I) (black solid lines) and the corresponding Kuramoto order parameter R (blue dash-dotted lines) for the three stationary solutions. Red dashed lines indicate the interval (r₁ − 3r₂, r₁ + 3r₂) relevant for the distribution of external inputs. (B) shows the continuous (black dots) and the discrete spectra (red crosses) for the stationary solutions: B₁ and B₂ are stable and unstable nodes, respectively, while B₃ is an unstable focus.

4.3 Stability and Bifurcation Analysis of Stationary Solutions

Given that Ott-Antonsen equation Eq. 11 contains both the global order parameter and its complex conjugate, stability and bifurcation analysis of the stationary solutions (Omel’chenko and Wolfrum, 2013; Klinshov and Franović, 2019) can be carried out by writing the local and global order parameters as z (I, t) = x (I, t) + iy (I, t), Z(t) = X(t) + iY(t) and separating for the real and imaginary parts. This results in the system

$\begin{array}{cc}\hfill \dot{x}& =F\left(x,y,X,Y\right)=\frac{1}{2}\left(1-x^2+y^2\right)-Iy+\frac{\sigma }{2}X-\frac{\sigma }{2}X\left(x^2-y^2\right)-\sigma xyY,\hfill \\ \hfill \dot{y}& =G\left(x,y,X,Y\right)=-xy+Ix+\frac{\sigma }{2}Y-\sigma xyX+\frac{\sigma }{2}Y\left(x^2-y^2\right),\hfill \end{array}(19)$

which can be linearized for variations ξ = (δx, δy)^T, Ξ = (δX, δY)^T around the stationary solution (x₀, y₀, X₀, Y₀), ultimately arriving at

$\frac{\mathrm{d}\xi }{\mathrm{d}t}=\widehat{P}\left(I\right)\xi \left(I,t\right)+\widehat{Q}\left(I\right)\mathrm{\Xi }\left(t\right),(20)$

where $\widehat{P}$ and $\widehat{Q}$ are the corresponding Jacobian matrices

$\widehat{P}=\left(\begin{array}{cc}\hfill \frac{\partial F}{\partial x}\hfill & \hfill \frac{\partial F}{\partial y}\hfill \\ \hfill \frac{\partial G}{\partial x}\hfill & \hfill \frac{\partial G}{\partial y}\hfill \end{array}\right),\widehat{Q}=\left(\begin{array}{cc}\hfill \frac{\partial F}{\partial X}\hfill & \hfill \frac{\partial F}{\partial Y}\hfill \\ \hfill \frac{\partial G}{\partial X}\hfill & \hfill \frac{\partial G}{\partial Y}\hfill \end{array}\right).$

Eq. 20 is augmented by the variational equation for Eq. 10:

$\mathrm{\Xi }\left(t\right)={\int }_{-\infty }^{\infty }g\left(I\right)\xi \left(I,t\right)\mathrm{d}I.(21)$

Assuming that the variations ξ (I, t) and Ξ (t) satisfy ξ (I, t) = ξ₀ (I)e^μt, Ξ (t) = Ξ₀e ^μt, systems Eq. 20 and Eq. 21 transform into

$\left(\widehat{P}\left(I\right)-\mu \widehat{\mathbb{I}}\right){\xi }_0\left(I\right)+Q\left(I\right){\mathrm{\Xi }}_0=0,{\mathrm{\Xi }}_0={\int }_{-\infty }^{\infty }g\left(I\right){\xi }_0\mathrm{d}I,(22)$

where $\widehat{\mathbb{I}}$ denotes the identity operator. From the general spectrum theory of linear operators (Omel’chenko and Wolfrum, 2013; Mirollo and Strogatz, 2007), it follows that the Lyapunov spectrum of Eq. 22 consists of a continuous and a discrete part. Here, the continuous spectrum turns out to be always stable or marginally stable, such that the stability of stationary solutions depends on the discrete spectrum. The latter can be determined by rewriting Eq. 22 in the form $\widehat{C}(\mu ){\mathrm{\Xi }}_0=0$, where

$\widehat{C}\left(\mu \right)=\widehat{\mathbb{I}}+{\int }_{-\infty }^{\infty }g\left(I\right){\left(\widehat{P}\left(I\right)-\mu \widehat{\mathbb{I}}\right)}^{-1}Q\left(I\right)\mathrm{d}I.(23)$

The discrete spectrum is then obtained by solving the characteristic equation $\det \widehat{C}(\mu )=0$ (Klinshov and Franović, 2019). An example of the discrete and continuous spectra calculated for the stationary states B₁, B₂ and B₃ at (r₁, r₂) = (0.9, 2) is provided in the right column of Figure 4.

4.4 Comparison Between Analysis and Numerics

The previous analysis allows for an analytic description of the existence and stability of stationary solutions in the limit of large populations (N → ∞). In particular, the bifurcation diagram for the Ott-Antonsen equation of the layer dynamics Eq. 11 in the (r₁, r₂) plane is organized around a co-dimension two cusp point, indicated by C in Figure 5 where the two branches of folds meet (black dashed lines). Both branches of folds are calculated by numerical continuation of the solutions of Eq. 17 using the software package BifurcationKit.jl (Veltz, 2020). The lower branch of folds which folds over for larger r₂ corresponds to annihilation and reemergence of a pair of equilibria, B₁ and B₂, whereby the former (latter) is always stable (unstable). For smaller r₂, crossing this branch either by enhancing r₁ or r₂ gives rise to long-period collective oscillations as a stable equilibrium B₁ vanishes by colliding with an unstable equilibrium B₂. The divergence of the oscillation period when approaching the curve indicates that it corresponds to a SNIPER bifurcation of the full system. For larger r₂, as the branch folds over, one observes the reappearance of a stable stationary state B₁, emerging in an inverse fold bifurcation together with an unstable equilibrium B₂. The upper branch of folds involves stationary states B₂ and B₃, such that they collide and disappear above the curve, where B₁ remains the only stable stationary state, cf. Figure 5.

FIGURE 5

FIGURE 5. Bifurcation diagram for the system of active rotators Eq. 1 in dependence on resource levels r₁ and r₂. For the simulation, we have chosen one set of random initial conditions for the phases and one set of parameters ν_k randomly drawn from a normalized Gaussian distribution $\mathcal{N}(\mathrm{0,1})$. Simulations comprise of 200 time units with activity averaged over the last 100 time units. Fold bifurcations involving stationary solutions B₁ and B₂ (lower branch) and B₂ and B₃ (upper branch) obtained from Eq. 17 are shown by black dashed lines that give rise to a cusp point marked with C. Existence of particular solutions and their stability, derived from the discrete spectrum of Eq. 23, are indicated by their corresponding letters and a circle, respectively, whereby the circle indicates a stable solution. Along the black dotted line, stationary solution B₃ changes its stability in a Hopf-like bifurcation. Other parameters: N = 5,000, σ = 5.

Note that apart from the fold bifurcation, the stability of B₃ is also affected by a Hopf-like bifurcation (black dotted line). Above the given curve, stability of B₃ is determined by a pair of complex conjugate eigenvalues which have the smallest negative real parts. However, crossing the curve, these eigenvalues merge with the imaginary axis and remain neutrally stable immediately below the curve, implying that the central manifold theorem associated to Hopf bifurcation cannot immediately be applied. Still, in close vicinity below the curve, starting from an initial condition corresponding to B₃ results in oscillations similar to a genuine scenario of Hopf bifurcation.

Using numerical continuation, we have verified that the described structure of bifurcation diagram for the layer dynamics remains qualitatively the same under variation of coupling strength σ. One only notes that for increasing σ, the branches of folds shift toward larger r₂, which corresponds to a higher diversity of external inputs.

Figure 5 further shows a comparison of the existence and stability conditions for the collective stationary states derived from Ott-Antonsen approach for the limit N → ∞ with simulations for a finite population of N = 5,000 active rotators with fixed resources r₁ and r₂. One observes that simulation results agree well with the fold bifurcation lines separating parameter regimes of low and high collective activity. The differences can be attributed to the finite size of assemblies considered in the simulations.

With Figure 6 we complement the analysis of the layer equation. In particular, we show how the dynamical regimes change in a wide range of parameters r₁ and r₂, and indicate the boundary (black dashed line) that separate parameter regions supporting stable stationary states from those admitting oscillatory states. We illustrate three different trajectories corresponding to qualitatively different collective regimes found by numerical analysis. For parameter pairs a and b, we observe the emergence of stationary states in accordance with the bifurcation analysis shown in Figure 5. In both cases, activity A (t) and the coherence measure R (t) settle to a constant value. We observe that with increasing r₂ (b to a) the activity level rises while the coherence level declines. For the parameter set c, there is no stable stationary state and we observe stable oscillations. The activity shows a regular, tonic-like spiking shape corresponding to an increase in the average activity. Meanwhile, variation of the order parameter causes its average value to decrease. In order to quantify the temporal variations of the order parameter, we also plot the difference max  R (t)−min R (t) for the considered average time interval. We observe that even though the activity level might be high, e.g. for r₁ > 1 close to r₂ = 0, the coherence within the population is not necessarily strongly varying. However, there are also regimes, e.g. for r₁ > 1 and r₂ > 2, where the order parameter varies strongly and covers almost the entire interval from 0 to 1. In this section, we have illustrated the stability regions of macroscopic stationary states in a heterogeneous population of active rotators. Numerically, we have also determined the values of resource parameters r₁ and r₂ where no stable stationary solutions exist. Using these insights, we are now able to qualitatively describe the phenomena in systems with a slow adaptation of the resources and the resource-dependent dynamics. The next section is devoted to explaining the emerging states of collective activity bursting.

FIGURE 6

FIGURE 6. Bifurcation diagram for the system of active rotators Eq. 1 in terms of resource levels (r₁, r₂). All simulations are carried out the same way as in Figure 5. Three diagrams at the (A) show the time-averaged values of activity A (t) (white-green), order parameter R(t) (yellow-red) and variations of the order parameter max R(t) − min R(t) (white-black). The black dashed line separates regions where the mean phase Θ of the complex order parameter Z features stationary or oscillating dynamics, respectively. The (B) show time traces of activity and order parameter for three parameter pairs (r₁, r₂): a—(0.9,2), b—(0.9,1), c—(1.1,2).

5 (Slow) Resource Dynamics and the Emergence of Multistability

The analysis of the layer dynamics in Section 4 provides insight on how the system evolves for constant resources r₁ and r₂. Due to the different timescales of the population (fast dynamics) and the pool of resources (slow dynamics), we can average (Sanders et al., 2007; Franović et al., 2020) the system Eqs 2, 3 as

$\dot{\mathit{r}}=f\left(\mathit{r}-\mathit{s},\mathrm{\lambda }\right),(24)$

$\dot{\mathrm{\lambda }}=-\mathrm{\lambda }+{\mathrm{\lambda }}_0+\rho ⟨A⟩,(25)$

where $⟨A⟩=\frac{1}{T}{\int }_0^{T}A(s)\mathrm{d}s$. Here we also assume ϵ′ = ϵ and rescale time t_new = ϵt_old. Note that the average activity shown in Figure 6 depends on the resource variables r, since the definition Eq. 4 implies A = r₁ − Im (Z). Hence, the system Eqs 24, 25 describes an effective three-dimensional coupled dynamics for the slow subsystem. A further analytical analysis of this system is beyond the scope of our study. However, we directly use the insight that the slow dynamics follows the average activity of the fast system to understand the emergence of collective activity bursting.

Figure 7 shows the trajectories of the resource variables r₁(t) and r₂ (t) for the collective bursting presented in Figure 6 along with the averaged values of population activity and the order parameter. We clearly see that the asymptotic orbit passes through both the regimes of an active and inactive population which explains the episodes of high and low activity in Figure 2A. Also the segments of increasing and decreasing average activity visible in Figure 2A can be explained by Figure 7. Here, the average activity shows the same pattern along the trajectory (r₁ (t), r₂ (t)). The same also holds for the average values of the order parameter and even its variations if we compare Figure 2B with Figure 7. Therefore, the splitting of the fast from the slow subsystem provides a very good qualitative explanation for the observed phenomenon. To understand the emergence of the full periodic orbit shown in Figure 7, we first note that without any population activity, i.e., ⟨A⟩ = 0 and λ₀ = − 0.05, the resource dynamics possess a stable focus close to the critical line (black dashed line in Figure 7) describing the transition from stationary to oscillatory dynamics of the mean phase Θ. During the stationary phase, r₁ and r₂ tend to s₁ = 0.97 and s₂ = 1.2, respectively. As in Figure 7, the trajectory (r₁ (t), r₂ (t)) may start in the active region, i.e., oscillatory mean phase dynamics. Due to the positive average value of activity, the variable λ (t) characterizing the resource activity increases according to (Eq. 25) and becomes positive, see Figure 2C. Hence, the resources become activated and (r₁(t), r₂(t)) follows the limit cycle solution of the resource dynamics revolving around (s₁, s₂). Note that the resources obey the Hopf normal form with a Hopf bifurcation at λ = 0, see Eq. 2. After passing the critical line, the average activity immediately drops to ⟨A⟩ ≈ 0 which causes λ to tend to λ₀, see Figure 2C. After λ falls below zero, the dynamics of the resources (r₁(t), r₂(t)) is described by a spiral towards (s₁, s₂). This spiral, however, enters the active region by passing the critical line which ultimately leads to the recurrent phenomenon observed in Figure 2. As we have seen, the emergence of collective activity bursting relies on the subtle interplay between activation and deactivation of resources and the population. Furthermore, the need for the spiraling dynamics towards a stable focus explains well the necessity for the resource basis levels (s₁, s₂) to be close to the critical line separating the population active and inactive regimes.

FIGURE 7

FIGURE 7. Bifurcation diagrams as in Figure 6 complemented with the trajectories of the resource variables r₁ (t) and r₂ (t) for the collective activity bursting shown in Figure 6.

With regards to the above description of the collective activity bursting, one might ask for the coexistence of a stable steady state in the system Eqs 1–3 as long as (s₁, s₂) lie in the inactive regime. This state might have a small basin of attraction such that the spiral towards the steady state cannot reach the active regime.

In order to get insights into the different stable states that exist in Eqs 1–3, we use the numerical method of adiabatic continuation. To do so, we fix the base level s₂ = 1.2 and gradually vary s₁ from 0.8 to 1.4 (sweep up) and from 1.4 to 0.8 (sweep down). For each value of s₁, we run the simulation starting from the final state of the previous simulation. In Figure 8, we show the results of both sweeps. We observe the existence of stable steady and stable oscillating states for various values of s₁. As expected, close to the boundary between active and inactive states of layer dynamics, we also find an interval of coexistence between collective activity bursting and stable steady states, see panels for (a) and (b) in Figure 8, respectively. For larger s₂, only the oscillatory state can be observed, which does not enter the inactive regime above a certain s₁, see panel (c) in Figure 8. Note that the character of the solution can be deduced from the maximal value of λ(t) on the averaging time interval. In particular, there is a stationary state only if max λ(t) < 0. In all other cases, there are time intervals where the trajectory of r diverges from the base level s and follows the periodic solution of Eq. 24.

FIGURE 8

FIGURE 8. Bifurcation diagram with respect to the resource base level s₁ for a system of active rotators with adaptive resource interaction Eqs 1–3. The (A) shows the results from two adiabatic continuations with step size Δs₁ = 0.002 from s₁ = 0.8 to s₁ = 1.4 (sweep up) and vice versa (sweep down). Sweeps up and down start from a stable stationary and a stable oscillatory state, respectively. For both sweeps are shown the average activity ⟨A(t)⟩ (green), average order parameter ⟨R(t)⟩ (red) and maximal resource activity λ (blue). The results were obtained by simulating (1)–(3) for 7,000 time units and taking the average over the last 5,000 time units. The branches corresponding to the two sweeps are marked by arrows. The black dashed lines indicate the value (s₁ ≈ 0.963) of the critical line shown in Figure 6. Three trajectories represented by the activity (green, first column), order parameter (red, middle column) and the resource activity (blue, last column) are shown in the (B) for (a,b) s₁ = 0.97 and (c) s₁ = 1.35. The panels in (a) and (b) represent the states found along the sweeps down and up, respectively. The simulations were performed using the same values of ν_k as in Figure 5. Parameters: N = 5,000, σ = 5, ϵ = 0.05, s₂ = 1.2, ω = 0.2, γ = 0.5.

From the arguments laid out in this section, we have seen that the mutual activation and deactivation between the neural population and the pool of resources close to criticality of layer dynamics induces a rich dynamical behavior. It is believed, particularly, that the human brain operates close to criticality (Chialvo, 2010; Haimovici et al., 2013; Yu et al., 2013; Cocchi et al., 2017; Wilting and Priesemann, 2019). Therefore, it is of major importance to understand the dynamics of neural populations in this regimes including the interaction with its environment. In the next section, we propose a simple mechanism which can induce a switch between coexisting macroscopic regimes.

6 Population Switching Dynamics Induced by Resource Activation and Inactivation

In the vicinity of the transition between the population inactivity and activity, we have observed collective activity bursting induced by an adaptive dynamical pool of resources. Moreover, this phenomenon emerges in a stable coexistence with a steady state. In this section, we consider two simple perturbation approaches that can induce a switch between these two functionally different states.

Figure 9 shows the results for two different perturbation approaches to system Eqs. 1–3. The first approach aims to induce a switch in population dynamics by an instantaneous resetting of the resource activity λ. In the second approach, we induce such a transition by maintaining the resource activity at a certain level for a certain period of time. The first approach works well for large resetting values of λ, see Figures 9A,C. Small values, however, would not be sufficient to induce the macroscopic regime switch. Furthermore, in case of an initial bursting state, eliciting the switch to a steady state depends on the moment at which the perturbation is applied. However, in our numerical simulations (not shown), we have always been able to induce a switching for sufficiently large resetting values of λ.

FIGURE 9

FIGURE 9. Two perturbation scenarios to induce a switch between an inactive steady state and collective activity bursting. Time traces of macroscopic activity A (t) and order parameter R (t) are shown green and red, respectively. In the panels (A) and (B) (C,D), we start from an initial steady state (bursting state). The first perturbation scenario is illustrated for the cases where the resource activity variable λ is set to λ = 20 [panel (A)] and λ = −5 [panel (C)] at t = 2000. The second perturbation scenario is demonstrated for the cases where the resource activity is kept fixed at λ = 1 [panel (B)] and λ = −0.5 [panel (D)] for a duration of 500 t. u. beginning at t = 2000 t. u. Simulations were run for s₁ = 0.97 and the remaining parameters fixed as in Figure 8.

Due to the functionally very different nature of the two stable states, there might be reasons to favor one over the other in light of potential applications in medicine. Therefore, it is of great interest to understand simple mechanisms that would induce a switch to the desired state. While the first perturbation approach provides such a mechanism, it still requires strong perturbations which might be undesirable for certain medical reasons, e.g. side effects. Therefore, we have proposed another perturbation approach that leads to a switch while keeping the reset level lower. For this approach, we have also been able to induce switches between a steady state and a bursting state in one or the other direction, see Figures 9B,D, with the advantage of having the resetting level of the resource activity much lower than for the first method.

In this section, we have proposed two simple perturbation approaches to induce a switch between the two functionally different macroscopic states of the full system which emerge near the transition in layer dynamics between the population activity and inactivity and due to an adaptive dynamical pool of resources. We note that the approaches we proposed are not the only way to induce macroscopic regime shifts. One might also think of perturbing the resource variables (r₁, r₂) or even the whole population. Thus, perturbation of the resource activity variable is perhaps the simplest but not the only approach possible.

7 Conclusion

We have investigated collective dynamics in a system of interacting excitable units coupled to a pool of resources with nontrivial dynamics. The feedback of the resources to the population of coupled excitable units has been realized by an adaptation of the individual units’ inputs, whereas in turn, the excitable population is capable of activating or deactivating the pool of resources depending on the population’s own activity. As a prototype of excitable local dynamics, we have considered active rotators. Following the ideas outlined by Roberts et al. (Roberts et al., 2014), we have assumed the processes at the pool of resources to occur much slower than the local dynamics of excitable units. As a consequence, we have ended up with a system featuring multiscale dynamics, allowing us to use the methods from singular perturbation theory (Desroches et al., 2012; Kuehn, 2015).

As our most important finding, we have reported on the phenomenon of collective activity bursting. The phenomenon is characterized by a recurrent switching between episodes of quiescence and episodes of activity bursts in the population of active rotators. To gain a better understanding of the emergence of collective activity bursting, we have made use of the explicit slow-fast timescale separation. In particular, we have divided the system dynamics into the fast layer dynamics of the population and the slow average dynamics of the resources.

Using the Ott-Antonsen approach, we have analyzed the stability and bifurcations of the stationary solutions of layer dynamics in the thermodynamic limit. For the population of active rotators with a heterogeneity given by a Gaussian distribution, we have derived a bifurcation diagram for the steady state solutions. The bifurcations of layer dynamics depending on the mean and the width of the Gaussian distribution have been corroborated by numerical simulations of a large ensemble of rotators. Doing so, we have determined the parameter regions admitting high or low (or even no) population activity and have obtained the critical lines separating these regions.

Taking the analysis of the layer problem into account, we have further analyzed how the slow averaged dynamics of the resources gives rise to a slow variation of the mean and width of the Gaussian distribution. We have observed the onset of collective activity bursting close to criticality where the population of active rotators undergoes a transition from an inactive to an active state. The emergence of collective bursting is due to a subtle interplay of co-activation and co-deactivation of the dynamical population of rotators and the pool of resources.

We have further found a region of bistability between collective activity bursting and an inactive steady state close to criticality of the layer dynamics. A similar observation has been also discussed in the context of collective bursting induced by synaptic short-term plasticity (Gast et al., 2020). Moreover, we have proposed two different perturbation methods that can trigger switches between coexisting macroscopic regimes. In particular, we have demonstrated that the regime shifts can be induced either by using instantaneous large perturbations or persistent perturbations of the resource activity.

In terms of theory, an important extension of our work could concern a further analytical study of the reduced slow-fast system governing the collective dynamics of the ensemble of excitable units and its interaction with the resources. For convenience, we summarize the reduced system here

$\begin{array}{ll}\hfill \dot{z}\left(I,t\right)& =\frac{1}{2}\left(1-z^2\left(I,t\right)\right)+iIz\left(I,t\right)+\frac{\sigma }{2}Z\left(t\right)-\frac{\sigma }{2}\bar{Z}\left(t\right)z^2\left(I,t\right),\hfill \\ \hfill \dot{\mathit{r}}\left(t\right)& =ϵf\left(\mathit{r}\left(t\right)-\mathit{s},\mathrm{\lambda }\left(t\right)\right),\hfill \\ \hfill \dot{\mathrm{\lambda }}\left(t\right)& =-ϵ\left[\mathrm{\lambda }\left(t\right)-{\mathrm{\lambda }}_0-\rho \left(r_1-\text{Im}\left(Z\left(t\right)\right)\right)\right],\hfill \end{array}$

with

$Z\left(t\right)=\int g\left(I\right)z\left(I,t\right)\mathrm{d}I.$

In a broader context, we have proposed a simple paradigmatic model to study the emergence of complex collective phenomena induced by a dynamically co-evolving pool of resources. The research on the impact of resource constraints on the dynamical regimes of populations of neurons or neuron-like units from the dynamical network perspective (Nicosia et al., 2017; Kroma-Wiley et al., 2021) has begun only recently. In our study, we have shown that even a simple model that includes nontrivial dynamical resources gives rise to the emergence of collective activity bursting close to criticality in a population of neuron-like excitable units. Our study underlines the potentially important role of resource constraints in the operating of the human brain that is often hypothesized to operate close to criticality. We have further shown that the collective activity bursting may stably coexist with a steady state. Either one of these regimes could be desirable or undesirable, which makes understanding of the control mechanisms to switch between the regimes highly important (Tang and Bassett, 2018). In this context, we have discussed two simple approaches that can successfully induce such regime shifts. Both approaches impose perturbations to the single activity variable of the resource pool and can thus be generalized to systems with even more complex dynamical resource pools.

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

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for submission.

Funding

The work of RB and SY was supported by the German Research Foundation DFG, Project Nos 411803875 and 440145547.

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.

Acknowledgments

IF acknowledges funding from the Institute of Physics Belgrade through the grant by the Ministry of Education, Science and Technological Development of the Republic of Serbia. We acknowledge support by the German Research Foundation (DFG) and the Open Access Publication Fund of Humboldt-Universität zu Berlin.

References

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

CrossRef Full Text | Google Scholar

Attwell, D., and Laughlin, S. B. (2001). An Energy Budget for Signaling in the Grey Matter of the Brain. J. Cereb. Blood Flow Metab. 21, 1133–1145. doi:10.1097/00004647-200110000-00001

PubMed Abstract | CrossRef Full Text | Google Scholar

Bačić, I., and Franović, I. (2020). Two Paradigmatic Scenarios for Inverse Stochastic Resonance. Chaos 30, 033123. doi:10.1063/1.5139628

PubMed Abstract | CrossRef Full Text | Google Scholar

Bačić, I., Klinshov, V. V., Nekorkin, V. I., Perc, M., and Franović, I. (2018a). Inverse Stochastic Resonance in a System of Excitable Active Rotators with Adaptive Coupling. EPL 124, 40004. doi:10.1209/0295-5075/124/40004

CrossRef Full Text | Google Scholar

Bačić, I., Yanchuk, S., Wolfrum, M., and Franović, I. (2018b). Noise-induced Switching in Two Adaptively Coupled Excitable Systems. Eur. Phys. J. Spec. Top. 227, 1077–1090. doi:10.1140/epjst/e2018-800084-6

CrossRef Full Text | Google Scholar

Berner, R., Sawicki, J., and Schöll, E. (2020). Birth and Stabilization of Phase Clusters by Multiplexing of Adaptive Networks. Phys. Rev. Lett. 124, 088301. doi:10.1103/PhysRevLett.124.088301

PubMed Abstract | CrossRef Full Text | Google Scholar

Berner, R., Vock, S., Schöll, E., and Yanchuk, S. (2021). Desynchronization Transitions in Adaptive Networks. Phys. Rev. Lett. 126, 028301. doi:10.1103/PhysRevLett.126.028301

PubMed Abstract | CrossRef Full Text | Google Scholar

Berner, R., Schöll, E., and Yanchuk, S. (2019a). Multiclusters in Networks of Adaptively Coupled Phase Oscillators. SIAM J. Appl. Dyn. Syst. 18, 2227–2266. doi:10.1137/18m1210150

CrossRef Full Text | Google Scholar

Berner, R., Fialkowski, J., Kasatkin, D., Nekorkin, V., Yanchuk, S., and Schöll, E. (2019b). Hierarchical Frequency Clusters in Adaptive Networks of Phase Oscillators. Chaos 29, 103134. doi:10.1063/1.5097835

PubMed Abstract | CrossRef Full Text | Google Scholar

Berner, R. (2021). Patterns of Synchrony in Complex Networks of Adaptively Coupled Oscillators (Cham: Springer). Springer Theses.

Bick, C., Goodfellow, M., Laing, C. R., and Martens, E. A. (2020). Understanding the Dynamics of Biological and Neural Oscillator Networks through Exact Mean-Field Reductions: a Review. J. Math. Neurosci. 10, 9. doi:10.1186/s13408-020-00086-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Boccaletti, S., Latora, V., Moreno, Y., Chavez, M., and Hwang, D. (2006). Complex Networks: Structure and Dynamics. Phys. Rep. 424, 175–308. doi:10.1016/j.physrep.2005.10.009

CrossRef Full Text | Google Scholar

Brandstetter, S., Dahlem, M. A., and Schöll, E. (2010). Interplay of Time-Delayed Feedback Control and Temporally Correlated Noise in Excitable Systems. Phil. Trans. R. Soc. A. 368, 391–421. doi:10.1098/rsta.2009.0233

PubMed Abstract | CrossRef Full Text | Google Scholar

Ceni, A., Ashwin, P., and Livi, L. (2019). Interpreting Recurrent Neural Networks Behaviour via Excitable Network Attractors. Cogn. Comput. 12, 330–356. doi:10.1007/s12559-019-09634-2

CrossRef Full Text | Google Scholar

Chialvo, D. R. (2010). Emergent Complex Neural Dynamics. Nat. Phys 6, 744–750. doi:10.1038/nphys1803

CrossRef Full Text | Google Scholar

Chigwada, T. R., Parmananda, P., and Showalter, K. (2006). Resonance Pacemakers in Excitable media. Phys. Rev. Lett. 96, 244101. doi:10.1103/physrevlett.96.244101

PubMed Abstract | CrossRef Full Text | Google Scholar

Cocchi, L., Gollo, L. L., Zalesky, A., and Breakspear, M. (2017). Criticality in the Brain: A Synthesis of Neurobiology, Models and Cognition. Prog. Neurobiol. 158, 132–152. doi:10.1016/j.pneurobio.2017.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

De Maesschalck, P., and Wechselberger, M. (2015). Neural Excitability and Singular Bifurcations. J. Math. Neurosci. 5, 16. doi:10.1186/s13408-015-0029-2

CrossRef Full Text | Google Scholar

Desroches, M., Guckenheimer, J., Krauskopf, B., Kuehn, C., Osinga, H. M., and Wechselberger, M. (2012). Mixed-mode Oscillations with Multiple Time Scales. SIAM Rev. 54, 211–288. doi:10.1137/100791233

CrossRef Full Text | Google Scholar

Dolmatova, A. V., Goldobin, D. S., and Pikovsky, A. (2017). Synchronization of Coupled Active Rotators by Common Noise. Phys. Rev. E 96, 062204. doi:10.1103/PhysRevE.96.062204

PubMed Abstract | CrossRef Full Text | Google Scholar

Ermentrout, G. B., and Kopell, N. (1986). Parabolic Bursting in an Excitable System Coupled with a Slow Oscillation. SIAM J. Appl. Math. 46, 233–253. doi:10.1137/0146017

CrossRef Full Text | Google Scholar

Feketa, P., Schaum, A., and Meurer, T. (2020). Synchronization and Multi-Cluster Capabilities of Oscillatory Networks with Adaptive Coupling. IEEE Trans. Autom. Control. 66, 3084.

Google Scholar

Fields, R. D. (2015). A New Mechanism of Nervous System Plasticity: Activity-dependent Myelination. Nat. Rev. Neurosci. 16, 756–767. doi:10.1038/nrn4023

PubMed Abstract | CrossRef Full Text | Google Scholar

Franović, I., Omel'chenko, O. E., and Wolfrum, M. (2018). Phase-sensitive Excitability of a Limit Cycle. Chaos 28, 071105. doi:10.1063/1.5045179

PubMed Abstract | CrossRef Full Text | Google Scholar

Franović, I., Yanchuk, S., Eydam, S., Bačić, I., and Wolfrum, M. (2020). Dynamics of a Stochastic Excitable System with Slowly Adapting Feedback. Chaos 30, 083109. doi:10.1063/1.5145176

PubMed Abstract | CrossRef Full Text | Google Scholar

Franović, I., Perc, M., Todorović, K., Kostić, S., and Burić, N. (2015). Activation Process in Excitable Systems with Multiple Noise Sources: Large Number of Units. Phys. Rev. E 92, 062912. doi:10.1103/physreve.92.062912

CrossRef Full Text | Google Scholar

Fuhrmann, G., Markram, H., and Tsodyks, M. (2002). Spike Frequency Adaptation and Neocortical Rhythms. J. Neurophysiol. 88, 761–770. doi:10.1152/jn.2002.88.2.761

PubMed Abstract | CrossRef Full Text | Google Scholar

Gast, R., Schmidt, H., and Knösche, T. R. (2020). A Mean-Field Description of Bursting Dynamics in Spiking Neural Networks with Short-Term Adaptation. Neural Comput. 32, 1615–1634. doi:10.1162/neco_a_01300

PubMed Abstract | CrossRef Full Text | Google Scholar

Gross, T., and Blasius, B. (2008). Adaptive Coevolutionary Networks: a Review. J. R. Soc. Interf. 5, 259–271. doi:10.1098/rsif.2007.1229

PubMed Abstract | CrossRef Full Text | Google Scholar

Gross, T., D’Lima, C. J. D., and Blasius, B. (2006). Epidemic Dynamics on an Adaptive Network. Phys. Rev. Lett. 96, 208701. doi:10.1103/physrevlett.96.208701

PubMed Abstract | CrossRef Full Text | Google Scholar

Gutiérrez, R., Amann, A., Assenza, S., Gómez-Gardeñes, J., Latora, V., and Boccaletti, S. (2011). Emerging Meso- and Macroscales from Synchronization of Adaptive Networks. Phys. Rev. Lett. 107, 234103. doi:10.1103/physrevlett.107.234103

PubMed Abstract | CrossRef Full Text | Google Scholar

Ha, G. E., and Cheong, E. (2017). Spike Frequency Adaptation in Neurons of the central Nervous System. Exp. Neurobiol. 26, 179–185. doi:10.5607/en.2017.26.4.179

PubMed Abstract | CrossRef Full Text | Google Scholar

Haimovici, A., Tagliazucchi, E., Balenzuela, P., and Chialvo, D. R. (2013). Brain Organization into Resting State Networks Emerges at Criticality on a Model of the Human Connectome. Phys. Rev. Lett. 110, 178101. doi:10.1103/physrevlett.110.178101

PubMed Abstract | CrossRef Full Text | Google Scholar

Horstmeyer, L., and Kuehn, C. (2020). Adaptive Voter Model on Simplicial Complexes. Phys. Rev. E 101, 022305. doi:10.1103/PhysRevE.101.022305

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanov, P. C. (2021). The New Field of Network Physiology: Building The Human Physiolome. Front. Net. Physiol. 1, 1. doi:10.3389/fnetp.2021.711778

CrossRef Full Text | Google Scholar

Izhikevich, E. M. (2007). Dynamical Systems in Neuroscience. Cambridge, MA: MIT Press.

Google Scholar

Jain, S., and Krishna, S. (2001). A Model for the Emergence of Cooperation, Interdependence, and Structure in Evolving Networks. Proc. Natl. Acad. Sci. 98, 543–547. doi:10.1073/pnas.98.2.543

CrossRef Full Text | Google Scholar

Kasatkin, D. V., Yanchuk, S., Schöll, E., and Nekorkin, V. I. (2017). Self-organized Emergence of Multilayer Structure and Chimera States in Dynamical Networks with Adaptive Couplings. Phys. Rev. E 96, 062211. doi:10.1103/PhysRevE.96.062211

PubMed Abstract | CrossRef Full Text | Google Scholar

Klinshov, V., and Franović, I. (2019). Two Scenarios for the Onset and Suppression of Collective Oscillations in Heterogeneous Populations of Active Rotators. Phys. Rev. E 100, 062211. doi:10.1103/PhysRevE.100.062211

PubMed Abstract | CrossRef Full Text | Google Scholar

Klinshov, V. V., Shchapin, D. S., Lücken, L., Yanchuk, S., and Nekorkin, V. I. (2016). Experimental Study of Jittering Chimeras in a Ring of Excitable Units. AIP Conf. Proc. 1738, 210007. doi:10.1063/1.4951990

CrossRef Full Text | Google Scholar

Klinshov, V. V., Zlobin, D. A., Maryshev, B. S., and Goldobin, D. S. (2021). Effect of Noise on the Collective Dynamics of a Heterogeneous Population of Active Rotators. Chaos 31, 043101. doi:10.1063/5.0030266

PubMed Abstract | CrossRef Full Text | Google Scholar

Kroma-Wiley, K. A., Mucha, P. J., and Bassett, D. S. (2021). Synchronization of Coupled Kuramoto Oscillators under Resource Constraints. Phys. Rev. E 104, 014211. doi:10.1103/PhysRevE.104.014211

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuehn, C. (2015). Multiple Time Scale Dynamics. Cham: Springer.

Google Scholar

Kuehn, C. (2019). Multiscale Dynamics of an Adaptive Catalytic Network. Math. Model. Nat. Phenom. 14, 402. doi:10.1051/mmnp/2019015

CrossRef Full Text | Google Scholar

Lafuerza, L. F., Colet, P., and Toral, R. (2010). Nonuniversal Results Induced by Diversity Distribution in Coupled Excitable Systems. Phys. Rev. Lett. 105, 084101. doi:10.1103/PhysRevLett.105.084101

PubMed Abstract | CrossRef Full Text | Google Scholar

Laing, C. R. (2014). Derivation Of a Neural Field Model from a Network of Theta Neurons. Phys. Rev. E Stat. Nonlin Soft Matter Phys. 90, 010901. doi:10.1103/PhysRevE.90.010901

PubMed Abstract | CrossRef Full Text | Google Scholar

Levina, A., Herrmann, J. M., and Geisel, T. (2007). Dynamical Synapses Causing Self-Organized Criticality in Neural Networks. Nat. Phys 3, 857–860. doi:10.1038/nphys758

CrossRef Full Text | Google Scholar

Lindner, B., García-Ojalvo, J., Neiman, A. B., and Schimansky-Geier, L. (2004). Effects of Noise in Excitable Systems. Phys. Rep. 392, 321–424. doi:10.1016/j.physrep.2003.10.015

CrossRef Full Text | Google Scholar

Lücken, L., Popovych, O. V., Tass, P. A., and Yanchuk, S. (2016). Noise-enhanced Coupling between Two Oscillators with Long-Term Plasticity. Phys. Rev. E 93, 032210. doi:10.1103/PhysRevE.93.032210

PubMed Abstract | CrossRef Full Text | Google Scholar

Lücken, L., Rosin, D. P., Worlitzer, V. M., and Yanchuk, S. (2017). Pattern Reverberation in Networks of Excitable Systems with Connection Delays. Chaos 27, 013114. doi:10.1063/1.4971971

PubMed Abstract | CrossRef Full Text | Google Scholar

Luke, T. B., Barreto, E., and So, P. (2013). Complete Classification of the Macroscopic Behavior of a Heterogeneous Network of Theta Neurons. Neural Comput. 25, 3207–3234. doi:10.1162/neco_a_00525

PubMed Abstract | CrossRef Full Text | Google Scholar

Markram, H., Gerstner, W., and Sjöström, P. J. (2011). A History of Spike-timing-dependent Plasticity. Front. Synaptic Neurosci. 3, 4. doi:10.3389/fnsyn.2011.00004

PubMed Abstract | CrossRef Full Text | Google Scholar

Meisel, C., and Gross, T. (2009). Adaptive Self-Organization in a Realistic Neural Network Model. Phys. Rev. E Stat. Nonlin Soft Matter Phys. 80, 061917. doi:10.1103/PhysRevE.80.061917

PubMed Abstract | CrossRef Full Text | Google Scholar

Mirollo, R., and Strogatz, S. H. (2007). The Spectrum of the Partially Locked State for the Kuramoto Model. J. Nonlinear Sci. 17, 309–347. doi:10.1007/s00332-006-0806-x

CrossRef Full Text | Google Scholar

Morris, C., and Lecar, H. (1981). Voltage Oscillations in the Barnacle Giant Muscle Fiber. Biophysical J. 35, 193–213. doi:10.1016/s0006-3495(81)84782-0

CrossRef Full Text | Google Scholar

Neiman, A., Schimansky-Geier, L., Cornell-Bell, A., and Moss, F. (1999). Noise-enhanced Phase Synchronization in Excitable media. Phys. Rev. Lett. 83, 4896–4899. doi:10.1103/physrevlett.83.4896

CrossRef Full Text | Google Scholar

Nicosia, V., Skardal, P. S., Arenas, A., and Latora, V. (2017). Collective Phenomena Emerging from the Interactions between Dynamical Processes in Multiplex Networks. Phys. Rev. Lett. 118, 138302. doi:10.1103/physrevlett.118.138302

PubMed Abstract | CrossRef Full Text | Google Scholar

Omel’chenko, O. E., and Wolfrum, M. (2013). Bifurcations in the Sakaguchi-Kuramoto Model. Physica D 263, 74. doi:10.1016/j.physd.2013.08.004

CrossRef Full Text | Google Scholar

Omel’chenko, O. E., and Wolfrum, M. (2012). Nonuniversal Transitions to Synchrony in the Sakaguchi-Kuramoto Model. Phys. Rev. Lett. 109, 164101.

PubMed Abstract | Google Scholar

Osipov, G. V., Kurths, J., and Zhou, C. (2007). Synchronization in Oscillatory Networks. Berlin, Heidelberg: Springer.

Google Scholar

Ott, E., and Antonsen, T. M. (2009). Long Time Evolution of Phase Oscillator Systems. Chaos 19, 023117. doi:10.1063/1.3136851

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 18, 037113. doi:10.1063/1.2930766

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, S. H., and Kim, S. (1996). Noise-induced Phase Transitions in Globally Coupled Active Rotators. Phys. Rev. E 53, 3425–3430. doi:10.1103/physreve.53.3425

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, S. H., and Lefebvre, J. (2020). Synchronization and Resilience in the Kuramoto white Matter Network Model with Adaptive State-dependent Delays. J. Math. Neurosci. 10, 16. doi:10.1186/s13408-020-00091-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Pikovsky, A. S., and Kurths, J. (1997). Coherence Resonance in a Noise-Driven Excitable System. Phys. Rev. Lett. 78, 775–778. doi:10.1103/physrevlett.78.775

CrossRef Full Text | Google Scholar

Popovych, O. V., Xenakis, M. N., and Tass, P. A. (2015). The Spacing Principle for Unlearning Abnormal Neuronal Synchrony. PLoS ONE 10, e0117205. doi:10.1371/journal.pone.0117205

PubMed Abstract | CrossRef Full Text | Google Scholar

Pototsky, A., and Janson, N. (2008). Excitable Systems with Noise and Delay, with Applications to Control: Renewal Theory Approach. Phys. Rev. E Stat. Nonlin Soft Matter Phys. 77, 031113. doi:10.1103/PhysRevE.77.031113

PubMed Abstract | CrossRef Full Text | Google Scholar

Roberts, J. A., Iyer, K. K., Vanhatalo, S., and Breakspear, M. (2014). Critical Role for Resource Constraints in Neural Models. Front. Syst. Neurosci. 8, 154. doi:10.3389/fnsys.2014.00154

PubMed Abstract | CrossRef Full Text | Google Scholar

Röhr, V., Berner, R., Lameu, E. L., Popovych, O. V., and Yanchuk, S. (2019). Frequency Cluster Formation and Slow Oscillations in Neural Populations with Plasticity. PLoS ONE 14, e0225094. doi:10.1371/journal.pone.0225094

PubMed Abstract | CrossRef Full Text | Google Scholar

Ronge, R., and Zaks, M. A. (2021). Emergence and Stability of Periodic Two-Cluster States for Ensembles of Excitable Units. Phys. Rev. E 103, 012206. doi:10.1103/PhysRevE.103.012206

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanders, J. A., Verhulst, F., and Murdock, J. (2007). Averaging Methods in Nonlinear Dynamical Systems. New York, NY: Springer.

Google Scholar

Scialla, S., Loppini, A., Patriarca, M., and Heinsalu, E. (2021). Hubs, Diversity, and Synchronization in FitzHugh-Nagumo Oscillator Networks: Resonance Effects and Biophysical Implications. Phys. Rev. E 103, 052211. doi:10.1103/PhysRevE.103.052211

PubMed Abstract | CrossRef Full Text | Google Scholar

Shinomoto, S., and Kuramoto, Y. (1986). Phase Transitions in Active Rotator Systems. Prog. Theor. Phys. 75, 1105–1110. doi:10.1143/ptp.75.1105

CrossRef Full Text | Google Scholar

Song, T., Kim, H., Son, S. W., and Jo, J. (2020). Synchronization of Active Rotators Interacting with Environment. Phys. Rev. E 101, 022613. doi:10.1103/PhysRevE.101.022613

PubMed Abstract | CrossRef Full Text | Google Scholar

Stoop, R., Blank, D., Kern, A., v.d. Vyver, J.-J., Christen, M., Lecchini, S., et al. (2002). Collective Bursting in Layer IV. Cogn. Brain Res. 13, 293–304. doi:10.1016/s0926-6410(01)00123-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Strogatz, S. H. (1994). Nonlinear Dynamics and Chaos. 1st ed. Cambridge, MA: Perseus Books.

Google Scholar

Strogatz, S. H. (2001). Exploring Complex Networks. Nature 410, 268–276. doi:10.1038/35065725

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, E., and Bassett, D. S. (2018). Colloquium: Control Of Dynamics in Brain Networks. Rev. Mod. Phys. 90, 031003. doi:10.1103/revmodphys.90.031003

CrossRef Full Text | Google Scholar

Taylor, D., Ott, E., and Restrepo, J. G. (2010). Spontaneous Synchronization of Coupled Oscillator Systems with Frequency Adaptation. Phys. Rev. E Stat. Nonlin Soft Matter Phys. 81, 046214. doi:10.1103/PhysRevE.81.046214

PubMed Abstract | CrossRef Full Text | Google Scholar

Terrien, S., Pammi, V. A., Krauskopf, B., Broderick, N. G. R., and Barbay, S. (2021). Pulse-timing Symmetry Breaking in an Excitable Optical System with Delay. Phys. Rev. E 103, 012210. doi:10.1103/PhysRevE.103.012210

PubMed Abstract | CrossRef Full Text | Google Scholar

Thamizharasan, S., Chandrasekar, V. K., Senthilvelan, M., Berner, R., Schöll, E., and Senthilkumar, D. V. (2021). Exotic States Induced by Co-evolving Connection Weights and Phases, arXiv:2111.09861

Google Scholar

Thiele, M., Berner, R., Tass, P. A., Schöll, E., and Yanchuk, S. (2022). Asymmetric Adaptivity Induces Recurrent Synchronization in Complex Networks, arXiv2112.08697. submitted

Google Scholar

Veltz, R. (2020). BifurcationKit.jl. URL https://hal.archives-ouvertes.fr/hal-02902346.

Google Scholar

Vergara, R. C., Jaramillo-Riveri, S., Luarte, A., Moënne-Loccoz, C., Fuentes, R., Couve, A., et al. (2019). The Energy Homeostasis Principle: Neuronal Energy Regulation Drives Local Network Dynamics Generating Behavior. Front. Comput. Neurosci. 13, 49. doi:10.3389/fncom.2019.00049

PubMed Abstract | CrossRef Full Text | Google Scholar

Virkar, Y. S., Shew, W. L., Restrepo, J. G., and Ott, E. (2016). Feedback Control Stabilization of Critical Dynamics via Resource Transport on Multilayer Networks: How Glia Enable Learning Dynamics in the Brain. Phys. Rev. E 94, 042310. doi:10.1103/PhysRevE.94.042310

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X.-J., and Buzsáki, G. (1996). Gamma Oscillation by Synaptic Inhibition in a Hippocampal Interneuronal Network Model. J. Neurosci. 16, 6402–6413. doi:10.1523/jneurosci.16-20-06402.1996

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilting, J., and Priesemann, V. (2019). 25 Years of Criticality in Neuroscience - Established Results, Open Controversies, Novel Concepts. Curr. Opin. Neurobiol. 58, 105–111. doi:10.1016/j.conb.2019.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Yanchuk, S., Ruschel, S., Sieber, J., and Wolfrum, M. (2019). Temporal Dissipative Solitons in Time-Delay Feedback Systems. Phys. Rev. Lett. 123, 053901. doi:10.1103/PhysRevLett.123.053901

PubMed Abstract | CrossRef Full Text | Google Scholar

Yanchuk, S., Roque, A. C., Macau, E. E. N., and Kurths, J. (2021). Dynamical Phenomena in Complex Networks: Fundamentals and Applications. Eur. Phys. J. Spec. Top. 230, 2711–2716. doi:10.1140/epjs/s11734-021-00282-y

CrossRef Full Text | Google Scholar

Yu, S., Yang, H., Shriki, O., and Plenz, D. (2013). Universal Organization of Resting Brain Activity at the Thermodynamic Critical point. Front. Syst. Neurosci. 7, 42. doi:10.3389/fnsys.2013.00042

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, C., and Pikovsky, A. (2018). Delay-induced Stochastic Bursting in Excitable Noisy Systems. Phys. Rev. E 98, 042148. doi:10.1103/physreve.98.042148

CrossRef Full Text | Google Scholar

Zierenberg, J., Wilting, J., and Priesemann, V. (2018). Homeostatic Plasticity and External Input Shape Neural Network Dynamics. Phys. Rev. X 8, 031018. doi:10.1103/physrevx.8.031018

CrossRef Full Text | Google Scholar