Amplitude death and restoration in networks of oscillators with random-walk diffusion

We study the death and restoration of collective oscillations in networks of oscillators coupled through random-walk diffusion. Differently than the usual diffusion coupling used to model chemical reactions, here the equilibria of the uncoupled unit is not a solution of the coupled ensemble. Instead, the connectivity modifies both, the original unstable fixed point and the stable limit-cycle, making them node-dependent. Using numerical simulations in random networks we show that, in some cases, this diffusion induced heterogeneity stabilizes the initially unstable fixed point via a Hopf bifurcation. Further increasing the coupling strength the oscillations can be restored as well. Upon numerical analysis of the stability properties we conclude that this is a novel case of amplitude death. Finally we use a heterogeneous mean-field reduction of the system in order to proof the robustness 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 a wide class of natural and man made systems that can be fruitfully represented in terms of complex networks [1]. Under this representation, in which nodes stand for units and edges for pairwise unit interactions, a natural framework arises that is capable to unify the functional and structural properties of a wide variety of different systems. Particularly important at this respect are those systems that are the substrate for some dynamic of transport process, whose properties can be strongly impacted by the topology 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 them inside 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 of this formalism can be found, for example, in the study of ecological species dispersal [4] or epidemic spreading [5]. In these studies, the network structure represents a metapopulation [6], defined as a group of populations or physical patches, joined by migration 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 the homogeneous equilibrium, preserved by the diffusion, can be established by means of the dispersion relation, which, in the complex network case, relies on the diagonalization of a Laplacian operator [7,9]. As a result, Turing patterns are observed to arise in networked substrates, characterized by a node-dependend pattern of species density [7]. The observation of Nakao and Mikhailov has spurred a flurry of activity in the field of reaction-diffusion processes on networks, leading to the study of the 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 chemical reactions and heat transport. In this case, the exchange rate of the physical quantities between two nodes is proportional to their density difference. A different choice for the diffusive coupling that can be considered on complex networks corresponds to random 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 diffusion is particularly relevant in the case of ecology dynamics where each node represents a metapopulation of a certain species in an ecosystem, which then might randomly migrate to the surrounding environments [15,16]. Random walks have been extensively studied 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, randomwalk diffusion does not accept, in general, a homogeneous equilibrium, hence the steady states are topology dependent. As direct implication of this situation, even if an equilibrium point of the network is known, its stability cannot be attained by means of a dispersion-relation. Overall, the emergent dynamical patterns that might arise in reaction systems with random-walk diffusion have been hitherto unexplored. In this paper we investigate one of such collective states which is forbidden in reaction systems with Fick's diffusion: the quenching of the oscillatory dynamics through the mechanism known as Amplitude Death (AD), as opposed 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 indeed used as a proxy to study many relevant phenomena such as chemical reactions [20,21], cardiac cells [22,23] neural dynam-ics [24,25], and ecological fluctuations [26,27]. Amplitude Death and Oscillation Death are the two main routes through which the coupling among the oscillatory units leads to a state of steadiness [19]. Although the two mechanisms have been prone 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 system falls into the steady state created by the coupling. Nevertheless, the (unstable) oscillatory solution and the inhomogeneous fixed point coexist [11,21]. On the other hand, AD occurs when different coupled oscillators pull each other out of the limit-cycle upon increasing the interaction strength. Thus, the amplitude of the oscillations diminish until it completely vanishes and the oscillators fall into the homogeneous fixed point of the system. Therefore, in this case, the oscillatory solution collapses into the steady state and there is no coexistence [28][29][30].
In oscillatory systems coupled with Fick's diffusion it is possible to compute the dispersion relation associated to the homogeneous time-varying solution, whose instability leads to Turing patterns. If these patterns are steady in time, then we are before a case of OD [11]. Nevertheless, as we review in detail in this paper, the AD phenomena is forbidden in networks with identical 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], delayed interactions [30] or dynamic coupling [31].
Choosing random-walk diffusion strongly modifies this scenario. Here we show that an increase of the diffusion strength diminishes the amplitude of the oscillations until they collapse into an inhomogeneous steady state. This phenomena differs from OD in the sense that there is no coexistence between the oscillatory solution and the fixed point. We show that the stationary solution 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. Here we extensively study this transition towards AD and the later restoration of the oscillations, as well as we show how suitable modifications of the network topology lead towards the disappearance of AD. We also perform heterogeneous mean-field analysis in 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 where D x and D y are the coupling (diffusion) coefficients, f and g are non-linear reactive terms, and the matrix ∆ ij is the discrete Laplacian operator specifying the diffusive transport of species between connected nodes. If diffusive transport is ruled by Fick's law, that is, by the sum of fluxes of incoming species at each node, where the flux is assumed to be proportional to 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 1 when nodes i and j are connected, and zero otherwise [1]. In this case the equations of motion ruling the dynamics of the ith node of the network read (1) For this gradient-driven diffusion scheme, any solution of the uncoupled system (D x = D y = 0) corresponds to a solution of the coupled system. Indeed, the diffusion term only depends on the density difference between connected nodes, so if all nodes evolve with the exact same dynamics, (x j (t), y j (t)) = (x(t), y(t)) for all j = 1, . . . , N , then the diffusive coupling vanishes and the solution is preserved. Such solutions are homogeneous, meaning that the dynamic evolution of the system is the same for 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) periodic solutions, also referred to as limit-cycles [33].
The addition of the coupling diffusive terms can spontaneously modify the stability of these homogeneous solutions. In principle, the stability of a homogeneous fixed point of system (1) boils down to study the eigenvalues λ i of a 2N × 2N Jacobian matrix. Nevertheless it is possible to relate the eigenvalues of such high-dimensional operator to the eigenvalues of the Laplacian matrix of the system by means of a dispersion relation which maps all the Laplacian matrix eigenvalues Λ j to the eigenvalues of the system Jacobian λ j . This corresponds to the extension 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 worked out 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 0 , y (0) ) as a function of the eigenvalues of the Laplacian Λ j for the case of an homogeneous Erdös-Rényi network. As it is well known, if the network has a single connected component, the associated Laplacian always contains a single zero eigenvalue corresponding to a uniform eigenvector, and the rest 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 case depicted in Figure 1(a), where the system eigenvalue corresponding to the uniform eigenvector has Re[λ 0 ] = 0.1. On the other hand, if the fixed point is originally stable, the rest of the eigenvalues associated to the strictly negative Laplacian eigenmodes might have a positive real part, thus destabilizing the homogeneous solution and generating spatio-temporal patterns. This is the well-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 be obtained by means of Floquet theory [11]. Figure 1(b) depicts the relation between the Floquet exponents µ j of the limit-cycle for the Brusselator system and the eigenvalues of the network Laplacian. The Floquet exponents corresponding to the uniform eigenvector are µ + 0 = 0 and µ − 0 −0.216, but there are many other Floquet exponents associated to different network Laplacian modes that are positive. Thus, in this case, the originally stable limit-cycle is being destabilized through a Turing instability, so an arbitrary perturbation of the homogeneous solution will develop heterogeneous patterns. If these heterogeneous patterns are stationary, then we would be before a case of OD [11]. Regardless of whether this is the case or not, the AD phenomena will always be forbidden here, since the instability of the original fixed point is preserved when the coupling is added, and the also unstable 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 represent the 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. Representing the pattern of connections between patches by a network, and assuming that animals move at random between patches, we have that population of node i diffuses uniformly through its k i neighbors, each receiving a flux proportional to 1/k i . The Laplacian matrix can thus be written as∆ = (∆ ij ), with∆ ij = a ij /k j − δ ij . The reaction-diffusion process then reads Here the solution of the uncoupled system (D x = D y = 0) is not a solution of the coupled system unless we have a regular network, 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 of a dispersion relation and the arguments exposed above still hold. However, for generic non-regular networks there is no theory that 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 system defined by Eq. (2) in an Erdös-Rényi (ER) network with average degree k = 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 generation of static Turing patterns in models with Fick's diffusion requires D x = D y [34]. In Figure 2 we show the evolution of the variable 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.02 does not strongly affect the behavior of the system. The limit cycle defined by the variables in each node oscillates with a very similar amplitude and period, this last one close to the estimate for an isolated Brusselator, T 2π/ √ a. The different oscillators fluctuate however out-of-phase. Upon increasing the coupling, (Figure 2(b), corresponding to D = 0.3), the oscillators evolve in periodic phase-synchrony, but each having a different amplitude of the limit-cycle, amplitude that is related to the degree of the corresponding node in the network. Increasing the coupling to D = 2, (Figure 2(c)), the system reaches a steady state, the heterogeneous oscillations being substituted by an heterogeneous fixed point, whose value is also related to the node's degree. Further increasing of the coupling up to D = 3.5 ( Figure 2(d)) restores the in-phase synchrony, with node dependent amplitude. These preliminary observations indicate the existence of an oscillation quenching phenomenon, and a later oscillation restoration. In the oscillatory regimes for D > 0, the nodes 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 whole network, and measure the temporal average y and standard deviation σ = σ(y), defined as 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 y as a function of the diffusion coefficient D. The associated error bars 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 for diffusion larger than D 1 0.3775. The system is frozen in a heterogeneous steady state until the oscillations are restored again for a diffusion D larger than D 2 3.095. We are therefore in front of a case of oscillation quenching occurring at D 1 , followed by a subsequent oscillation restoration when D crosses D 2 . Whereas the inhomogeneity of the fixed point would support the idea that the phenomenology corresponds to OD, the two transitions between oscillatory and steady regimes strongly resemble the AD phenomena, where the periodic solutions collapse into the fixed point, which in this case happens to be heterogeneous as opposed to the common cases of AD [19,35]. In the folowing we study the two bifurcations through which the oscillations vanish and restore in order to unveil the mechanism responsible for the cessation of the oscillations.  5 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000

A. Heterogeneous fixed points
In order to investigate the nature of these transitions, we start considering the underlying fixed points, corresponding to the solutions of the nonlinear system 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 > 0, we obtain an inhomogeneous solution, 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 values which correspond to nodes with the same degree. It is worth noticing that we find only one solution of Eq. (4) for each value of D. Overall, it looks like the random-walk diffusion induces the hetoregeneity of the equilibria, which becomes homogeneous only for D = 0.

B. Dispersion of the fixed point
The previous numerical analysis points out that the heterogeneous equilibrium of the network is linked with the equilibrium of the uncoupled system for small diffusion values, thus indicating that the fixed point of the network corresponds to a coupling induced modification of the original steady state. We validate this assumption by exploring the effect that small modifications of the coupling strength D have on the fixed point.
) 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, > 0. Expanding up to first order terms in equation (4) one obtains (see Appendix C for a detailed derivation) the system of equations for i = 1, . . . , N , where J(x, y) is the 2 × 2 Jacobian matrix of the (uncoupled) reactive field, (f (x, y), g(x, y)). Such equations form 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 equation using Euler's method, starting from the uncoupled system equilibrium, leads to the average activity reported in black crosses in Figure 3(a), whereas black squares in Figure 3 words, the steady state is not a completely new state induced by the coupling, but a smooth transformation of the original equilibrium of the system, which turns out to be node-dependent as soon as D > 0. Since the steady state associated with OD corresponds to new states created by the coupling, the observation that here the fixed point is not new is key to classify the observed 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 by numerically computing the eigenvalues and eigenvectors of the full system 2N ×2N Jacobian. Figure 3(c) shows the dependence of the real part of the maximum eigenvalue as the diffusion D is tuned and Figure 3(d) shows the full spectra for three different values of the diffusion. As it happens with the fixed point, for small D all the eigenvalues are located close to the eigenvalues of the uncoupled system λ ± = 0.1 ± 0.7i. Thus, there is at least a pair of complex conjugate eigenvalues with positive real part. As D increases, most of the eigenvalues rapidly spread and are pushed to the left of the complex plane. Nevertheless, a pair of complex conjugate eigenvalues remain isolated from the rest and do not cross the imaginary axis until D 1 (see isolated symbols close to the imaginary axis in figure 3(d)). Upon further increasing the diffusion, the same isolated pair of eigenvalues crosses again the imaginary axis at D 2 , signaling the restoration of the oscillations (see Fig.3(c)). According to these results it is clear that both transitions, D 1 and D 2 , are Hopf supercritical bifurcations. Indeed, as shown in figure 4, the amplitude of the oscillations vanishes as σ |D − D 1 |, as expected for a supercritical Hopf bifurcation, and is restored then again at D 2 with the same exponent. Also, the limit cycle centroid of each node approaches the fixed point solution as D approaches the bifurcation 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 a saddle-node bifurcation [19]. On the other hand, OD is associated to a symmetry breaking pitchfork bifurcation, which allows for the coexistence of the oscillatory solution and the steady state [19]. Here there is no such symmetry breaking, the inhomogeneity of the fixed point being induced by the heterogeneous coupling structure, that is not arbitrary. Moreover, the limit-cycle perishes at the bifurcation points, and thus no coexistence between the two solutions is possible. This set of observations altogether indicates that the oscillation quenching mechanism presented here falls into the category of Amplitude Death, with the novelty that 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 same phenomenology using the Holling-Tanner predator-prey system [36].

A. Network density
In order to clarify the impact that irregular network topologies have on the stabilization of the fixed point and the later restoration of the oscillations we have performed numerical simulations of the Brusselator system Eq. (B1) with random-walk diffusion in different classes of complex networks (in the Appendix F we report results corresponding to different reaction terms).
Since an heterogeneous network topology is key for the emergence of the amplitude death, we first consider ER networks with different average degree k . In order to avoid fluctuations due to the independent generations of the topology, we construct networks starting form an initial ER configuration with given size and k = 6, and progressively add at random new connections in order to generate denser topologies with largest average degree. Figures 5(a,b) show the amplitude of the mean-field oscillations for different D using ER networks with increasing average degree. In Figure 5(a) we can observe that AD takes place only in sparse networks, as increasing the connectivity favors the oscillation amplitude. In Figure 5(b) the black region corresponding to AD clearly vanishes smoothly with increasing k . In fact, for the smaller values of the average degree, oscillations are not restored even for very large values of D. As more connections are added to the topology, the stable steady state shrinks, until it completely vanishes for k 45. 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 and restoration. Indeed, for regular networks (where all nodes have the same degree), no quenching can emerge, whereas irregularity seems to induce AD. To unveil this situation, we repeat the same analysis on small-world networks generated according to the Watts-Strogatz (WS) algorithm [37] (see Methods). In this case the network topology depends on a rewiring parameter p ∈ [0, 1]. For p = 0 the generated architecture corresponds to a regular network and the stability of the limit-cycle can be studied through the dispersion relation. In this case, for p = 0, no Turing-pattern is triggered and the fully synchronized solution is the only stable attractor. Instead, p = 1 corresponds to a random network similar to the ER model, and the behavior of the system does not differ much from the results reported in figure 3(a). Therefore, the behavior observed for intermediate values of p should reveal 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 the topology generation algorithm, we construct the networks by progressively rewiring a larger fraction of random edges in the same network. We use 1 − p as topological order parameter, measuring the network departure from randomness , i.e., random for 1 − p = 0 and regular for 1 − p = 1. As expected, for 1 − p = 0 the situation is similar to the ER topology with k = 20: there is a first bifurcation towards steadiness for small D, and a second bifurcation where oscillations are restored. For larger regularity, the two bifurcation points come closer together until they collide and vanish for p 0.36. Overall the scenario is analogous to that seen in ER networks with increasing density, but now with the disorder parameter 1 − p playing the role of control parameter, instead of the average degree. Both scenarios are also reminiscent of the amplitude death bifurcation observed in systems of oscillators with quenched heterogeneity [29,38]. In our case, however, the heterogeneity resides in the connectivity structure 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) shows the outcome of numerical simulations for the same ER network with k = 20, keeping fixed the parameter a = 0.5 and varying b. For b < a + 1 = 1.5 the reactive part does not display oscillations and the fixed point is stable. For b > a + 1 = 1.5 the uncoupled system undergoes a Hopf bifurcation and the stable limit-cycle emerges. Here we investigate the dynamics of the coupled system for b varying from 1.52 to 2.1. Not surprisingly, for values of b close to the bifurcation of the uncoupled system the amplitude death regime is quickly attained at small values of D, and the restoration of the oscillations is not triggered even for large coupling, i.e., the fixed point stabilizes quickly with D. As b increases, the limit-cycle has a larger amplitude and the diffusion-induced stabilization of the fixed point requires larger diffusion values. Indeed, the region of amplitude death becomes smaller as they move away from the single-node bifurcation point, until it completely vanishes for b 1.9 in the same manner as 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 apply the standard tool of heterogeneous mean-field (HMF) theory [3,39], specialized to reaction-diffusion processes [14]. The basis of HMF consists in the annealed network approximation [3,40], that replaces the static adjacency matrix of a real network by an average over degree classesā ij that, in the case of uncorrelated networks, takes the form Introducing this expression into Eq. (2), we obtain the HMF dynamical equations wherek x j and y = 1 N N j=1 y j are the average (mean-field) activities of x and y variables. Within this framework, the dynamics of each node depends only on its own degree, the mean-field values x and y, and the average connectivity of the network. In fact, one can assume then that all nodes with the same degree behave identically, thus formally reducing the 2N -dimensional system (6) to a 2n-dimensional system where n is the number of different degrees in the network [3]. An additional interest of this approach lies in the fact that it allows to consider network sizes much larger than those permitted in a direct numerical integration, since in general n N .
System (6) can be solved semi-analytically assuming that x and y are parameters that are to be determined self-consistently a posteriori (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 black squares 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 accuracy that of the actual system (red continuous curve). Moreover, the stability analysis of the fixed point also coincides with that of the original system, with a cloud of eigenvalues lying far on the stable part, and a pair of complex conjugate eigenvalues crossing back and forth the imaginary axis upon modifying D. Thus, the mean-field reduction does not only provide a good proxy to study the fixed point, but also allows to study the amplitude death phenomena in terms of stability.

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 to larger systems. Figure 6(d) shows the largest eigenvalue real part for different values of D as obtained from the heterogeneous mean-field analysis. We analyzed 4 sets of ER networks of size N = 10 3 , 10 4 , 10 5 , and 10 6 , with different average degrees k = 20, 30, 40, and 50. For a specific value of k , we observe that networks with different sizes produce qualitatively similar results, converging to a well defined limit for sufficiently large N . For instance, for k = 20 the upper red curve, corresponding to N = 10 3 , shows a smaller region of amplitude death, whereas the red curves below, corresponding to larger systems, converge nicely upon increasing N . The same situation is repeated for the other values of k , thus confirming that the effects of the average degree on the amplitude death are not finite-size but rather robust. Indeed, the overall conclusion of this analysis 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 the Holling-Tanner system.

VI. DISCUSSION
Reaction-diffusion processes are a powerful formalism to represent general dynamical processes on networks, in which particles or species interact inside nodes while moving diffusively between pairs of connected nodes. By analogy with chemical reactions, a gradient-driven diffusion term, given by Fick's law, is usually assumed. However, in certain circumstances, a random walk diffusion term, in which particles jump at random along edges, might be more realistic. Here we have shown how the nature of the diffusion term can alter the behavior of the limit cycles in oscillatory reaction-diffusion processes when driven by 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 amplitude death quenching characterized by a inhomogeneous steady state induced by the structure of the diffusion operator. In this case, a further increase in diffusion restores the oscillations, in the form of a set of limit cycles with different amplitude for each 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 the Brusselator systems in random networks, are confirmed by a linear stability analysis of the eigenvalues of the diffusion operator and by a seminumerical heterogeneous mean-field approximation, which allows us to check the phenomenology in very large networks. The role of random walk diffusion in the amplitude quenching of oscillations, which is observed for different classes of homogeneous networks, different sets of parameters and even different reaction terms, is thus a robust phenomenology, which could be relevant in processes such as epidemic spreading or ecological dispersion. Future promising work along these lines would 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 check Murray's book [34]. Here we summarize the results of Nakao et al. [7] on reaction-diffusion systems in complex networks. Let us consider a reaction-diffusion system in a complex network with a gradient-driven diffusion term with a Laplacian matrix ∆ ij = a ij − k i δ ij . This Laplacian matrix is semi-definite negative, and its eigenvalues are all real and non-positive. Let (x (0) , y (0) ) be a fixed point of the uncoupled system (D x = D y = 0), and hence the homogeneous solution of the coupled system. Inserting an arbitrary perturbation {(δx i (t), δy i (t))} N i=1 of the homogeneous solution to the equations and retaining up to first order terms one gets that the evolution of the perturbation behaves as where J is the Jacobian of the homogeneous system, evaluated at (x (0) , y (0) ). Let Φ (α) = (φ We can express the perturbation of the homogeneous solution in terms of the basis of the eigenvectors as Applying this change of coordinates on Eq. (A1) and making use of the relation (A2) and of the linear independence of the eigenvectors {Φ (α) } N α=1 one obtains that the evolution of (u (α) ,v (α) ) T becomes independent for each α = 1, . . . , N through the relation ∂f ∂y x (0) , y (0) ∂g ∂y (A3) Thus the stability of the homogeneous solution simplifies to the study of the eigenvalues (and eigenvectors) of the previous 2 × 2 matrix, which are obtained as functions of the Laplacian eigenvalues Λ α . See below for an explicit application to the Brusselator system. The extension of this analysis to limit-cycles makes use of Floquet theory, see [11].
Dividing then both sides of these equations by one finally can write the differential system that rules the dependence of (x where J(x, y) is the 2 × 2 Jacobian matrix of the (uncoupled) reactive field, (f (x, y), g(x, y)). The previous equation is thus an implicit linear non-autonomous system of differential equations that can be solved by means of usual numerical integrators using as initial condition for D = 0 the solution of the uncoupled oscillators.
Appendix D: Network models

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, 1]. The average degree of the network is then pN and the degree distribution follows a binomial distribution with parameters N − 1 and p.

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 halfway between random and regular topologies. The model starts with N nodes distributed on a ring, each node being connected to its k nearest neighbors. Then, every edge of the network is rewired with probability p ∈ [0, 1]. The small-world network class is represented for small values of the rewiring probability, whereas for p close to 1 the topology becomes closer to that of a ER network.

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 obtains that (dropping the subindices for simplicity), we obtain from where one can obtain a cubic equation for the value of x: − 9a 2 D(D + 1)(b + 1 + D)(1 + D(x + y)) + 27a 2 (D + 1)D(1 + Dx) .
The Brusselator system has the property that, in the fixed point and even without using any kind of approximation, x = 1. However, in order to obtain y self-consistently one needs to rely on numerics. In practice, one starts with an educated guess for the mean-field y and determines then all x j and y j using equations (E2) and (E1). Using a bisection method one can reduce the error 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 the Holling-Tanner predator-prey system [36]. In this model x ≥ 0 and y ≥ 0 respectively represent the population fraction of a prey and predator species in an ecosystem. On its dimensionless form, the corresponding reactive terms read where a, b, and d are system parameters [42].
The system has a fixed point at and the corresponding Jacobian is The resulting eigenvalues have positive real part if and only if Fixing a = 1 and d = 0.1, the fixed point becomes unstable through a supercritical Hopf bifurcation at b c 0.2625..., so that for b < b c the single node dynamics exhibit periodic oscillations. Unless otherwise stated, we consider here b = 0.25.

Amplitude death and restoration
Now we consider a network of Holling-Tanner oscillators coupled with random-walk diffusion. Let A = (a ij ) be the adjacency matrix of the (symmetric and fully connected) network, and let k j be the degree of node j. The random-walk Laplacian matrix is then defined as∆ := (∆ ij ) where∆ ij = a ij /k j − δ ij . With this notation, the dynamical rule governing the evolution of the ith node of the system reads where D is the diffusion or coupling strength. (see green crosses, blue pluses and red circles in Figure 7(b)). The analysis of the dispersion of the fixed point, explained in detail in the paper, unveils again that the steady state corresponds to a diffusion-induced modification of the single-node solution, (x (0) , y (0) ) (0.27, 0.27). Indeed, the black crosses in Figure 7(a) and black squares in 7(b) indicate the goodness of the dispersion equation to obtain the underlying fixed point. Finally, stability analysis reveals again a Hopf bifurcation towards the 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 terms for both systems. Figure 8 shows the same instances reported in the paper where topology and parameter modification vary the scenario depicted in the previous section. In particular, an increase of the density in the network (see Figures 8(a,b)) reduces the region of steadiness, until both critical points collapse. The same occurs for the WS networks upon increasing the regularity of the networks (see Figures 8(c,d)). Finally, the same results are obtained by keeping the topology fixed and tuning the system parameter b. Closer to the single-node bifurcation point b c , the coupling rapidly tends to kill the periodic behavior, whereas increasing b finally vanishes the AD phenomena (see Figures 8(e,f)).

Heterogeneous mean-field approach
The heterogeneous mean-field approach also provides a good agreement with the original simulations, although here the self-consistency of the problem needs further numerical effort to be solved, since the nonlinearities in the formulation of the system limit the analytical results. Let us recover system (1) but we now substitute the elements of the adjacency matrix a ij by the annealed network approximation, After applying this transformation in the coupling terms of Eq. (1) and doing some straightforward calculations one obtains the mean-field system wherek i := k i / N j=1 k j , and x = 1 N N j=1 x j and y = 1 N N j=1 y j are the mean-field activities of x and y variables. Here, in order to solve for the fixed point,ẋ i = 0 andẏ i = 0 one needs to find the roots of a 5th order polynomial, thus, unlike the Brusselator model, the expressions of x i and y i cannot be explicitly obtained. Therefore, given a value of x and y we obtain the fixed points x  Figure 9 shows the results corresponding to this analysis. Again, the heterogeneous mean-field approach provides a good representation of the underlying scenario with ER. Therefore, we implement this dimensionality reduction approach to study the amplitude death and restoration in large networks. Again, as it happens in the Brusselator system, increasing the network size does not effect the observed phenomenology. Therefore, we finally repeat the analysis depicted in Figure 8 but using networks with N = 10 6 nodes analyzed through HMF, using the maximum of the largest eigenvalue of the system, as a proxy to identify whether there oscillations or not. Results in Figure 10.  the peak is only slightly shifted around a similar value, thus the global dynamics corresponds to irregular fluctuations evolving at similar frequencies. For D = 0.3 and D = 3.5 the spectra shows a sharp peak, indicating that the frequency entrainment is 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 1 . The centroid are computed as the time-averaged value over a long enough time series, y j . Figure 11(c) shows the decrease of the amplitude of each node's oscillations until it vanishes at D 1 . Finally, figure 12 displays the results of the HMF analysis for networks with N = 10 6 analogous to figure 5 of the main paper.