Noise and Correlations in a Spatial Population Model with Cyclic Competition
Noise and spatial degrees of freedom characterize most ecosystems. Some aspects of their influence on the coevolution of populations with cyclic interspecies competition have been demonstrated in recent experiments [e.g. B. Kerr et al., Nature 418, 171 (2002)]. To reach a better theoretical understanding of these phenomena, we consider a paradigmatic spatial model where three species exhibit cyclic dominance. Using an individual-based description, as well as stochastic partial differential and deterministic reaction-diffusion equations, we account for stochastic fluctuations and spatial diffusion at different levels, and show how fascinating patterns of entangled spirals emerge. We rationalize our analysis by computing the spatio-temporal correlation functions and provide analytical expressions for the front velocity and the wavelength of the propagating spiral waves.
Understanding the combined influence of spatial degrees of freedom and noise on biodiversity is an important issue in theoretical biology and ecology. This implies to face the challenging problem of studying complex nonequilibrium structures, which form in the course of nonlinear evolution Turing (1952); May (1974); Murray (2002); Kerr et al. (2002); Tainaka (1989); Reichenbach et al. (2007a). More generally, self-organized nonequilibrium patterns and traveling waves are ubiquitous in nature and appear, for instance, in chemical reactions, biological systems, as well as in epidemic outbreaks Kap (1995). Among the most studied types of patterns are spiral waves, which are relevant to autocatalytic chemical reactions, aggregating slime-mold cells and cardiac muscle tissue Zaikin and Zhabotinskii (1970). In all these nonequilibrium and nonlinear processes, as well as in population dynamics models Turing (1952); Murray (2002); Tainaka (1989), pattern formation is driven by diffusion which, together with internal noise, act as mechanisms allowing for stabilization and coevolution of the reactants. In this work, we consider a paradigmatic spatially-extended species population system with cyclic competition, which can be regarded as a simple food-chain model Hofbauer and Sigmund (1998). In fact, such a system is inspired by recent experiments on the coevolution of species of bacteria in cyclic competition Kerr et al. (2002). Using methods of statistical physics, we study the influence of spatial degrees of freedom and internal noise on the coevolution of the species and on the emerging spiral patterns. In particular, we compute the correlation functions and provide analytical expressions for the spreading speed and wavelength of the propagating fronts. To underpin the role of internal noise, the results of the stochastic description are compared with those of the deterministic equations.
In this Letter, we investigate a stochastic spatial variant of the rock-paper-scissors game Hofbauer and Sigmund (1998) (also referred to as cyclic Lotka-Volterra model). These kinds of systems have been studied both from a game-theoretic perspective, see e.g. Hauert and Szabó (2005); Szabó and Fath (2007) and references therein, and within the framework of chemical reactions Tainaka (1989); Frachebourg et al. (1996), revealing rich spatio-temporal behaviors (e.g. emergence of rotating spirals). While our methods have a broad range of applicability, they are illustrated for a prototypical model introduced by May and Leonard May and Leonard (1975) where species, and undergo a cyclic competition (codominance with rate ) and reproduction (with rate ), according to the reactions
Hence, an individual of species will consume one of species () with rate and will reproduce with rate if an empty spot, denoted , is available (, i.e. there is a finite carrying capacity). In addition, to mimic the possibility of migration, it is realistic to endow the individuals with a form of mobility. For the sake of simplicity, we consider a simple exchange process, with rate , among any nearest-neighbor pairs of agents: . If one ignores the spatial structure and assumes the system to be well-mixed (with an infinite number of individuals), the population’s mobility plays no role and the dynamics is aptly described by the deterministic rate equations (RE) for the densities of species and , respectively. Introducing , the RE read:
where the index is taken modulo 3 and is the total density. As shown by May and Leonard May and Leonard (1975) (see also Durrett and Levin (1998)), these equations possess 4 absorbing fixed points, corresponding to a system filled with only one species and to an empty system. In addition, there is a reactive fixed point , corresponding to a total density . A linear stability analysis shows that is unstable. The absorbing steady states , and are heteroclinic points. The existence of a Lyapunov function allows to prove that, within the realm of the above RE, the phase portrait is characterized by flows spiraling outward from , with frequency in its vicinity. Approaching the boundaries of the phase portrait, the trajectories form (heteroclinic) cycles indefinitely close to the edges (without ever reaching them), with densities approaching in turn the value one. Despite its mathematical elegance, this behavior has been recognized to be unrealistic May and Leonard (1975); Durrett and Levin (1998). In fact, for finite populations, fluctuations arise and always cause the extinction of two species in finite time (see e.g. Ref Reichenbach et al. (2006)).
In this work, considering the spatial version of the above model in the presence of internal noise, we show that a robust (and, arguably, more realistic) scenario for the evolution arises. The reaction schemes (2) and the exchange events are considered to occur on a dimensional regular lattice of sites, labeled . Each lattice site has neighbors at a distance (e.g. and for hypercubic lattices of linear size ) and is either empty or occupied by at most one individual. On the lattice, the binary reactions (2) and exchanges only occur among pairs of nearest-neighbors. In the situation of large system sizes, the continuum limit reveals that for the exchange process to be an efficient driving mechanism, the rate has to scale as , with and . In fact, if the system is dominated by the local reactions (2) among neighboring individuals; while effective diffusion renders locality irrelevant when . Only when , there is an effective competition between the stirring process and the local reactions (2). It is therefore useful to introduce the effective diffusion constant . Because of the discreteness of the number of individuals involved in the reactions, internal fluctuations arise in the system. The latter originate from (i) the interspecies reactions (2) and (ii) the exchange processes. In the continuum limit, where with (and finite ), there is a separation of time scales and the pair exchanges occur much faster than the reactions (). Actually, the fluctuations associated with (2) and the agents’ mobility scale respectively as and , with the former dominating over the latter and being the only relevant contribution. This result is revealed by a system size, also called Kramers-Moyal (see e.g. Ref. Gardiner (1983), Chap. 8), expansion (SZE) of the master equation underlying the exchange processes and the reactions (2) Reichenbach et al. (2007d). Furthermore, the SZE yields a proper Fokker-Planck equation, which is equivalent to a set of (Ito) stochastic partial differential equations (SPDE) with white noise. The derivation, obtained in the continuum limit from the master equation, is outlined in the supplementary EPAPS document Reichenbach et al. (2007d) and will be detailed elsewhere Reichenbach et al. (2007b). Here, we quote the expression of the SPDE:
where is the Laplacian operator; , and
Again, the indices are taken modulo and now . As explained in Gardiner (1983); Reichenbach et al. (2007b), these SPDE have to be interpreted in the sense of Ito calculus. While Eqs. (4) and our approach are valid in any dimension Reichenbach et al. (2007d), for specificity, we now analyze the spatio-temporal properties of the system in two dimensions with periodic boundary conditions. On the one hand we have solved numerically the SPDE (4) using the open software from the XMDS project xmd (). On the other hand, we have carried out individual-based simulations of the reactions (2) for mobile (exchange process) particles on lattices of size , with . This allows to carefully study the convergence towards the continuum limit, where the description in terms of (4) is expected to be accurate.
As other spatially-extended dynamical systems Tainaka (1989); Hauert and Szabó (2005); Szabó and Fath (2007); Frachebourg et al. (1996), the model under consideration displays fascinating nonequilibrium patterns emerging in the course of the evolution. In Fig. 1 (a) and (b), we report typical long-time snapshots of the system for low (a) and high (b) exchange rates (but keeping fixed), as obtained from lattice simulations. In both cases we notice that intriguing patterns form. For slow exchange rate, the system displays non-geometrical patches, similarly to what happens in systems with self-organized criticality Schenk et al. (2000). When the exchange rate is raised, the patterns display spiral structures. In fact, starting from a spatially homogeneous initial condition, , the system is randomly perturbed by the internal noise and the resulting spatial inhomogeneities grow and form wavefronts moving through the system. The emergence of spiral patterns is a feature shared by other excitable systems (see e.g. Kap (1995); Zaikin and Zhabotinskii (1970)) and corresponds to the ability of the system to sustain the propagation of oscillating waves. For sufficiently large , one observes a striking resemblance between the size and structure of the patterns obtained from the lattice simulations [Fig. 1 (b)] and those from the SPDE (4) [Fig. 1 (c)]. To further compare the predictions of the SPDE (4) with the lattice simulations, and to gain additional information on the structure of the the emerging patterns, we have computed the correlation functions, in two dimensions. In Fig. 2 (red and blue curves), we report the results for in the steady state and notice an excellent agreement between the results of the lattice simulations and the predictions of the SPDE (4). The inset of Fig. 2, displays the correlation length 111The correlation length denotes the characteristic length at which the spatial correlations decay by a factor from their maximal value. as function of ( is kept fixed, varies) obtained in the lattice simulations, which is found to coincide with the prediction of the SPDE already for . We have also computed the autocorrelation function and found, both in the lattice simulations and from the solutions of the SPDE, an oscillating behavior with a similar characteristic frequency, markedly different from Reichenbach et al. (2007b). This confirms that, even for finite exchange rates, the solution of the SPDE (4) provides an excellent approximation of the lattice simulations of the system. This is rather surprising since Eqs. (4) have been derived in the continuum limit, where and . A comparable influence of finite exchange rate in a predator-prey system has been reported recently Mobilia et al. (2006).
According to the SPDE (4), scales as , so that by raising the diffusion one increases the size of the spirals. As we have shown in Ref. Reichenbach et al. (2007a), this happens up to a critical value (e.g. for ): above that threshold, the spiral structures outgrow the system size and only one species survives, corresponding to an absorbing steady state predicted by Eqs. (3).
As the properties of the lattice simulations are well captured by the SPDE (4), where the strength of the noise scales as , with , it is natural to investigate the actual influence of this internal noise on the steady state of the system. To address this issue, we have solved numerically (in 2D, with periodic boundary conditions) the deterministic reaction-diffusion equation (DRDE) obtained from (4) by dropping the noise terms, i.e. Of course, to obtain a nontrivial steady state for the DRDE one has to assume spatially inhomogeneous initial conditions. In Fig. 1 (d), we have reported a snapshot of the long-time behavior predicted by the DRDE starting from . In this case, the dynamics evolves towards a reactive steady state which also exhibits spiral waves. However, the latter do not form entangled structures, but ordered geometrical patterns. As an example, only four spirals cover the system in Fig. 1 (d) [while noise leads to 106 entangled spirals in Fig. 1 (c)]. The correlation functions associated with the DRDE therefore exhibit only weakly damped spatial oscillations (see Fig. 2, green triangles). By analyzing typical snapshots like those of Fig. 1 (d), we have noted that in the deterministic and stochastic [i.e. lattice simulations with “large” and solutions of Eqs. (4)] descriptions, the spiral waves share the same propagation velocity, frequency and wavelength. However, a major difference between these descriptions lies on the crucial dependence of the DRDE on initial conditions, which determine the overall number of spirals and their size. On the contrary, because the internal noise acts a random source of spatial inhomogeneities, the lattice stochastic system and the SPDE display robust features. In particular, we have found noise to induce a universal spiral density of about per square wavelength.
Analytical expressions for the spreading velocity and the wavelength of the propagating fronts of the DRDE can be obtained by considering the dynamics on the invariant manifold of the RE Wiggins (1990), given by , with . On , up to order, the DRDE can be recast in the form of a forced complex Ginzburg-Landau equation (CGLE) Saarloos (2003); Cross and Hohenberg (1993). By performing the nonlinear transformation and , upon ignoring nonlinearities like , one is left with the following CGLE in the variable Reichenbach et al. (2007b):
with , and . The general theory of front propagation Saarloos (2003); Cross and Hohenberg (1993) predicts that Eq. (7) always admits traveling waves as stable solutions (i.e. no Benjamin-Feir or Eckhaus instabilities occur). We have determined such periodic solutions by computing, from the dispersion relation of (7), the spreading velocity and the spirals’ wavelength (details will be given in Reichenbach et al. (2007b)):
In the stochastic version of the model, the wavelength and velocity of the wavefronts have been found to agree with those of the deterministic treatment. Hence, the expressions (8) also apply (for large , with ) to the results of lattice simulations (rescaled by a factor ) and to the solution of the SPDE (4). For instance, on a square grid with , lattice simulations and Eqs. (4) yield , in good agreement with the prediction of (8): . For the spirals’ wavelength, numerical results (lattice simulations and SPDE) yield as predicted by (8). In Fig. 3, the analytical prediction (8) for is compared with the values obtained from the SPDE (4), yielding a remarkable agreement for the functional dependence on the parameter . Yet, as Eq. (7) does not account for all nonlinearities, the analytical and numerical values differ by a prefactor (considered in Fig. 3) Reichenbach et al. (2007b). It can still be noted that (7) and the predictions (8) are valid in all dimensions Reichenbach et al. (2007d, b).
Motivated by recent experiments Kerr et al. (2002), we have considered a spatially-extended model with three species in cyclic competition, and focused on the spatial and stochastic effects. The local character of the reactions and internal noise allow mobile populations to coexist and lead to pattern formation. We have shown that already for finite mobility the lattice model can be described by SPDE. With the latter and lattice simulations, we have studied how entanglement of spirals form and obtained expressions for their spreading velocity and wavelength. The size of the patterns crucially depends on the diffusivity: above a certain threshold the system is covered by one species Reichenbach et al. (2007a). In the absence of noise, the equations still predict the formation of spiral waves, but their spatial arrangement depends on the initial conditions.
Support of the German Excellence Initiative via the program “Nanosystems
Initiative Munich” is gratefully acknowledged.
M. M. is grateful to the Humboldt Foundation
for support through the grant IV-SCZ/1119205.
*Currently at: Mathematics Institute & Centre for Complexity Science, University of Warwick, Coventry CV4 7AL, U.K.
- Turing (1952) A. M. Turing, Phil. Trans. R. Soc. B 237, 37 (1952).
- May (1974) R. M. May, Stability and Complexity in Model Ecosystems (Cambridge University Press, 1974); S. J. Maynard, Models in Ecology (Cambridge University Press, 1974).
- Murray (2002) J. D. Murray, Mathematical Biology (Springer, 2002).
- Kerr et al. (2002) B. Kerr et al., Nature 418, 171 (2002).
- Tainaka (1989) K. Tainaka, Phys. Rev. Lett. 63, 2688 (1989); Phys. Rev. E 50, 3401 (1994).
- Reichenbach et al. (2007a) T. Reichenbach, M. Mobilia, and E. Frey, Nature 448, 1046 (2007a).
- Kap (1995) in Chemical Waves and Patterns, Part One, edited by R. Kapral and K. Showalter (Kluwer Academic Publishers, Dordrecht, 1995); J. Lechleiter et al., Science 252, 123 (1991); B. T. Grenfell, O. N. Bjornstad, and J. Kappey, Nature 414, 716 (2001).
- Zaikin and Zhabotinskii (1970) A. N. Zaikin and A. M. Zhabotinskii, Nature 225, 535 (1970)
- Hofbauer and Sigmund (1998) J. Hofbauer and K. Sigmund, Evolutionary Games and Population Dynamics (Cambridge University Press, 1998).
- Hauert and Szabó (2005) C. Hauert and G. Szabó, Am. J. Phys. 73, 405 (2005).
- Szabó and Fath (2007) G. Szabó and G. Fáth, Phys. Rep. 446, 97 (2007).
- Frachebourg et al. (1996) L. Frachebourg, P. L. Krapivsky, and E. Ben-Naim, Phys. Rev. Lett. 77, 2125 (1996).
- May and Leonard (1975) R. M. May and W. J. Leonard, SIAM J. Appl. Math 29, 243 (1975).
- Durrett and Levin (1998) R. Durrett and S. Levin, Theor. Pop. Biol. 53, 30 (1998).
- Reichenbach et al. (2006) T. Reichenbach, M. Mobilia, and E. Frey, Phys. Rev. E 74, 051907 (2006).
- Gardiner (1983) C. W. Gardiner, Handbook of Stochastic Methods (Springer, 1983), 1st ed.
- Reichenbach et al. (2007d) See the Supplementary Material (EPAPS Doc.), pages 5-6, for the derivation of Eqs. (4) from the master equation and a discussion of the general dimensional situation.
- Reichenbach et al. (2007b) T. Reichenbach, M. Mobilia, and E. Frey, in preparation.
- (20) URL http://www.xmds.org.
- Schenk et al. (2000) K. Schenk et al., Eur. Phys. J. B 15, 177 (2000).
- Mobilia et al. (2006) M. Mobilia, I. T. Georgiev, and U. C. Täuber, Phys. Rev. E 73, 040903(R) (2006).
- Wiggins (1990) S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos (Springer, 1990), 1st ed.
- Saarloos (2003) W. Saarloos, Phys. Rep. 386, 29 (2003).
- Cross and Hohenberg (1993) M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. 65, 851 (1993).
Supplementary Material (EPAPS Document)
Appendix A From the master equation to stochastic partial differential equations
In this Supplementary Notes, we want to explicitly outline how, starting from the master equation associated with the stochastic May-Leonard model [defined by the reactions (1)], the set of stochastic partial differential equations (3) can be obtained via system size expansion.
For the sake of illustration, here we focus on the role of internal noise stemming from the reactions (1), and ignore spatial degrees of freedom. As detailed in a forthcoming publication , and following the reasoning presented in ,
in a proper continuum limit, the spatial dispersal of individuals is accounted by diffusive terms in the SPDE (3).
As in the main text, the overall number of individuals is denoted and stands for the frequencies (or densities) of the species , , and in the population (i.e. and ). The master equation giving the time-evolution of the probability of finding the system in the state at time then reads
where denotes the transition probability from state to the state within one time step (loss term), is the analogous gain term, and the summation extends over all possible changes . As an example, the relevant changes in the density resulting from the basic reactions (1) are in the third reaction, in the fourth, and zero in all others. We also choose the unit of time such that, on average, every individual reacts once per time step. The transition rates resulting from the reactions (1) then read for the reaction (the prefactor of enters due to our choice of time scale, where reactions occur in one unit of time) and for . Transition probabilities associated with all other reactions (1) follow similarly.
The system size, or Kramers-Moyal, expansion (SZE)  of the Master equation is an expansion in the increment , which is proportional to . Therefore, the SZE may be understood as an expansion in the inverse system size . To second order in , it yields the (generic) Fokker-Planck equation :
For the system under consideration, in the above the indices and the summation convention in (10) implies sums carried over them. According to the Kramers-Moyal expansion (or SZE), the quantities and are given by 
Note that is symmetric. For the sake of clarity, we outline the calculation of : The relevant changes result from the third and fourth reactions in (1), as described above. The corresponding rates respectively read and , resulting in . The other quantities are computed analogously. All explicit expressions for and will be derived and given in detail in . The well-known correspondence between Fokker-Planck equations and Ito calculus  implies that (10) is equivalent to the following set of Ito stochastic differential equations (with the above summation convention):
where the matrix is defined from via the relation , and the ’s denote (uncorrelated) Gaussian white noise terms. Note that for the model under consideration is diagonal and one can therefore always choose a diagonal matrix for [see Eq. (5)], with only ’s contributing to the right-hand side of Eqs. (11).
In Ref., we demonstrate that (in a proper continuum limit) spatial degrees of freedom and exchange processes simply yield additional diffusive terms in (15). This leads to the SPDE (3), given and discussed in the main text, in which the s still denote Gaussian white noise contributions.
Appendix B Generalization in arbitrary dimensions
The SPDE (3) and the CGLE (6), as well as the analytical predictions (7), are valid in all spatial dimensions. While in the text, for the sake of specificity, we mainly focus on the two-dimensional situation, here we comment on some features of the one-dimensional version of the system, as well as on some properties of the general situation in higher dimensions.
In one dimension and in the absence of exchange processes (mixing), our model is expected to exhibit coarsening like the cyclic Lotka-Volterra model (see e.g. ). The same phenomenon still occurs in the presence of very slow mixing rate (for the model under consideration, this means very low values of the exchange rate ). On the other hand, and as shown in the main text, in the presence of (finite) mixing the systems’ behavior is aptly described by the SPDE (3). The underlying CGLE (6) predicts the propagation of traveling waves, with velocity and wavelength still given by the analytical expression (7). Similarly to what has been found in the two-dimensional system , if the exchange rate is “moderate”, i.e. below a certain mobility threshold but still finite [6,15], the system confirms these predictions and the species coexist. However, due to stochastic events and the presence of absorbing boundaries, some domains will occasionally merge, resulting in growing domains. This coarsening phenomenon will happen on a much longer time-scale than in the absence of the mixing processes. Above the threshold value for the diffusivity, the system is well-mixed, the underlying spatial structure plays no role, and the description in terms of the rate equations (2) is valid. In this mean-field scenario, the system approaches an absorbing steady state and no patterns emerge.
In higher dimensions (i.e. in dimensions ), the description of the stochastic spatial system in terms of the SPDE (3) and of the CGLE (7) is again valid for “moderate” (or intermediate) mixing (i.e. for a finite mobility rate which is below a certain critical threshold, see ). In this situation, and as discussed in the main text, the description in terms of the CGLE (6) is qualitatively valid and predicts the emergence of moving spiral waves in two dimensions (which is the case discussed in detail in the main text, see also [6,15]), and to “scroll waves”, i.e. vortex filaments, in three dimensions .