Stasis in heterogeneous networks of coupled oscillators: discontinuous transition with hysteresis

We consider a heterogeneous ensemble of dynamical systems in R4 that individually are either attracted to fixed points (and are termed inactive) or to limit cycles (in which case they are termed active). These distinct states are separated by bifurcations that are controlled by a single parameter. Upon coupling them globally, we find a discontinuous transition to global inactivity (or stasis) when the proportion of inactive components in the ensemble exceeds a threshold: there is a first–order phase transition from a globally oscillatory state to global oscillation death. There is hysteresis associated with these phase transitions. Numerical results for a representative system are supported by analysis using a system-reduction technique and different dynamical regimes can be rationalised through the corresponding bifurcation diagrams of the reduced set of equations.


Introduction
In 2004, Daido and Nakanishi [1] drew attention to the fact that in heterogeneous ensembles of coupled dynamical systems, the collective behaviour could overtake the dynamics of any of the individual units. They considered systems that underwent a Hopf bifurcation so that by varying a parameter one has either limit cycles or damped oscillations. By changing the proportion of systems that were active (namely limit cycle dynamics) and inactive (fixed point behaviour) the collective dynamics of the mixture could either be globally sustained oscillations or global oscillation death. This latter behaviour they denoted ageing, but since this term has been used in earlier in a different sense in the context of the glass transition and other phenomena [2], we prefer to describe the transition to a global loss of oscillation also as a transition to stasis [3].
In an effort to understand the origin of this transition-which does not appear to be the consequence of a simple majority rule-a number of variations of the initial model [1] have been studied, by incorporating distributed time delays [4,5], by varying the network topology [6,7], or by introducing heterogeneity in the parameters [8,9] as well as by including noise [10]. The nature of the coupling between the oscillators-whether through similar or dissimilar variables [11,12]-has also been explored in some detail over the past few years. One conclusion that can be drawn from this large variety of studies is that the essential feature that is required in an ensemble that exhibits the transition to stasis or ageing is heterogeneity or diversity in the components or in the interactions: some form of dynamical frustration seems necessary. The transition itself has typically been observed to be continuous in earlier studies [1] although an abnormal route to stasis, with the phase space dynamics becoming unbounded prior to the continuous transition has also been noted in a recent study of coupled van der Pol oscillators [13].
Here we consider dynamical systems with more than two state variables so that the dynamics in the active regime can potentially be richer. Specifically we consider ensembles of nonlinear systems in R 4 wherein the transition from a static state to periodic oscillations is via two saddle-node (SN) bifurcations. We study heterogeneous ensembles, namely mixtures of active and inactive oscillators, and find that as the proportion of inactive oscillators is increased, there is a transition to stasis that is abrupt, namely first-order. In contrast to the continuous transitions to stasis that have typically been observed so far [1], this is also different from the related phenomenon of explosive death that occurs in a homogeneous ensemble of coupled oscillators. Such an abrupt transition, from collective oscillations to a global fixed point, has been seen in networks of Stuart-Landau oscillators for a specific coupling and network topology [14,15], as well as in ensembles of identical limit-cycle and chaotic oscillators coupled via mean-field diffusion [16,17]. While the eventual dynamics is similar, the underlying mechanism appears to be different. The basic oscillator system that we study here is the Stuart-Landau dimer [18], a dynamical system in 4-dimensions. This system of a pair of coupled Stuart-Landau oscillators has been studied earlier [19] and is known to undergo bifurcations from a fixed point to periodic motion; there are thus well-defined parameter regimes of active and inactive dynamics. The general framework of the present study is described in the next section. In section 3 we present a detailed bifurcation analysis of the Stuart-Landau dimer. Analytic and numerical results for mixed oscillator ensembles of Stuart-Landau dimers are presented in section 4. The consequences of a first-order transition to globally static solutions are discussed along with a summary of our results in section 5.

The model
Consider an ensemble of N globally coupled m−dimensional nonlinear systems. The dynamics is schematically represented asẊ where the superscript j = 1, 2, . . . , N is the oscillator index. The state variables of an individual system are For the jth oscillator, the bifurcation parameter is µ j , and in the absence of coupling, the flow of each oscillator is governed by the functions The coupling between oscillator j and k is specified by the function G(X j , X k ), and the coupling strength can be tuned via the parameter K.
In the uncoupled systems, there is a bifurcation in the dynamics at parameter value µ = µ * , with active systems having µ j < µ * and inactive ones having µ j > µ * . The ageing transition is studied by taking a fraction Np of the oscillators to be inactive; this set is denoted S I . The remaining oscillators are active, and this is the set S A . For convenience, all the oscillators in set S I have parameter µ I = µ * + α, while those in set S A have parameter µ A = µ * − α, symmetric about the bifurcation point.
When p = 1 all the systems are inactive for all values of the coupling K; each oscillator is at a fixed point. When all systems in the ensemble are active, namely when p = 0, for suitable coupling strength the individual oscillators synchronize [1,20]. For small p, the active elements in S A effectively drive the initially inactive ones in S I away from the fixed points. Within each set the respective oscillators have coherent dynamics as the coupling is increased, although global synchronization between S A and S I typically does not occur. With increasing p the proportion of inactive oscillators increases, and there is a transition to global inactivity (or stasis) at p = p c which depends on the strength of the coupling K.
The order parameter [1] where ⟨·⟩ is the time average, and i and j are the variable and oscillator indices respectively is useful in studying the transition as a function of p. The normalised order parameter Q is where M(0) corresponds to the value of the order parameter when all oscillators are active. The equations of motion were integrated numerically using the standard fourth-order Runge-Kutta method [21] with a sufficiently small step size of 0.01 up to 100 000 steps. The last 20 000 steps were considered for computing the time averages of various quantities after transients were reliably removed. The initial conditions for the state variables X were chosen from a uniform random distribution in the range (−2, 2) using a pseudo-random number generator, and the final state of the first computation is treated as the initial condition for the next one when p or K is changed to a nearest value. On a backward sweep, similarly, we start with the same attractor we left as the final state. We vary p and K adiabatically in steps of 0.05 or 0.01 and 0.1 respectively to obtain the phase diagrams.
In the next section we study the bifurcations of a Stuart-Landau dimer and the consequences of an ensemble of such dimers numerically and analytically in a mean-field approximation.

The Stuart-Landau dimer
The Stuart-Landau system is specified by the equation of motion, where z is a complex variable, ω is the frequency of oscillation. The bifurcation parameter a is such that the oscillators settle in their own limit cycle with radius √ a for a > 0, and to the origin when a < 0. The Stuart-Landau dimer is a pair of such oscillators coupled so as to give a system in R 4 [18,19] [ The internal coupling for the Stuart-Landau dimer is controlled by the parameter ε, which is the relevant bifurcation parameter for the system. The fixed points for the flow are the origin (0,0,0,0) and in addition the four points with x * 3 = −x * 1 and x * 4 = −x * 2 that exist when ε > ω. These fixed points are created via two simultaneous SN bifurcations at ε = ω; see figure 1(a). The limit cycles for the dimer can be obtained by switching to plane polar coordinates, replacing (z 1 , z 2 ) by (r 1 exp(ιθ 1 ), r 2 exp(ιθ 2 )).
The radii can be shown to be given by A detailed bifurcation diagram for the dimer (equation (5)) is shown in figure 1(a) for the case of ω = 2 and a = 1, so that the origin is always unstable. The saddle node bifurcations occur at ε = 2, below which there are two stable limit cycles, the time series of the variable x 1 for these limit cycles is shown in figure 1(b). They are distinguished as having a lower-amplitude, marked LC 1 (green) and a higher amplitude, LC 2 (blue) respectively. The solid green and blue lines in figure 1(a) trace the extrema of the variable x 1 in these two cases.
At ε = ε S = 2, two simultaneous SN bifurcations occur and LC 2 vanishes, leaving behind two pairs of stable and unstable fixed points. The stable fixed points are denoted S + and S − , and are indicated by the (black) solid lines, while the unstable fixed points are indicated by the (red) dotted lines. The low amplitude limit cycle LC 1 continues to exist beyond ε S , becoming unstable at ε = 2.5 when a Pitchfork-like (PF) bifurcation occurs [22]. At this point, LC 1 becomes unstable but the two stable fixed points S + and S − continue to exist.
To summarize, the SN and PF bifurcations divide the ε parameter range as follows. Below ε S there are two stable limit cycles and one unstable fixed point. In the range ε S < ε < ε P (shaded region), there is one limit cycle LC 1 coexisting with two stable equilibria and three unstable equilibria. Above ε P there are only two stable attractors in the system: S + and S − corresponding to an inhomogeneous steady state (IHSS) [23] as well as one unstable equilibria and an unstable limit cycle attractor. Their basins of attraction have smooth boundaries in all these regimes, and are shown in figures 1(c)-(e).  (5)) with ω = 2, as a function of the coupling strength ε [19]. The solid (green/blue) and dashed (green) lines indicate the extrema of the stable and unstable limit cycle oscillations respectively while the dotted lines (red) indicate the unstable fixed points. The locus of the stable fixed points is shown as a solid line (black) marked as S + and S − . The shaded region corresponds to the coexistence of oscillatory/active state and oscillation death. Two saddle-node bifurcations and a Pitchfork-like bifurcation occur at εS and εP respectively. (b) Time series of x1 at ε = 1.5 shows the existence of two limit cycles having different frequencies, LC2 (blue) and LC1 (green). The basins of attraction for the fixed points S + (black) and S − (grey), and limit cycle attractors LC2 (blue), and LC1 (green) are shown in (c) at ε = 1.5, (d) at ε = 2.25, (e) at ε = 2.75. These basin images are numerically generated with 6400 random initial conditions in the (x1, x2) in the region shown, keeping the initial values x3 = x4 = 0.

Ensemble of Stuart-Landau dimers
We study an ensemble of N Stuart-Landau dimers as discussed above, with G taken to be simple pairwise diffusive coupling, namely (X j − X k ). The set of initially active oscillators S A in the ensemble all have the bifurcation parameter ε A = ε S − α while the set S I are initially inactive and have ε I = ε S + α. Since there are two limit cycles at ε A the set of initial conditions also plays a crucial role in determining the collective state of the ensemble. The emergent collective behaviour of the ensemble results in all oscillators either being active or all of them being inactive.

Simulation results
Shown in figure 2(a), is a schematic phase diagram at α = 1, obtained from simulations with N = 50 Stuart-Landau dimers. In the weak coupling limit (when K ≈ 0.5) and low enough p, the set S A oscillate at one amplitude, while the set S I have small amplitude oscillations about the inhomogeneous steady states S + & S − : we term this region SL in figure 2(a). The time series of variable x j 1 corresponding these sets are shown in figure 2(b) in red and blue respectively. In this parameter window there are different limit cycles in the set S A as a result of dynamical frustration. When the values of p or K are increased further, the dynamics of the S A oscillators becomes more homogenous and synchronised. The spread σ [24] in the amplitudes of the limit cycles is related to the number of distinct limit cycles in the set S A , and the dependence on K for a range of different sets of initial conditions at a given p has a characteristic behaviour, an example of which is shown in figure 3. The value of K after which the spread is minimum (here roughly K ∼ 1.6) marks the change from mixed oscillations to synchrony, namely the region marked LL in figure 2(c). Note that the sets S A and S I are The dynamics is examined in detail as a function of p at three different values of K along the lines marked A, B, and C in figure 2(a). The normalised order-parameter Q is 1 when there is global activity in the ensemble, and 0 when there is stasis or global oscillation death. The transition between these limits is  . Along A and B, the transition to stasis is continuous whereas along C the transition is abrupt. The values of pc 1 and pc 2 remain 1 along A, whereas along B, pc 2 is 0.58 but pc 1 is initial conditions dependent and varies between 0.58 & 1. There is power law scaling along A & B with β ≈ 1 & 0.6 respectively. Along C, the abrupt jump to stasis occur at pc 1 ≈ 0.7 and is always initial condition dependent, whereas in backward sweep an abrupt transition from stasis to active occurs at pc 2 ≈ 0.48 which remains constant regardless of initial conditions. continuous along the lines A and B in figure 2(a), and as can be seen in figure 4(a), the order parameter obeys a power law scaling [1,11] as p → p c , β being the scaling exponent. This continuous transition to stasis occurs at low K values; the phase transition curves in figures 4(a) and (b) correspond to increasing/decreasing p along the lines A and B in figure 2(a) for K = 0.5 and 1 respectively. In the forward sweep, p is varied from 0 to 1 following the attractor generated at p = 0 and at a certain p c = p c1 the system attains stasis. Similarly, in the backward sweep p is varied from 1 to 0, the system being in either of the IHSS: S + or S − at p = 0 switches to an oscillatory state at p c = p c2 . In the weak coupling regime (K below 0.5, say) both p c1 and p c2 are the same, namely 1 but above K = 0.5, p c1 and p c2 have distinct values. p c1 remains constant at 0.88 while p c2 is initial condition dependent and varies between 0.58 & 1. There are hysteresis loops associated with Q due to the occurrence of multiple limit cycle attractors in set S A . The scaling exponent β ≈ 1 ± 0.009 for K = 0.5 ( figure 4(a)), decreasing to β ≈ 0.6 ± 0.048 for K = 1 ( figure 4(b)). At K = K c (see figure 2(a)) however, the power-law breaks down and the transition to stasis becomes discontinuous. Shown in figure 4(c) for K = 4, namely along the line marked C in figure 2(a) is the variation in Q: below p = p c1 , there are large amplitude global oscillations in the ensemble and Q is nearly 1 but at p c1 = 0.7 there is complete cessation of oscillation in the system with a sudden jump to stasis. Sweeping backwards from p = 1, activity resumes at a lower value of p, i.e. p c2 ≈ 0.48, the hysteresis loop being quite prominent. K c plays a crucial role in forming a boundary between first order and second order phase transitions in the system and its dependence on α (i.e. the distance from ε S ) is briefly discussed in the appendix B (see figure B.1).

Mean-field analysis
In order to understand the discontinuous transition we apply the system reduction technique [1] to equation (1). Note that all the oscillators in S A and S I are identical, so it is convenient to define the system through reduced mean-field active variable A and inactive variable I respectively defined as where i = 1, 2, 3, 4. For the active set, and correspondingly, for the inactive set: This reduces the set of 4N equations in equation (1) to the above eight equations. Although obtaining an analytical expression for the fixed points is difficult, bifurcation diagrams can be plotted for the reduced system as shown in figure 5. For K = 4 below p ≈ 0.5 the origin is the only fixed point. Since this is unstable there are oscillations about it (figure 5(e)) and at p = p S ≈ 0.48, two simultaneous SN bifurcations occur, leading to a pair of stable fixed points as well as a pair of unstable fixed points, shown as solid and dotted lines, respectively in , there are three coexisting stable states, namely the system has multistability. In addition to periodic motion there are two stable equilibrium states, see figure 5(f). At p = p P ≈ 0.64 the two unstable fixed points merge into the origin as a consequence of Pitchfork-like bifurcation, the periodic oscillation about the origin disappears, and the origin remains as an unstable fixed point with no periodic oscillation about it. Thus the bifurcation, which is similar to subcritical Pitchfork case marks the end of the multistable region. Since the stable fixed points continue to exist, in this region we observe that the system undergoes oscillation death. Details of the stability analysis about the corresponding fixed points are given in appendix A.
The resulting phase diagram based on the above analysis is presented in figure 5(g) and is in fairly good agreement with the numerically obtained figure 2(a). In the vicinity of the Hopf bifurcation, namely the red dashed curve separating the regions marked SL and ST in figure 5(g) the system goes through a continuous phase transition. These points can be identified as the boundary between regions of small amplitude oscillation and stasis in the corresponding numerically obtained result, figure 2(a).
At the SN bifurcation points, the nearly vertical dashed curve in figure 5(g) there is a discontinuous abrupt transition to stasis along the boundary between the state of large amplitude oscillations and stasis in figure 2(a). The multistable region marked as LST in the former figure corresponds to the shaded region in figure 5(b). The window of multistability gradually expands as α is decreased. Note that some of the features of figure 2(a), i.e. SST, do not appear in figure 5(g) as the later is a simplified model.
The single Stuart Landau system (equation (4)) has a Hopf bifurcation [1], while when two of these are coupled in an Stuart-Landau dimer (equation (5)), the Hopf bifurcation disappears and there are two other bifurcations, namely the saddle-node and Pitchfork-like bifurcations in the system. In an ensemble of dimers with global coupling, it is actually possible to get all the three bifurcations, depending on the parameter regimes as in figure 5(g). The critical point above which one has SN and below which there is Hopf bifurcation is identified as K c , are indicated in figures 2(a) and 5(g), and are both in good agreement.

Summary and discussion
Studies of the dynamics of heterogeneous ensembles of dynamical systems have practical relevance to a variety of problems in the physical and biological sciences. The models that have been studied so far in this context [1,11,13,20,26] have typically used a model oscillator with a single parameter controlling the dynamics: by varying the parameter one has either active (oscillatory) or inactive (fixed point) dynamics. When coupled, there are either global oscillations or global oscillation death and the transition to this latter state has been described as an ageing process which is also termed a transition to stasis [11].
There can be several different bifurcations that separate the active and inactive states of an individual system, and in earlier work systems undergoing the super or subcritical Hopf bifurcation [1,27], super or subcritical Pitchfork bifurcation [28] or the saddle-node bifurcation has been studied in the context of ageing [1,13,20]. While the type of bifurcation that separates the active and inactive regimes is clearly important in determining the characteristics of the transition from a global oscillatory state to global death, the correlation is not straightforward. For example, the study [20] of an ensemble of Morris-Lecar oscillators that have two time-scales corresponding to fast and slow manifolds [29] and in which the loss of oscillatory behaviour is via an SN homoclinic bifurcation [30] also has a continuous transition to stasis as in the case of the Hopf bifurcation [1] while complex networks of van der Pol oscillators also show a discontinuous transition to stasis [27].
Our aim in the present work has been to examine heterogeneous ensembles of higher dimensional oscillating units which have richer bifurcations as well as dynamics. By coupling two limit cycle oscillators, we obtained a dynamical system in R 4 , and found that ensembles of such dimers show a discontinuous, first-order ageing transition to stasis. In the weak coupling regime, the dynamical frustration makes the attractor set more diverse. Consequently, one has different hysteresis loops. When the coupling is sufficiently strong, the collective dynamics of the ensemble resembles that of an uncoupled oscillating unit and hence the jumps. The hysteresis in the strong coupling region, however, is due to the coexistence of steady states as well as oscillatory states. A mean-field analysis along the lines proposed earlier by Daido et al was carried out, and although only a numerical solution of the mean-field equations can be obtained, we find excellent agreement with numerical simulations.
It will be interesting to see whether this transition persists in the case when the parameters are chosen from a distribution so that the active and inactive elements are also different, and also in the situation when there is diversity in the types of oscillators & a variety of coupling scenarios other than the linear diffusive one considered here in the interacting ensemble. The increase in the dimensionality of the constituent units does not, in the present instance, bring new dynamical features into play, but it offers the possibility of new bifurcations, and new active dynamical states-such as chaotic attractors-that could potentially exhibit new collective behaviour. This is currently under investigation [31].

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).  (10) and (11) using the mean field analysis, namely Kc ≈ 2.8 for α = 1 figure 5(g). The data points obtained are within numerical error but the trend is similar.