Network resilience in the presence of non-equilibrium dynamics
Subhendu Bhandary, Taranjot Kaur, Tanmoy Banerjee, Partha Sharathi Dutta
aa r X i v : . [ n li n . AO ] A ug Network resilience in the presence of non-equilibrium dynamics
Subhendu Bhandary , Taranjot Kaur , Tanmoy Banerjee , ∗ and Partha Sharathi Dutta † Department of Mathematics, Indian Institute of Technology Ropar, Rupnagar 140 001, Punjab, India. Chaos and Complex Systems Research Laboratory, Department of Physics,University of Burdwan, Burdwan 713 104, West Bengal, India. (Dated: September 1, 2020; Received :to be included by reviewer)Many complex networks are known to exhibit sudden transitions between alternative steady stateswith contrasting properties. Such a sudden transition demonstrates a network’s resilience, whichis the ability of a system to persist in the face of perturbations. Most of the research on networkresilience has focused on the transition from one equilibrium state to an alternative equilibriumstate. Although the presence of non-equilibrium dynamics in some nodes may advance or delaysudden transitions in networks and give early warning signals of an impending collapse, it has notbeen studied much in the context of network resilience. Here we bridge this gap by studying a neu-ronal network model with diverse topologies, in which non-equilibrium dynamics may appear in thenetwork even before the transition to a resting state from an active state in response to environmen-tal stress deteriorating their external conditions. We find that the percentage of uncoupled nodesexhibiting non-equilibrium dynamics plays a vital role in determining the network’s transition type.We show that a higher proportion of nodes with non-equilibrium dynamics can delay the tipping andincrease networks’ resilience against environmental stress, irrespective of their topology. Further,predictability of an upcoming transition weakens, as the network topology moves from regular todisordered.
I. INTRODUCTION
Biological systems are often composed of nonlinear in-teractions and may exhibit complex multistable dynam-ics in response to various stressors, such as biochemicalchanges, environmental fluctuations, etc [1]. Empiricaland theoretical evidences suggest that many multistablesystems; including lakes [2], tropical forests [3], climate[4, 5], cellular systems [6, 7], global finance [8] can un-dergo sudden transitions to alternative stable states [9].These transitions occur when a system bypasses a criti-cal value, known as tipping points, due to the loss of re-silience to perturbations [10–14]. Determining such sud-den transitions is relatively easy when a species or anisolated system governs the state of a system. However,for complex networked systems where many nodes inter-act, e.g., genetic networks, neural networks, and ecolog-ical networks of interacting species – sudden transitionsare highly unpredictable [15, 16]. This is because currentanalytical frameworks of measuring resilience are mostlylimited to low-dimensional systems and are underdevel-oped for high-dimensional systems, where multiple com-ponents simultaneously interact through a complex net-work. Understanding, either the loss or the maintenanceof network resilience in the face of internal disturbancesor external environmental changes, is an important prob-lem with broad interest.Given a complex network, variations in the networktopology can act as internal stress, which may per-turb the system’s resilience when exposed to stress de-teriorating their external conditions (i.e., environmental ∗ [email protected] † Corresponding author; [email protected] stress) [17]. Investigations on how system parameters,such as growth rate and interaction strength, affect net-work resilience under changing environmental conditions,are primarily limited to mutualistic interactions [18, 19].However, real-world interactions are not restricted to mu-tualism but reveal diverse network architecture. Varioussocial interactions, biochemical pathways, and protein-protein interactions primarily exhibit non-homogeneousproperties that often vary in degree distribution with thepower-law formation or the scale-free interactions. Apartfrom diversity in network structure, the stability of com-plex networks is also accounted for by damages in thenetwork parameters such as changes in the number oflinks or nodes, coupling strength, average degree, etc .For example, [20] studied the robustness of communi-cation networks to loss of nodes in the disintegrationof the network from a well-functioning state to a dis-turbed state. [21] investigated the resilience to damagethe networks in graphs with degree correlations. Further,[22, 23] suggested that strength up to which a system canadapt changes to external conditions can increase by in-troducing substantial heterogeneity, noise, or weak cou-pling amongst the connected components. However, [24]reported consideration of perfectly regular interaction toinvestigate network resilience [23] leads to bias in the out-comes. Nevertheless, the investigations mentioned abovemostly neglect the role of complex local dynamics, whileinvestigating the interplay between external stress andnetwork topology in the occurrence of tipping points.While investigating the resilience patterns of complexnetworks, for mathematical tractability, it is appealing toconsider either one-dimensional or two-dimensional non-linear systems in each node with equilibrium steady statedynamics. However, this may overlook the concurrent ef-fects of non-equilibrium dynamics preceding a transition[25]. Such non-equilibrium/cyclic dynamics may gener-ate trade-off amongst other nodes, which can further sup-press or fuel the occurrence of tipping points in the inter-action network. In dynamical systems, cyclic dynamicsmostly occur due to negative feedback, and positive feed-back generates hysteresis loops, which creates the possi-bility of sudden transition of an equilibrium steady state.It would be interesting to study the resilience of a het-erogeneous network, where a fraction of its nodes canexhibit equilibrium steady state, and a fraction of nodescan show cyclic dynamics. This new direction bear rel-evance for the management of tipping points in manynetworks, like neuronal and ecological networks, whereoscillations are ubiquitous.Here we consider a network of FitzHugh-Nagumo(FHN) neurons [26, 27]. The FHN model representsa paradigmatic neuronal model and have widely beenused to study the collective behaviors of an ensembleof excitable as well as self excitatory spiking neurons.Several concepts of synchronization [28], coherence res-onance [29], and chimera states [30] have been discov-ered using networks of FHN neurons. Some of the emer-gent dynamics have been believed to have connectionswith observable physiological and cognitive behaviors.For example, the so-called bump state was proposed asthe mechanism behind the visual orientation tuning inthe rat’s brain [31]. The unihemispheric slow-wave sleepin aquatic mammals and some migratory birds is con-nected with the coherent-incoherent dynamics of coupledoscillating neurons [32]. Recently, it has been shown ex-plicitly that FHN oscillators on complex networks mimicepileptic-seizure-related synchronization phenomena [33].In this work, one of our main interests is to study theinfluence of diverse network topologies and environmen-tal stress on the FHN network resilience in the presenceof non-equilibrium dynamics [25]. We find that in a regu-lar FHN network with heterogeneity in nodes, sometimestransitions from one stable state to another is gradual(i.e., continuous/second-order phase transition). At thesame time, it can be sudden (i.e., discontinuous/first-order phase transition) for disordered networks. Themanagement of internal parameters, such as couplingstrength or average degree, can significantly influence thetype of transitions. When the network is dominated bycyclic dynamics in proportion of nodes, an equilibriumnetwork state precedes non-equilibrium dynamics. Wefind that a precise classification of transition type is notalways possible. We observe fluctuations before a tran-sition in both regular and disordered networks, underexternal stress, which indicates that the presence of non-equilibrium dynamics in a proportion of nodes plays asignificant role in determining network resilience. Fur-ther, we find that the efficacy of covariance-based earlywarning signals weakens, as the network topology movesfrom regular to disordered.
II. MATHEMATICAL MODELA. A single FHN neuron
We consider the following model of the FHN neuron: dxdt = 1 ǫ ( x − x − y ) , (1a) dydt = gx − y + b , (1b)where x is the membrane potential, and y is the recoveryvariable [34]. ǫ controls the time scale separation betweenthe membrane potential ( x ) and the recovery variables( y ) (0 < ǫ < b and g are thecontrol parameter of the neuron: depending upon themeach isolated FHN neuron can undergo several transitionsbetween equilibrium and cyclic (or oscillating) dynamics.This transitions are governed by saddle-node bifurcationor/and Andronov-Hopf bifurcation.For a clear understanding of a single neuron’s un-derlying dynamics, bifurcation diagrams with b for dif-ferent values of g are shown in Fig. 1 (for ǫ = 0 . g = 0 .
2, with variations in b the model, showsa bistable region and forms a hysteresis loop governedby two saddle-node bifurcation points (Fig. 1A). For anincrease in g , i.e., g = 0 .
6, the stable branches losetheir stability via a subcritical Andronov-Hopf bifurca-tion (Fig. 1B). With a further increase in g , there is asudden transition from a steady state to a large ampli-tude limit cycle oscillations via a subcritical Andronov-Hopf bifurcation and saddle-node bifurcation of the limitcycle (Figs. 1C and 1D) with changes in b . The systemexperience sudden transitions from an equilibrium steadystate to a non-equilibrium state (i.e., limit cycle) and viceversa. B. Networks of FHN neuron
Next, we consider a network of FHN neurons wheredynamics of the i -th neuron in a network of N neuronsare governed as follows: dx i dt = 1 ǫ ( x i − x i − y i ) + dk i N X j =1 A ij ( x j − x i ) , (2a) dy i dt = g i x i − y i + b, (2b)where i = 1 , , . . . , N denotes the node index and d isthe interaction strength between neurons. A ij is the ele-ment of the adjacency matrix ( A ) of the network, whichis equal to 1 if nodes i and j are connected and 0 other-wise. k i denotes the number of connections correspond-ing to the i -th node ( k i = P Nj =1 A ij ). g i is the param-eter through which heterogeneity among the nodes canbe introduced in the network, and b is considered as the -1 -0.5 0 0.5 1-2-1012 -1 -0.5 0 0.5 1-2-1012-1 -0.5 0 0.5 1-2-1012 -1 -0.5 0 0.5 1-2-1012 AC DB
FIG. 1. One-parameter bifurcation diagram of the FHN equa-tions (1) along changing the parameter b : for (A) g = 0 .
2, (B) g = 0 .
6, (C) g = 0 .
8, and (D) g = 1. Solid lines and dashedlines represent stable and unstable steady states, respectively.Whereas yellow and blue filled circles represent stable and un-stable limit cycles, respectively. Here, SN represents a saddle-node bifurcation point, HB represents a subcritical Andronov-Hopf bifurcation point, and SNLC represents a saddle-nodebifurcation of the limit cycle. Other parameter value: ǫ = 0 . external stress that drives the system to the transitionpoint.By simple algebraic manipulation we can rewrite (2)in terms of the Laplacian matrix ( L ) as follows: dx i dt = 1 ǫ ( x i − x i − y i ) − dk i N X j =1 L ij x j , (3a) dy i dt = g i x i − y i + b, (3b)where L ij is the element of the Laplacian matrix, suchthat L ii = P Nj =1 A ij and L ij = − A ij for i = j . Here weconsider three different network architectures; the Watts-Strogatz (WS) model [35], the Erd˝os-R´enyi (ER) model[36], and the Barab´asi-Albert (BA) model [37]. For theWS model, the interaction networks are generated usinga rewiring probability ( p ). This probability determinesthe randomness of connectivity in a network. For ex-ample, the network is regular for p = 0 with each nodecorresponds to same degree distribution, completely ran-dom when p = 1, and is small-world for 0 < p <
1. TheER and BA models also generate random networks, how-ever one follows a Poisson distribution and another scale-free distribution, respectively. For all the aforementionedinteraction networks, without any loss of generality, theaverage degree is kept at k = 1 N P Ni =1 k i = 4, with k i be-ing degree of the i -th node; also, we consider N = 1000unless mentioned otherwise. III. RESULTS
For numerical integration, we use the 4th order Runge-Kutta method with adaptive step size and initial con-ditions ( x i , y i ) are randomly selected from [ − . , . × [0 , ǫ = 0 . X = N P Ni =1 x i . When X >
X <
A. Role of network topology on determining thenature of transitions
We identify the role of network topology in responseto the external stress ( b ) by obtaining the stationary net-work state X for the WS, ER, and BA models. The sta-tionary state X is calculated by increasing the stress b from − . . δb = 0 .
01 and g i is randomly chosen from [0 . ,
1] with a uniform distribu-tion. We take the coupling strength d = 1 .
5. Figure 2 de-picts the response of the network with different rewiringprobability ( p ) of the WS model, varying from a perfectlyregular network to a completely random network via thesmall-world networks. For low stresses, each network ex-hibits an active state while moving to a resting state as b increases beyond a critical threshold. However, the na-ture of the transition from an active to a resting statedepends upon the network topology. Regular networksreveal gradual response to the stress (Fig. 2A), which is ahallmark of second-order phase transition [22, 23]. How-ever, disordered networks (WS networks with p ≥ . X >
X < b is de-creased from b = 0 . b = − . X cor-responds to increasing the stress b from − . b = 0 . b = 0 . b = 0 . b = 0 . b = 0 .
15, and b = 0 .
17. Since all thepathways of transition from one state to another in boththe scenarios are comparable, it suggests no strategy tomitigate the hysteresis in networks by controlling the ex-ternal parameter b in the WS networks (regular to disor-dered behavior). Thus, we find the transition thresholdand changing environmental stress associated with tip-ping points in disordered networks, which can lead to acatastrophic transition. However, for regular networks, -202-202-202-0.2 -0.1 0 0.1 0.2-202 -1-0.500.51 -0.2 -0.1 0 0.1 0.2-1-0.500.51 -4 -4 BACD FE
FIG. 2. Stationary state ( X ) of the FHN network with variations in b : for (A) regular ( p = 0) and small world ( p = 0 . p = 0 . .
5, WS, BA, and ER networks. (C-D) Each dashed line representssingle realisation obtained by increasing the stress b from − . b once it reaches b = 0 . b = 0 . b = 0 . b = 0 . b = 0 .
15, and b = 0 .
17, for the WS network with p = 0 and p = 0 .
5. (E-F) Stationary state distribution ofregular network ( p = 0), and small-world network ( p = 0 . b . To estimate the stationary state distribution of a network we generate 5 × stationary states, where b is randomly selected from [ − . , .
2] following uniform distribution, similarly g i is randomly chosen from uniform distribution[0 . , x i , y i ) are randomly selected from [ − . , . × [0 , ǫ = 0 .
2, and d = 1 . a transition from active to resting-state or vice-versa isnot characterized by a tipping point.To elucidate the nature of the transition thresholds,we generate stationary state distributions of (3), withdifferent network topologies. A large number of sta-tionary states are estimated to construct the station-ary state distributions. We generate 5 × station-ary states ( X ’s), where values of b are randomly selectedfrom the interval [ − . , .
2] following the uniform distri-bution, initial conditions ( x i , y i ) are randomly selectedfrom [ − . , . × [0 , g i is also randomly selectedfrom [0 . ,
1] following the uniform distribution. The sta-tionary state distribution of the regular network differsqualitatively to the transition type viz a viz disorderednetworks (see Figs. 2E–2F). A higher probability of sta-tionary states is distributed in the active (
X >
0) andresting (
X <
0) states for both the regular and the dis-ordered networks. Within the transition regions betweentwo states, stationary states are dense for regular net-works (Fig. 2E), while it becomes scattered for the dis-ordered networks (Fig. 2F). This implies that the tran-sition in regular networks is continuous. Simultaneously,the stationary state distribution depicts a sudden jumpfrom an upper to a lower state as randomness in networktopology increases. Thus, regular networks have a higher degree of resilience to environmental perturbations thanthat of disordered networks. Overall, network topologysignificantly decides the nature of a system’s transitionin virtue of external stress. The state curve depicts grad-ual response for regular networks and abrupt change fordisordered ( p ≥ . B. Effects of node heterogeneity
In this section, we study effects of different ranges ofheterogeneity g i on the network resilience for diverse net-work topologies, when exposed to external stress b . Fig-ure 3 depicts the state curves for regular networks ( p = 0)considering different ranges of g i . As the range of g i getsnarrow (i.e., weak heterogeneity), we observe the signif-icant influence of the limit cycle oscillations present inan isolated FHN neuron (see Figs. 1C–1D), during thetransition in the network. When g i ∈ [0 . ,
1] (Fig. 3A)and g i ∈ [0 . ,
1] (Fig. 3B), there exists continuous tran-sition like a transition from
X >
X <
0; however,the transition looks relatively flattered for the latter case.Further decrease in the heterogeneity range ( g i ∈ [0 . , g i ∈ [0 . , -202-202-202-0.3 -0.15 0 0.15 0.3 External stress ( b ) -202 -2-1012 × -3 -0.3 -0.15 0 0.15 0.3 External stress ( b ) -2-1012 DCB EF
XXX A X FIG. 3. Stationary states ( X ) of the regular network ( p =0) with changes in b , where g i ’s are chosen from differentintervals. g i ’s are randomly chosen, following the uniformdistribution, from the intervals: (A) g i ∈ [0 . , g i ∈ [0 . , g i ∈ [0 . , g i ∈ [0 . , g i ∈ [0 . , g i ∈ [0 . , to decide the nature of the transition. Since the dynam-ics of a few uncoupled neurons have limit cycle oscilla-tions, to identify the significance of the cyclic dynamicson the network, we calculate the distribution of the sta-tionary states when g i ∈ [0 . ,
1] (Fig. 3E) and g i ∈ [0 . , g i from a narrowinterval). This suggests that the networks with weakheterogeneity in g i are more resilient as sudden transi-tions seem impossible. Therefore, the presence of limitcycle oscillations influences the network’s resilience by in-ducing intermediate fluctuating stationary states as thenetwork transits from an active to a resting state. Fig-ure 4 depicts the proportion of nodes which are initially( d = 0) either in upper steady state, or with cyclic dy-namics, or in lower steady state with variations in thestress b and the interval in which g i lies. The intermedi-ate ranges of b , the 1000 nodes are in different states dueto the variations in g i . It is evident that when the net-work reaches its stationary state, the proportion of nodedynamics alters from its initial configuration.Next, we analyze the impact of cyclic dynamics for thedisordered small-world network ( p = 0 . g i ∈ [0 . ,
1] to g i ∈ [0 . ,
1] can sup-press the sudden transition in the network and increasethe resilience to the stress (see Figs. 5A–5B). Moreover,the cyclic dynamics’ influence becomes more evident inthe small-world network than the regular network. As b moves from a negative to a positive value, we initiallysee that the state X ( >
0) slowly approaching towards thetransition point; however, the states are largely scatteredbefore the transition to
X < g i can serve as asignificant factor in determining the type of transition inthe network. C. The role of network parameters
This section focuses on investigating the impact of net-work parameters, such as interaction strength ( d ) and av-erage degree ( k ) on the resilience of the network (3) to thestress b . We observe the influence of the network topol-ogy and the cyclic dynamics on the resilience patterns, inSubsections III Aand III B. Here, we calculate stationarystate curves for different interaction strength ( d ), each forthe regular and the disordered (small-world) networks aswell as for different choice of the heterogeneity param-eter g i (Fig. 6). In particular, to isolate the impact ofcyclic dynamics and d on the resilience of the network,we analyse the stationary state curves when g i ∈ [0 . , g i ∈ [0 . ,
1] (Figs. 6D-6F).We find that the resilience of regular and small-worldnetworks with substantial heterogeneity ( g i ∈ [0 . , d increases. For the weak interaction strength d = 0 .
5, the type of transition is robust to networktopology choice. In particular, we see continuous orsecond-order phase transitions (Fig. 6A). However, thenetwork tends toward an abrupt transition as d increases(Figs. 6A–6C). Therefore, in case of substantial hetero-geneity in g i , an increase in d can reduce adaptive re-sponses to stress and trigger sudden collapse to the rest-ing state, irrespective of the network topology. Further,regular networks experience a relatively smooth transi-tion to external stress for g i ∈ [0 . , d . Therefore,when g i ∈ [0 . ,
1] for large enough d , it is not possibleto determine the transition type of the stationary state.They are neither first-order nor second-order transition.Similarly, an increase in the network’s average degree k leads to more stronger abrupt transition from onestate to another for the system with strong heterogene-ity ( g i ∈ [0 . , p = 0 and p = 0 . b for high aver-age degrees and strong heterogeneity in the system. For g i ∈ [0 . , k for both the network topologiesleads to strong fluctuations prior to a transition to a rest- Upper steady state Limit cycle Lower steady state
71% 29% 100%100% 39% 100% 39% 61% g i ∈ [ . , ] g i ∈ [ . , ] g i ∈ [ . , ] g i ∈ [ . , ] b = − . b = 0 b = − . b = 0 . b = 0 . FIG. 4. Pie diagrams, when d = 0, representing the percentage of nodes exhibiting steady states and cyclic dynamics, fordifferent values of b and range of heterogeneity g i . We consider an ensemble average of 100 realizations for each of the piediagrams. -202-202-202-0.3 -0.15 0 0.15 0.3 External stress ( b ) -202 -2-1012-0.3 -0.15 0 0.15 0.3 External stress ( b ) -2-1012 × -3 ABCD
XXXX EF FIG. 5. Stationary states ( X ) of a disordered network( p = 0 .
5) with changes in b , where g i ’s are randomly cho-sen, following the uniform distribution, from the intervals:(A) g i ∈ [0 . , g i ∈ [0 . , g i ∈ [0 . , g i ∈ [0 . , g i ∈ [0 . , g i ∈ [0 . , ing state (Figs. 7D–7F). Moreover, the fluctuations ob-served here for p = 0 . p = 0 ( Figs. 6D–6F). Therefore, an increasein k and incorporating weak heterogeneity makes it dif- ficult to determine the transition type. Overall, strongheterogeneity in g can improve networks’ resilience with alow average degree and low coupling strength, while weakheterogeneity increases networks’ resilience with strong d and k .To check the stationary state of each node x i , i =1 , , . . . , d andexternal stress b , we calculate the percentage of nodesthose exhibit either steady state or cyclic dynamics (seeFig. 8). When the stress is close to the transition thresh-old ( b = 0), we observe that for each isolated node( d = 0), the system is dominated by steady state dy-namics (Fig. 8A). An increase in d leads to an oscillatorystate in a large number of nodes (Figs. 8B–8C), which fur-ther shows only oscillatory dynamics in all the nodes as d becomes very high (Fig. 8D). Hence, when d is low, somenodes exhibit steady state behavior while others fluctuatebetween active and resting states. The above results canlead to two exciting outcomes: First, if the steady stateof the nodes represents their presence at an upper attrac-tor, a large number of nodes exhibit functional stability.An increase in d decreases the network’s functional sta-bility as all the nodes have changed their state from anupper attractor to an oscillatory state between the twoattractors. Similarly, suppose the steady state for d = 0represents the nodes are at resting state. In that case,an increase in d will enhance the functional stability asmany nodes have shifted from a completely resting state -0.2 -0.1 0 0.1 0.2-2-1012-0.2 -0.1 0 0.1 0.2-2-1012 -0.2 -0.1 0 0.1 0.2-2-1012-0.2 -0.1 0 0.1 0.2-2-1012 -0.2 -0.1 0 0.1 0.2-2-1012-0.2 -0.1 0 0.1 0.2-2-1012 AD EB CF
FIG. 6. The effect of different interaction strengths d on network resilience for both the regular and disordered networkswith changes in b : for (A-C) g i ∈ [0 . , g i ∈ [0 . , g i becomes weak. The initial conditions and other parameters are same as in Fig. 2. -2-1012-0.3 -0.15 0 0.15 0.3-2-1012 -2-1012-0.3 -0.15 0 0.15 0.3-2-1012 -2-1012-0.3 -0.15 0 0.15 0.3-2-10 A B CE FD
FIG. 7. The effect of different average degree k on network resilience for both the regular and disordered networks with changesin b : for (A-C) g i ∈ [0 . , g i ∈ [0 . , g i becomes weak. The initial conditions and other parameters are same as in Fig. 2. to fluctuate between active and resting states.The outcome is different for a negative value of thestress ( b = − . d = 0), steady states aredominating a majority of the nodes. As d increases, thepercentage of cyclic states reduce, and finally, for large enough d dynamics in all the nodes become a steady state(Figs. 8E–8G).
59% 36%5% Upper Steady state Limit cycle Lower Steady state29% 71% 20% 79%1% 100%85% 15% 99% 1% 100% 100%
HGFBA DCE
FIG. 8. Pie diagrams representing the percentage of nodesexhibiting steady state and cyclic dynamics, for different val-ues of d in a regular network: for (A-D) b = 0, and (E-H) b = − .
15. The initial conditions ( x i , y i ) are randomly cho-sen following uniform distribution from [ − . , . × [0 ,
1] with d = 0, and then we calculate the proportion of nodes exhibit-ing steady state and oscillatory dynamics after removing thetransients. For increasing values of d , the initial conditionsused are same as for the system with no coupling ( d = 0). Weconsider ensemble average of 100 realisations for each of thepie diagram. D. Role of network topology on predictability ofstate transition
Analyses of diverse network topologies have shown dif-ferent system responses while exhibiting transition froman active to a resting state. The regular network revealsa gradual transition of the FHN neurons network froman active state, transmitting signals to different parts ofcells, to a resting state. However, this transition is some-times discontinuous, giving rise to the tipping point insmall-world interactions. Also, we observe fluctuatingdynamics in the system with weak heterogeneity beforean actual transition takes place.In many systems, a transition from a stable state to analternative state preceded by a decrease in the equilib-rium state’s resilience. Studies suggest various statisticalmethods due to critical slowing down to anticipate tran-sitions that occur in natural systems. In spatial models,critical slowing down is preceded by an increase in thestate variable’s spatial correlation and can serve as anearly warning signal. Hence, it leads to an increase influctuations in the state variable, influencing fluctuationsin the nearby state variables, leading to a more signifi-cant correlation in neighbouring components. Here, weinvestigate the impact of network architecture; in otherwords, different transition types on the early warning in-dicator of spatial systems.Consider a general form of a deterministic frameworkfollowed by stochastic perturbations [38, 39]: d z ( t ) = f ( z ( t ) , b ) dt + DdW ( t ) , (4)where deterministic part of the system is expressed as avector function f = { f , f , f , . . . , f N } and state of the system z ( t ) = { z ( t ) , z ( t ) , . . . , z N ( t ) } is controlled bythe dynamics of (4). W ( t ) is a Wiener process with D asthe noise intensity. Consider a vector of small deviationin the state u = z − z ∗ , where z ∗ is the equilibrium pointin the absence of any perturbation. Local stability of thefixed point z ∗ is calculated from the linearisation of theperturbed differential equation: d z dt = d u dt = f ( z ) = f ( z ∗ + u ) ≈ f ( z ∗ ) + J u = J u , (5)where J is the Jacobian matrix that determines the ex-pected trajectory of the perturbed state as it returns toits equilibrium state z ∗ . Consider a probability distri-bution function ρ ( u , t ) to be a solution of (4) which canbe approximated as a solution of the below linear FokkerPlanck equation [40]: ∂ρ ( u , t ) ∂t = − N X i,j =1 J ij ∂ ( ρu i ) ∂u i + 12 N X i,j =1 D ij ∂ ρ∂ u i . (6)The correlation between two state variables can be de-termined by a covariance matrix, where the diagonalelements represent fluctuations at each node and off-diagonal entries represent fluctuations between pair ofnodes. The covariance matrix C can be written as [41]: C = − J − ( D + Q ) / , (7)where Q is an antisymmetric matrix that can be esti-mated as: JQ + QJ T = JD − DJ T , (8)where the superscript T denotes the transpose of a ma-trix.Near the bifurcation/tipping threshold, the dominanteigenvalue of the Jacobian matrix J tends to zero, andthe eigenvector corresponding to the dominant eigenvaluewill become slower. The rate at which the covariancematrix’s dominant eigenvalue converges to zero is rela-tively faster than those of the other eigenvalues. Thus,the dominant eigenvalue and its ratio over the Euclideannorm of a vector consisting of all eigenvalues of the co-variance matrix serve as an early warning indicator inspatial interactions [42].We begin our investigation by considering f as the de-terministic FHN network (3) with N nodes. Each nodehas z = { x , y , x , y , . . . x N , y N } as the pair-wise poten-tial and recovery variables, respectively. Therefore, (4)can be expressed as: dz i − = dx i = ǫ ( x i − x i − y i ) − dk i N X j =1 L ij x j dt + σdW i , (9a) dz i = dy i = ( g i x i − y i + b ) dt + σdW i , (9b)where, σI = D , and σ = 0 .
01 represents the noise in-tensity of each node. We generate the stochastic time
FIG. 9. Stochastic time series of the FHN neuron model (Eq. 9) for different network topology and heterogeneity in g . (A-B)Each line represents the time series of the potential variable at each node, when the network follows strong heterogeneity in g and has perfectly regular and small-world topology ( p = 0 . p = 0, p = 0 . p = 1. (D)-(E) stochastic time series of each node when the heterogeneity become weak, that is, g i ischosen randomly from uniform distribution [0 . , Index E i gen v a l ue o f t he c o v a r i an c e m a t r i x FIG. 10. Spectrum of the covariance matrix as the systemapproaches close to the transition threshold. Each line rep-resents the eigenvalues of the covariance matrix changes fordifferent parametric values of the stress b . The Index i repre-sents the index of i th eigenvalue corresponding to the potentialvariable of the network, when each eigenvalue is arranged indescending order. series of the FHN network along changing external stress b , having both strong (Figs. 9A–9C) as well as weak het-erogeneity (Figs. 9D–9F) in g i , using (9). The stochas-tic time series along changing external stress, for eachneuron ( x i ) (Figs. 9A–9B) reveals transition from activeto resting state. Furthermore, the stationary state ( X )stochastic time series reveals that this transition is rela-tively gradual for low rewiring probability of the network(Fig. 9C). Investigating the stochastic time series for net-work with weak heterogeneity ( g i ∈ [0 . , BDCA
FIG. 11. Largest eigenvalue λ max of the covariance matrixalong with changes in the parameter b for two different net-work structure depending on rewiring probability. (A-B) g i is randomly selected from uniform distribution [0 . , g i is randomly selected from uniform distribution [0 . , hibit high amplitudes and demonstrate large fluctuationsbetween alternative stable states for the small-world tocompletely random network. However, smaller fluctua-tions are identified when the network is perfectly regular(Fig. 9F). Further, we calculate eigenvalues of the co-variance matrix (using Eq. 7 and Eq. 8) of the stochastic0system, where the Jacobian matrix J is given by: J i − , j − = ǫ (1 − x i ) − dk i L ij | z = z ∗ , if i = j, dk i , if i = j, and i, j are connected,0 , if i = j, and i, j are not connected, J i − , j = ( − , if i = j, , if i = j,J i, j − = ( g i , if i = j, , if i = j, and J i, j = ( − , if i = j, , if i = j . Let λ max = λ ≤ λ ≤ . . . ≤ λ = λ min be the eigen-values corresponding to the potential variable of the co-variance matrix C . Index is the integer varying from 1to 1000 such that Index-1 denotes the index of largesteigenvalue, Index-2 corresponds to the index of secondlargest eigenvalue, and so on. We observe that the dom-inant eigenvalue (Index-1) of C increases at a larger ratethan the other eigenvalues. Also, the difference betweenthe dominant eigenvalue and other eigenvalues increasesas the system approaches close to the transition thresh-old (Fig. 10). Thus, the dominant eigenvalue serves as anearly warning indicator for the network. To investigatethe role of different network topologies in predictabilityof the transition, we analyse the trend in the dominanteigenvalue of the covariance matrix along with the in-creasing/decreasing the stress b (Fig. 11). We observethat as the network approaches in the vicinity of transi-tion point, the dominant eigenvalue shows an increasingtrend (Figs. 11A–11B). In case of weak heterogeneity, anincrease in the largest eigenvalue is observed at relativelylow value of b (Figs. 11C–11D). This is due to fluctuationspreceded by the state of the system before the transitionhappens. Importantly, we notice that the increase is rel-atively higher for the regular network than that for thesmall world interaction. This suggests that the strengthof predictability of an upcoming transition depends uponthe network topology. IV. CONCLUSIONS AND DISCUSSION
Dynamical transitions from one stable state to anotheralternative stable state are hallmarks of loss of resiliencein many complex systems [10, 11, 43]. Studies providea myriad of mechanisms that lead to the occurrence ofsuch transitions. These mechanisms propel the transi-tion to an alternative stable state when an input con-dition by-pass a critical point [3, 44, 45]. For instance, change in population growth due to reduced populationsize, changes in interaction strength between spatiallyconnected components, etc . However, it has been chal-lenging to understand how mechanisms altering the staticdynamics of complex networks affect the occurrence ofcatastrophic and non-catastrophic transitions. Here, weaim to uncover, together with the presence of equilib-rium dynamics, how non-equilibrium/cyclic dynamics ina few nodes of a network influence its resilience. We findthat the proportion of uncoupled nodes exhibiting non-equilibrium dynamics plays a crucial role in determiningthe transition type of the network. A higher proportion ofnodes with non-equilibrium dynamics can delay the tip-ping by inducing flickering between alternative networkstates against environmental stress, irrespective of theirtopology. Moreover, the predictability of a forthcomingtransition weakens, as the network topology moves fromregular to disordered.Real-world disordered networks undergo deterioratedfunctional stability, often exhibiting discontinuous orfirst-order phase transition to an alternative state, asthe input conditions change. Such sudden transitionsin a network occur due to the inability of many nodesto cope with internal or external stresses. On the otherhand, regular networks can buffer the collapse of a net-work state in response to increased environmental stress,overall supporting the network beyond a tipping point.Notably, the network’s internal factors, such as averagedegree and coupling strength, can also regulate an abrupttransition. For weak interaction strengths, the change innetwork topology does not impact the system’s resilience,while it can delay or remove the discontinuous transitionin the system’s state. The strong coupling can decreasethe resilience of both regular as well as disordered net-works. In contrast, disordered networks are robust tochange in the network’s average degree, as they reveal asudden transition to an alternative state at low and highnetwork degrees. We observe that increasing network de-gree decreases the resilience of the regular networks, andtransitions become abrupt.Furthermore, we find that each node’s dynamics cansignificantly influence the response of a network’s statewith variations in external environmental stress. Thedominance of limit cycle oscillations in the isolated com-ponents can shape the interaction network’s response tochanging stress. When cyclic dynamics dominate, in-creased coupling strength buffers the sudden transitionby the emergence of fluctuations before transitioning toan alternate state. The propensity of fluctuations, how-ever, is observed to be more for disordered networks.Thus, cyclic dynamics significantly delay the onset ofa system’s sudden transition to a collapsed state, andthe resilience of disordered networks can be enhanced bycontrolling the coupling parameter and the network’s av-erage degree. We find that entirely random, ER, andBA networks share similar outcomes observed for disor-dered networks with small-world topology. Overall, thepresence of non-equilibrium dynamics in nodes can im-1prove network resilience by delaying sudden transitionsand flickering between alternative states in the networkpreceding a transition can also work as an early warningsignal [46].While some real networks are highly resilient to en-vironmental stress, others show vulnerable response to-wards it. The nature and scale of sudden transitionsin a network due to changing environmental stress cansignificantly owe to the diversity of network topology, in-ternal perturbations, and external heterogeneity. Hence,we analyze early warning signals [47, 48] to foresee suchvarying nature of transitions. We test the predictabilityof network-based early warning indicators for both theregular and the small-world topology [42] by investigat-ing the covariance matrix of our system in the presenceof additive Gaussian white noise. We observe that in the vicinity of the transition phase, the variance of thefluctuations increases. This increase is more substan-tial for regular than that for the small-world networks.While many real-world networks may undergo an abrupttransition from one stable state to an alternative sta-ble state, such transitions’ predictability depends uponthe network topology. Finally, understanding networkresilience and forecasting sudden transitions in the faceof rapidly evolving stress can be considered a challengingfuture problem.
ACKNOWLEDGMENTS
P.S.D. acknowledges financial support from the Science& Engineering Research Board (SERB), Govt. of India[Grant No.: CRG/2019/002402]. [1] U. Feudel, A. N. Pisarchik, and K. Showalter, “Multi-stability and tipping: From mathematics and physics toclimate and brain — minireview and preface to the focusissue,”
Chaos , vol. 28, no. 3, p. 033501, 2018.[2] M. Scheffer, S. Carpenter, J. A. Foley, C. Folke, andB. Walker, “Catastrophic shifts in ecosystems,”
Nature ,vol. 413, no. 6856, pp. 591–596, 2001.[3] M. Hirota, M. Holmgren, E. H. Van Nes, and M. Scheffer,“Global resilience of tropical forest and savanna to criti-cal transitions,”
Science , vol. 334, no. 6053, pp. 232–235,2011.[4] T. M. Lenton, H. Held, E. Kriegler, J. W. Hall, W. Lucht,S. Rahmstorf, and H. J. Schellnhuber, “Tipping ele-ments in the earth’s climate system,”
Proceedings ofthe National Academy of Sciences USA , vol. 105, no. 6,pp. 1786–1793, 2008.[5] T. M. Lenton, J. Rockstr¨om, O. Gaffney, S. Rahmstorf,K. Richardson, W. Steffen, and H. J. Schellnhuber, “Cli-mate tipping points – too risky to bet against,”
Nature ,vol. 575, pp. 592–595, 2019.[6] Y. Sharma and P. S. Dutta, “Regime shifts driven bydynamic correlations in gene expression noise,”
PhysicalReview E , vol. 96, no. 2, p. 022409, 2017.[7] S. Sarkar, S. K. Sinha, H. Levine, M. K. Jolly, and P. S.Dutta, “Anticipating critical transitions in epithelial–hybrid-mesenchymal cell-fate determination,”
Proceed-ings of the National Academy of Sciences USA , vol. 116,no. 52, pp. 26343–26352, 2019.[8] R. M. May, S. A. Levin, and G. Sugihara, “Complexsystems: Ecology for bankers,”
Nature , vol. 451, pp. 893–895, 2008.[9] M. Scheffer,
Critical transitions in nature and society ,vol. 16. Princeton University Press, 2009.[10] B. Walker, C. S. Holling, S. R. Carpenter, and A. Kinzig,“Resilience, adaptability and transformability in social–ecological systems,”
Ecology and Society , vol. 9, no. 2,2004.[11] C. Folke, S. R. Carpenter, B. Walker, M. Scheffer,T. Chapin, and J. Rockstr¨om, “Resilience thinking: in-tegrating resilience, adaptability and transformability,”
Ecology and Society , vol. 15, no. 4, 2010. [12] S. Saavedra, R. P. Rohr, V. Dakos, and J. Bascompte,“Estimating the tolerance of species to the effects ofglobal environmental change,”
Nature Communications ,vol. 4, p. 2350, 2013.[13] J. J. Lever, E. H. van Nes, M. Scheffer, and J. Bascompte,“The sudden collapse of pollinator communities,”
EcologyLetters , vol. 17, no. 3, pp. 350–359, 2014.[14] A.-L. Barab´asi et al. , Network science . Cambridge Uni-versity Press, 2016.[15] J. Gao, B. Barzel, and A.-L. Barab´asi, “Universal re-silience patterns in complex networks,”
Nature , vol. 530,no. 7590, pp. 307–312, 2016.[16] X. Liu, D. Li, M. Ma, B. K. Szymanski, H. E. Stan-ley, and J. Gao, “Network resilience,” arXiv preprintarXiv:2007.14464 , 2020.[17] Y.-H. Eom, “Resilience of networks to environmentalstress: From regular to random networks,”
Physical Re-view E , vol. 97, no. 4, p. 042313, 2018.[18] S. Saavedra, D. B. Stouffer, B. Uzzi, and J. Bascompte,“Strong contributors to network persistence are the mostvulnerable to extinction,”
Nature , vol. 478, no. 7368,pp. 233–235, 2011.[19] J. J. Lever, E. H. van Nes, M. Scheffer, and J. Bascompte,“The sudden collapse of pollinator communities,”
EcologyLetters , vol. 17, no. 3, pp. 350–359, 2014.[20] R. Albert, H. Jeong, and A.-L. Barab´asi, “Error and at-tack tolerance of complex networks,”
Nature , vol. 406,no. 6794, pp. 378–382, 2000.[21] A. V´azquez and Y. Moreno, “Resilience to damage ofgraphs with degree correlations,”
Physical Review E ,vol. 67, no. 1, p. 015101, 2003.[22] E. H. van Nes and M. Scheffer, “Implications of spatialheterogeneity for catastrophic regime shifts in ecosys-tems,”
Ecology , vol. 86, no. 7, pp. 1797–1807, 2005.[23] P. V. Mart´ın, J. A. Bonachela, S. A. Levin, and M. A.Mu˜noz, “Eluding catastrophic shifts,”
Proceedings of theNational Academy of Sciences USA , vol. 112, no. 15,pp. E1828–E1836, 2015.[24] H. Weissmann, N. M. Shnerb, and D. A. Kessler, “Simu-lation of spatial systems with demographic noise,”
Phys-ical Review E , vol. 98, no. 2, p. 022131, 2018. [25] J. M. Tylianakis and C. Coux, “Tipping points in eco-logical networks,” Trends in Plant Science , vol. 19, no. 5,pp. 281–283, 2014.[26] R. FitzHugh, “Impulses and physiological states in theo-retical models of nerve membrane,”
Biophysical Journal ,vol. 1, no. 6, p. 445, 1961.[27] J. Nagumo, S. Arimoto, and S. Yoshizawa, “An activepulse transmission line simulating nerve axon,”
Proc.IRE , vol. 50, pp. 2061–2070, 1962.[28] A. S. Pikovsky, M. G. Rosenblum, and J. Kurths,
Syn-chronization. A Universal Concept in Nonlinear Sciences .Cambridge, UK: Cambridge University Press, 2001.[29] N. Semenova, A. Zakharova, V. Anishchenko, andE. Sch¨oll, “Coherence-resonance chimeras in a networkof excitable elements,”
Physical Review Letters , vol. 117,p. 014102, 2016.[30] A. Zakharova,
Chimera Patterns in Networks: Inter-play between Dynamics, Structure, Noise, and Delay.
Springer, 2020.[31] G. B. Ermentrout and D. H. Terman,
MathematicalFoundations of Neuroscience . Springer-Verlag New York,2010.[32] L. Ramlow, J. Sawicki, A. Zakharova, J. Hlinka, J. C.Claussen, , and E. Sch¨oll, “Partial synchronization inempirical brain networks as a model for unihemisphericsleep,”
Europhysics Letters , vol. 126, pp. 2061–2070,50007.[33] M. Gerster, R. Berner, J. Sawicki, A. Zakharova,A. ˇSkoch, J. Hlinka, K. Lehnertz, and E. Sch¨oll,“Fitzhugh-nagumo oscillators on complex networksmimic epileptic-seizure-related synchronization phenom-ena.” arXiv:2007.05497 [nlin.AO].[34] I. Shepelev, D. Shamshin, G. Strelkova, and T. Vadi-vasova, “Bifurcations of spatiotemporal structures in amedium of fitzhugh–nagumo neurons with diffusive cou-pling,”
Chaos, Solitons & Fractals , vol. 104, pp. 153–160,2017.[35] D. J. Watts and S. H. Strogatz, “Collective dynamicsof ?small-world?networks,”
Nature , vol. 393, no. 6684,pp. 440–442, 1998.[36] P. Erd˝os and A. R´enyi, “On the evolution of randomgraphs,”
Publ. Math. Inst. Hung. Acad. Sci , vol. 5, no. 1,pp. 17–60, 1960.[37] A.-L. Barab´asi and R. Albert, “Emergence of scaling in random networks,”
Science , vol. 286, no. 5439, pp. 509–512, 1999.[38] D. J. Wilkinson, “Stochastic modelling for quantitativedescription of heterogeneous biological systems,”
NatureReviews Genetics , vol. 10, no. 2, pp. 122–133, 2009.[39] L. Gammaitoni, “Stochastic resonance and the ditheringeffect in threshold physical systems,”
Physical Review E ,vol. 52, no. 5, p. 4691, 1995.[40] H. Risken, “Fokker-planck equation,” in
The Fokker-Planck Equation , pp. 63–95, Springer, 1996.[41] C. Kwon, P. Ao, and D. J. Thouless, “Structure ofstochastic dynamics near fixed points,”
Proceedings of theNational Academy of Sciences U.S.A. , vol. 102, no. 37,pp. 13029–13033, 2005.[42] S. Chen, E. B. O?Dea, J. M. Drake, and B. I. Epureanu,“Eigenvalues of the covariance matrix as early warningsignals for critical transitions in ecological systems,”
Sci-entific Reports , vol. 9, no. 1, pp. 1–14, 2019.[43] C. Perrings, “Resilience in the dynamics of economy-environment systems,”
Environmental and Resource Eco-nomics , vol. 11, no. 3-4, pp. 503–520, 1998.[44] R. M. May, “Thresholds and breakpoints in ecosystemswith a multiplicity of stable states,”
Nature , vol. 269,no. 5628, pp. 471–477, 1977.[45] M. Rietkerk, S. C. Dekker, P. C. De Ruiter, andJ. van de Koppel, “Self-organized patchiness and catas-trophic shifts in ecosystems,”
Science , vol. 305, no. 5692,pp. 1926–1929, 2004.[46] R. Wang, J. A. Dearing, P. G. Langdon, E. Zhang,X. Yang, V. Dakos, and M. Scheffer, “Flickering givesearly warning signals of a critical transition to a eutrophiclake state,”
Nature , vol. 492, no. 7429, p. 419, 2012.[47] V. Dakos, S. R. Carpenter, W. A. Brock, A. M. Ellison,V. Guttal, A. R. Ives, S. K´efi, V. Livina, D. A. Seekell,E. H. van Nes, and M. Scheffer, “Methods for DetectingEarly Warnings of Critical Transitions in Time Series Il-lustrated Using Simulated Ecological Data,”
PLoS One ,vol. 7, p. e41010, 2012.[48] M. Scheffer, S. R. Carpenter, T. M. Lenton, J. Bas-compte, W. A. Brock, V. Dakos, J. van de Koppel, I. A.van de Leemput, S. A. Levin, E. H. van Nes, M. Pascual,and J. Vandermeer, “Anticipating critical transitions,”