Amplitude death and restoration in networks of oscillators with random-walk diffusion
AAmplitude death and restoration in networks of oscillators with random-walk diffusion
Pau Clusella, M. Carmen Miguel,
2, 3 and Romualdo Pastor-Satorras Departament de Física, Universitat Politècnica de Catalunya, Campus Nord B4, 08034 Barcelona, Spain Departament de Física de la Matèria Condensada, Universitat de Barcelona, Martí i Franquès 1, 08028 Barcelona, Spain Universitat de Barcelona Institute of Complex Systems (UBICS), Universitat de Barcelona, Barcelona, Spain
We study the death and restoration of collective oscillations in networks of oscillators coupled throughrandom-walk diffusion. Differently than the usual diffusion coupling used to model chemical reactions, herethe equilibria of the uncoupled unit is not a solution of the coupled ensemble. Instead, the connectivity modifiesboth, the original unstable fixed point and the stable limit-cycle, making them node-dependent. Using numericalsimulations in random networks we show that, in some cases, this diffusion induced heterogeneity stabilizes theinitially unstable fixed point via a Hopf bifurcation. Further increasing the coupling strength the oscillationscan be restored as well. Upon numerical analysis of the stability properties we conclude that this is a novel caseof amplitude death. Finally we use a heterogeneous mean-field reduction of the system in order to proof therobustness of this phenomena upon increasing the system size.
I. INTRODUCTION
Complex systems composed by many units interrelated through a heterogeneous topological pattern of interactions form awide class of natural and man made systems that can be fruitfully represented in terms of complex networks [1]. Under thisrepresentation, in which nodes stand for units and edges for pairwise unit interactions, a natural framework arises that is capableto unify the functional and structural properties of a wide variety of different systems. Particularly important at this respect arethose systems that are the substrate for some dynamic of transport process, whose properties can be strongly impacted by thetopology of interaction [2, 3].A versatile formalism to describe dynamical processes on networks is given by the theory of reaction-diffusion processes,described in terms of different kinds of particles or species that diffuse along the edges of the network and interact among theminside the nodes. Reaction-diffusion systems find a natural representation in terms of sets of nonlinear differential equations,representing in-node interactions, coupled by a diffusion term representing transport between nodes. Successful applications ofthis formalism can be found, for example, in the study of ecological species dispersal [4] or epidemic spreading [5]. In thesestudies, the network structure represents a metapopulation [6], defined as a group of populations or physical patches, joined bymigration or mobility paths.In the context of reaction-diffusion systems on networks, it is particularly noteworthy the work of Nakao and Mikhailov [7],where it was outlined the relation between the usual understanding of pattern formation in reaction-diffusion systems in lattices,established by the seminal work of Turing [8], and its counterpart in irregular topologies. In both cases, the stability of thehomogeneous equilibrium, preserved by the diffusion, can be established by means of the dispersion relation, which, in thecomplex network case, relies on the diagonalization of a Laplacian operator [7, 9]. As a result, Turing patterns are observedto arise in networked substrates, characterized by a node-dependend pattern of species density [7]. The observation of Nakaoand Mikhailov has spurred a flurry of activity in the field of reaction-diffusion processes on networks, leading to the study ofthe effects of directedness in the network edges [9], the competition in predator-prey models [10], effects on limit cycles [11],collective synchronization [12] or irregular Turing patterns [13].Most previous works, however, are based on models where the coupling follows Fick’s diffusion law, used to model chemicalreactions and heat transport. In this case, the exchange rate of the physical quantities between two nodes is proportional to theirdensity difference. A different choice for the diffusive coupling that can be considered on complex networks corresponds torandom walk diffusion [5, 14]. In this case the exchange rate between two nodes is proportional to the inverse of their degree,thus corresponding to particles diffusing by jumping between randomly chosen nearest neighbor sites. This version of diffusionis particularly relevant in the case of ecology dynamics where each node represents a metapopulation of a certain species in anecosystem, which then might randomly migrate to the surrounding environments [15, 16]. Random walks have been extensivelystudied in complex networks [17], but their application in the context of reaction-diffusion systems is rather limited [5, 14, 18].The differences between the two coupling schemes significantly affect the nature of the system. As starting point, random-walk diffusion does not accept, in general, a homogeneous equilibrium, hence the steady states are topology dependent. As directimplication of this situation, even if an equilibrium point of the network is known, its stability cannot be attained by means of adispersion-relation. Overall, the emergent dynamical patterns that might arise in reaction systems with random-walk diffusionhave been hitherto unexplored. In this paper we investigate one of such collective states which is forbidden in reaction systemswith Fick’s diffusion: the quenching of the oscillatory dynamics through the mechanism known as Amplitude Death (AD), asopposed to the Oscillation Death (OD) mechanism that has been indeed studied in systems with Fick’s diffusion [19].Oscillatory systems represent an important class for the choice of the reaction terms. Complex oscillatory systems are indeedused as a proxy to study many relevant phenomena such as chemical reactions [20, 21], cardiac cells [22, 23] neural dynam- a r X i v : . [ n li n . AO ] S e p ics [24, 25], and ecological fluctuations [26, 27]. Amplitude Death and Oscillation Death are the two main routes throughwhich the coupling among the oscillatory units leads to a state of steadiness [19]. Although the two mechanisms have beenprone to confusion, they correspond to two different dynamical phenomena, with different implications on their applicability.The emergence of OD occurs when the coupling among units induces the creation of new inhomogeneous stationary solutions.Upon modifying the coupling strength, an originally oscillatory solution, such as a limit-cycle, is destabilized and the systemfalls into the steady state created by the coupling. Nevertheless, the (unstable) oscillatory solution and the inhomogeneous fixedpoint coexist [11, 21]. On the other hand, AD occurs when different coupled oscillators pull each other out of the limit-cycleupon increasing the interaction strength. Thus, the amplitude of the oscillations diminish until it completely vanishes and theoscillators fall into the homogeneous fixed point of the system. Therefore, in this case, the oscillatory solution collapses into thesteady state and there is no coexistence [28–30].In oscillatory systems coupled with Fick’s diffusion it is possible to compute the dispersion relation associated to the homo-geneous time-varying solution, whose instability leads to Turing patterns. If these patterns are steady in time, then we are beforea case of OD [11]. Nevertheless, as we review in detail in this paper, the AD phenomena is forbidden in networks with identi-cal reaction terms and Fick’s diffusion, since the instability of the homogeneous fixed point is preserved by the coupling [28].Therefore one needs to invoke further complexity on the description of the model, such as distributed frequencies [29], delayedinteractions [30] or dynamic coupling [31].Choosing random-walk diffusion strongly modifies this scenario. Here we show that an increase of the diffusion strengthdiminishes the amplitude of the oscillations until they collapse into an inhomogeneous steady state. This phenomena differs fromOD in the sense that there is no coexistence between the oscillatory solution and the fixed point. We show that the stationarysolution corresponds to the uncoupled local equilibria of each node that has been modified anisotropically by the coupling.Therefore, the cease of the oscillations corresponds to a novel case of AD where the stabilized solution is inhomogeneous. Herewe extensively study this transition towards AD and the later restoration of the oscillations, as well as we show how suitablemodifications of the network topology lead towards the disappearance of AD. We also perform heterogeneous mean-field analysisin order to validate the generality of this phenomena for large systems. II. GRADIENT-DRIVEN DIFFUSION
General two-species reaction-diffusion processes on a network can be represented by the set of equations ˙ x i = f ( x i , y i ) + D x N (cid:88) j =1 ∆ ij x j ˙ y i = g ( x i , y i ) + D y N (cid:88) j =1 ∆ ij y j , where D x and D y are the coupling (diffusion) coefficients, f and g are non-linear reactive terms, and the matrix ∆ ij is thediscrete Laplacian operator specifying the diffusive transport of species between connected nodes. If diffusive transport is ruledby Fick’s law, that is, by the sum of fluxes of incoming species at each node, where the flux is assumed to be proportionalto the concentration [32], the discrete Laplacian can be written as ∆ ij = a ij − k i δ ij , where k i is the degree of node i , δ ij is the Kronecker symbol and a ij is the network adjacency matrix, taking value when nodes i and j are connected, and zerootherwise [1]. In this case the equations of motion ruling the dynamics of the i th node of the network read ˙ x i = f ( x i , y i ) + D x N (cid:88) j =1 a ij ( x j − x i )˙ y i = g ( x i , y i ) + D y N (cid:88) j =1 a ij ( y j − y i ) . (1)For this gradient-driven diffusion scheme, any solution of the uncoupled system ( D x = D y = 0 ) corresponds to a solution ofthe coupled system. Indeed, the diffusion term only depends on the density difference between connected nodes, so if all nodesevolve with the exact same dynamics, ( x j ( t ) , y j ( t )) = ( x ( t ) , y ( t )) for all j = 1 , . . . , N , then the diffusive coupling vanishesand the solution is preserved. Such solutions are homogeneous , meaning that the dynamic evolution of the system is the samefor all nodes independently of their topological properties. Since we consider that the dynamics of each node has 2 dimensions,only two types of attractors are possible here: (i) fixed points of the original uncoupled system, ( x (0) , y (0) ) , and (ii) periodicsolutions, also referred to as limit-cycles [33].The addition of the coupling diffusive terms can spontaneously modify the stability of these homogeneous solutions. Inprinciple, the stability of a homogeneous fixed point of system (1) boils down to study the eigenvalues λ i of a N × N -5-4-3-2-101 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 (a) -1-0.8-0.6-0.4-0.200.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 (b) R e [ λ j ] − Λ j µ j − Λ j FIG. 1. Dispersion relation for the Brusselator system coupled with Fick’s diffusion on an Erdös–Rényi network with N = 1000 and averagedegree (cid:104) k (cid:105) = 20 . (a) Real part of the eigenvalues λ j controlling the stability of the homogeneous fixed point as a function of the Laplacianeigenvalues Λ j . System parameters are a = 1 . , b = 2 . , D x = 0 . , and D y = 5 . (b) Real part of the Floquet exponents µ j controllingthe stability of the homogeneous limit-cycle solution as a function of the Laplacian eigenvalues Λ j . Continuous lines show the continuousdispersion relation, associated to a lattice at the thermodynamic limit (see Appendix A). Notice that with our definition of the Laplacian, theeigenvalues Λ j are non-positive, thus we display − Λ j in the x -axis. Other works might use a different definition. Jacobian matrix. Nevertheless it is possible to relate the eigenvalues of such high-dimensional operator to the eigenvalues of theLaplacian matrix of the system by means of a dispersion relation λ = F (Λ) which maps all the Laplacian matrix eigenvalues Λ j to the eigenvalues of the system Jacobian λ j . This corresponds to the exten-sion of discrete dispersion-relations in lattices to the complex network case [7, 34] (see Appendix A for a detailed explanation).In order to illustrate this situation, we consider as an example the Brusselator model, whose dispersion relation can be workedout analytically (see Appendix B). The resulting dispersion relation for specific system parameters is depicted in Figure 1(a),where we plot the real part of the eigenvalues controlling the stability of ( x , y (0) ) as a function of the eigenvalues of theLaplacian Λ j for the case of an homogeneous Erdös–Rényi network. As it is well known, if the network has a single connectedcomponent, the associated Laplacian always contains a single zero eigenvalue corresponding to a uniform eigenvector, and therest are all negative. In the dispersion relation, such zero eigenvalue always returns the eigenvalues of the uncoupled Jacobian(see Eq. (A3)). Thus, if the fixed point of the original system is unstable, it will remain so despite the coupling. This is the casedepicted in Figure 1(a), where the system eigenvalue corresponding to the uniform eigenvector has Re[ λ ] = 0 . . On the otherhand, if the fixed point is originally stable, the rest of the eigenvalues associated to the strictly negative Laplacian eigenmodesmight have a positive real part, thus destabilizing the homogeneous solution and generating spatio-temporal patterns. This is thewell-known Turing instability, which triggers the so-called Turing patterns [7, 8, 34].This argument holds also for limit-cycles, where a numerical (and, in some cases, analytical), dispersion relation can beobtained by means of Floquet theory [11]. Figure 1(b) depicts the relation between the Floquet exponents µ j of the limit-cyclefor the Brusselator system and the eigenvalues of the network Laplacian. The Floquet exponents corresponding to the uniformeigenvector are µ +0 = 0 and µ − (cid:39) − . , but there are many other Floquet exponents associated to different network Laplacianmodes that are positive. Thus, in this case, the originally stable limit-cycle is being destabilized through a Turing instability, soan arbitrary perturbation of the homogeneous solution will develop heterogeneous patterns. If these heterogeneous patterns arestationary, then we would be before a case of OD [11]. Regardless of whether this is the case or not, the AD phenomena willalways be forbidden here, since the instability of the original fixed point is preserved when the coupling is added, and the alsounstable oscillatory regime still exists. III. RANDOM WALK DIFFUSION
Let us move now to a different coupling, the so called random-walk diffusion [14, 18]. Here variables x j and y j might representthe population density of animals that coexist in an ecosystem divided in different discrete patches (a metapopulation [6]).Animals interact inside patches following some non-linear dynamics and can migrate between connected patches. Representingthe pattern of connections between patches by a network, and assuming that animals move at random between patches, we havethat population of node i diffuses uniformly through its k i neighbors, each receiving a flux proportional to /k i . The Laplacianmatrix can thus be written as ˜∆ = ( ˜∆ ij ) , with ˜∆ ij = a ij /k j − δ ij . The reaction-diffusion process then reads ˙ x i = f ( x i , y i ) + D x N (cid:88) j =1 a ij x j k j − x i ˙ y i = g ( x i , y i ) + D y N (cid:88) j =1 a ij y j k j − y i . (2)Here the solution of the uncoupled system ( D x = D y = 0 ) is not a solution of the coupled system unless we have a regularnetwork, i.e. all nodes have the same degree k i = k ∀ i . In this case, one can again study the stability of the system by means ofa dispersion relation and the arguments exposed above still hold. However, for generic non-regular networks there is no theorythat applies, and thus, the collective effects induced by the coupling are unknown.To study the effects of random walk diffusion, we consider the Brusselator dynamics and integrate numerically the systemdefined by Eq. (2) in an Erdös–Rényi (ER) network with average degree (cid:104) k (cid:105) = 20 . For simplicity we set D x = D y = D throughout the rest of the paper. Unequal diffusion coefficients can lead to much more complex dynamics; in fact, the generationof static Turing patterns in models with Fick’s diffusion requires D x (cid:54) = D y [34]. In Figure 2 we show the evolution of thevariable y i resulting from numerical integration for different values of the diffusion coefficient D .In Figure 2(a) we can see that a very small diffusion D = 0 . does not strongly affect the behavior of the system. Thelimit cycle defined by the variables in each node oscillates with a very similar amplitude and period, this last one close to theestimate for an isolated Brusselator, T (cid:39) π/ √ a . The different oscillators fluctuate however out-of-phase. Upon increasingthe coupling, (Figure 2(b), corresponding to D = 0 . ), the oscillators evolve in periodic phase-synchrony, but each having adifferent amplitude of the limit-cycle, amplitude that is related to the degree of the corresponding node in the network. Increasingthe coupling to D = 2 , (Figure 2(c)), the system reaches a steady state, the heterogeneous oscillations being substituted by anheterogeneous fixed point, whose value is also related to the node’s degree. Further increasing of the coupling up to D = 3 . (Figure 2(d)) restores the in-phase synchrony, with node dependent amplitude. These preliminary observations indicate theexistence of an oscillation quenching phenomenon, and a later oscillation restoration. In the oscillatory regimes for D > , thenodes appear to have the same period (except in the vicinity of the first oscillation quenching transition), period that depends on D , as a power spectrum analysis reveals (see Appendix F 3). Whether it corresponds to OD or AD is not yet clear.In order to characterize the oscillation quenching behavior of the system we consider the average of variable y over the wholenetwork, y ( t ) = 1 N N (cid:88) j =1 y j ( t ) , and measure the temporal average (cid:104) y (cid:105) and standard deviation σ = σ ( y ) , defined as (cid:104) y (cid:105) = lim T →∞ (cid:90) T y ( t (cid:48) ) dt (cid:48) , σ ( y ) = lim T →∞ (cid:90) T y ( t (cid:48) ) dt (cid:48) − (cid:104) y (cid:105) , (3)that works as a measure the average amplitude of oscillations, with σ = 0 indicating a steady state.In Figure. 3(a) we show in blue circles the value of (cid:104) y (cid:105) as a function of the diffusion coefficient D . The associated errorbars indicate the corresponding standard deviation σ ( y ) . For small values of D the network displays oscillatory dynamics,as indicated by the non-zero amplitude σ ( y ) . Upon increasing D , such oscillations diminish until they vanish completely fordiffusion larger than D (cid:39) . . The system is frozen in a heterogeneous steady state until the oscillations are restored againfor a diffusion D larger than D (cid:39) . . We are therefore in front of a case of oscillation quenching occurring at D , followedby a subsequent oscillation restoration when D crosses D . Whereas the inhomogeneity of the fixed point would support theidea that the phenomenology corresponds to OD, the two transitions between oscillatory and steady regimes strongly resemblethe AD phenomena, where the periodic solutions collapse into the fixed point, which in this case happens to be heterogeneousas opposed to the common cases of AD [19, 35]. In the folowing we study the two bifurcations through which the oscillationsvanish and restore in order to unveil the mechanism responsible for the cessation of the oscillations. y i ( t ) (a) D = 0 . y i ( t ) (b) D = 0 . y i ( t ) (c) D = 2 y i ( t ) t (d) D = 3 . FIG. 2. Time series of the variable y i corresponding to 15 randomly chosen nodes (same nodes for each panel) obtained from the numericalintegration of the Brusselator model Eq. (B1) with random walk diffusion Eq. (2) on an Erdös–Rényi networks of size N = 1000 and averagedegree (cid:104) k (cid:105) = 20 . System parameters are a = 0 . , b = 1 . . The values of diffusion coefficient are (a) D = 0 . , (b) D = 0 . ,(c) D = 2 , and(d) D = 3 . . A. Heterogeneous fixed points
In order to investigate the nature of these transitions, we start considering the underlying fixed points, corresponding to thesolutions of the nonlinear system f ( x i , y i ) + D N (cid:88) j =1 ˜∆ ij x j = 0 ,g ( x i , y i ) + D N (cid:88) j =1 ˜∆ ij y j = 0 , (4)that we solve numerically using a standard multidimensional root solver based on the Newton-Raphson method. Figure 3(b)shows the results obtained for different values of D in the phase space. For any value of D > , we obtain an inhomogeneoussolution, i.e., depending on the node i , that for small values of D is close to the fixed point of the uncoupled system ( x (0) =1 , y (0) = b/a ) (see green points). Upon increasing D , the network equilibria spreads, covering a wider area of the phase space.For large diffusion (see for instance red dots, corresponding to D = 2 ) the solution shows clusters of nodes with similar valueswhich correspond to nodes with the same degree. It is worth noticing that we find only one solution of Eq. (4) for each valueof D . Overall, it looks like the random-walk diffusion induces the hetoregeneity of the equilibria, which becomes homogeneousonly for D = 0 . (a) D D -0.100.1 0 0.5 1 1.5 2 2.5 3 3.5 4 (c) (b) -3-2-10123-5 -4 -3 -2 -1 0 (d) h y i ± σ ( y ) m a x ( R e [ λ j ] ) D y ( ) j x (0) j Dis. sol. D = 2 . D = 2 . D = 0 . D = 0 . I m [ λ j ] Re[ λ j ] D = 0 . D = 0 . D = 2 . FIG. 3. Results of simulations and numerical analysis from the Brusselator model coupled through an ER network with N = 1000 and averagedegree (cid:104) k (cid:105) = 20 . System parameters are a = 0 . and b = 1 . . (a) Blue circles indicate the time-averaged variable y , (cid:104) y (cid:105) obtained fromnumerical simulations for different values of diffusion D . Error bars indicate the standard deviation σ ( y ) . Red continuous line shows themean-field corresponding to the heterogeneous fixed point obtained solving numerically the system of equations (4). (b) Heterogeneous fixedpoint. Green plusses, blue crosses, and red points correspond to numerical solutions of system (4) for different values of D . Black squarescorrespond to the solution obtained by integrating Eq.(5) up to D = 2 .The crossing of the two black dashed lines indicate the equilibria of theuncoupled system, ( x (0) , y (0) ) = (1 , b/a ) . (c) Largest eigenvalue’s real part for different values of D . (d) Eigenvalue spectra in the complexplane for different values of D (same symbols as in (b). B. Dispersion of the fixed point
The previous numerical analysis points out that the heterogeneous equilibrium of the network is linked with the equilibriumof the uncoupled system for small diffusion values, thus indicating that the fixed point of the network corresponds to a couplinginduced modification of the original steady state. We validate this assumption by exploring the effect that small modificationsof the coupling strength D have on the fixed point. Let ( x (0) j ( D ) , y (0) j ( D )) be the solution of Eqs. (4) for the diffusion value D . Assuming that the dependence on D is smooth, we consider a small increment of the diffusion, (cid:15) > . Expanding up to firstorder terms in equation (4) one obtains (see Appendix C for a detailed derivation) the system of equations J (cid:16) x (0) i , y (0) i (cid:17) dx (0) i dDdy (0) i dD + D N (cid:88) j =1 ˜∆ ij dx (0) j dDdy (0) j dD = − N (cid:88) j =1 ˜∆ ij x (0) j y (0) j (5)for i = 1 , . . . , N , where J ( x, y ) is the × Jacobian matrix of the (uncoupled) reactive field, ( f ( x, y ) , g ( x, y )) . Such equationsform an implicit linear non-autonomous system of differential equations with the diffusion strength D as independent variable.An analytical solution of system in Eq. (5) is generally unfeasible, but a simple numerical integration of such equationusing Euler’s method, starting from the uncoupled system equilibrium, leads to the average activity reported in black crossesin Figure 3(a), whereas black squares in Figure 3(b) correspond to the fixed point obtained with this method for D = 2 . Bothresults match the solutions obtained by directly solving Eqs. (4), hence confirming that the heterogeneous equilibria of thesystem coupled through random-walk diffusion corresponds to a modification of the solution of the uncoupled system. In other − − − − D ’ . D ’ . σ | D − D c | FIG. 4. Dependence of the oscillation amplitude σ on | D − D | (red circles) and | D − D | (blue squares). Black dashed lines correspond topower laws with exponent β = 0 . . words, the steady state is not a completely new state induced by the coupling, but a smooth transformation of the originalequilibrium of the system, which turns out to be node-dependent as soon as D > . Since the steady state associated with ODcorresponds to new states created by the coupling, the observation that here the fixed point is not new is key to classify theobserved oscillation quenching mechanism as AD rather than OD. C. Stability analysis
With the numerical solution of the system obtained by directly solving Eq. (4) one can study the stability of the system bynumerically computing the eigenvalues and eigenvectors of the full system N × N Jacobian. Figure 3(c) shows the dependenceof the real part of the maximum eigenvalue as the diffusion D is tuned and Figure 3(d) shows the full spectra for three differentvalues of the diffusion. As it happens with the fixed point, for small D all the eigenvalues are located close to the eigenvaluesof the uncoupled system λ ± = 0 . ± . i . Thus, there is at least a pair of complex conjugate eigenvalues with positive realpart. As D increases, most of the eigenvalues rapidly spread and are pushed to the left of the complex plane. Nevertheless, apair of complex conjugate eigenvalues remain isolated from the rest and do not cross the imaginary axis until D (see isolatedsymbols close to the imaginary axis in figure 3(d)). Upon further increasing the diffusion, the same isolated pair of eigenvaluescrosses again the imaginary axis at D , signaling the restoration of the oscillations (see Fig.3(c)). According to these resultsit is clear that both transitions, D and D , are Hopf supercritical bifurcations. Indeed, as shown in figure 4, the amplitude ofthe oscillations vanishes as σ (cid:39) (cid:112) | D − D | , as expected for a supercritical Hopf bifurcation, and is restored then again at D with the same exponent. Also, the limit cycle centroid of each node approaches the fixed point solution as D approaches thebifurcation point (see Appendix F 3).A supercritical Hopf bifurcation induced by the coupling is one of the main routes to Amplitude Death, the other being asaddle-node bifurcation [19]. On the other hand, OD is associated to a symmetry breaking pitchfork bifurcation, which allows forthe coexistence of the oscillatory solution and the steady state [19]. Here there is no such symmetry breaking, the inhomogeneityof the fixed point being induced by the heterogeneous coupling structure, that is not arbitrary. Moreover, the limit-cycle perishesat the bifurcation points, and thus no coexistence between the two solutions is possible. This set of observations altogetherindicates that the oscillation quenching mechanism presented here falls into the category of Amplitude Death, with the noveltythat the underlying fixed point is heterogeneous due to the combination of random-walk diffusion and irregular network structure.Although all the results reported so far have been obtained using the Brusselator model, in Appendix F we show the exact samephenomenology using the Holling-Tanner predator-prey system [36]. IV. SENSITIVITY TO TOPOLOGY AND PARAMETER MODIFICATIONA. Network density
In order to clarify the impact that irregular network topologies have on the stabilization of the fixed point and the laterrestoration of the oscillations we have performed numerical simulations of the Brusselator system Eq. (B1) with random-walkdiffusion in different classes of complex networks (in the Appendix F we report results corresponding to different reactionterms).Since an heterogeneous network topology is key for the emergence of the amplitude death, we first consider ER networks withdifferent average degree (cid:104) k (cid:105) . In order to avoid fluctuations due to the independent generations of the topology, we construct net-works starting form an initial ER configuration with given size and (cid:104) k (cid:105) = 6 , and progressively add at random new connections inorder to generate denser topologies with largest average degree. Figures 5(a,b) show the amplitude of the mean-field oscillationsfor different D using ER networks with increasing average degree. In Figure 5(a) we can observe that AD takes place only insparse networks, as increasing the connectivity favors the oscillation amplitude. In Figure 5(b) the black region correspondingto AD clearly vanishes smoothly with increasing (cid:104) k (cid:105) . In fact, for the smaller values of the average degree, oscillations are notrestored even for very large values of D . As more connections are added to the topology, the stable steady state shrinks, until itcompletely vanishes for (cid:104) k (cid:105) (cid:39) . The density of the connections thus is an important factor to determine the existence of AD. B. Small-world network
Apart from the overall connectivity density, the inner topological structure of the network also plays a role in the AD andrestoration. Indeed, for regular networks (where all nodes have the same degree), no quenching can emerge, whereas irregularityseems to induce AD. To unveil this situation, we repeat the same analysis on small-world networks generated according to theWatts-Strogatz (WS) algorithm [37] (see Methods). In this case the network topology depends on a rewiring parameter p ∈ [0 , .For p = 0 the generated architecture corresponds to a regular network and the stability of the limit-cycle can be studied throughthe dispersion relation. In this case, for p = 0 , no Turing-pattern is triggered and the fully synchronized solution is the onlystable attractor. Instead, p = 1 corresponds to a random network similar to the ER model, and the behavior of the system doesnot differ much from the results reported in figure 3(a). Therefore, the behavior observed for intermediate values of p shouldreveal a topologically induced transition from synchronization to amplitude death.Figures 5(c,d) show the outcome of such simulations in WS networks for different values of p . To avoid fluctuations on thetopology generation algorithm, we construct the networks by progressively rewiring a larger fraction of random edges in thesame network. We use − p as topological order parameter, measuring the network departure from randomness , i.e., randomfor − p = 0 and regular for − p = 1 . As expected, for − p = 0 the situation is similar to the ER topology with (cid:104) k (cid:105) = 20 :there is a first bifurcation towards steadiness for small D , and a second bifurcation where oscillations are restored. For largerregularity, the two bifurcation points come closer together until they collide and vanish for p (cid:39) . . Overall the scenario isanalogous to that seen in ER networks with increasing density, but now with the disorder parameter − p playing the role ofcontrol parameter, instead of the average degree. Both scenarios are also reminiscent of the amplitude death bifurcation observedin systems of oscillators with quenched heterogeneity [29, 38]. In our case, however, the heterogeneity resides in the connectivitystructure rather than in the reactive terms, which remain homogeneous. C. System parameter
Finally, the choice of the dynamical parameters also plays a role on whether amplitude death arises or not. Figure 5(e,f) showsthe outcome of numerical simulations for the same ER network with (cid:104) k (cid:105) = 20 , keeping fixed the parameter a = 0 . and varying b . For b < a + 1 = 1 . the reactive part does not display oscillations and the fixed point is stable. For b > a + 1 = 1 . theuncoupled system undergoes a Hopf bifurcation and the stable limit-cycle emerges. Here we investigate the dynamics of thecoupled system for b varying from . to . . Not surprisingly, for values of b close to the bifurcation of the uncoupled systemthe amplitude death regime is quickly attained at small values of D , and the restoration of the oscillations is not triggered evenfor large coupling, i.e., the fixed point stabilizes quickly with D . As b increases, the limit-cycle has a larger amplitude and thediffusion-induced stabilization of the fixed point requires larger diffusion values. Indeed, the region of amplitude death becomessmaller as they move away from the single-node bifurcation point, until it completely vanishes for b (cid:38) . in the same manneras it does when the modifications are done in the topology of the network. V. HETEROGENEOUS MEAN-FIELD ANALYSIS
In order to gain some insight on the behavior of amplitude death an restoration induced by random walk diffusion, we applythe standard tool of heterogeneous mean-field (HMF) theory [3, 39], specialized to reaction-diffusion processes [14]. The basisof HMF consists in the annealed network approximation [3, 40], that replaces the static adjacency matrix of a real network by (a) (b) (c) (d) (e) (f) σ h k ih k i D σ σ − p − p D σ σ bb D σ FIG. 5. Amplitude death and restoration dependence on different topological and system parameters. Top panels show the dependence of σ on D where each line corresponds to a specific topological or system parameter value. Bottom panels show heatmaps of the amplitude upontuning the diffusion and the corresponding parameter. Results obtained with simulations of the Brusselator model in networks with N = 1000 and system parameters a = 0 . and b = 1 . unless otherwise stated. (a,b) Dependence of the oscillation amplitude σ on the diffusion D ina range of ER networks with different average degree k . (c,d) Dependence of the oscillation amplitude σ on the diffusion D for a range ofWatts-Strogatz networks with average degree (cid:104) k (cid:105) = 20 and rewiring probability p ∈ [0 , . (d). (e,f) Dependence of the oscillation amplitude σ on the diffusion D for different values of the system parameter b with fixed . in a single ER networks with (cid:104) k (cid:105) = 20 . an average over degree classes ¯ a ij that, in the case of uncorrelated networks, takes the form ¯ a ij = k i k j (cid:104) k (cid:105) N .
Introducing this expression into Eq. (2), we obtain the HMF dynamical equations ˙ x i = f ( x i , y i ) + D (cid:16) ˜ k i x − x i (cid:17) ˙ y i = g ( x i , y i ) + D (cid:16) ˜ k i y − y i (cid:17) (6)where ˜ k i = k i / (cid:104) k (cid:105) , x = N (cid:80) Nj =1 x j and y = N (cid:80) Nj =1 y j are the average (mean-field) activities of x and y variables. Withinthis framework, the dynamics of each node depends only on its own degree, the mean-field values x and y , and the averageconnectivity of the network. In fact, one can assume then that all nodes with the same degree behave identically, thus formallyreducing the N -dimensional system (6) to a n -dimensional system where n is the number of different degrees in the net-work [3]. An additional interest of this approach lies in the fact that it allows to consider network sizes much larger than thosepermitted in a direct numerical integration, since in general n (cid:28) N .System (6) can be solved semi-analytically assuming that x and y are parameters that are to be determined self-consistently aposteriori (see Appendix E for details). Figure 6 shows the results of the heterogeneous mean-field approximation. In Figure 6(a)we show the fixed point obtained for the reduced system (black squares) compared to the actual fixed point obtained solving Eq.(4) (red circles) for D = 2 . The solution of the reduced system matches the overall spreading of the solution, with the blacksquares fitting nicely the center of the bands displayed by the actual fixed point, corresponding to nodes with the same degree.In Figure 6(b) the mean-field activity of y determined by the HFM equations (dashed black curve) matches with good accuracythat of the actual system (red continuous curve). Moreover, the stability analysis of the fixed point also coincides with that of theoriginal system, with a cloud of eigenvalues lying far on the stable part, and a pair of complex conjugate eigenvalues crossingback and forth the imaginary axis upon modifying D . Thus, the mean-field reduction does not only provide a good proxy tostudy the fixed point, but also allows to study the amplitude death phenomena in terms of stability.0 (a) (c) -3-2-10123-5 -4 -3 -2 -1 0 (b) -0.100.1 0 1 2 3 4 h k i = 20 h k i = 30 h k i = 40 h k i = 50 (d) y ( ) j x (0) j h y i D I m [ λ j ] Re[ λ j ] m a x ( R e [ λ j ] ) D FIG. 6. Mean-field analysis of the Brusselator model. (a) Red circles correspond to the fixed point for D = 2 as obtained from solvingsystem (4). Black squares correspond to the mean-field result for the same D . (b) Red circles correspond to the eigenvalues of the originalfixed point for D = 2 , whereas open black square indicate the spectra resulting from the mean-field reduction. (c) Average activity of the fixedpoint obtained from directly solving system (4) (red continuous curve), and from the mean-field reduction (black dashed curve). (d) Largesteigenvalue’s real part for the mean-field solution with ER networks with average degree (cid:104) k (cid:105) = 30 (red), 30 (blue), 40 (green), and 50 (purple).For each set of networks each line denotes a different network size. From top to bottom, N = 10 , , , and . A. Large network size limit
The good agreement of the mean-field theory with the numerical results of the ER networks allows to extend our analysis tolarger systems. Figure 6(d) shows the largest eigenvalue real part for different values of D as obtained from the heterogeneousmean-field analysis. We analyzed 4 sets of ER networks of size N = 10 , , , and , with different average degrees (cid:104) k (cid:105) = 20 , , , and . For a specific value of (cid:104) k (cid:105) , we observe that networks with different sizes produce qualitativelysimilar results, converging to a well defined limit for sufficiently large N . For instance, for (cid:104) k (cid:105) = 20 the upper red curve,corresponding to N = 10 , shows a smaller region of amplitude death, whereas the red curves below, corresponding to largersystems, converge nicely upon increasing N . The same situation is repeated for the other values of (cid:104) k (cid:105) , thus confirming that theeffects of the average degree on the amplitude death are not finite-size but rather robust. Indeed, the overall conclusion of thisanalysis is that the network behavior depends strongly on the degree of each node, but is only mildly dependent on system size.In the Appendix F 3 we extend this analysis to the other topological parameters displayed in Figure 5, and in Appendix F to theHolling-Tanner system. VI. DISCUSSION
Reaction-diffusion processes are a powerful formalism to represent general dynamical processes on networks, in which par-ticles or species interact inside nodes while moving diffusively between pairs of connected nodes. By analogy with chemicalreactions, a gradient-driven diffusion term, given by Fick’s law, is usually assumed. However, in certain circumstances, a ran-dom walk diffusion term, in which particles jump at random along edges, might be more realistic. Here we have shown howthe nature of the diffusion term can alter the behavior of the limit cycles in oscillatory reaction-diffusion processes when drivenby the diffusion term. Thus, while gradient-diffusion can quench the oscillations by means of an oscillation death mechanism,in which the stability of the original fixed point is preserved by the diffusion operator, random walk diffusion generates an am-plitude death quenching characterized by a inhomogeneous steady state induced by the structure of the diffusion operator. In1this case, a further increase in diffusion restores the oscillations, in the form of a set of limit cycles with different amplitude foreach node. The transitions to amplitude quenching and restoration are observed to correspond to Hopf supercritical bifurcations,in agreement with an amplitude death mechanism at work. Our observations, backed up by direct numerical integration of theBrusselator systems in random networks, are confirmed by a linear stability analysis of the eigenvalues of the diffusion operatorand by a seminumerical heterogeneous mean-field approximation, which allows us to check the phenomenology in very largenetworks. The role of random walk diffusion in the amplitude quenching of oscillations, which is observed for different classesof homogeneous networks, different sets of parameters and even different reaction terms, is thus a robust phenomenology, whichcould be relevant in processes such as epidemic spreading or ecological dispersion. Future promising work along these lineswould include studying the effects of an heterogeneous, scale-free topology [41], as observed in many natural networks.
Appendix A: Dispersion relation in reaction-diffusion processes
For a detailed introduction on Turing patterns and how to compute dispersion-relations on regular lattices the reader can checkMurray’s book [34]. Here we summarize the results of Nakao et al. [7] on reaction-diffusion systems in complex networks. Letus consider a reaction-diffusion system in a complex network with a gradient-driven diffusion term ˙ x i = f ( x i , y i ) + D x N (cid:88) j =1 ∆ ij x j , ˙ y i = g ( x i , y i ) + D y N (cid:88) j =1 ∆ ij y j , with a Laplacian matrix ∆ ij = a ij − k i δ ij . This Laplacian matrix is semi-definite negative, and its eigenvalues are all real andnon-positive. Let ( x (0) , y (0) ) be a fixed point of the uncoupled system ( D x = D y = 0 ), and hence the homogeneous solution ofthe coupled system. Inserting an arbitrary perturbation { ( δx i ( t ) , δy i ( t )) } Ni =1 of the homogeneous solution to the equations andretaining up to first order terms one gets that the evolution of the perturbation behaves as ˙ δx i ˙ δy i = J (cid:16) x (0) , y (0) (cid:17) (cid:32) δx i δy i (cid:33) + N (cid:88) j =1 ∆ ij (cid:32) D x δx j D y δy j (cid:33) , (A1)where J is the Jacobian of the homogeneous system, evaluated at ( x (0) , y (0) ) . Let Φ ( α ) = ( φ ( α )1 , . . . , φ ( α ) N ) T be the Laplaciannormalized eigenvector associated the eigenvalue Λ α for α = 1 , . . . , N , such that (cid:88) m ∆ jm φ ( α ) m = Λ α φ ( α ) j . (A2)We can express the perturbation of the homogeneous solution in terms of the basis of the eigenvectors as (cid:32) δx j δy j (cid:33) = N (cid:88) α =1 u ( α ) v ( α ) φ ( α ) j . Applying this change of coordinates on Eq. (A1) and making use of the relation (A2) and of the linear independence of theeigenvectors { Φ ( α ) } Nα =1 one obtains that the evolution of ( ˙ u ( α ) , ˙ v ( α ) ) T becomes independent for each α = 1 , . . . , N throughthe relation ˙ u ( α ) ˙ v ( α ) = J (cid:16) x (0) , y (0) (cid:17) u ( α ) v ( α ) + Λ α D x u ( α ) D y v ( α ) = ∂f∂x (cid:16) x (0) , y (0) (cid:17) + D x Λ α ∂g∂x (cid:16) x (0) , y (0) (cid:17) ∂f∂y (cid:16) x (0) , y (0) (cid:17) ∂g∂y (cid:16) x (0) , y (0) (cid:17) + D y Λ α u ( α ) v ( α ) . (A3)Thus the stability of the homogeneous solution simplifies to the study of the eigenvalues (and eigenvectors) of the previous × matrix, which are obtained as functions of the Laplacian eigenvalues Λ α . See below for an explicit application to the Brusselatorsystem. The extension of this analysis to limit-cycles makes use of Floquet theory, see [11].2 Appendix B: The Brusselator
The Brusselator [21] is a prototypical model of autocatalytic chemical reaction with two species showing oscillatory behaviorin the form of a limit cycle. In its dimensionless form, in the absence of diffusion, it is defined by the reaction terms f ( x, y ) = 1 − x ( b + 1) + ax y, (B1) g ( x, y ) = bx − ax y, where a, b > are model parameters. The system has a single fixed point at ( x (0) , y (0) ) = (1 , b/a ) , whose stability is ruled bythe Jacobian J (cid:16) x (0) , y (0) (cid:17) = (cid:18) b − a − b − a (cid:19) . The fixed point is therefore stable for b < a + 1 . Otherwise, the system exhibits periodic behavior, with a period, for b close to a , approximately equal to T = 2 π/ √ a . The transition from the fixed point to the oscillatory regime is through a supercriticalHopf bifurcation. Unless otherwise stated, the system parameter values used in the paper are a = 0 . and b = 1 . .When considering the Brusselator reaction terms (B1) with Fick’s diffusion (1) one can obtain the dispersion relation of thehomogeneous fixed point. Following Eq. (A3), one only needs to to compute the eigenvalues of (cid:32) b − D x Λ α a − b − a + D y Λ α (cid:33) , which yields λ ± (Λ α ) = 12 (cid:18) b − a − − Λ α ( D x + D y ) ± (cid:113) ( D y − D x )Λ α + b − a − ab (cid:19) . This is the continuous curve displayed in Figure 1(a).
Appendix C: Dispersion of the fixed point
Let (cid:16) x (0) i ( D ) , y (0) i ( D ) (cid:17) be the fixed point solution of the reaction-diffusion system (2) with diffusion value D , so that f (cid:16) x (0) i ( D ) , y (0) i ( D ) (cid:17) + D N (cid:88) j =1 ˜∆ ij x (0) j ( D ) = 0 g (cid:16) x (0) i ( D ) , y (0) i ( D ) (cid:17) + D N (cid:88) j =1 ˜∆ ij y (0) j ( D ) = 0 . (C1)for i = 1 , . . . , N . Assuming a smooth dependence of ( x (0) j ( D ) , y (0) j ( D )) on D for j = 1 , . . . , N , we consider a small incrementof the diffusion (cid:15) > . The aim is to solve then f (cid:16) x (0) i ( D + (cid:15) ) , y (0) i ( D + (cid:15) ) (cid:17) + ( D + (cid:15) ) N (cid:88) j =1 ˜∆ ij x (0) j ( D + (cid:15) ) = 0 and g (cid:16) x (0) i ( D + (cid:15) ) , y (0) i ( D + (cid:15) ) (cid:17) + ( D + (cid:15) ) N (cid:88) j =1 ˜∆ ij y (0) j ( D + (cid:15) ) = 0 . Expanding such equations around ( x (0) j ( D ) , y (0) j ( D )) and retaining the first-order terms one obtains, (cid:15) f x (cid:16) x (0) i ( D ) , y (0) i ( D ) (cid:17) dx (0) i ( D ) dD + (cid:15) f y (cid:16) x (0) i ( D ) , y (0) i ( D ) (cid:17) dy (0) i ( D ) dD + (cid:15) N (cid:88) j =1 ˜∆ ij (cid:32) D dx (0) j ( D ) dD + x (0) j ( D ) (cid:33) = 0 and (cid:15) g x (cid:16) x (0) i ( D ) , y (0) i ( D ) (cid:17) dx (0) i ( D ) dD + (cid:15) g y (cid:16) x (0) i ( D ) , y (0) i ( D ) (cid:17) dy (0) i ( D ) dD + (cid:15) N (cid:88) j =1 ˜∆ ij (cid:32) D dy (0) j ( D ) dD + y (0) j ( D ) (cid:33) . (cid:15) one finally can write the differential system that rules the dependence of ( x (0) i , y (0) i ) on D , J (cid:16) x (0) i , y (0) i (cid:17) dx (0) i dDdy (0) i dD + D N (cid:88) j =1 ˜∆ ij dx (0) j dDdy (0) j dD = − N (cid:88) j =1 ˜∆ ij x (0) j y (0) j for i = 1 , . . . , N where J ( x, y ) is the × Jacobian matrix of the (uncoupled) reactive field, ( f ( x, y ) , g ( x, y )) . The previous equation is thusan implicit linear non-autonomous system of differential equations that can be solved by means of usual numerical integratorsusing as initial condition for D = 0 the solution of the uncoupled oscillators. Appendix D: Network models1. Erdös–Rényi networks
The Erdös–Rényi (ER) model for network generation provides random network topologies [1]. In particular we use the G ( N, p ) model, where a each pair of nodes is connected with probability p ∈ [0 , . The average degree of the network is then pN and the degree distribution follows a binomial distribution with parameters N − and p .
2. Watts-Strogatz model
The Watts-Strogatz (WS) model generates networks with small-world connectivity, i.e., with small density and diameter,while still showing a large degree of transitivity or clustering [37]. We use this model as way to generate networks halfwaybetween random and regular topologies. The model starts with N nodes distributed on a ring, each node being connected to its (cid:104) k (cid:105) nearest neighbors. Then, every edge of the network is rewired with probability p ∈ [0 , . The small-world network classis represented for small values of the rewiring probability, whereas for p close to 1 the topology becomes closer to that of a ERnetwork. Appendix E: Heterogeneous mean-field analysis
We consider the HMF equations Eq. (6) specialized for the Brusselator reaction terms. Solving for the fixed point one obtainsthat (dropping the subindices for simplicity), we obtain y = 1 − x ( b + 1) + D (˜ kx − x ) − ax = bx + Dyax + D , (E1)from where one can obtain a cubic equation for the value of x : − a ( D + 1) x + a (1 + D ( x + y )) x − D ( b + 1 + D ) x + D (1 + Dx ) = 0 . Thus x = − a (cid:18) b + C + ∆ C (cid:19) (E2)where C = (cid:32) ∆ ± (cid:112) ∆ − (cid:33) / ∆ = a (1 + D ( x + y )) − aD ( D + 1)( b + 1 + D )∆ = 2 a (1 + D ( x + y )) − a D ( D + 1)( b + 1 + D )(1 + D ( x + y ))+ 27 a ( D + 1) D (1 + Dx ) . x = 1 .However, in order to obtain y self-consistently one needs to rely on numerics. In practice, one starts with an educated guess forthe mean-field y and determines then all x j and y j using equations (E2) and (E1). Using a bisection method one can reduce theerror between the new estimated values x and y and the initial ones to the desired accuracy. Appendix F: The Holling-Tanner predator-prey system
Here we revisit the results exposed on the main manuscript using a different system for the single-node dynamics. We use theHolling-Tanner predator-prey system [36]. In this model x ≥ and y ≥ respectively represent the population fraction of aprey and predator species in an ecosystem. On its dimensionless form, the corresponding reactive terms read f ( x, y ) = x (1 − x ) − axyx + dg ( x, y ) = yb (cid:16) − yx (cid:17) , where a, b , and d are system parameters [42].The system has a fixed point at x (0) = 12 (cid:16) (1 − a − d ) + (cid:112) (1 − a − d ) + 4 d (cid:17) , y (0) = x (0) and the corresponding Jacobian is J (cid:16) x (0) , y (0) (cid:17) = x (0) (cid:18) ax (0) ( x (0) + d ) − (cid:19) − ax (0) x (0) + db − b . The resulting eigenvalues have positive real part if and only if x (0) (cid:18) ax (0) ( x (0) + d ) − (cid:19) < b and au (0) + d − ax (0) ( x (0) + d ) > . Fixing a = 1 and d = 0 . , the fixed point becomes unstable through a supercritical Hopf bifurcation at b c (cid:39) . ... , so thatfor b < b c the single node dynamics exhibit periodic oscillations. Unless otherwise stated, we consider here b = 0 . .
1. Amplitude death and restoration
Now we consider a network of Holling-Tanner oscillators coupled with random-walk diffusion. Let A = ( a ij ) be the adja-cency matrix of the (symmetric and fully connected) network, and let k j be the degree of node j . The random-walk Laplacianmatrix is then defined as ˜∆ := ( ˜∆ ij ) where ˜∆ ij = a ij /k j − δ ij . With this notation, the dynamical rule governing the evolutionof the i th node of the system reads ˙ x i = f ( x i , y i ) + D N (cid:88) j =1 ˜∆ ij x j ˙ y i = g ( x i , y i ) + D N (cid:88) j =1 ˜∆ ij y j (F1)where D is the diffusion or coupling strength.Figure 7(a) shows numerical results obtained using this system, analogous to Figure 2 of the main manuscript. Overall, thescenario is akin to the one exposed with the Brusselator model. The error bars of the blue points in Figure 7(a) indicate theamplitude of the mean-field oscillations, which decrease until become steady at D (cid:39) . . Further increasing the coupling,the oscillations restore around D (cid:39) . . Again, solving for the fixed point of the system with a Newton-Raphson methodone obtains the steady state solution towards which the oscillations die (see red curve in Figure 7(a)). As in the Brusselatormodel, single snapshots of the fixed point already show that the heterogeneity of the steady state increases with the diffusion5 (a) D D -0.006-0.004-0.00200.0020.0040.006 0 0.5 1 1.5 2 (c) h y i ± σ ( y ) m a x ( R e [ λ j ] ) D y ( ) j x (0) j Disp. D = 0 . D = 2 . D = 0 . D = 0 . (b) I m [ λ j ] Re[ λ j ] (d) FIG. 7. Results of simulations and numerical analysis from the Holling-Tanner model coupled through an Erd˝os–Rényi network with N = 1000 and average degree (cid:104) k (cid:105) = 20 . System parameters are a = 1 , b = 0 . , and d = 0 . . (a) Blue circles indicate the time-averagedmean-field activity of variable y , (cid:104) y (cid:105) obtained from numerical simulations for different values of diffusion D . Error bars indicate the standarddeviation σ ( y ) . Red continuous line shows the mean-field corresponding to the heterogeneous fixed point obtained solving numerically thecorresponding system of nonlinear equations. Black crosses indicate the same mean-field, but as obtained from the dispersion equation. (b)Heterogeneous fixed point as obtained numerically solving the system of nonlinear equations for different values of D (see legend). Blacksquares mark the result obtained by numerical integration of the dispersion equation. (c) Largest eigenvalue’s real part for different values of D . (d) Eigenvalue spectra in the complex plane for different values of D (symbols as in (b)). (see green crosses, blue pluses and red circles in Figure 7(b)). The analysis of the dispersion of the fixed point, explainedin detail in the paper, unveils again that the steady state corresponds to a diffusion-induced modification of the single-nodesolution, ( x (0) , y (0) ) (cid:39) (0 . , . . Indeed, the black crosses in Figure 7(a) and black squares in 7(b) indicate the goodness ofthe dispersion equation to obtain the underlying fixed point. Finally, stability analysis reveals again a Hopf bifurcation towardsthe amplitude death and its posterior restoration, in a similar fashion as for the Brusselator system (see Figure 7(c) and 7(d)).Overall, the observed phenomenology proves ubiquitous despite the strong differences on the formulation of the reactive termsfor both systems.
2. Topology and parameter modification
Figure 8 shows the same instances reported in the paper where topology and parameter modification vary the scenario depictedin the previous section. In particular, an increase of the density in the network (see Figures 8(a,b)) reduces the region ofsteadiness, until both critical points collapse. The same occurs for the WS networks upon increasing the regularity of thenetworks (see Figures 8(c,d)). Finally, the same results are obtained by keeping the topology fixed and tuning the systemparameter b . Closer to the single-node bifurcation point b c , the coupling rapidly tends to kill the periodic behavior, whereasincreasing b finally vanishes the AD phenomena (see Figures 8(e,f)).
3. Heterogeneous mean-field approach
The heterogeneous mean-field approach also provides a good agreement with the original simulations, although here theself-consistency of the problem needs further numerical effort to be solved, since the nonlinearities in the formulation of thesystem limit the analytical results.6 (a) (b) (c) (d) (e) (f) σ h k ih k i D σσ D − p − p D σσ b c − bb c − b D σ FIG. 8. Amplitude death and restoration dependence on different topological and system parameters. Top panels show the dependence of σ on D where each line corresponds to a specific parameter value. Bottom panels show heatmaps of the amplitude upon tuning the diffusionand the parameter. System parameters are a = 1 , d = 0 . , and b = 0 . , except in panels (e,f) where b is modified as indicated. (a,b)Dependence of the oscillation amplitude σ on the diffusion D in a range of ER networks with different average degree (cid:104) k (cid:105) . (c,d) Dependenceof the oscillation amplitude σ on the diffusion D for a range of Watts-Strogatz networkss with rewiring probability p ∈ [0 , . (e,f) Dependenceof the oscillation amplitude σ on the diffusion D for different values of the system parameter b in a single ER network with (cid:104) k (cid:105) = 20 . Let us recover system (1) but we now substitute the elements of the adjacency matrix a ij by the annealed network approxi-mation, a ij = k i k j (cid:104) k (cid:105) N .
After applying this transformation in the coupling terms of Eq. (1) and doing some straightforward calculations one obtains themean-field system ˙ x i = f ( x i , y i ) + D (cid:16) ˜ k i x − x i (cid:17) ˙ y i = g ( x i , y i ) + D (cid:16) ˜ k i y − y i (cid:17) (F2)where ˜ k i := k i / (cid:80) Nj =1 k j , and x = N (cid:80) Nj =1 x j and y = N (cid:80) Nj =1 y j are the mean-field activities of x and y variables. Here,in order to solve for the fixed point, ˙ x i = 0 and ˙ y i = 0 one needs to find the roots of a 5th order polynomial, thus, unlike theBrusselator model, the expressions of x i and y i cannot be explicitly obtained. Therefore, given a value of x and y we obtain thefixed points (cid:16) x (0) i , y (0) i ) (cid:17) for i = 1 , . . . , N numerically. Proceeding this way, one can approximate the values of x and y up tothe desired accuracy.Figure 9 shows the results corresponding to this analysis. Again, the heterogeneous mean-field approach provides a goodrepresentation of the underlying scenario with ER. Therefore, we implement this dimensionality reduction approach to study theamplitude death and restoration in large networks. Again, as it happens in the Brusselator system, increasing the network sizedoes not effect the observed phenomenology. Therefore, we finally repeat the analysis depicted in Figure 8 but using networkswith N = 10 nodes analyzed through HMF, using the maximum of the largest eigenvalue of the system, ρ := max( Re [ λ j ]) as a proxy to identify whether there oscillations or not. Results in Figure 10. Appendix G: Supplementary material for the Brusselator system
Figure 11(a) shows the Power Spectral Density corresponding to time series depicted in Figure 2 of the main paper fordiffusion values D = 0 . , D = 0 . , and D = 3 . . For D = 0 . the spectra corresponding to different nodes are irregular, but7 (a) (c) -0.6-0.4-0.200.20.40.6-0.8 -0.6 -0.4 -0.2 0 (b) -0.01-0.00500.005 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 (d) y ( ) j x (0) j h y i D I m [ λ j ] Re[ λ j ] m a x ( R e [ λ j ] ) D FIG. 9. Mean-field analysis of the Holling-Tanner model with a = 1 , b = 0 . , and d = 0 . . (a) Red circles correspond to the fixed point for D = 2 as obtained from the numerical solution for the fixed point on the original system. Black squares correspond to the mean-field resultfor the same value of D . (b) Red circles correspond to the eigenvalues of the fixed point for D = 2 , whereas open black square indicate thespectra resulting from the mean-field reduction. (c) Average activity of the fixed point obtained from directly solving the fixed point system ofequations. (red continuous curve), and from the mean-field reduction (black dashed curve). (d) Largest eigenvalue’s real part for the mean-fieldsolution with ER networks with average degree (cid:104) k (cid:105) = 30 (red), 30 (blue), 40 (green), and 50 (purple). For each set of networks each linedenotes a different network size. From top to bottom, N = 10 , , , and . the peak is only slightly shifted around a similar value, thus the global dynamics corresponds to irregular fluctuations evolvingat similar frequencies. For D = 0 . and D = 3 . the spectra shows a sharp peak, indicating that the frequency entrainmentis total. In Figure 11(b) we show how the centroid of each node’s limit-cycle approaches the (unstable) fixed point y (0) j as D increases up to D . The centroid are computed as the time-averaged value over a long enough time series, (cid:104) y j (cid:105) . Figure 11(c)shows the decrease of the amplitude of each node’s oscillations until it vanishes at D . Finally, figure 12 displays the results ofthe HMF analysis for networks with N = 10 analogous to figure 5 of the main paper. [1] Mark Newman, Networks: An Introduction (Oxford University Press, Inc., New York, NY, USA, 2010).[2] Alain Barrat, Marc Barthelemy, and Alessandro Vespignani,
Dynamical processes on complex networks (Cambridge university press,2008).[3] S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, “Critical phenomena in complex networks,” Rev. Mod. Phys. , 1275–1335(2008).[4] Shigefumi Hata, Hiroya Nakao, and Alexander S. Mikhailov, “Dispersal-induced destabilization of metapopulations and oscillatoryturing patterns in ecological networks,” Scientific Reports , 3585 (2014).[5] V. Colizza, R. Pastor-Satorras, and A. Vespignani, “Reaction-diffusion processes and metapopulation models in heterogeneous networks,”Nature Physics , 276–282 (2007).[6] I.A. Hanski and O.E. Gaggiotti, Ecology, Genetics and Evolution of Metapopulations (Elsevier Science, Princeton, 2004).[7] Hiroya Nakao and Alexander S. Mikhailov, “Turing patterns in network-organized activator-inhibitor systems,” Nature Physics , 544–550 (2010).[8] Alan Mathison Turing, “The chemical basis of morphogenesis,” Philosophical Transactions of the Royal Society of London B: BiologicalSciences , 37–72 (1952).[9] Malbor Asllani, Joseph D. Challenger, Francesco Saverio Pavone, Leonardo Sacconi, and Duccio Fanelli, “The theory of pattern forma-tion on directed networks,” Nature Communications , 4517 (2014). -0.015-0.01-0.00500.0050.01 0 0.5 1 1.5 2 (a) (b) -0.01-0.00500.0050.01 0 0.5 1 1.5 2 (c) (d) -0.02-0.0100.010.020.03 0 0.5 1 1.5 2 (e) (f) ρ h k ih k i D ρρ − p − p D ρρ b c − bb c − b D ρ FIG. 10. Results of the HMF analysis for the Holling-Tanner model with a = 1 , b = 0 . , and d = 0 . . Top panels show the dependence of ρ on D where each line corresponds to a specific parameter value. Bottom panels show heatmaps of the amplitude upon tuning the diffusion andthe parameter. System parameters are a = 1 , d = 0 . , and b = 0 . , except in panels (e,f) where b is modified as indicated. (a,b) Dependenceof ρ on the diffusion D in a range of ER networks with different average degree k . (c,d) Dependence of the oscillation amplitude ρ on thediffusion D for a range of Watts-Strogatz networkss with rewiring probability p ∈ [0 , . (e,f) Dependence of the oscillation amplitude ρ onthe diffusion D for different values of the system parameter b in a single ER network with (cid:104) k (cid:105) = 20 . − − − (a) (c)(b) P S D Frequency σ ( y i ) D | h y i i − y ( ) i | FIG. 11. (a) Power Spectral Density of the time series corresponding to figure 2(a) of the main paper. Red curves correspond to D = 0 . ,blue corresponds to D = 0 . , and purple to . . Results corresponding to 15 different nodes. For each value of the diffusion the PSD of the15 nodes depicted in figure 2(a) of the paper are displayed. (b) Distance between the focus of each node limit-cycle, (cid:104) y j (cid:105) and the underlyingfixed point y (0) j for different values of D . Only results of 100 nodes shown. (c) Amplitude of each node evolution σ ( y j ) for different valuesof D . Only results of 100 shown.[10] Lucas D. Fernandes and M. A. M. de Aguiar, “Turing patterns and apparent competition in predator-prey food webs on networks,”Physical Review E , 056203 (2012), arXiv:1207.3424.[11] Joseph D. Challenger, Raffaella Burioni, and Duccio Fanelli, “Turing-like instabilities from a limit cycle,” Phys. Rev. E , 022818(2015).[12] Giulia Cencetti, Franco Bagnoli, Giorgio Battistelli, Luigi Chisci, Francesca Di Patti, and Duccio Fanelli, “Topological stabilization forsynchronized dynamics on networks,” The European Physical Journal B , 9 (2017).[13] Giulia Cencetti, Pau Clusella, and Duccio Fanelli, “Pattern invariance for reaction-diffusion systems on complex networks,” ScientificReports , 16226 (2018).[14] Andrea Baronchelli, Michele Catanzaro, and Romualdo Pastor-Satorras, “Bosonic reaction-diffusion processes on scale-free networks,”Phys. Rev. E , 016111 (2008).[15] P. Turchin, Quantitative Analysis of Movement: Measuring and Modeling Population Redistribution in Animals and Plants , Weimar andNow; 13 (Sinauer, 1998). -0.15-0.1-0.0500.050.1 0 0.5 1 1.5 2 2.5 3 3.5 4 (a) (b) -0.1-0.0500.050.1 0 0.5 1 1.5 2 2.5 3 3.5 4 (c) (d) -0.100.10.20.3 0 0.5 1 1.5 2 2.5 3 3.5 4 (e) (f) ρ h k ih k i D ρρ − p − p D ρρ bb D ρ FIG. 12. Results of the HMF analysis for the Brusselator model. Top panels show the dependence of ρ on D where each line corresponds toa specific parameter value. Bottom panels show heatmaps of the amplitude upon tuning the diffusion and the parameter. System parametersare a = 1 , d = 0 . , and b = 0 . , except in panels (e,f) where b is modified as indicated. Recall that ρ := max j ( Re [ λ j ]) . (a,b) Dependenceof ρ on the diffusion D in a range of ER networks with different average degree k . (c,d) Dependence of the oscillation amplitude ρ on thediffusion D for a range of Watts-Strogatz networks with rewiring probability p ∈ [0 , . . (e,f) Dependence of the oscillation amplitude ρ onthe diffusion D for different values of the system parameter b in a single ER network with (cid:104) k (cid:105) = 20 .[16] A. Okubo and S.A. Levin, Diffusion and Ecological Problems: Modern Perspectives , Interdisciplinary Applied Mathematics (SpringerNew York, 2013).[17] Naoki Masuda, Mason A. Porter, and Renaud Lambiotte, “Random walks and diffusion on networks,” Physics Reports , 1 – 58(2017), random walks and diffusion on networks.[18] Giulia Cencetti, Federico Battiston, Duccio Fanelli, and Vito Latora, “Reactive random walkers on complex networks,” Phys. Rev. E ,052302 (2018).[19] Aneta Koseska, Evgeny Volkov, and Jürgen Kurths, “Oscillation quenching mechanisms: Amplitude vs. oscillation death,” PhysicsReports , 173 – 199 (2013), oscillation quenching mechanisms: Amplitude vs. oscillation death.[20] Y. Kuramoto, Chemical Oscillations, Waves and Turbulence (Springer, Berlin, 1984).[21] I. Prigogine and R. Lefever, “Symmetry breaking instabilities in dissipative systems. ii,” The Journal of Chemical Physics , 1695–1700(1968), https://doi.org/10.1063/1.1668896.[22] MR Guevara, L Glass, and A Shrier, “Phase locking, period-doubling bifurcations, and irregular dynamics in periodically stimulatedcardiac cells,” Science , 1350–1353 (1981), https://science.sciencemag.org/content/214/4527/1350.full.pdf.[23] D C Michaels, E P Matyas, and J Jalife, “Mechanisms of sinoatrial pacemaker synchronization: a new hypothesis.” Circulation Research Rhythms of the Brain (Oxford University Press, 2006).[25] Eugene M. Izhikevich and Gerald M. Edelman, “Large-scale model of mammalian thalamocortical systems,” Proceedings of the NationalAcademy of Sciences , 129 – 139 (2006).[27] Martin Baurmann, Thilo Gross, and Ulrike Feudel, “Instabilities in spatially extended predator–prey systems: Spatio-temporal patternsin the neighborhood of turing–hopf bifurcations,” Journal of Theoretical Biology , 220 – 229 (2007).[28] D.G. Aronson, G.B. Ermentrout, and N. Kopell, “Amplitude response of coupled oscillators,” Physica D: Nonlinear Phenomena , 403– 449 (1990).[29] Renato E. Mirollo and Steven H. Strogatz, “Amplitude death in an array of limit-cycle oscillators,” Journal of Statistical Physics ,245–262 (1990).[30] D. V. Ramana Reddy, A. Sen, and G. L. Johnston, “Time delay induced death in coupled limit cycle oscillators,” Phys. Rev. Lett. ,5109–5112 (1998).[31] Keiji Konishi, “Amplitude death induced by dynamic coupling,” Phys. Rev. E , 067202 (2003).[32] A. N. Samukhin, S. N. Dorogovtsev, and J. F. F. Mendes, “Laplacian spectra of, and random walks on, complex networks: Are scale-freearchitectures really important?” Physical Review E , 036115 (2008).[33] S.H. Strogatz, Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering (CRC Press, 2018).[34] J.D. Murray,
Mathematical Biology II: Spatial Models and Biomedical Applications , Interdisciplinary Applied Mathematics (SpringerNew York, 2006). [35] Aneta Koseska, Evgenii Volkov, and Jürgen Kurths, “Transition from amplitude to oscillation death via turing bifurcation,” Phys. Rev.Lett. , 024103 (2013).[36] James T. Tanner, “The stability and the intrinsic growth rates of prey and predator populations,” Ecology , 855–867 (1975).[37] Duncan J. Watts and Steven H. Strogatz, “Collective dynamics of ’small-world’ networks,” Nature , 440–442 (1998).[38] Garima Saxena, Awadhesh Prasad, and Ram Ramaswamy, “Amplitude death: The emergence of stationarity in coupled nonlinearsystems,” Physics Reports , 205–228 (2012).[39] R. Pastor-Satorras and A. Vespignani, “Epidemic spreading in scale-free networks,” Phys. Rev. Lett. , 3200–3203 (2001).[40] M. Boguñá, C. Castellano, and R. Pastor-Satorras, “Langevin approach for the dynamics of the contact process on annealed scale-freenetworks.” Phys. Rev. E , 036110 (2009).[41] Albert-László Barabási and R. Albert, “Emergence of scaling in random networks,” Science , 509–512 (1999).[42] J.D. Murray, Mathematical Biology: I. An Introduction , Interdisciplinary Applied Mathematics (Springer New York, 2011).