Extinction in neutrally stable stochastic Lotka-Volterra models

# Extinction in neutrally stable stochastic Lotka-Volterra models

Alexander Dobrinevski Arnold Sommerfeld Center for Theoretical Physics and Center for NanoScience, Department of Physics, Ludwig-Maximilians-Universität München, Theresienstrasse 37, D-80333 München, Germany CNRS-Laboratoire de Physique Théorique de l’Ecole Normale Supérieure, 24 rue Lhomond, 75005 Paris Cedex-France thanks: LPTENS is a Unité Propre du C.N.R.S. associée à l’Ecole Normale Supérieure et à l’Université Paris Sud    Erwin Frey Arnold Sommerfeld Center for Theoretical Physics and Center for NanoScience, Department of Physics, Ludwig-Maximilians-Universität München, Theresienstrasse 37, D-80333 München, Germany
###### Abstract

Populations of competing biological species exhibit a fascinating interplay between the nonlinear dynamics of evolutionary selection forces and random fluctuations arising from the stochastic nature of the interactions. The processes leading to extinction of species, whose understanding is a key component in the study of evolution and biodiversity, are influenced by both of these factors. Here, we investigate a class of stochastic population dynamics models based on generalized Lotka-Volterra systems. In the case of neutral stability of the underlying deterministic model, the impact of intrinsic noise on the survival of species is dramatic: it destroys coexistence of interacting species on a time scale proportional to the population size. We introduce a new method based on stochastic averaging which allows one to understand this extinction process quantitatively by reduction to a lower-dimensional effective dynamics. This is performed analytically for two highly symmetrical models and can be generalized numerically to more complex situations. The extinction probability distributions and other quantities of interest we obtain show excellent agreement with simulations.

###### pacs:
05.40.-a, 05.10.Gg, 02.50.Ey, 02.50.Fz, 87.23.Cc

## I Introduction

Interactions between biological species are known to lead to very diverse and intricate behaviour of a population. This includes, just to give a few examples, coexistence of a surprisingly high number of competing species in the same ecological niche Hutchinson1961 (), oscillating population cycles Gilg2003 () and chaos Turchin2000 (). The question which (if any) of the species in a web of interactions survive, and for how long, is thus very nontrivial but central for the understanding of evolution and biodiversity Frey2010 ().

A classical and long-established model for the interaction of species in a well-mixed habitat is the Lotka-Volterra model Lotka1920 (); Volterra1931 (). Since its introduction, it has also been successfully applied in many different contexts outside of population dynamics: among others, neural networks Rabinovich2008 (), game theory Hofbauer1998 () and physiology McCarley1975 (). This model attempts to describe the interaction between species through a set of coupled ordinary differential equations of the form

 ∂txi(t)=xi(t)(bi+S∑j=1Aijxj(t)). (1)

The abundance of each species is given by a continuous, real-valued variable with . are constant source terms describing the growth (or decline) of each species in the absence of the others, and is a constant matrix modelling the interactions between the species. Within this model, survival or extinction of species is purely deterministic: Any fixed initial condition determines unambiguously which, if any, of the species survive. Technically, the main underpinning for this is the stability or instability of certain stationary solutions of the differential equations (1). Some rather precise criteria for determining the persistence of species directly from the vector , the matrix and the initial conditions have been obtained in literature Goh1977 (); Zeeman1995 ().

However, in a real biological situation, the population is made up from a large but still finite number of individuals. Hence, the abundances of each species can only change in discrete steps and not continuously. Furthermore, the interactions between them, as well as birth and death processes, have – to some extent – a stochastic nature. All these features cannot be modelled by the deterministic equations (1). In fact, such effects of finite system size and fluctuations due to some intrinsic randomness (or, likewise, due to external noise) have recently been recognized to be very important for extinction processes, especially in the case when the deterministic solutions exhibit neutral stability (see, for example, Reichenbach2006 (); Traulsen2006 (); Traulsen2005 (); Cremer2009 (); Nowak (); Antal2006 (); Nowak2004a (); Taylor2004 (); Fogel1998 (); Ficici2000 () and many others).

In this paper, we propose a new method based on the idea of stochastic averaging which allows one to gain a quantitative understanding of the stochastic extinction process in the case when the deterministic limit of the model is neutrally stable. The idea of stochastic averaging was first introduced by Khasminskii in Ref. Khasminskii1968 (). Later on, it was rigorously justified in Ref. Freidlin2004 () for two-dimensional systems possessing a conservation law. So far, however, it has not gained a lot of popularity in physical literature.

In section II we will formulate a stochastic model of population dynamics based on a graph of interactions between species, whose dynamics in the absence of noise reduces to a Lotka-Volterra model of Eq. (1).

In sections III and IV we will treat two pedagogical examples of such models, the three-species and four-species systems with cyclic dominance. Their deterministic dynamics will be shown to be neutrally stable, i.e. to lead to perpetual coexistence of all species with periodically oscillating abundances. However, we will see that taking into account fluctuations due to the stochastic nature of the interactions introduces a finite mean extinction time proportional to the population size. Using stochastic averaging, we will characterize the extinction process by an effective stochastic process in the deterministically conserved quantities. This will allow us to obtain quantitative results on extinction times and their dependence on the initial conditions.

The generalization to more complex models will be discussed in section V.

## Ii Stochastic Lotka-Volterra Models

Let us consider a well-mixed population with species, and interactions between them defined by a graph with vertices and (directed) edges

 E={(Xi,Xj)|i≠j;i,j=1,...,S}.

An edge from to is denoted by an arrow in figure 1, and is taken to indicate that the species dominates over . We allow at most one edge between each pair of species. From this graph, we can now write down a set of reaction equations implementing the interactions of the species. For every edge we formulate an interaction between and in the formalism of chemical reaction equations:

 Xi+Xj\lx@stackrelkij⟶Xi+Xi, (2)

where the reaction rate is the probability for this reaction to occur per (infinitesimal) time unit and per possible pair of individuals.

Note that the model described by the reactions in Eq. (2) provides an individual-based, discrete description of the interactions, and includes stochastic fluctuations since the reaction rates are interpreted probabilitistically.

Since the reactions in Eq. (2) keep the total number of individuals fixed, we can assume a constant system size . The system state is then described by the -tuple of species counts , with and . We define the “basis vectors” by where the is on the -th position. The reactions in Eq. (2) can then be translated into a master equation giving the evolution of the occupation probabilities for each state :

 ∂tP→n(t) =∑(Xi,Xj)∈Ekij [(ni−1)(nj+1)P→n−→ei+→ej(t) (3) −ninjP→n(t)].

For biological applications, one is mostly interested in large populations, i.e. in the limit of large . The relative abundances of each species, , can then be assumed to be real-valued variables in the interval . Using a standard Kramers-Moyal expansion VanKampen (), the master equation (3) can then be approximated by a Fokker-Planck equation in the intensive variables :

 ∂tP({xk},t) = −S∑i=1∂i[αiP({xk},t)] (4) +12NS∑i,j=1∂i∂j[BijP({xk},t)].

The conservation of the total population size gives rise to the normalization condition . The drift and noise terms in Eq. (4) are given by:

 αi = xiS∑j=1Aijxj, (5) Bij = {∑Sk=1|Aik|xixkfori=j−∣∣Aij∣∣xixjfori≠j. (6)

The entries of the interaction matrix are:

 Aij=⎧⎪⎨⎪⎩kijif(Xi,Xj)∈E−kijif(Xj,Xi)∈E0otherwise. (7)

As is well known Gardiner (), the Fokker-Planck equation (4) can be reformulated as a set of stochastic differential equations, or Langevin equations:

 dxi=αidt+1√NS∑j=1CijdWj. (8)

Throughout this paper, we take all stochastic differential equations to be in the Itô interpretation. is a matrix satisfying , with defined by Eq. (6). Certainly, this does not fix uniquely, but the precise choice has no influence on the stochastic process Gardiner (). are independent Wiener processes, or Brownian motions, with zero mean and unit variance.

From equation (8) we see that the deterministic, noiseless limit of the general model of Eq. (2) is given by the following set of coupled ordinary differential equations (so-called rate equations):

 ∂txi=xiS∑j=1Aijxj. (9)

This can now be immediately identified as a generalized Lotka-Volterra model of the form of Eq.(1). The interaction matrix is given by Eq. (7), and the source terms vanish, i.e. , simply since all reactions in Eq. (2) have exactly two reactants. A model with non-vanishing source terms can be built by considering reactions of the form of death processes, , and birth (branching) processes, , in addition to Eq. (2).

The noise terms in Eq. (8), proportional to , encapsulate the fluctuations due to the discreteness of the individuals and the stochastic nature of the reactions in Eq. (2).

In total, we have given a procedure allowing us to obtain a stochastic model (in terms of Fokker-Planck or Langevin equations) of a system with a large population of individuals from a general interaction graph of the species. This prescription is certainly not unique, however, it has the nice property that the deterministic dynamics of the resulting model is in one-to-one correspondence with a Lotka-Volterra model whose interaction matrix is the adjacency matrix of the original graph.

Considering the widespread use and the importance of Lotka-Volterra models, it seems worthwhile to study models of the form of Eq. (8) and the effects of stochasticity in them.

In general, the deterministic rate equations (9) possess extinction fixed points (where some species are extinct, i.e. there are some with ) and coexistence fixed points (where all species are present, for all ). From Eq. (9) we see immediately that the coexistence fixed points form a linear subspace of the phase space of the system which is given by the kernel of the matrix . When the full stochastic model in Eq. (8) is considered, fluctuations cause the system to touch an extinction hyperplane where for some sooner or later. Since a species which has died out cannot be re-introduced (this is apparent from the reaction equations (2)), this means that in the stochastic system, extinction always occurs eventually.

However, the time scale on which this process occurs can vary greatly. A classification of the possible scenarios, characterized by the scaling of the mean extinction time with the population size was discussed in Ref. Antal2006 () and further developed in Ref. Reichenbach2007 (); Cremer2009 ():

1. Stable coexistence for , occurring when the deterministic dynamics has a stable attractor in the coexistence region. Here, extinction is driven by rare large deviations and hence the extinction times for large populations are extremely long.

2. Unstable coexistence for . This occurs when the flow of the deterministic dynamics approaches one of the extinction hyperplanes for large times, and weak fluctuations are already sufficient to make one of the species go extinct.

3. Neutrally stable coexistence for a power-law dependence, . This occurs when the deterministic dynamics possesses a family of neutrally stable, closed orbits, corresponding to the existence of a conservation law.

For simple models, these criteria correspond to (linearly) stable, unstable, or neutrally stable coexistence fixed points in the rate equations (2).

While here we are only considering well-mixed populations, a similar classification is possible for models which include spatial degrees of freedom Reichenbach2007 (). For a review on population dynamics in spatially extended systems see e.g. Szabo2007 (); Frey2010 (). There is yet another interesting connection to extinction times close to absorbing-state phase transitions; see e.g. Refs. Odor2004 (); Henkel2009 (); Kuhr2011 ().

Observe that the effects of stochasticity are most dramatic in a neutrally stable model: while the deterministic dynamics predicts perpetual coexistence far away from the extinction planes, inclusion of fluctuations introduces a finite mean extinction time which only scales as a power law with the population size. In the following two sections, we will now analyze the stochastic extinction process for two pedagogical example models of this kind.

## Iii Cyclic Three-Species Model: The Rock-Paper-Scissors Game

### iii.1 Introduction

The first example that shall be considered in detail is a three-species model with cyclic dominance, whose interaction graph is shown in figure 2. One of the most popular areas where such cyclic, intransitive relationships between three entities arise is in game theory as a so-called rock-paper-scissors game Hofbauer1998 (). In a more biological context, they have been observed between strains of E. coli bacteria Kerr2002 () and between lizard morphs Sinervo1996 (). Another rather different application is to forest fire models Clar1996 (), where the three states trees, fire and ash obey a similar relatioship.

The reaction equations for this model, according to the general treatment in section II, read:

 A+B →A+A, B+C →B+B, C+A →C+C. (10)

The interaction matrix is, accordingly:

 (Aij)=⎛⎜⎝01−1−1011−10⎞⎟⎠. (11)

In order to simplify the calculations, we have set all reaction rates to be equal in Eq. (11). By rescaling time we can then set them to without loss of generality. According to the general treatment in section II, the stochastic model in the large- limit is then described by the Fokker-Planck equation (4) (or, equivalently, the Langevin equations (8)) with . The drift term in Eq. (5) and the noise term in Eq. (6) evaluate explicitely to the following expressions:

 →α = ⎛⎜⎝a(b−c)b(c−a)c(a−b)⎞⎟⎠, (12) B = ⎛⎜⎝a(b+c)−ab−ac−abb(c+a)−bc−ac−bcc(a+b)⎞⎟⎠. (13)

Note that here and in the following we shall use the variable names , , and , , interchangeably. A qualitative treatment of precisely this model was given in Ref. Reichenbach2006 (). In the following, we shall briefly summarize the previous results relevant for our considerations.

The deterministic model, obtained by dropping the noise terms from Eq. (8), is given by the rate equations with as in Eq. (12). Due to the normalization condition , its phase space can be viewed as the -simplex, i.e. an equilateral triangle. Its corners correspond to complete extinction (i.e. only one species of the three species is present), and its edges to states where one species is extinct and two are still present. The dynamics of the rate equations yields oscillations along closed, periodic orbits around a coexistence fixed point at . Close to the fixed point, the orbits are almost circular, whereas further away, they approach the triangular shape of the simplex boundaries (see figure 3). These orbits are neutrally stable due to the existence of a conserved quantity

 ρ=abc. (14)

assumes its maximum value at the coexistence fixed point in the center of the phase space triangle, and its minimum value on the edges, corresponding to extinction of at least one species.

In Ref. Reichenbach2006 (), it was also shown that with the inclusion of noise in the full stochastic model in Eq. (4), the evolution of the population no longer takes place deterministically along one closed orbit, but can fluctuate randomly between different orbits, cf. Fig. 3. By means of a linearization around the coexistence fixed point, it was derived that eventually, the stochastic trajectory will hit one of the simplex boundaries, and from there move towards one of the absorbing corners of the triangle. This process means that two of the three species go extinct when stochasticity is included. It was also motivated that the mean extinction time scales as .

Here, instead of linearizing the stochastic model in Eq. (4), we will perform a stochastic averaging procedure over the deterministic orbits. This will remove the fast, oscillatory degrees of freedom (taking into account all non-linearities and the precise geometry of the phase space) and produce an effective one-dimensional stochastic differential equation for . Through this, we will obtain an exact description of the extinction process and quantitatively correct results for mean extinction times.

### iii.2 Stochastic Averaging

Let us start with the formulation of the stochastic model using the Itô stochastic differential equations (8). Since the deterministic drift terms in Eq. (8) keep conserved, this quantity changes only due to the noise terms , i.e. much more slowly than the oscillations along an orbit with constant . Furthermore, is a measure for closeness to extinction in the sense that the time when becomes for the first time is exactly the time when the first of the three species goes extinct. Thus, a description in terms of allows us to separate the deterministic dynamics (i.e. the rapid oscillations along the closed orbits), which does not contribute to extinction, from the stochastic fluctuations which lead to movement between different orbits and ultimately cause one of the species to die out.

To determine the dynamics of quantitatively, we use Eq. (8) and apply the Itô chain rule Gardiner (), giving:

 dρ = (3∑i=1αi∂iρ+12N3∑i,j=1Bij∂i∂jρ)dt (15) +1√N3∑i,j=1(Cij∂iρ)dWj.

The first term is zero since is conserved by the rate equations. Eq. (15) then implies that , i.e. that changes on a slow time scale and that coexistence in our model is neutrally stable. Actually, there is a more general relationship between the existence of conserved quantities and neutral stability of species 111Any quantity conserved by the rate equations (9) in a system of species with arbitrary reaction rates also satisfies . Thus, as in Eq. (15), the evolution of such a conserved quantity in the stochastic system occurs on a slow time scale . Suppose now there exists a conserved quantity which is a non-trivial function of the abundances of a subset of species , and that denotes a manifold in phase space where at least one of the species in the subset goes extinct. Starting from an initial condition with distinct from , the time to reach this manifold has to scale . This implies that the coexistence of the species in subset is neutrally stable in the classification scheme of section II..

The second term in Eq. (15) is a “stochastic drift” term arising from the fluctuations, and evaluates to:

 12N3∑i,j=1Bij∂i∂jρ=−3Nρ. (16)

Khasminskii’s stochastic averaging theorem Khasminskii1968 () now states that to leading order in , the evolution of is exactly described by the stochastic differential equation

 dρ=−3Nρdt+1√N√¯¯¯¯¯D(ρ)dV, (17)

where is a Wiener process with zero mean and unit variance. is an averaged diffusion coefficient given by

 ¯¯¯¯¯D(ρ) := 1T(ρ)∫T(ρ)0D(a(t),b(t),c(t))dt, (18) D(a,b,c) = [CT(∇ρ)]T[CT(∇ρ)] = (∇ρ)TB(∇ρ)=ρ2(−9+1a+1b+1c).

Note that is the total variance of the noise terms in the Langevin equation for , Eq. (15). Khasminskii’s theorem thus tells us that due to time scale separation, on the slow time scale these noise terms may be treated as independent and replaced by a single effective noise source 222 Khasminskii’s work Khasminskii1968 () actually applies to general Markovian processes with an arbitrary number of slow and fast variables. It shows that on a fixed time interval, when the time scale separation becomes stronger the stochastic process in the slow variables of the full model converges (in some weak sense) to the solution of the stochastically averaged SDE for the slow variables. We refer the interested reader to Ref. Khasminskii1968 () for details on the mathematical formulation, including the precise continuity requirements on the drift and noise terms and the precise statement of the convergence results.. Its variance, given in Eq. (18), is the time-average of over the closed orbit of the deterministic rate equations corresponding to a fixed value of . is the period of this orbit, and , , parametrize this deterministic orbit in terms of the time .

Note that Eq. (15) is not a closed equation. It describes the dynamics of only as a function depending on the dynamics of the variables , , , since the prefactors of the noise terms depended on all three variables individually. In contrast, Eq. (17) is a closed equation uniquely defining a stochastic process , now viewed as a stochastic variable evolving on the interval .

Intuitively, this averaging procedure is justified by the time scale separation described above: The deterministic drift along the orbit with constant takes place on a time scale of fixed by the rate equations, and the movement between different orbits – i.e. the changes in due to the noise terms in the stochastic differential equations (8) – occur on a time scale of . This allows one to average over the fast deterministic evolution and noise in the phase variable, leaving an effective slow process given by Eq. (17).

Eq. (17) can be reformulated equivalently as a Fokker-Planck equation for the probability distribution :

 ∂tP(ρ,t)=1N{∂ρ[3ρP(ρ,t)]+12∂2ρ[¯¯¯¯¯D(ρ)P(ρ,t)]}. (19)

In this form it is most apparent that the time scale of the extinction process is .

In this simple model with its high degree of symmetry, the integral in Eq. (18) can be performed analytically to give a closed expression for in terms of . A sketch of the computation is given in appendix A, and the resulting formula is:

 (20)

Here, and are complete elliptic integrals of the first and second kind, respectively. The elliptic modulus is given by

 k2=(amax−amin)a1(a1−amin)amax, (21)

and , and are the three real roots of the polynomial

 a(1−a)2=4ρ. (22)

As is well-known, these roots can be written down explicitely in terms of . Graphically, is shown in Fig. 4. At the boundary , i.e. near extinction, one obtains the asymptotic form

 ¯¯¯¯¯D(ρ)=−ρlnρ+O(ρ2,ρ2lnρ,ρ2(lnρ)2).

On the other hand, at the boundary , i.e. with nearly equal concentrations of the three species, one obtains the asymptotics

 ¯¯¯¯¯D(ρ)=29(127−ρ)+O(127−ρ)2. (23)

We can now make contact with the results of Ref. Reichenbach2006 (), which were obtained by linearizing in the vicinity of the coexistence fixed point . The radial variable in Eq. (12) of Ref. Reichenbach2006 (), which gives the distance to the coexistence fixed point, is related to our by

 R2=3(127−ρ).

Our SDE (17), independently of the precise form of , predicts an exponential decay for :

 ¯¯¯¯¯¯¯¯¯ρ(t)=ρ(0)e−3t.

Starting from the coexistence fixed point , we obtain for the mean :

 ¯¯¯¯¯¯¯¯¯¯¯¯¯R2(t)=19(1−e−3t). (24)

For small times, this gives which is just the expression in Eq. (30) of Ref. Reichenbach2006 (). The full form of (24) is its correct extension to long times.

We can now go beyond this and compute the fluctuations of near the coexistence fixed point. Rewriting our SDE (17) in terms of the variable , we get

 d(R2)=−3dρ=(13−3R2)dt+√¯¯¯¯¯D(127−13R2)dV.

Near the coexistence fixed point , taking the leading order of the drift term and the asymptotics (23) of , we obtain

 d(R2)=13dt+√23R2dV.

This shows that for small times, the distribution of is exponential:

 P(R2,t)=3te−3R2t.

In particular, we obtain the variance

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯R4−¯¯¯¯¯¯¯R22=t9.

In total, with Eq. (17) and Eq. (19) we have provided a description of the extinction process in the rock-paper-scissors game as an effective one-dimensional stochastic process on the space of the deterministically conserved quantity . This process has a linear drift term and a complicated multiplicative noise which we computed exactly and which is given by Eq. (20). The asymptotics of our results near the coexistence fixed point reproduce the known results in Reichenbach2006 ().

### iii.3 Constant Noise Approximation

In order to avoid the stochastic averaging procedure and the long computation leading to the complicated expression for the noise coefficient in Eq. (20), one might be tempted to take Eq. (17) and simplify it by replacing the multiplicative noise by a constant, additive noise as a rough approximation.

In this section, we will perform this constant noise approximation and compute some observables analytically. In the next section III.4, we will see that close to the boundaries of phase space, numerical results deviate significantly from such a constant noise approximation. This shows that the computation of the nontrivial form of the diffusion coefficient in Eq. (20) is essential in order to obtain quantitatively correct results, especially close to the phase space boundaries.

Replacing the complicated function in Eq. (20) by a constant , the extinction process in Eq.(17) is reduced to a standard Ornstein-Uhlenbeck process:

 dρ=−3Nρdt+1√ND0dV. (25)

In this approximation, the dependence of the mean extinction time on the starting value of can be computed analytically in terms of the generalized hypergeometric function :

 Text(ρ)N = C6√3πD0erfi(√3D0ρ) (26) −ρ2D02F2(1,1;32,2;3ρ2D0).

is a constant fixed by the appropriate boundary conditions. From the singular nature of the boundary at in Eq. (19) and Eq. (20), one can derive that the mean extinction time must satisfy the boundary condition . In the constant noise approximation, this fixes the constant to be:

 C=1ρmaxe−3ρ2maxD0(13+3D02)+√π3D0erf(√3D0ρmax).

The full extinction probability distribution depending on time, starting from a fixed initial condition , can also be written down explicitely if the boundary at is neglected. It then reads:

 Pext(t,ρ0)=erfc⎛⎝ρ0√3D0(e6t−1)⎞⎠. (27)

Asymptotically, this probability distribution possesses an exponential tail, , independent of the precise value of . The exponent coincides very well with previous numerical results obtained in Ref. Ifti2003 ().

### iii.4 Comparison to simulations

To verify the accuracy of the stochastic averaging procedure and the precise noise structure in Eq. (17), we simulated the mean extinction times and the extinction probability distribution of the original reaction system in Eq. (10) using an efficient algorithm due to Gillespie Gillespie1976 (); Gillespie1977 (). The results are shown as crosses in figure 5.

This can then be compared to the predictions of the effective Fokker-Planck equation (19) with the full form of in Eq. (18) obtained by stochastic averaging. Although this effective Fokker-Planck equation cannot be solved analytically, the mean extinction times and survival probabilities can easily be determined numerically from the corresponding backward Fokker-Planck equation Gardiner (). The results are shown as blue lines in Fig. 5. One can observe that the agreement to the simulation of the original reaction system in Eq. (10) is excellent.

Furthermore, Fig. 5 shows the purely analytical results of the constant noise approximation in the previous section (namely, Eqs. (26) and (27)) for . Close to the extinction boundary , there is qualitative but no quantitative agreement between simulations and the constant noise approximation (which is not surprising, considering the shape of close to ). This shows that close to the boundary, the precise form of the multiplicative noise plays a significant role, and provides further evidence for the correctness of Eq. (18). However, we find it interesting that when starting from the coexistence fixed point , the constant noise approximation is in good quantitative agreement with simulations of both the mean extinction time and the extinction probability. This is surprising, since before going extinct the system will have to pass near the boundary , where the approximation fails.

Having given an extensive treatment of the cyclic three-species model, we will now increase the number of species by one and consider the four-species model with cyclic dominance.

## Iv Cyclic Four-Species Model

In this section, we shall apply the formalism developed above to another example. We will consider the cyclic four-species Lotka-Volterra model, which is a natural object to study after the three-species rock-paper-scissors game discussed in the previous section. The reaction equations are given by:

 A+B →A+A, B+C →B+B, C+D →C+C, D+A →D+D. (28)

For the case of equal reaction rates, the drift and noise terms for the Fokker-Planck equation (4) are obtained from Eq. (5) and Eq. (6) as:

 →α=⎛⎜ ⎜ ⎜ ⎜⎝a(b−d)b(c−a)c(d−b)d(a−c)⎞⎟ ⎟ ⎟ ⎟⎠, (29) B=⎛⎜ ⎜ ⎜ ⎜⎝a(b+d)−ab0−ad−abb(c+a)−bc00−bcc(d+b)−cd−ad0−cdd(a+c)⎞⎟ ⎟ ⎟ ⎟⎠. (30)

### iv.1 Rate Equations

The rate equations for this model can be written down directly from the drift term in Eq. (29), and read:

 ∂txi=αi⇔⎛⎜ ⎜ ⎜⎝∂ta∂tb∂tc∂td⎞⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜ ⎜⎝a(b−d)b(c−a)c(d−b)d(a−c)⎞⎟ ⎟ ⎟ ⎟⎠. (31)

As expected, the equations (31) keep the normalization condition invariant. The phase space is now a three-dimensional simplex (i.e. a regular tetrahedron). Again, vertices correspond to extinction of all but one species, edges correspond to extinction of two out of the four species, and faces to extinction of one species.

The fixed points of the rate equations (31) form three lines in the phase space simplex:

• One line of coexistence fixed points given by and , i.e. parametrized by

 (a,b,c,d)=(t,12−t,t,12−t)witht∈[0;12]. (32)
• The edge AC, i.e. all states with coexistence of A and C only, parametrized by

 (a,b,c,d)=(t,0,1−t,0)witht∈[0;1]. (33)
• The edge BD, i.e. all states with coexistence of B and D only, parametrized by

 (a,b,c,d)=(0,t,0,1−t)witht∈[0;1]. (34)

A graphical representation of the structure of the fixed points is shown in figure 6.

It is straightforward to check that the trajectories solving the equations (31) now exhibit two conserved quantities:

 τ1 = ac, τ2 = bd. (35)

The curves given by , are closed, neutrally stable orbits around the line of coexistence fixed points. Close to this line, they are almost circular, while further away they approach the shape of the simplex boundaries. A few exemplary orbits are shown in figure 6.

In total, the deterministic dynamics for the four-species cyclic Lotka-Volterra model is quite similar to the behaviour in the three-species case. The relative abundances of the species oscillate indefinitely along a fixed, closed orbit in phase space, and no extinction takes place.

### iv.2 Stochastic Extinction Process

We would now like to investigate the behaviour of the four-species model when stochastic fluctuations, modelled by the noise terms in Eq. (8), are included. According to the experience from the three-species rock-paper-scissors model, we again expect to see extinction on a time scale since the deterministic orbits are neutrally stable.

However, as already apparent from the description of the rate equation dynamics in the previous section, the set of fixed points is now much larger than in the rock-paper-scissors game. Inserting the parametrizations of the fixed lines in equations (32), (33), (34) into (30), we see that the noise matrix for the Fokker-Planck equation vanishes on the edges AC and BD, but not on the line of coexistence fixed points.

Hence, the absorbing states for the stochastic process are precisely the edges AC and BD of the phase space simplex, parametrized by the equations (33) and (34), corresponding to states with a mixture of non-interacting species (A and C or B and D). This is consistent with the picture obtained directly from the reaction equations (28), where it is clear that these are exactly the states where no more reactions can occur.

In order to analyze the stochastic model quantitatively, we now pursue the same approach used for the rock-paper-scissors model and investigate the behaviour of the deterministically conserved quantities when stochasticity is included. Applying the Itô formula and using the Langevin equations (8), we obtain the following stochastic differential equation for the conserved quantities , :

 dτμ=1√N∑i,j(Cij∂iτμ)dWj. (36)

Note that in contrast to the corresponding calculation in III.2, there are no stochastic drift terms due to the specific form of in Eq. (30). Just as for the noise in the 3-species case, the prefactor of the noises in Eq. (36) still depends on the individual species concentrations , , , which vary along a trajectory with fixed , .

Stochastic averaging now again permits us to obtain a closed system of equations describing the two-dimensional slow process :

 dτμ=1√N2∑ν=1DμνdVν, (37)

As previously, and are now treated as free variables. , are two independent Gaussian noises. The region of - space on which the process given by Eq. (37) occurs is bounded by the absorbing boundaries and , as well as a reflecting boundary at

 √τ1+√τ2=12. (38)

The matrix of the effective diffusion coefficients is defined by , where is the orbit average of the correlation matrix in Eq. (36):

 ¯¯¯¯¯¯Bτ = N(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯⟨dτ1dτ1⟩¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯⟨dτ1dτ2⟩¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯⟨dτ2dτ1⟩¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯⟨dτ2dτ2⟩) (39) = (¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(∇τ1)B(∇τ1)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(∇τ1)B(∇τ2)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(∇τ2)B(∇τ1)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(∇τ2)B(∇τ2)) = (τ1¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(a+c)(b+d)−4τ1τ2−4τ1τ2τ2¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(a+c)(b+d)) =: (τ1f(τ1,τ2)−4τ1τ2−4τ1τ2τ2f(τ1,τ2)).

The high degree of symmetry and the simplicity of the model allow one to give an analytic expression for the function , providing an exact expression for Eq. (39):

 f(τ1,τ2):=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(a+c)(b+d)=h(τ1,τ2)E(k(τ1,τ2))K(k(τ1,τ2)) (40)

Again, and are complete elliptic integrals of the first and second kind, respectively. The elliptic modulus and the helper function are given by:

 k2(τ1,τ2)=4√σ21−σ2h(τ1,τ2), (41) h(τ1,τ2)=12(σ1+√σ21−σ2), (42)

with

 σ1=1−4τ1−4τ2σ2=64τ1τ2.

For details of the calculation leading to these expressions, see appendix B.

Equivalently to Eq. (37), we can write the effective stochastic process as a Fokker-Planck equation:

 ∂tP(τ1,τ2,t)=1N2∑μ,ν=1∂μ∂ν[(¯¯¯¯¯¯Bτ)μνP(τ1,τ2,t)]. (43)

Note how this calculation provides a natural generalization of the analysis performed for the rock-paper-scissors game (section III.2). As is apparent from Eq. (37) (or its Fokker-Planck equivalent Eq. (43)), we again obtain a dynamics for the variables which occurs on a time scale .

In total, we have obtained a complete description of the extinction process in the four-species cyclic model as a two-dimensional diffusion process with varying diffusion coefficients.

### iv.3 Comparison to simulations

As for the rock-paper-scissors game, we can now verify the accuracy with which various quantities of interest for the extinction process can be predicted by the effective stochastic process in Eq. (37). Since this is now a two-dimensional stochastic process in a region with a complicated shape and mixed boundary conditions, it is much harder to treat than the one-dimensional effective process in the three-species case.

Determining the mean extinction times and extinction probabilities from the effective Fokker-Planck equation (37) as was done in the three-species case is not feasible here, since it would require solving elliptic second-order PDE’s over a domain bounded by Eq. (38). Instead, we obtained mean extinction times and extinction probabilities from stochastic simulations of the effective Langevin equations (37) using the XmdS package Hope ().

The results are shown in figure 7. It can be observed that the predictions of the effective Langevin equations (37) compare very well to the results of direct simulation of the original reaction system in Eq. (28). We attribute the slight discrepancies in mean extinction times close to the boundary given by Eq. (38) to the general difficulty of simulating a stochastic process near a curved reflecting boundary. The stochastic averaging method becomes exact at this boundary, since the deterministic orbits are the individual coexistence fixed points.

## V Conclusion

We have analyzed the extinction process in two stochastic Lotka-Volterra models which are neutrally stable in the deterministic formulation. We have seen that when fluctuations are included, the deterministically conserved quantities change slowly (on a time scale proportional to the population size) and drive extinction. After some finite time, only non-interacting species remain.

The separation of time scales between the rapid oscillations described by the deterministic rate equations and the slow movement between different orbits due to noise allowed us to apply the method of stochastic averaging. By doing this, we removed the fast, oscillatory degrees of freedom and gave a quantitative description of the extinction process using effective stochastic differential equations on the space of deterministically conserved quantities.

We have obtained various quantities of interest for the extinction process from these effective equations, and observed that they agree very well with direct simulations.

The stochastic averaging procedure required computation of certain integrals over the closed deterministic orbits, which we were able to perform analytically in the toy models we considered. In more complicated models with less symmetry (e.g. different reaction rates), this may not be possible anymore. However, even for complicated and asymmetric closed deterministic orbits, the averaging required to determine the effective drift and diffusion coefficients can easily be performed numerically. Thus, we think that our approach should be just as helpful for elucidating the impact of noise in more general models.

As we saw, a considerable advantage of the stochastic averaging method (especially in comparison to the treatment in Ref. Reichenbach2006 ()) is that neither the drift nor the noise terms need to be linearized. The full, nonlinear dynamics of the model and the multiplicative noise structure, as well as the complex geometry of the phase space can be taken into account. Furthermore, it is not necessary to write the dynamics explicitely in terms of a radial and a phase variable. This is useful since e.g. in the rock-paper-scissors model in section III, there is no obvious choice for a canonical phase variable.

The idea of describing the extinction process by the evolution of a deterministically conserved variable was also utilized by Parker and Kamenev in Ref. Parker2009 (). They applied it in a semiclassical approximation and obtained the asymptotics of the extinction probability distribution in the standard two-species Lotka-Volterra model. However, our approach is quite different technically and allows a straightforward generalization to more complex models containing more than two species (as in the example in section IV).

We also expect that it should be possible to extend this treatment away from the borderline case of neutral stability to weakly stable or unstable models. Heuristically, this would give rise to a deterministic drift term in equation (17) (or its analogues) which is independent of but controlled by some other small expansion parameter. Investigating how such models can be constructed in a natural way and the details of this generalization requires further research. Indeed, for a recent study along these lines see Refs. Case2010 (); Durney2011 ().

It is also important to understand if and how the present method can be extended to models with spatial degrees of freedom, as discussed in Refs. He2010 (); Tauber2011 (). They exhibit much more complex phenomena (see e.g. Reichenbach2007 (); Szabo2007 (); Perc2007 (); Frey2010 (); Tauber2011 ()) and are more realistic than the well-mixed models we discuss here.

Another approach to investigating the effects of stochasticity on similar models was pursued in Refs. Alonso2007 (); Dauxois2009 (); McKane2005 (); Boland2009 (), where the effects of the stochasticity on the phase variable and the spectral distribution of its oscillation were investigated. This is complementary to our treatment, where we focus on the radial variables instead. This allows us to capture the dynamics of the extinction process.

In a more general context, the present discussion gives an illustration of how stochasticity may change the behaviour in a nonlinear dynamical system qualitatively by adding a stochastic drift to a deterministically conserved quantity. It also provides some ideas for treating the effects of complicated, multiplicative noise on such systems analytically. This may be of considerable interest for non-equilibrium statistical physics in general.

###### Acknowledgements.
We would like to thank T. Reichenbach and J. Cremer for helpful discussions. Financial support of the German Research Foundation via the SFB TR12 “Symmetries and Universalities in Mesoscopic Systems” is gratefully acknowledged.

## Appendix A Computation of noise term for three-species cyclic model

In this appendix, we shall sketch the computation of the stochastic averaging integral in Eq. (18).

We will parametrize an orbit with fixed by one of the species’ concentrations, e.g. . As is easily verified using the method of Lagrange multipliers, the extremal values and which assumes on such an orbit are real roots of the polynomial

 a(1−a)2=4ρ. (44)

The third root of this polynomial is then also real, and will be denoted by . As is well known, explicit expressions for all three roots exist.

We hence write down the factorization

 a(1−a)2−4ρ=(a−amin)(a−amax)(a−a1). (45)

Now, let us give an explicit parametrization of each orbit. We will choose and as the independent variables, with and . Then and are given by

 b1,2 = a(1−a)±√a2(1−a)2−4aρ2a, (46) c1,2 = a(1−a)∓√a2(1−a)2−4aρ2a. (47)

In each case, the + respectively - signs correspond to the two branches of the orbit for a fixed value of .

Inserting Eq. (46) and Eq. (47) into the rate equation we get:

 ∂ta=±√a2(1−a)2−4aρ. (48)

Now, let us calculate the period of an orbit, . By a simple substitution we have

 T(ρ)=∫T(ρ)01dt=2∫amaxaminda∂ta. (49)

The factor arises since each orbit has two symmetric branches, when parametrized by e.g. . Inserting Eq. (48) and Eq. (45), we get

 T(ρ)=2∫amaxaminda√a(a−amin)(a−amax)(a−a1). (50)

This is a standard integral that can be expressed in terms of , the complete elliptic integral of the first kind (see e.g. Byrd1971 ()):

 T(ρ)=4K(k)√(a1−amin)amax. (51)

The elliptic modulus is given by Eq. (21).

The last remaining piece we need is the average of over a deterministic orbit. Again applying a substitution and using Eq. (48) and Eq. (45), we have:

 ∫T(ρ)0dta=2∫amaxamindaa∂ta =2∫amaxamindaa√a(a−amin)(a−amax)(a−a1).

This, too, is a standard integral that can be expressed in terms of complete elliptic integrals (see e.g. Byrd1971 ()), giving:

 ∫T(ρ)0dta = 4[(a1−amax)Π(k2,k)+amaxK(k)]a1amax√(a1−amin)amax (52) = 4√(a1−amin)amaxa1amaxaminE(k) +4K(k)a1√(a1−amin)amax.

Here, is the complete elliptic integral of the second kind, is the complete elliptic integral of the third kind, and the elliptic modulus is again given by Eq. (21). The second line in Eq. (52) was obtained by applying the relation .

Combining Eq. (52) and Eq. (51) we get the result used in Eq. (20).

## Appendix B Computation of noise term for four-species cyclic model

The computation of the average of the correlation matrix for the four-species model, Eq. (39), works along the same lines as the three-species case in appendix A.

We parametrize the deterministic orbit with fixed , in terms of . The extremal values of are given as:

 amax,min=12(1−2√τ2±√(1−2√τ2)2−4τ1). (53)

The other variables are expressed in terms of , and as:

 b1,2 = 12(1−a−τ1a±√(1−a−τ1a)2−4τ2), d1,2 = 12(1−a−τ1a∓√(1−a−τ1a)2−4τ2), c = τ1a. (54)

As for the orbits seen in the three-species case, the plus and minus signs correspond to the two (symmetrical) branches of an orbit for each . From Eq. (54) we get:

 ∂ta=a(b−d)=±√(a(1−a)−τ1)2−4a2τ2. (55)

With all of these results, we can now calculate the period of an orbit with fixed , :

 T(τ1,τ2) = 2∫amaxaminda√(a(1−a)−τ1)2−4a2τ2. (56)

Again, the factor arises from the two symmetric branches of each orbit. The denominator of the integrand can now be factored as:

 √(a(1−a)−τ1)2−4a2τ2 =√(amax−a)(a−amin)(a1−a)(a−a2), (57)

where

 a1,2=12(1+2√τ2±√(1+2√τ2)2−4τ1). (58)

Note that and . Furthermore, Vieta’s theorem gives following relations between , , , :

 a1a2=aminamax=τ1, a1+a2=1+2√τ2, amin+amax=1−2√τ2. (59)

These will be very useful for simplifying some expressions later on.

Upon inserting Eq. (57) into the integral in Eq. (56), we get a standard elliptic integral

 T(τ1,τ2)=4√h(τ1,τ2)K(k(τ1,τ2)). (60)

Here, is a helper function defined by

 h(τ1,τ2):=(a1−amin)(amax−a2). (61)

is the complete elliptic integral of the first kind, with the elliptic modulus given by

 k2=(amax−amin)(a1−a2)(a1−amin)(amax−a2). (62)

By using Eq. (59), Eq. (53) and Eq. (58) we can express the helper function and the elliptic modulus explicitly in terms of and , giving the expressions in Eq. (41) and Eq. (42). Observe that in Eq. (41) and Eq. (42), the symmetry in and (which is required by the cyclic symmetry of the model, but was lost when we chose to parametrize the orbit explicitely using the variable ) is again manifest.

To calculate the average of the noise matrix