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.

We consider a system of N coupled active rotators (Strogatz, 1994) with a Kuramoto-type coupling given by, ϕ ̇ k = I k r − sin ϕ k + σ N ∑ j = 1 N sin ϕ j − ϕ k , (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 ν k ∈ N ( 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 ) = 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 r ̇ = ϵ f r − s , λ , (2) λ ̇ = − ϵ ′ λ − λ 0 − γ A t , (3) with activity A t ≔ 1 N ∑ j = 1 N ϕ ̇ 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 ϕ t = 1 N ∑ j = 1 N e i ϕ j t = R ϕ t e i Θ ϕ t , (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.

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

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 ϕ ̇ i > 0 or inactivity ϕ ̇ 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 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.

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 ϕ k ̇ = r 1 + r 2 ν k − sin ϕ k + σ N ∑ j = 1 N sin ϕ j − ϕ k . (6)

In the absence of adaptation of the resource variables r ₁ and r ₂, the local dynamics ϕ ̇ k = I k − sin ϕ 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.

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 r ̇ = f r − s , λ , (24) λ ̇ = − λ + λ 0 + ρ ⟨ A ⟩ , (25) where ⟨ A ⟩ = 1 T ∫ 0 T A ( s ) 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.

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.

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.

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 λ.

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.

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 z ̇ I , t = 1 2 1 − z 2 I , t + i I z I , t + σ 2 Z t − σ 2 Z ̄ t z 2 I , t , r ̇ t = ϵ f r t − s , λ t , λ ̇ t = − ϵ λ t − λ 0 − ρ r 1 − Im Z t , with Z t = ∫ g I z I , t 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.