Dispersal-induced destabilization of metapopulations and oscillatory Turing patterns in ecological networks

As proposed by Alan Turing in 1952 as a ubiquitous mechanism for nonequilibrium pattern formation, diffusional effects may destabilize uniform distributions of reacting chemical species and lead to both spatially and temporally heterogeneous patterns. While stationary Turing patterns are broadly known, the oscillatory instability, leading to traveling waves in continuous media and also called the wave bifurcation, is rare for chemical systems. Here, we extend the analysis by Turing to general networks and apply it to ecological metapopulations of biological species with dispersal connections between habitats. Remarkably, the oscillatory Turing instability does not lead to wave patterns in networks, but to spontaneous development of heterogeneous oscillations and possible extinction of some species, even though they are absent for isolated populations. Furthermore, our theoretical analysis reveals that this instability is more common in ecological metapopulations than in chemical reactions. Indeed, we find the instabilities for all possible food webs with three predator or prey species, under various assumptions about the mobility of individual species and nonlinear interactions between them. Therefore, we suggest that the oscillatory Turing instability is generic and must play a fundamental role in metapopulation dynamics, providing a common mechanism for dispersal-induced destabilization of ecosystems.


I. INTRODUCTION
A rich variety of nonequilibrium pattern formation is supported by reaction-diffusion processes.One of the universal mechanisms of such pattern formation is provided by the Turing instability [1]; a diffusion-induced instability of the homogenous state which leads to spontaneous development of self-organized patterns.The Turing instability can play an important role in biological morphogenesis and it has been extensively studied in various applications, including biological [2][3][4][5][6], chemical [7,8] and physical systems [9].The classical Turing instability leads to the establishment of stationary spatial patterns.However, the oscillatory analogue of this instability is possible and it has also been discovered by Alan Turing [1].This oscillatory intability produces traveling or standing waves and therefore it is often called "the wave instability" (see [10]).He has shown that at least three species are needed for the oscillatory instability, while the stationary instability is possible already with two species.The stationary Turing instability has been extensively studied both theoretically [2][3][4]9]and experimentally [5][6][7][8], whereas the oscillatory instability is more rare and it was found only for special chemical systems [11][12][13].Note that complex spatio-temporal patterns can also emerge as a result of interactions between the stationary Turing instability and other bifurcations.For example, near the Turing-Hopf bifurcation point, complex mixed modes leading to standing waves and spatio-temporal chaos can exist [14].However, their mechanism is different from that of the oscillatory Turing instability (or the wave bifurcation).
Reaction-diffusion processes are also characteristic to ecological systems.The reactions correspond in this case to the predator-prey and other interactions between the species.Both passive diffusion and active random migration are possible in ecological populations.Moreover, there are situations when a population occupies a large habitat and therefore can be considered as an extended spatial system.The classical Turing instability is possible in ecosystems.It has been shown that such instability should be generic for the two-species predator-prey models [15] (see also [16]).Complex spatio-temporal ecological dynamics related to the Turing-Hopf [15,17] and Turing-Takens-Bogdanov [15] bifurcations has been discussed.Stationary Turing patterns have been found in realistic models describing plantparasite [18], plankton-fish [19] and plant-insect [20] interactions.The oscillatory Turing instability is also possible in ecology.In a study of a three-species plant-parasite-hyperparasite system, such instability leading to standing waves, has been previously considered [18].
While some ecological systems can be described by reaction-diffusion equations for continuous media, there are also many ecosystems that are spatially fragmented and represent networks [21][22][23].Such networks are formed by individual habitats which are linked by dispersal connections.Ecological species populate the habitats and diffusively migrate over a network.Such network-organized ecosystems are known as metapopulations [24][25][26].The metapopulation concept has been applied to describe and investigate real ecological systems (see, e.g.[27][28][29]).It has also been used in the context of the epidemic research [30][31][32][33].In the framework of the metapopulation concept, the role of dispersal connections in the stability enhancement of an ecosystem (the rescue effect) has been discussed [25] (see also [28]).The theoretical results have been tested in the experiments with specially prepared metapopulations [27,29].
Ecological metapopulations provide examples of reaction-diffusion systems with a network structure.Such systems can however be also found in other research fields.For instance, a biological embryo can be viewed as a network of cells with the chemicals diffusing over the pattern of intercellular connections [34][35][36].Networks formed by coupled reactors can also be considered [37,38].Theoretical studies of reaction-diffusion processes on networks have already attracted much attention [39].Effects of infection spreading over the networks have been discussed in detail [30][31][32][33].The role of network topology on the phase diagrams of nonequilibrium phase transitions on networks has been considered [40].Traveling and pinned fronts in networks of diffusively coupled bistable elements were analyzed [41] and control of front propagation by global feedback has been considered [42].There is large literature on the networks formed by diffusively coupled oscillators [43,44] (see also [45]).The role of dispersal connections in the synchronization effects in ecological networks has also been discussed [46].
The stationary Turing instability for networks has been first analyzed in 1971 by Othmer and Scriven [34].The authors have introduced a general mathematical description of the classical Turing instability in two-component reaction-diffusion networks and have applied their theory for regular lattices (see also [35]).Stationary Turing patterns in small networks of coupled chemical reactors have subsequently been discussed [37].The properties of such instability and of the final established stationary patterns in large random networks of diffusively coupled activator-inhibitor elements have been investigated and the mean-field theory of Turing patterns in such network systems has been constructed [9].The global feedback control of the stationary Turing patterns in networks has been studied [47].A detailed mathematical analysis of the hysteresis phenomena related to the network Turing bifurcation has recently been performed [48].
In the present article, we provide, for the first time, a general mathematical theory of the oscillatory Turing instability (the analogue of the wave bifurcation) in network-organized reaction-diffusion systems and apply the theory to large ecological networks.Our numerical simulations are performed for the metapopulations with various threespecies food webs under different assumptions about the nonlinearities of the population dynamics.The instability could be found for all such systems and, as we therefore believe, it should be common in ecology.In contrast to the wave instability in continuous media, traveling or standing waves do not develop and oscillations, localized on a subset of network nodes, are instead observed.The bifurcation is always supercritical and therefore the final pattern is usually well described by the first critical mode.This diffusion-induced instability leads to destabilization of metapopulations and the extinction of some species may be its result.
This article is organized as follows: In section II, the addressed problem is mathematically formulated.General descriptions for ecological networks are introduced and specific mathematical models which we employ are explained.In section III, analytical investigations, including the linear stability analysis for general reaction-diffusion networks with three reacting species, is undertaken.Sufficient conditions for the oscillatory Turing instability in ecological systems are constructed.Numerical investigations for various ecological models are performed in section IV.First, examples of oscillatory and stationary patterns in networks are demonstrated.After that, a detailed analysis of the simulation results is performed.As we find, the oscillatory bifurcation is always supercritical and the small-amplitude patterns are well described by the critical modes which correspond to certain Laplacian eigenvectors.The localization of the Laplacian eigenvectors for networks explains the observed localization of the oscillatory Turing patterns on a subset of network nodes.An example of a secondary instability of the oscillatory patterns, leading to the extinction of one of the species, is given.Finally, we summarizes the results and discuss them in Section V.In Appendices, numerical simulations for additional ecological models and for two chemical network models are reported.
FIG. 1: Food web diagrams.Each arrow goes from the consuming to the consumed species.In the last food web (d), no persistence with fixed values of (ū, v, w) can be achieved; therefore this web is excluded from our analysis.

II. FORMULATION OF THE PROBLEM
We consider ecological networks formed by individual populations which occupy separate habitats and are coupled by dispersal connections.Our attention is focused on the populations which consist of three interacting species.All possible food webs with three different species are displayed in Fig. 1.Generally, such ecological networks are described by equations for i, j = 1, • • • , N , where population densities of species on node i are denoted as are the differences of reproduction (Q) and death (R) rates for each species, and σ u,v,w are the mobilities of the three species; the common parameter ǫ is introduced for convenience, so that the mobility of all species can be varied without changing relative mobilities.The Laplacian matrix L has elements L ij = A ij − j A ij δ ij where A ij is the matrix of connections between the habitats.For simplicity, we assume that all connections, if present, have the same strength.Therefore, we have A ij = 1, if there is a connection between habitats i and j, and A ij = 0 otherwise.We assume that, in absence of diffusive coupling, a stable stationary state (u 0 , v 0 , w 0 ) exists which is determined by F (u 0 , v 0 , w 0 ) = G(u 0 , v 0 , w 0 ) = H(u 0 , v 0 , w 0 ) = 0 where u 0 > 0, v 0 > 0 and w 0 > 0. We examine three different food webs shown in Fig. 1(a)-(c).Below, we provide explicit expressions for the reproduction rates and death rates in each model.While different nonlinear dependence can be considered for predator-prey interactions, we will employ the Holling type II functions [3,49] (See Supporting Information for other dependences).Note that the food web shown in Fig. 1(d) is excluded from our analysis.In such system, two competing species V and W cannot coexist in a steady uniform state and Turing instabilities with three reacting species are impossible.Model A. In the food chain shown in Fig. 1(a), top predator W is feeding on intermediate species V which is in turn a predator for prey U .The collective dynamics of such metapopulation is described by Eqs.(1) with In the food web shown in Fig. 1(b), both species V and W play a role of the predators for prey U while V is also a prey for W .Such system is modelled as follows: Model C. In the food web shown in Fig. 1(c), species U and V are the prey for predator W . Thus we have For comparison, the continuous analog of the model (Eqs.(1)) will be also considered.Then, node variables u i , v i and w i are replaced by space-dependent densities u, v and w and differential diffusion terms like ǫ ′ σ u ∇ 2 u are introduced instead of the last terms in Eqs.(1).

III. ANALYTICAL INVESTIGATIONS A. Linear stability analysis
Below we give the precise definition of the oscillatory Turing instability in ecological networks.The stability of the uniform state of the model can be analyzed by the linear stability analysis.Small perturbations (δu i , δv i , δw i ) are introduced to the steady state as (u i , v i , w i ) = (u 0 , v 0 , w 0 ) + (δu i , δv i , δw i ).Substituting this into Eqs.(1), the following linearized differential equations are obtained: where F u = ∂F/∂u| (u0,v0,w0) , F v = ∂F/∂v| (u0,v0,w0) , F w = ∂F/∂w| (u0,v0,w0) , ... are partial derivatives at the uniform steady state (u 0 , v 0 , w 0 ).The Laplacian eigenvectors { φ (α) } are introduced to decompose the perturbations.They are defined as where Λ (α) is the Laplacian eigenvalue of the α th mode (α = 1, • • • , N ).The mode indices {α} are sorted in the increasing order of the Laplacian eigenvalues The perturbations (δu i , δv i , δw i ) are expanded over the set of the Laplacian eigenvectors as where γ (α) = λ (α) + iω (α) is a complex growth (or decay) rate of the α th eigenmode.Substituting these expressions into Eqs.( 5), the following equations are obtained for each eigenmode: Thus, the decay rate of each eigenmode is determined by the characteristic equation The uniform steady state is stable if λ (α) < 0 for all α.The Turing instability takes place if λ (α) becomes positive at some α = α c which represents the critical mode for the instability.The critical modes are stationary, φ , if ω (αc) = 0. On the other hand, the critical modes can also be oscillatory, φ (αc) i e iω (αc ) t , if ω (αc) = 0.As noticed already by Turing [1], oscillatory instabilities are possible only if the number of species is at least three.
In the continuous case, the diffusion operator is ∇ 2 and its eigenvectors are plane waves, since ∇ 2 e ikx = −k 2 e ikx .Thus, small perturbations are decomposed over plane waves as to obtain the characteristic equation The instability occurs if λ (k) becomes positive at some critical wave number k c .The critical mode corresponds to a traveling wave φ (kc) i e iω (kc ) t if ω (kc) = 0 or to a stationary wave e ikcx if ω (kc) = 0.For continuous media, the oscillatory Turing instability is usually called the wave instability, since the first unstable modes are the traveling waves.In networks, traveling waves do not appear and therefore it is not appropriate to talk about a wave instability in this case.

B. Sufficient conditions for the oscillatory Turing instability in ecological networks
In our recent work [50], we considered general continuous reaction-diffusion systems with three species and derived sufficient conditions for the oscillatory Turing instability in continuous media.The derivation was based on the linear stability analysis.This analysis can be straightforwardly extended to the case of discrete networks.Below, we only reformulate the previously derived sufficient conditions for the networks and apply them to the considered ecological models.
Suppose that the mobilities σ u and σ v of the species U and V are fixed (and at least one of them is non-vanishing) and we gradually increase the mobility σ w of species W .It follows from our previous analysis [50] that the oscillatory Turing instability will always be found at sufficiently high mobility σ w if the following three sufficient conditions are satisfied: Note that one of the species U or V may be immobile, i.e. either σ u = 0 or σ v = 0.
If the mobilities σ v and σ w of the species V and W are fixed (and at least one of them is non-vanishing) and we gradually increase the mobility σ u of species U , the oscillatory Turing instability will always be found at sufficiently high mobility σ u if the following three sufficient conditions are satisfied: Here, one of the species V or W may be immobile so that either σ v = 0 or σ w = 0.When the mobilities σ w and σ u of the species W and U are instead fixed (and at least one of them is non-vanishing) and we gradually increase the mobility σ v of species V , the oscillatory Turing instability will always be found at sufficiently high mobility σ v if the following three sufficient conditions are satisfied: Det Here, one of the species W or U may be immobile so that either σ w = 0 or σ u = 0. Thus, if conditions ( 12), ( 13) or ( 14) are satisfied, the oscillatory Turing instability will be found in the respective ecological networks as the mobilities σ u , σ v or σ w are increased.
Note that the prey death rates R should be increasing functions of the predator densities.Hence, in the food chain (Fig. 1(a)), we have ∂R u /∂v > 0 and ∂R v /∂w > 0 and therefore F v = ∂F/∂v < 0 and G w = ∂G/∂w < 0.Moreover, the predator reproduction rates should generally increase with the prey densities.This implies that ∂Q v /∂u > 0 and ∂Q w /∂v > 0 so that we have G u = ∂G/∂u > 0 and H v = ∂H/∂v > 0. Thus, the following inequalities are satisfied.
In the food web shown in Fig. 1(b), both species V and W play a role of the predators for prey U while V is also a prey for W .Therefore, . Hence, we have ∂R u /∂v > 0, ∂R u /∂w > 0, ∂R v /∂w > 0, ∂Q v /∂u > 0, ∂Q w /∂u > 0 and ∂Q w /∂v > 0. This leads to the conditions In the food web shown in Fig. 1(c), species U and V are the prey for predator W . Now, we have ) and R w = R w (w).Therefore ∂R u /∂w > 0, ∂R v /∂w > 0, ∂Q w /∂u > 0 and ∂Q w /∂v > 0. This implies Thus, in all considered food webs, the inequalities hold and at least one of the inequalities holds strictly.We notice that they imply that the conditions (12b), (13b), and (14b) are all satisfied.Therefore, we can conclude that, in contrast to general reaction-diffusion models, at least one of the three sufficient conditions will be always satisfied for ecological networks with three species.Hence, the oscillatory Turing instability should be more common in ecology, as compared with chemical systems.

IV. NUMERICAL INVESTIGATIONS A. Setup of numerical simulation
Below we numerically investigate the oscillatory Turing instability in ecological networks.The 4th-order Runge-Kutta scheme with time step ∆t = 10 −3 was employed in numerical integration.The simulations were started from the uniform steady state (u 0 , v 0 , w 0 ) with random small perturbations with standard deviations (u 0 , v 0 , w 0 ) × 10 −3 .As a typical example of network architecture, we employed scale-free networks generated by the Barabási-Albert preferential attachment algorithm [51].The network size was N = 50 and the mean degree was k = 8.For convenience, the nodes i were sorted in the decreasing order of their degrees k i = j A ij so that we had To quantify emergent patterns, variation amplitudes for individual nodes i at time t were used; where are the mean network quantities.Time-averaged amplitudes were computed as after discarding the transients with the averaging time T = 100000.

B. Turing instabilities in ecological networks
We first show the results for Model A. The parameters in Eq. ( 2) are fixed at yielding a uniform stationary state (u 0 , v 0 , w 0 ) = (1/2, 1/2, 1/4).Jacobian matrix at the steady state is which satisfies the sufficient conditions (13) irrespective of the mobilities σ v and σ w of the two species V and W .Then, starting from equal mobilities σ u = σ v = σ w , we gradually increase the mobility σ u of the bottom prey U and find the oscillatory Turing instability.The stationary steady state becomes unstable above a certain threshold and non-uniform oscillations develop.The emerging oscillatory Turing pattern is displayed in Fig. 2(a).One can notice that oscillations in the pattern are localized on a subset of nodes.Network nodes with relatively high degrees (hubs) exhibit oscillations while the other nodes remain staying near the uniform steady state.We will show later that such localized oscillations are typical for the oscillatory Turing patterns in networks.Remarkably, the stationary Turing instability can also be found in the same model when the mobility σ v of the intermediate species V is decreased (Fig. 2(b)).
Numerical simulations of the continuous model yield the behavior shown in Figs.2(c)(d), where periodic boundary conditions are employed.The stationary pattern (Fig. 2(b)) becomes transformed to a periodic stationary Turing pattern (Fig. 2(d)), whereas the localized oscillatory pattern (Fig. 2(a)) gives rise to a traveling wave (Fig. 2(c)).While patterns in networks and continuous media may appear different, they indeed correspond to the same, stationary or oscillatory, Turing bifurcations.
The oscillatory Turing instability is observed also in other food webs, Models B (Eqs.  4), yielding a uniform steady state (u 0 , v 0 , w 0 ) ≃ (0.295, 2.330, 3.975).As well as in Model A, the instability takes place in each system when the mobility of one species is increased to exceed a certain threshold.In Model B, the oscillatory Turing instability occurs when the mobility σ u of the bottom prey U is increased.The emerging pattern is shown in the left panel of Fig. 3(a).In Model C, the instability occurs as the mobility σ v of a prey V is increased up (Left panel in Fig. 3(b)).Thus, the oscillatory Turing instability is observed in all possible food webs with three species shown in Fig. 1(a)-(c).Furthermore, the oscillatory Turing instablity is not sensitive to the choice of nonlinear functions in Eqs.(1).As shown in Supporting Information, different nonlinear functions can be used to describe predator-prey interactions.We can observe the oscillatory Turing instability under various assumptions [3,49] about the functional form of the reproduction and death rates.Moreover, the instability can be observed even if one of the three species was immobile (see Supporting Information).
Therefore, we found the oscillatory Turing instability in a wide range of ecological models and conclude that the oscillatory Turing instability is generic for ecosystems.This discovery agrees with the statement of our analytical investigation (Sec.III B).Examining emergent patterns, one can notice that developing oscillations in all considered systems are localized on a subset of network nodes with close degrees.Although the localizing nodes are different depending on a system, localized oscillations were always found.As we discuss below, localization is a common

C. Subcritical vs. supercritical bifurcations
Above the instability threshold, nonlinear effects become important.They can lead to the saturation of growth and the establishment of a final pattern.Generally, they also determine whether a bifurcation is subcritical or supercritical.When the bifurcation is subcritical, the pattern with large magnitude becomes immediately established once the instability threshold is exceeded.Such bifurcations are characterized by hysteresis, so that the pattern persists even below the instability threshold.In contrast to this, a supercritical bifurcation does not show a hysteresis and the magnitudes of established patterns are small close to the threshold.In this case, the final patterns do not differ much from the first critical modes near the bifurcation point.
To investigate the bifurcation behavior of the oscillatory and stationary Turing instabilities shown in Fig. 2, amplitudes of developing patterns can be examined.The amplitude Ā of a final pattern on a network can be defined as where the network averages u(t) , v(t) and w(t) were defined in Eqs.(20).Note that this definition holds both for the oscillatory and stationary patterns.Figure 4 displays the results of the linear stability analysis and the global amplitude Ā for the oscillatory (a)(b) and stationary (c)(d) Turing instabilities in Model A. Increasing the mobility σ u of the prey U , we calculated the linear growth rate γ = λ + iω to find the instability at a threshold σ u,crit. .The critical eigenmode has a complex linear growth rate (Fig. 4(a)) and therefore the instability is oscillatory.After passing the instability, we gradually increase the mobility σ u further away from the threshold and calculate the global amplitude Ā (Fig. 4(b)).As can be clearly seen in the figure, the oscillatory Turing instability corresponds to a supercritical bifurcation.We have checked that near the threshold Ā ∝ (σ u − σ u,crit.) 1/2 holds.Thus, small amplitude patterns can be established near the instability threshold.In contrast, the stationary Turing instability, which corresponds to a growth rate of a real Due to its supercritical bifurcation, oscillatory Turing patterns near the instability threshold are well described only by the first critical eigenmode α c , so that we have Therefore, in this case, final patterns can be predicted by means of the critical Laplacian eigenvectors φ (αc) .
In Fig. 5, we demonstrate the prediction for Model A. Figure 5(a) shows the results of the linear stability analysis at two different values of overall mobility ǫ in Model A. The respective critical Laplacian eigenvectors are displayed in Figs.5(c)(d) and the actual oscillatory Turing patterns in Figs.5(e)(f).As seen in Fig. 5, critical Laplacian eigenvectors and developing oscillatory patterns are localized on subsets of network nodes with close degrees.The localization for the considered scale-free network of size N = 50 is demonstrated in Fig. 5(b), where magnitudes φ (α) i of the components of Laplacian eigenvectors φ (α) are displayed as a function of α.As we explained in Sec.III A, the mode indices {α} are sorted in the increasing order of the Laplacian eigenvalues {Λ (α) } so that Λ (1) holds.This localization is consistent with the previous discovery for large random networks [9,52].The critical eigenmode α c depend on the overall mobility ǫ shown in Fig. 5(a).As discussed previously [9], the Laplacian eigenvalue Λ (α) appears in the characteristic function (9) multiplied by ǫ.Then, if we change the overall diffusion mobility ǫ, the critical eigenmode α c shifts so that the product ǫΛ (αc) is kept constant.Thus, the characteristic localization of the critical Laplacian eigenvector changes depending on the overall mobility ǫ shown in Fig. 5(c)(d).Correspondingly, in emerging oscillatory Turing patterns, localizing nodes shift from hubs (Fig. 5(e)) to peripheral nodes (Fig. 5(f)).Thus, examining Fig. 5, one can predict the emerging patterns by means of respective critical Laplacian eigenvectors.The oscillatory Turing instabilities found in other models (Models B and C, and other models shown in Supporting Information) are also identified with supercritical bifurcations shown in Figs. 3, 7 and 8. Thus, the instability generically leads to pattern formations of small amplitudes, which are described by the Laplacian eigenvectors of critical modes.This indicates that the localization property of the oscillatory patterns originates in the localizing Laplacian eigenvectors.
In contrast to the oscillatory instability, the stationary Turing instability corresponds to a subcritical bifurcation.This gives a striking difference between patterns arising from both the instabilities.It has been previously found in two-component reaction-diffusion networks that the stationary Turing instability always corresponds to a subcritical bifurcation [9].This holds also in three-component systems, that is, ecological networks with three species.The bifurcation is subcritical and the stable nonlinear state is far from the uniform steady state.Although differentiation starts from the nodes with the characteristic degree of the critical eigenvector, it takes place in surrounding nodes sequentially and spreads over the entire network.The amplitude of the developed final pattern becomes large.Thus, as we have seen in Fig. 2, although both instabilities are induced by the diffusion effect, the uniform steady state is destabilized in distinctly different fashions.

E. Secondary instability -Extinction of ecological species
We have studied the dynamical behavior induced by the oscillatory Turing instability far from the threshold and found a potential secondary instability, leading to the extinction of ecological species.Note that the dynamical behavior far from the instability threshold largely depends on the considered model.We have identified such extinction in Models A and B. As an illustration, we here show the results for Model A. Figure 6 displays numerical results far from the first instability threshold.Now, all eigenmodes are unstable except for the zero eigenmode α = 0 (Hopf mode).In panel (a), the network average w(t) for the top predator W is plotted as a function of time.After an oscillatory transient, the average w(t) tends to zero and therefore the top predator vanishes.Two remaining species, the prey U and the intermediate predator V , exhibit uniform oscillations after the extinction of the top predator W (Fig. 6(b)).Thus, the secondary instability may be inherent in the oscillatory Turing instability, which destabilizes the considered ecosystem and leads to the extinction of ecological species.

V. DISCUSSION AND CONCLUSIONS
The mathematical description for the classical (stationary) Turing instability in networks has been proposed already in 1971 [34].However, it has been first applied only to regular lattices [34,35] and small networks [37,38].Recently, such instability in large random networks has been investigated and characteristic properties of stationary Turing patterns in large networks have been discussed [9,47,48].In contrast to this, the mathematical theory of the oscillatory Turing bifurcation (the analogue of the wave bifurcation) in networks has been missing and our work is the first report where such theory is constructed.
Similar to continuous reaction-diffusion systems, the oscillatory instability needs three reacting species [1] and it occurs when diffusion mobilities of the species are largely different.However, wave patterns, which are typical for such instability in continuous media, do not emerge in networks.Therefore, we prefer not to use the term "wave instability" in the present study.As we find, heterogeneous oscillations spontaneously develop and they are localized on a subset of network nodes with similar degrees.The localization of developing oscillations could be explained by taking into account the known statistical properties of Laplacian eigenvectors in large random networks [9,52].Such eigenvectors correspond to the critical modes of the oscillatory Turing instability and hence they determine the properties of developing patterns.
Previously, it has been shown that the classical (stationary) Turing bifurcation in networks is always subcritical; it is characterized by strong hysteresis and the properties of final stationary patterns in networks are largely different from those of the first critical modes.In contrast to this, we have found that the oscillatory Turing bifurcation in networks is typically supercritical.Thus, the hysteresis is absent and the amplitude of the Turing patterns gradually grows as the control parameter is increased.Near the bifurcation point, oscillatory Turing patterns with small amplitudes are observed and they agree well with the first critical modes.This means that by considering Laplacian eigenvectors the properties of final small-amplitude patterns can be analyzed.
While the proposed mathematical theory is general and applicable to systems with various origins, our detailed numerical simulations have been performed for the models which correspond to ecological networks also known as metapopulations [26].We have considered metapopulations with the graph structure of scale-free networks.The habitats, representing network nodes, were occupied by local three-species ecosystems forming food webs.All possible food webs with three predator or prey species were considered.We have only excluded the food web with two predators and one prey shown in Fig. 1 (d) because steady coexistence of all three species is not possible in this case.Moreover, simulations under various assumptions for nonlinear predator-prey interactions have been carried out.
The oscillatory Turing instability could be observed for all considered metapopulation models.It has been found as the network mobility (dispersal rate) of one of the species was gradually increased.Note that it could be the mobility of a prey or a predator, depending on the food web and the mathematical model applied.The instability has been observed also in the three-species metapopulations where one of the species was immobile.The results of our numerical simulations are supported by the general sufficient conditions for the oscillatory bifurcation in threecomponent ecological systems which we have derived.While our analysis has been performed only for systems with three species, it can be straightforwardly extended to the systems with a larger number of components and hence for the metapopulations with more complex food webs.
We conclude that the oscillatory Turing instability may be generic for ecosystems.It should be generally expected whenever metapopulations with at least three species and sufficiently large differences in the mobilities of the species are investigated.This result may be of principal importance.Previously, the discussion has been focused on the role of dispersal connections in enhancing the stability of uniform steady states [25,27].We find however that dispersal connections would often destabilize the uniform steady state and lead to the development of oscillations on a subset of network nodes.We would like to stress that such oscillations are a consequence of the differential dispersal mobilities of species and they were always absent for isolated populations in individual habitats.
While the developing oscillations are weak near the instability threshold, their amplitude increases with the distance from the critical point and large-amplitude oscillations are also possible.We have shown that the secondary instabilities of such oscillations may take place and, in the considered example Fig. 6, they have resulted in the extinction of one of the species and, therefore, in the degradation of an ecosystem.So far, only the global extinction through the development of uniform oscillations via a Hopf bifurcation has been discussed.Our work provides a different scenario for the extinction of species through the development of dispersal-induced hegerogeneous oscillations in ecological networks.
It would be interesting to perform experimental or field studies of the oscillatory Turing instability and the resulting Turing patterns in real ecological systems.For this purpose, it may be beneficial to work first with the artificially constructed metapopulations as it has been done before in the experimental studies on the role of dispersal connections [27,29,53].To prove the presence of oscillatory Turing patterns in natural ecosystems, the development of such patterns and their responses to local perturbations should be analyzed, similar to what has been done to demostrate the existence of classical stationary Turing patterns in biological organisms [6].

VI. APPENDICES A. Numerical simulations of other ecological models
Here we show the results for different ecological models.Different dependence of reproduction and death rates on the densities of the species in predator-prey models are possible [3].In the main text, Holling type II dependences have been used as a typical example.Below we give models with other dependences including Holling type III [3,49] and the functions introduced by Murray [3].As the food web architecture, the food chain shown in Fig. 1(a) is employed.Model D. The food chain (Fig. 1(a)) with the Holling type III dependence for both prey and predator: Results are shown in Fig. 7(a).Parameters in Eqs. ( 25) are fixed at a 1, a w = 0.75, d w = 1 and µ = ν = 0.125.A uniform steady state is found at (u 0 , v 0 , w 0 ) ≃ (0.509, 0.612, 0.497).The oscillatory Turing instability is observed as σ u is increased.
Model E. The food chain (Fig. 1(a)) with the Holling type II dependence for prey and a linear function for predator: Results are shown in Fig. 7(b).Parameters in Eqs. ( 26) are fixed at 1, a w = 4, d w = 0.25 and µ = ν = 0.25.A uniform steady state is found at (u 0 , v 0 , w 0 ) ≃ (1.084, 2.557, 10.23).The oscillatory Turing instability takes place in this system as σ w is increased.Model F. The food chain (Fig. 1(a)) with the Holling type III dependence for prey and a linear function for predator: A uniform steady state is (u 0 , v 0 , w 0 ) ≃ (1.477, 1.672, 4.179).The oscillatory Turing instability is found as σ w is increased.Note that the constants C v,w enter both into the reproduction and death rates of species V and W in Eqs. ( 26) and (27).Therefore, they become cancelled in the expressions for G and H. Their numerical values are not relevant and are not provided.Model A with one immobile species.As implied by the sufficient conditions ( 12)-( 14), the oscillatory Turing instability is possible even if one of the species is immobile.Figure 8 shows results for the case that the intermediate species V is immobile, i.e. σ v = 0 (Fig. 8(a)) and that the top predator W is immobile, i.e. σ w = 0 (Fig. 8(b)).Parameters in Eqs.(2) are the same as those used in the main text.The oscillatory Turing instability is observed as σ u is increased in both cases.
Thus, the oscillatory Turing instability could be observed for different food web architectures (Fig. 4 and 3) and for different nonlinearities in the reproduction and death rates (Fig. 7).Moreover, the oscillatory Turing patterns were observed even when one of the three species was immobile (Fig. 8).In all considered systems, oscillatory Turing bifurcations were supercritical.

B. Numerical simulations of chemical reaction-diffusion networks
In chemical systems, autocatalytic reactions and diffusion may induce the Turing instability.The chemical Brusselator [54] and the Oregonator [55] are typical mathematical models to describe autocatalytic reactions and their numerical investigations are extensively used for exploring pattern formations in chemical systems.The Brusselator was originaly proposed to examine chemical oscillatory dynamics.It is a hypothetical model, which gives a rich variety of chemical dissipative structures.In contrast to the Brusselator, the Oregonator has been proposed to describe the real Belousov-Zhabotinsky reaction.
Both chemical models have been previously extended to explore the oscillatory Turing instability (wave instability) in continuous media [11][12][13].It was suggested that a chemical reaction in which the activator is reversibly transformed into an unreactive chemical species can explain wave patterns observed in the BZ aerosol OT system [12].Therefore, an additional third component taking such reversible transformation was added to the Brusselator.For the Oregonator, which originally contains three components, an additional reversal reaction was added.Thus, extended Brusselator and Oregonator were constructed [13].
Here, we consider the behaviors of such extended models in network-organized systems.The extended network Brusselator is described by equations We fix parameters as a = 1, b = 2.9, c = 1 and d = 1, yielding a steady state (u 0 , v 0 , w 0 ) = (1, 2.9, 1).The extended

FIG. 6 :
FIG. 6: Secondary instability in Model A. (a) Time-dependence of the network average w(t) far from the instability threshold.(b) Local densities of three species after a transient.The same dynamics are observed in all network nodes (Uniform oscillations).Dispersal mobilities are ǫ = 0.1, σu = 6.75, σv = σw = 0.01.