Basin stability measure of different steady states in coupled oscillators

In this report, we investigate the stabilization of saddle fixed points in coupled oscillators where individual oscillators exhibit the saddle fixed points. The coupled oscillators may have two structurally different types of suppressed states, namely amplitude death and oscillation death. The stabilization of saddle equilibrium point refers to the amplitude death state where oscillations are ceased and all the oscillators converge to the single stable steady state via inverse pitchfork bifurcation. Due to multistability features of oscillation death states, linear stability theory fails to analyze the stability of such states analytically, so we quantify all the states by basin stability measurement which is an universal nonlocal nonlinear concept and it interplays with the volume of basins of attractions. We also observe multi-clustered oscillation death states in a random network and measure them using basin stability framework. To explore such phenomena we choose a network of coupled Duffing-Holmes and Lorenz oscillators which are interacting through mean-field coupling. We investigate how basin stability for different steady states depends on mean-field density and coupling strength. We also analytically derive stability conditions for different steady states and confirm by rigorous bifurcation analysis.


Results
We start with a network of coupled oscillators with the isolated dynamics of each node of the network is given by =  X F X ( ), where X is a m-dimensional vector of the dynamical variables and F(X) is the vector field. The general framework of coupled network is given by the following equation: where N is the total number of nodes in the network, ε is the coupling strength, C ij are the elements of connectivity matrix and H(X i , X j ) is the coupling function between i-th and j-th node.
Duffing-Holmes oscillator. We first consider a two-dimensional physical example, namely Duffing-Holmes (DH) oscillator 43 : 3 The oscillator has three steady states, namely two symmetrical stable steady states (± 1, 0) (which are spiral or node depending on the damping coefficient b > 0) and a saddle point at (0, 0) irrespective of the values of the parameter b. For b < 0 each individual DH oscillator exhibit oscillatory state. Recently, Tamaševičiūtė et al. 44 discussed the stabilization of saddle fixed points of an uncoupled DH oscillator using modified unstable filter 45 method. The proposed technique is applicable only for b > 0 where the DH oscillator is either stable node or spiral. But they did not discuss the stabilization of saddle point in coupled oscillators. Here we study the stabilization of saddle point of coupled systems by taking all values of damping parameter b. In this context, detection and Scientific RepoRts | 7:45909 | DOI: 10.1038/srep45909 controlling both saddle and nonsaddle types of unstable steady states in high-dimensional nonlinear dynamical systems based on fast-slow manifold separation and Markov chain theory is articulated in ref. 46.
We consider the coupled network through mean-field in the following form: for i = 1, … , N. Here ε is the mean-field coupling strength, d(i) is the degree of the i-th node and ≤ < Q Q (0 1) is the mean-field density parameter. This mean-field density parameter Q gives an additional free parameter that control the mean-field dynamics while Q → 0 represents self-feedback case and Q → 1 indicates the maximum mean-field density. The elements of the connectivity matrix A ij = 1 if i-th and j-th nodes are connected and zero otherwise. At first we consider a minimal network of two (N = 2) coupled Duffing-Holmes oscillators with mean field coupling and identify the parameter region for stabilized saddle point at origin. The coupled DH oscillator has a trivial steady state E 0 = (0,0,0,0) which is the HSS solution of the system and the other four coupling dependent steady states: non-trivial homogeneous steady state (NHSS) and inhomogeneous steady state (IHSS) and that stabilization of saddle point occurred through inverse pitchfork bifurcation. By performing the stability analysis we analytically obtain the inverse pitchfork bifurcation (IPB) point at the coupling strength From linear stability analysis we also analytically derive the Hopf bifurcation (HB) point at where up to this critical value of the coupling strength, coupled systems exhibit oscillatory states ( Fig. 1(a)). Further increment of the coupling strength leads to co-existence of IHSS and NHSS up to a certain threshold of interaction strength and after ε PB , IHSS are completely eliminated and only NHSS sustained up to ε IPB . So, using linear stability analysis and combining the above results, structurally different dynamical states occur: AD exist for ε ε > IPB , IHSS and NHSS (OD) coexist for ε ε ε < < HB PB and only NHSS exist for ε ε ε < < PB IPB . For numerical simulation, we choose the damping coefficient b = − 0.01 for which an isolated oscillator exhibits oscillatory dynamics. At lower value of coupling strength ε, four coupling dependent fixed points (i.e. NHSS and IHSS) that arise through Hopf bifurcation at ε = ε HB , are stable. But as ε increases, two of these stable steady states E 3,4 become unstable at ε PB and E 12 remains stable for the value of ε upto ε IPB . At ε IPB , saddle point (0,0,0,0) turns stable through IPB and remains stable for ε > ε IPB . The corresponding bifurcation diagram (using XPPAUT 47 ) is shown in Fig. 1(a). Figure 1(b) shows the bifurcation diagram with respect to coupling strength ε when b = 0.5, for which an isolated oscillator approaches either to the steady state (1, 0) or to (− 1, 0) where for negative values of b different coupling dependent stable steady states appear through oscillatory states as shown in Fig. 1(a). Here again, due to the introduction of coupling ε, above mentioned four fixed points E 1,2 and E 3,4 become stable but E 3,4 remain stable only upto ε PB . Further increment in the value of ε makes the saddle point (0, 0, 0, 0) stable through an inverse pitchfork bifurcation at ε IPB . Figure 1(c) shows how the BS of the steady states E 1,2 and E 3,4 change for different values of ε. As can be seen, initially after the occurrence of Hopf bifurcation at ε ε = HB , all the fixed points (E 1,2,3,4 ) are equally probable although the probabilistic dominance of E 3,4 are shrinking gradually whereas E 1,2 acquire more and more space in the basin volume. Such changes on BS measure of E 3,4 gives a hint of the annihilation of E 3,4 which finally occurs at ε ε = PB where E 1 and E 2 share the basin with equal probability. But at ε ε = IPB , BS of these steady states abruptly decrease to zero without any presage and further increase of ε, the basin volume is fully covered by this HSS shown by cyan color with maximum BS i.e., 1. From this figure, we conclude that the BS for multi stable states (i.e. OD and NHSS) change with the variation of mean-field coupling strength ε while the BS for monostable state i.e. AD state remains unchanged with the variation of ε. Therefore, the trend in the changes of the percentage of the basin volume gives us a clear idea how the different steady states are evolving in a coupled system and which states will dominate the system and which will disappear early. We also obtain similar results on stabilization of saddle point in two coupled DH oscillators when they are coupled through cross mean-field type configuration (see Supplementary Information section I). ..
, and (for even number oscillators) where α β γ , , and δ are same as above. The fixed points E 1,2 are same for any choice of N whereas E 3,4 are same only for even number of N.
The distinct eigenvalues and critical bifurcation points of globally connected network (3) are same as for two coupled oscillators.
First, we consider N = 4 i.e. four globally coupled DH oscillators via mean-field coupling and the results are shown in Fig. 2. Analytically it is not easy to calculate all the coupling dependent fixed points (i.e. IHSS and NHSS), using bifurcation diagram (performed in XPPAUT 47 ) we identify all the fixed points and by BS measurement we measure the amount of their stability for different values of coupling strength. In Fig. 2(a), bifurcation diagram for the variables x i with respect to the coupling strength ε is plotted. For small values of ε, through Hopf bifurcation at ε = , eight coupling dependent fixed points are stable. But as ε is increased, firstly six of these stable fixed points lose their stability, among them two lose their stability through PB at ε = and remaining four lose their stability at earlier of ∈ PB and only two retain their stability. Even more increment in ε makes the two fixed points unstable and the saddle point (i.e., the origin) becomes stable through IPB. The process of stabilization and destabilization of all the coupling dependent fixed points are clarified in terms of their BS which validates the whole mechanism in global scale. Figure 2(b) shows the variation of BS for different steady states by varying the mean-field coupling strength ε. As mentioned earlier, the blue and magenta color in Fig. 2(b) belong to class NHSS and they acquire more and more space in the basin if we increase the coupling strength continuously. On the other hand, the other cluster belonging to IHSS (six states in three symmetric groups) losing their stability and finally all of them vanish at ε ε = PB point. Further changes in ε makes those two fixed points (NHSS) equally probable in the basin i.e. each of them acquires half of the whole basin and they become unstable at the point ε ε = IPB . Then the saddle point (i.e. the origin) becomes stable for all points in the basin of attraction i.e. the basin volume is fully covered by this HSS.
Next we will verify numerically whether the stabilization of saddle and all the coupling dependent steady states using the proposed coupling scheme is working in a large network. For our case, we choose N = 1000 globally coupled DH oscillators via mean-field and the analyzed results are illustrated in Fig. 3. Figure 3(a) shows time series of x-components of all the 1000 oscillators with ε = 0.3 that depicts the stabilization of the IHSS resulting in OD. For larger value of ε (ε = 2.2), all the oscillators populate to a single steady state, that is, the saddle point (the origin) gets stabilized, time series are shown in Fig. 3(b). The corresponding space-time plots are shown in Fig. 3(c) and (d) respectively. The parameter space in ε − Q plane for global network is same as in Fig. 1(d) for two coupled DH oscillators, as the distinct eigenvalues of the characteristic equation (4) are identical with two coupled case but with different multiplicity.

Random network of Duffing-Holmes oscillators.
In this section we investigate the stabilization of saddle point in Erdös-Rényi random networks of DH oscillators and the results are shown in Fig. 4 where the probability of existence of an edge between any two vertices of the random network is taken as p = 0.01. Figure 4   Lorenz oscillators. For quantifying the different stable steady states using BS measure, we extend our investigation on coupled paradigmatic chaotic Lorenz oscillator 48 . We consider N Lorenz oscillators interacting through mean-field diffusive coupling. The mathematical equations of the coupled systems are described as: In absence of coupling term, each oscillators oscillate chaotically for σ = = r 10, 28 and = b 8 3 and the individual systems have a saddle fixed point at origin and two unstable fixed point at , 1). Here ε and Q are the coupling strength and mean-field density parameter respectively.
For N = 2, the fixed points are The trivial fixed point E 0 is stable through inverse pitchfork bifurcation at ε = .
For fixed values of the above system parameters, from eigenvalue analysis the NHSS points E 1,2 becomes stable The results are shown in Fig. 5. Figures 5(a) and (b) show bifurcation diagrams with respect to ε for N = 2 and N = 4 respectively with Q = 0.5 fixed. As in Fig. 5(a), due to the presence , six coupling dependent stable fixed points (comprising of both IHSS and NHSS states) emerge together with six unstable fixed points. But among them, the fixed points except the NHSS E 1,2 lose their stability soon and only E 1,2 remain stable for higher values of ε. Similarly as before, through IPB at ε ε = IPB , E 1 and E 2 collides and E 0 turns stable. Figure 5(c) shows the bifurcation diagram against ε for a network of N = 4 randomly connected nodes where the appearance of six coupling dependent stable fixed points along with many other unstable fixed points can be seen. Similarly as in the previous cases, here also E 1,2 retain their stability for higher values of ε than the others and lose it stability at ε IPB and further higher coupling strength promotes the entire systems to the AD state. Figures 5(d) and (e) measure all the stable steady states that appear for N = 2 and N = 4 respectively in terms of their BS. Figure 5(d) shows that the BS of both E 1 and E 2 are non-zero and more or less the same for all values of ε upto ε IPB . As ε increases further, BS of E 1 and E 2 turns into zero and BS of E 0 becomes unity. On the other hand, soon after the Hopf bifurcation all the six coupling dependent stable fixed points get non-zero BS but E 1 and E 2 have larger BS than the others, as in Fig. 5(e) (left part). Increasing ε, BS of the other fixed points become zero and E 1 and E 2 shares almost the same BS value upto ε ε = IPB . After that BS of both E 1,2 becomes zero and that of E 0 appears to be 1 (right part in Fig. 5(e)).
Finally, Fig. 5(f) depicts the parameter region in ε − Q plane for globally coupled N = 4 oscillators. Here blue, yellow, cyan and red regions signify oscillatory state, co-existence of OD and NHSS states, stable NHSS state and AD state (i.e., the stabilization of saddle E 0 ) respectively. The oscillatory state (blue region) and coexistence of OD and NHSS (yellow region) or stable NHSS (cyan region) are separated by the Hopf bifurcation curve ε =  Networks of Lorenz oscillators. Next we will explore the proposed coupling scheme is applicable for large number of chaotic oscillators. To quantify the stability of different steady states using BS measure in global and random network. The characteristic equation at E 0 of network (5) is Taking a network of N = 1000 globally coupled Lorenz oscillators with Q = 0.5, the numerical results are shown in Fig. 6. Figure 6(a) shows time evolution of the − x components of all the 1000 oscillators with ε = 5.5 that represents the stabilization of IHSS resembling OD. Whereas for ε = 30, the saddle point (origin) appears to be stable, time series shown in Fig. 6(b). Figure 6(c) and (d) depict the corresponding space-time plots respectively.
Results regarding MCOD state and saddle stabilization in Erdös-Rényi random networks of coupled Lorenz systems are given in Fig. 7. Time series of the − x components of all the N = 1000 oscillators revealing MCOD state for ε = 6.0 and Q = 0.5 are shown in Fig. 7(a). In Fig. 7(b) stable AD state ensuing after a concise transient window and Fig. 7(c) shows the corresponding space-time plots representing stabilization of the saddle point (the origin) reflecting AD for ε = .
30 0 and = . Q 0 5. The dependence of the BS of the MCOD and NHSS states on coupling strength ε for the random network of Lorenz systems is characterized in Fig. 7(d). Here, after the Hopf bifurcation at ε .  0 9, the MCOD state retains BS almost 1 and the BS of NHSS is very small upto ε .  4 3. In fact, for ε . < < . 4 3 1 3 5 the states MCOD and NHSS co-exist but then BS of NHSS develops with tantamount and that of the MCOD state vanishes at ε .  13 5 and NHSS remains stable further upto ε .  23 7 from where the saddle point becomes stable abruptly without any pre-warning and carries BS unity further.

Discussion
In this work we have studied basin stability (BS) measure to quantify the stability of different stable steady states of coupled dynamical systems interacting through mean-field coupling. BS is an universal concept to quantify the stability of governing dynamical systems under a non-uniform distribution of perturbations. Using mean-field coupling configuration, we have obtained a homogeneous stable steady state (i.e. AD state) which is inherently saddle equilibrium point of the individual oscillator and also showed that the transition from inhomogeneous steady states (resembling OD) to homogeneous steady state (i.e. AD state) via stabilization of NHSS state. We identify the underlying mechanism to stabilize the saddle fixed points in a network of coupled oscillatory systems. The transition routes between different states of coupled systems are discussed through rigorous bifurcation analysis and confirmed with the obtained analytical results. We also map the different steady states in the wide parameter space by varying the mean-field coupling strength ε and mean-field density parameter Q. All the steady states are quantified by the value of BS. In contrast to this, we found that the BS of OD states gradually decreases as coupling strength increases. After annihilation of the BS of multi-stable OD states, the BS of NHSS states become prevalent with almost equal ratio. But further increasing of coupling strength NHSS states become unstable without any presage and immediately AD state is stabilized. In the context of oscillations suppression studies, all the previous works have been done by considering the specific initial conditions in the phase space and no one examined the whole basin volume therefore ignoring the multistability nature of the steady states. As multistable character is ubiquitous in natural systems so we clearly elucidate a global stability measure by means of basin stability. To validate the BS measure, we have considered a large number of initial states following 16 . All these phenomena and measures are performed using smaller size of networks (for N = 2 and N = 4) as well as network of bigger size (N = 1000). We test our proposition and statistical measure not only in complete graph but also in random network. For both cases, our analytical and numerical simulations give proper insight to track the multistablity features present in the systems. The models considered here cover the characteristics of limit cycle (Duffing-Holmes oscillator) or chaotic attractor (Lorenz system) having hyperbolic fixed points. There are many real systems such as laser 49 and geomagnetic 50 which are modeled like Lorenz systems or mimic of Lorenz systems after some transformations and the results of our approach can be easily implemented.
Our considered mean-field coupling is one of the most natural coupling scheme which is previously extensively applied to different branches of science and engineering. This strongly means that our approach is not limited to a particular situation or for some particular systems, rather this mechanism is applicable in wide range of systems throughout all these disciplines. Also multistable feature is omnipresent in nature and widespread phenomenon in dynamical systems that appears in diverse fields ranging from physics, chemistry, biology to social systems 51 . There are numerous systems in which multistability originates that include the human brain, semiconductor materials, chemical reactions, metabolic system, arrays of coupled lasers, hydrodynamical systems, various ecological systems, artificial and living neural systems etc. We believe that this study will broaden our understanding of stabilization of saddle points in multistable dynamical networks where units are connected via mean-field. Further we have shown that the critical mean-field coupling strength is independent of the size of the network but only depends on the largest real part of the eigenvalue of individual oscillator (refers to Linear Stability Theorem in Method Section). is an asymptotically stable equilibrium point of the given system. Now let ⊂ B I be the basin of attraction of the stable state X k (i.e, the solution of the system starting from any ∈ z B asymptotically converges to X k as t → ∞ ).
We numerically integrate the given system for V points which are drawn uniformly at random (sufficiently large) from I. Let V k be the count of the initial conditions that finally arrives at the stable steady state X k . Then the BS for the fixed point X k is estimated as V , where λ* is the maximum real part of eigenvalues of the isolated system at the equilibrium point O and ≤ < Q Q (0 1) is the mean-field density parameter. Proof: Consider N identical systems interacting through global mean-field diffusive coupling as follows: where f(X i ) be the evolution equation of the i th system, X i denotes − m dimensional state vector, k be the mean-field coupling strength, Q is the mean-field density parameter and = ∑ .