Mean extinction times in cyclic coevolutionary rock-paper-scissors dynamics
aa r X i v : . [ q - b i o . P E ] M a r Mean extinction times in cyclic coevolutionaryrock-paper-scissors dynamics
Markus Sch ¨utt and Jens Christian Claussen ∗ Institut f ¨ur Neuro- und BioinformatikUniversit¨at zu L ¨ubeck, D-23562 L ¨ubeck, GermanyFebruary 24, 2010
Abstract
Dynamical mechanisms that can stabilize the coexistence or diversity in biology are generallyof fundamental interest. In contrast to many two-strategy evolutionary games, games with threestrategies and cyclic dominance like the rock-paper-scissors game (RPS) stabilize coexistenceand thus preserve biodiversity in this system. In the limit of infinite populations, resembling thetraditional picture of evolutionary game theory, replicator equations predict the existence of afixed point in the interior of the phase space. But in finite populations, strategy frequencies willrun out of the fixed point because of stochastic fluctuations, and strategies can even go extinct.For three different processes and for zero-sum and non-zero-sum RPS as well, we present resultsof extensive simulations for the mean extinction time (MET), depending on the number of agents N , and we introduce two analytical approaches for the derivation of the MET. There are several examples of cyclic dominance in biology: Males of the californian lizards
Utastansburiana are known to inherit three different mating strategies which cyclically dominateeach other [22, 24]. Another example is given by three strains of
Escherichia coli [14], in tropicalmarine ecosystems [4] or in high arctic vertebrates [10]. Cyclic dominance is also importantin some theoretical models like the susceptible-infected-recovered-susceptible (SIRS) model inepidemiology [2], in cyclic extensions of the Lotka-Volterra model [31, 32] or the Public GoodsGame [13, 11]. Here, cyclic dominance is a way to preserve the coexistence of strategies (what isoften called ‘biodiversity’ in a biological context) in the models.A straightforward model system for cyclic dominance is the rock-paper-scissors game (RPS)well known from schoolyards. It contains the three strategies ‘rock’, ‘paper’, and ‘scissors’ witha simple domination rule: Paper wraps rock, scissors cut paper, rock crushes scissors. The payoffa dominating strategy gets from the dominated one is set to for normalization, the payoff for atie is , and a dominated strategy gets a payoff − s , where we assume − s < . Hence the payoffmatrix of this game reads P = − s
11 0 − s − s . (1) ∗ Corresponding author. or the standard choice s = 1 we have a zero-sum RPS, for all other values the game is nonzero-sum. One can simply intuit the impact of s [5]: For huge s it is important for a player toavoid loosing, so it is successful playing the same strategy as the majority; hence an equilibriumthat includes all three strategies is unstable. On the other hand, for small s it is more important towin occasionally so that the mixed equilibrium becomes stable. Experiments indicate s < forthe lizard example [36] and s > for E. coli [14].For a long time, the resulting dynamics have been described by a meanfield approximation inthe limit of infinite populations. The traditional way to describe such evolutionary games was thestandard replicator equation, an equation of motion for the density of an arbitrary strategy α [12], ˙ x α = ψx α ( π α ( ~x ) − h π ( ~x ) i ) , (2)where ψ is a constant prefactor which can be absorbed in the time scale by a constant rescaling oftime. π α ( ~x ) = x P α, + x P α, + x P α, and h π ( ~x ) i are the payoff of strategy α and the averagepayoff of the whole population, respectively. The adjusted replicator equation [23] ˙ x α = x α Γ + h π ( ~x ) i ( π α ( ~x ) − h π ( ~x ) i ) , (3)has been used less frequently. The prefactor / (Γ + h π ( ~x ) i ) can be absorbed in the time scale by adynamical rescaling of time for symmetric conflicts like the RPS, preserving stability properties.For asymmetric conflicts, like Dawkins Battle of the Sexes [7], this is not possible, as the averagepayoffs in each population do not coincide; hence the adjusted replicator equation can lead toqualitatively modified dynamics [26, 6]. The derivation of the replicator equations and Fokker-Planck-equations (comprising the first-order corrections) for the intrinsic noise are commonlytaken from the master equation of the underlying discrete stochastic process [27, 28, 26, 31, 29].Both replicator equations predict the existence of a fixed point at ( x R , x P , x S ) = (cid:0) , , (cid:1) .The FP is neutrally stable for s = 1 , asymptotically stable for s < , and unstable for s > . Inthe case of finite populations of size N , the population can drive out of the fixed point becauseof stochastic fluctuations, and after some time one of the strategies will go extinct. Once thishas happened, a second strategy will die out soon because of the lack of cyclic dominance, andthe third strategy will therefore survive forever. Several efforts have been made to overbear thisshortcoming especially of the zero-sum RPS, for example spatial discretization of populations (fora review see [25]), mobility [21], the introduction of intelligent update rules (best response [3]),the introduction of the parameter s as mentioned above [5], or the computation of a critical systemsize above which coexistence of strategies is likely [26], but it is still an open question how long ittakes on average until the first strategy has gone extinct. In this paper we investigate the scaling ofthe mean time to the extinction of the first strategy (mean extinction time, MET), depending on thesystem size N and the parameter s , for well mixed populations and three evolutionary processes(with two of them having the standard and accordingly the adjusted replicator equation as limits N → ∞ ), and present two analytical approaches which give theoretical insight in the scaling ofthe MET. Contrary to replicator equations describing the dynamics of relative abundance densities, real pop-ulations are finite (and discrete), appropriate modeling thus bases on discrete stochastic processes(of birth and death). In the classical Moran process [16] an individual a is chosen with probabil-ity proportional to its fitness. It reproduces, and the offspring replaces another randomly chosenindividual b . In the frequency-dependent Moran process [18] which extends the classical Moranprocess [16] by considering coevolution, the fitness is not static but depends on the frequencies of he strategies. For better comparison with the processes mentioned beneath, in each time step wechoose an individual a at random which reproduces with probability φ α ( a ) M ( n α , n β , n γ ) = 12 1 − ω + ωπ α ( a ) − ω + ω h π i , (4)and the offspring replaces another randomly chosen individual b . Greek lettes α, β, γ can assumeeach of the strategies R, S, P , respectively. Here, ω is the selection strength, the n i are the numberof agents playing strategy i , π α = ( n P α, + n P α, + n P α, ) / ( N − the payoff for an individualof strategy α , and h π i the average payoff of the whole population. In the limit N → ∞ , the Moranprocess leads to the so-called adjusted replicator equation [26]. Note that a factor is introducedhere for a better comparability with the two processes mentioned beneath.The Moran process is a well-established stochastic process for evolutionary birth-death (forgrowing population sizes, see e.g. [33]) dynamics with overlapping generations and therefore(in its frequency-dependent extension) serves as a standard model of evolutionary game theory.However, the Moran process requires perfect global information about the whole population, anassumption that can be unrealistic and undesirable. Because of that, two local processes havebeen mentioned. Again we choose an individual a at random for reproduction. Another randomlychosen individual b changes its strategy to the strategy of a with probability φ β ( b ) → α ( a ) lu ( n α , n β , n γ ) = 12 + ω ( π α ( a ) − π β ( b ) )2∆ π max (5)for the local update [26] and φ β ( b ) → α ( a ) F ( n α , n β , n γ ) = 11 + e − ω ( π α ( a ) − π β ( b ) ) (6)for the Fermi process [34], respectively, where ∆ π max is the maximum of the possible payoffdifference. In the limit N → ∞ , these processes lead to the standard replicator equation (localupdate) [26], and a nonlinear replicator equation (Fermi process) [30] with similar properties,respectively. The common approach for the derivation of the equations of motion and the first-order corrections for demographic stochasticity are to derive a Fokker-Planck equation from themaster equation for the respective stochastic process [27, 28, 26, 6, 29].Here we focus on the mean extinction time for the Moran and Local Update processes. Foreach process, the probability of increasing the population strength of the strategy α by in a singletime step and decreasing the population strength of strategy β by , is given by T process ( n α , n β , n γ → n α + 1 , n β − , n γ ) = n α N n β N φ β → αprocess ( n α , n β , n γ ) . (7)This quantity is known as hopping rate. Note that the sum over all hopping rates is = 1 becausereactions are possible that do not lead to changes in strategy frequencies. The time scale is chosenso that one reaction takes place every unit time step, including empty steps with no strategychange. Compared to the time scale a mutant occurs – re-introducing a strategy – the time scale charac-terizing the survival properties of a genetic strain is defined by the mean extinction time. For twostrategies, and fixed N , one has a onedimensional Markov process for which a closed expres-sion allows to proceed analytically as has been demonstrated by Antal and Scheuring [1]. Forhigher-dimensional cases unfortunately exact general solutions of the mean extinction time (or a ææææææææææææææææææææææææææææààààààààààààààààààààààààààààààààà ììììììììììììììììììììììììììììììììì òòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôççççççççççççççççççççççççççççççççç
100 200 300 400 N N (a) ææææææææææææææææææææ æææææææææààààààààààààààààààààààààààààààààà ììììììììììììììììììììììììììììììììì òòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôççççççççççççççççççççççççççççççççç
100 200 300 400 N N (b) æææææææææææææææææææææææææææææààààààààààààààààààààààààààààààààà ììììììììììììììììììììììììììììììììì òòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôççççççççççççççççççççççççççççççççç
100 200 300 400 N N (c) Figure 1: Semi-logarithmic plots of the simulation data for the METs for the local update process(a), the Moran process (b) and the Fermi process (c). The METs were divided by N . Each pointis an average over the extinction times of 1000 simulations. Red squares: s = 1 . , ω = 0 . , blackrhombi: s = 1 . , ω = 0 . , orange upwards triangles: s = 1 . , ω = 0 . , yellow downwardstriangles: s = 1 . , ω = 0 . , green open circles: s = 0 . , ω = 0 . , blue filled circles: s = 0 . , ω = 0 . . For the Moran process ω = 0 . has been used instead of ω = 0 . to keep the transitionprobabilities in the interval [0 , . mean-first-passage time [19]) problems are not known; in this case the situation is furthermorehampered by the location-dependent dynamics within a simplex boundary. So in most cases onehas to rely on systematic numerical investigations. We have carried out extensive simulations toquantify the mean times until one of the three strategies has gone extinct. We have analyzed themean extinction times for all three described processes depending on the number of agents N andthe parameters s and ω . In general, we observe the following properties (see Fig. 1): For all threeprocesses the MET becomes independent of the selection strength ω for s = 1 . . Further on, theMET for s = 1 . is identical for all three processes and proportional to N , h t ext i ( s = 1) ≈ (0 . ± . N . (8)As expected, for s < the MET is larger and for s > smaller than for s = 1 , but for weakselection the difference is small. In a logarithmic plot of the MET divided by N one can easilysee that h t ext ( s = 1 . i ∝ N , but also h t ext ( s < i ∝ N e f ( s,ω ) N with f ( s, ω ) > andindependent of N . For the Moran process the stabilizing effect for s < is larger than for theother processes.On the first view, for s > the MET seems to be h t ext ( s > i ∝ N e − f ( s,ω ) N . Butsuch a proportionality would imply a disappearance of the MET for large N . This is unrealisticbecause even the shortest way would require N/ steps. For that a dependence h t ext ( s > i = c N + c N e − f ( s,ω ) N is more likely (with f ( s, ω ) having the same properties as for the s < case), as it is predicted by one of our analytical approaches.In summary, it can be said that the MET switches from polynomial to exponential scaling at s = 1 . This agrees with the drift reversal picture from [26, 5]. So the (global) Moran process andthe (local) local update and Fermi process behave similar, but the impact of the parameters s and ω is much stronger for the Moran process. An exact analytical solution for two- or higher dimensional Markov chains is impossible, and a di-rect computation of the MET of such a three-strategy evolutionary game is not feasible. Althoughsome efforts have been made, for example the work of Reichenbach et al. [20], wich have ana-lyzed mean extinction properties in an urn model of a three-strategy game, employing the usual pproach of applying van Kampen’s linear noise approximation and deriving a Fokker-Planckequation [35, 27, 28, 26].Adapting this approach to the zero-sum RPS for the local update and the Moran process, wefind a predicted MET proportional to N and independent of ω , which is in accordance with thesimulation data. Unfornately this approach has some shortcomings: Although the overall scalingis correct, the numerical value (prefactor) of the predicted MET is too great by a factor of ≈ , andit is not possible to use this approach for the non zero-sum RPS or in the Fermi process. For thatwe have developed two analytical approaches which do better. We will present both approachesfor the local update. Following the same schemes for the other two processes is, although possible,a bit more sophisticated because not all integrals can be done analytically here in general. As wehave seen from the numerical investigations, the s = 1 payoffs override differences betweenthe underlying processes, so that it is warranted to concentrate on one analytically more feasibleprocess in the remainder. To compute – with help of approximations – the MET for the general case, we need an appropriatecoordinate system. For the standard replicator equation and s = 1 , H = − x R x P x S (9)is a constant of motion which assumes the value H = − / in the fixed point and H = 0 on theedge of the phase space, the simplex S . Here, x R , x P and x S are the frequencies of the strategies R , P and S , respectively. For s < , H is a Lyapunov function of the replicator equations with ˙ H < , and so the inner fixed point ist asymptotically stable [5]. For s > the fixed point isunstable, and the attractor of the system moves to the edge of the simplex. Via a transformationof the fixed point to the origin of the phase space, and by inserting x R + x P + x S = 1 , whichguarantees the conservation of the total density, we find for H : H = − ( x + 1 / y + 1 / / − x − y ) (10)For large absolute values of H , the curves with H = const. resemble slightly deformed circles,but for H → they approach the triangle form of the simplex. In the following we will use H asan observable to measure the effective distance to the fixed point. But obviously we need a secondcoordinate to describe a point in the phase space, so we will use a system with the variables x and H by eliminating the y -coordinate. This gives us two solutions for y ( x, H ) , an upper branchabove the fixed point, and a corresponding lower branch below. We need both because it is notpossible to describe all points of the phase space ( x, y ) by only one solution (see figure 2(b)). Sowe have y u = − x − x − √ H + 12 x + 324 Hx − x − x + 81 x x ) (11)and y d = − x − x + √ H + 12 x + 324 Hx − x − x + 81 x x ) . (12)Here the indices ‘ u ’ and ‘ d ’ stand for ‘up’ and ‘down’, respectively, because by re-inserting of x and H in the first (second) equation one will get only values of y that lie above (under) a separator.In two points, x min and x max , both solutions are identical so that the curve H = const. is closed.We can compute these from the condition y u = y d and solving for x , x R - x P (a) æææææææææææææææææææææææææææææææææààààààààààààààààààààààààààààààààà ììììììììììììììììììììììììììììììììì òòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòò ææàà ìì òòôô ççáá íí óó
100 200 300 400 N €€€€€€€€ tN (b) Figure 2: (a) A plot of curves for which H = const. holds, for different values of H . Differentcolours on each curve mark the different solution branches y u und y d . The black line connects thevalues x min and x max , respectively. (b) Plots of the simulation data for the MET from the localupdate (symbols) for different pairs of s and ω depending on N and the corresponding results fromthe analytical approach presented in section 4.1 (lines, with a not filled version of the correspondingsymbol at the end). Red circles: s = 1 . , ω = 0 . , orange squares: s = 1 . , ω = 0 . , yellow triangles: s = 1 . , ω = 0 . , green rhombi: s = 0 . , ω = 0 . . x min = 13 (cid:18) − cos (cid:18)
13 arg (cid:16) − H + 6 √ p H (27 H + 1) − (cid:17)(cid:19) (13) −√ (cid:18)
13 arg (cid:16) − H + 6 √ p H (27 H + 1) − (cid:17)(cid:19) + 1 (cid:19) and x max = 13 (cid:18) − cos (cid:18)
13 arg (cid:16) − H + 6 √ p H (27 H + 1) − (cid:17)(cid:19) (14) + √ (cid:18)
13 arg (cid:16) − H + 6 √ p H (27 H + 1) − (cid:17)(cid:19) + 1 (cid:19) and a third (formal) solution which here is of no use because for typical values of H it producesvalues outside the simplex. arg( z ) is the argument of the complex number z .Now we compute the expectation value of the change in H for a single time step for the localupdate process. For time t , the system may be in an arbitrary allowed point ( x, H ) . For that wewill need the hopping rates in the ( x, H ) coordinates. For example, for the change δ = N , δ = − N the hopping rate reads after inserting y = y u , T u (cid:18) N , − N (cid:19) = − (cid:0) − x + 3 x + f + 2 (cid:1) x + 1) (cid:0) ψx + 9( ψ − x + (2 s + 1) f ψ − (cid:1) (15)with f = p (3 x + 1) ((3 x + 1)(2 − x ) + 108 H ) . (16)Note that we will use all variables as if they were continous, although in fact they are not, andwe use T ( δ , δ ) as an abbreviation for T ( x , x , x → x + δ , x + δ , x + δ ) due to theconservation of total density (where δ und δ can take all values of {− /N, , /N } as long asthey are not identical). - - H - N ∆ H mean (a) - - - H - ´ - - ´ - - ´ - ´ - ´ - ∆ H mean (b) Figure 3: (a) A plot of the expectation values of the changes in H , δH ( H ) , multiplied with N , forthe local update process. The dashed lines show this quantity for N = 798 , the continous lines for N = 396 . red : s = 1 . , ω = 0 . , orange: s = 1 . , ω = 0 . , yellow: s = 1 . , ω = 0 . , green: s = 0 . , ω = 0 . , blue: s = 0 . , ω = 0 . . (b) Additional illustration of the meaning of δH ( H ) ( s = 0 . , ω = 0 . ). If the observable h has a certain value H t at a time t , H is expected to make astep in the direction of the arrow belonging to this specific value of H . The lengths of the arrows areover-emphasized by a factor of . Further on we will need the change in H that emerges from a change of the ( x, y ) -values inthe ( x, H ) -coordinates as well. To derive this, we write H ( x + δ , y + δ ) of the changed valuesof x and y , subduct the initial value H ( x, y ) , and transform the result into the ( x, H ) -coordinates.For y = y u and the same change as in the example above we get ∆ u (cid:18) N , − N (cid:19) = (cid:0) − N x − N + 2) x + N g − (cid:1) xN + N ) (cid:0) x − x + f − (cid:1) (17)with g = p (3 x + 1) (27 x − x + 108 H + 4) (18)By multiplying all hopping rates T ( δ , δ ) with the associated changes in H , ∆ ( δ , δ ) , andsumming over all terms, we get a relatively simple equation for the expectation value of thechange in H in a single time step, δH ( x, H ) = X δ ,δ T ( δ , δ ) ∆ ( δ , δ ) (19) = − H (9 + 27 x + (1 + 27 H ) N ( − s ) ω ∆ π max + 27 N ( − s ) x ω ∆ π max )3 N (1 + 3 x ) , with δH u ( x, H ) = δH d ( x, H ) =: δH ( x, H ) . Averaging over x by integrating over x from x min to x max and normalizing by the length of whole interval we derive δH ( H ) = (cid:0) a + (6 b − a + 6 b − b + 2 (cid:1) N ( s − ψ + 186 N − H ( s − ψ log (cid:16) b +13 a +1 (cid:17) ( a − b ) N (20)with a = x min , b = x max . The simplest way to approximate the mean extinction time for asystem starting in the fixed point of the replicator equations is then to define a map H i +1 = H i + δH ( H i ) (21) nd to iterate it until H h t i ≈ is reached. The MET is then simply the number of necessary steps.As one can see in the right hand side of Figure 2(a), this result is in good agreement with thesimulation data, for the dependence of s and ω , as well as for the numerical values. For s = 1 theMET is predicted to be proportional to N and independent of ω because δH ( H ) is independentof ω for s = 1 . It is in good agreement with the simulations that the MET for s = 1 and weakselection differs only slightly from the MET for s = 1 , while for strong selection and s > it isclearly smaller than for s = 1 .Unfornately, for s < this approach works only for small selection and not to great N ,because for s < non-positive expectation values of the changes in H are possible, so if we startthe mentioned iterated map in the fixed point and repeat the iteration again and again, we wouldarrive at a value of H where the expected change is at some time, and the predicted MET wouldbe infinite. The part of the interval where the expectation value is negative is bigger for larger ω than for smaller ω , and it grows with increasing N , so in this cases this approach is not sufficient. For the local update process, a Fokker-Planck equation (FPE) for the time evolution of the proba-bility density is given by (summation over double indices) ∂ t P ( ~x, t ) = − ∂ i [ α i ( ~x ) P ( ~x, t )] + 12 ∂ ij [ B ij ( ~x ) P ( ~x, t )] , (22)where the coefficients are α i = X δ~x δx i T ( ~x → ~x + δ~x ) (23)and B ij ( ~x ) = X δ~x δx i δx j T ( ~x → ~x + δ~x ) , (24)which gives the following results α R = − (2 x P + x R )(1 + 3 x R ) ψ N (25) α P = (1 + 3 x P )( x P + 2 x R ) ψ N (26) B RR = 2 + 3 x R − x R N (27) B RP = B P R = − (1 + 3 x P )(1 + 3 x R ) N (28) B P P = 2 + 3 x P − x P N . (29) ransforming this into polar coordinates x R = r cos φ , x P = r sin φ , we find for the coefficientsof the FPE α R = − rψ (3 r cos( φ ) + 1)(2 cos( φ ) + 2( s + 1) sin( φ ) + r ( s − φ ) + 2))6 N (30) α P = − rψ (3 r sin( φ ) + 1)( − s + 1) cos( φ ) − s sin( φ ) + r ( s − φ ) + 2))6 N (31) B RR = − r cos ( φ ) + 3 r cos( φ ) + 29 N (32) B RP = B P R = − (3 r cos( φ ) + 1)(3 r sin( φ ) + 1)9 N (33) B P P = − r sin ( φ ) + 3 r sin( φ ) + 29 N , (34)where ψ = ω/ ∆ π max , and hence for the FPE itself P (0 , , = 6 N ( s − ψ sin(2 φ ) r + N (cid:0) r − (cid:1) ( s − ψ − N P (0 , , (35) + (cid:18) rψ ( − s + 2) cos( φ ) + (7 s + 2) cos(3 φ ) − s + (2 s + 7) cos(2 φ ) + 5) sin( φ ))12 N − ( s + 1) ψ (sin(2 φ ) + 2)6 N + cos( φ ) + 3 cos(3 φ ) − sin( φ ) + 3 sin(3 φ )12 N r + cos(2 φ )9 N r (cid:19) P (0 , , + (cid:18) r ψ ((2 s + 1) cos( φ ) + (2 s + 7) cos(3 φ ) − ( s + 2) sin( φ ) + (7 s + 2) sin(3 φ ))12 N + r ( s − ψ (cos( φ ) sin( φ ) + 1) N + r ( − N ( s − ψ + N ( s + 1) cos(2 φ ) ψ − N + cos( φ ) − cos(3 φ ) + sin( φ ) + sin(3 φ )8 N + cos( φ ) sin( φ ) + 19 N r (cid:19) P (1 , , − (cos( φ ) − sin( φ ))(9 cos( φ ) sin( φ ) r + 3 r + cos( φ ) + sin( φ ))9 N r P (1 , , + 9 r cos( φ ) − r cos(3 φ ) + 9 r sin( φ ) + 4 sin(2 φ ) + 9 r sin(3 φ ) + 872 N r P (0 , , + − r + 3 cos( φ ) r + 9 cos(3 φ ) r + 3 sin( φ ) r − φ ) r − φ ) + 872 N P (2 , , , with P ( λ,µ,ν ) := ∂ λ ∂r λ ∂ µ ∂φ µ ∂ ν ∂t ν P ( r, φ, t ) .If we now assume the probability density to be independent of the angle φ (which is agood approximation at least for small r ), all terms which contain derivatives with respectto φ will drop out. By averaging over φ afterwards by integrating from to π over φ anddividing the result by π we get a probability density which is only independent of r and t . Hence the FPE reads ∂ t P ( r, t ) = r ( N (12 r −
1) ( s − ψ − N √ r P ( r, t ) (36) + 3 ( N (6 r −
1) ( s − ψ − r + 218 N √ r ∂ r P ( r, t )+ r (2 − r )18 N √ r ∂ r P ( r, t )=: a ( r ) P ( r, t ) − v ( r ) ∂ r P ( r, t ) + D ( r ) ∂ r P ( r, t ) . (37)9 æàà ìì òò ôô r - ´ - - ´ - - ´ - ´ - ´ - ´ - v H r L æææææææææææææææææææææææææææææææææà à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à ììììììììììììììììììììììììììììììììì òòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòôôôôôôôôôôôôôôôôôôôôôôôôôôôôô ææàà ìì òòôô ççááíí óóõõ
100 200 300 400 N €€€€€€€€ tN Figure 4: Left: A plot of the convection velocity v ( r ) for the local update with N = 3960 . As onecan see, v ( r ) is ≈ for s = 1 . This corresponds to the neutral stability of the fixed point in thelimit N → ∞ for s = 1 . For s < v ( r ) is negative, corresponding to the asymptotical stable FP,and for s > v ( r ) is positive corresponding to the unstable FP. For r → v ( r ) diverges for allparameter values to −∞ , which results from the simplifications in combination with the existence ofa FP at r = 0 . Right: Semilogarithmic plot of simulation data of the MET in the local update process,divided by N (symbols) for several pairs of s and ω , and the corresponding values predicted by theapproach in section 4.2 (lines, with a not filled version of the corresponding symbol at the end). Redcircles: s = 1 . , ω = 0 . , orange squares: s = 1 . , ω = 0 . , yellow rhombi: s = 1 . , ω = 0 . ,green upward triangles: s = 0 . , ω = 0 . , blue downward triangles: s = 0 . , ω = 0 . Let us now approximate the diffusion by its value at r = 0 , this gives us D = 19 N . (38) r can theoretically vary between and / for some values of φ , but most extinctionevents will happen around those point of the edge of the simplex that are closest to theinner fixed point. These points have the distance r & from the inner FP. Because ofthat we search for the solution of the one-dimensional random walk descripted by eq. 37in the interval [0 , L ] = [0 , ] . Next, let us approximate the convection velocity v ( r ) by itsvalue for mediate r , e.g. r = = L , and we get v = − . Nψ − . N , s = 0 . − N , s = 1 . . Nψ − . N , s = 1 . (39)(We cannot approximate it by its value at r = 0 as we have done for the diffusion because v (0) is −∞ corresponding to the FP.) For a last simplification, we neglect the term thatincludes the probability density itself. Hence, the FPE reads ∂ t P ( r, t ) = − v∂ r P ( r, t ) + D∂ r P ( r, t ) . (40)The computation of the MET for such an equation with one reflecting border at r = 0 and one absorbing border at r = L is a standard problem, for example, see [19]. One can10nd the MET by introducing a time integrated probability density, P ( r ) = Z ∞ P ( r, t ) dt. (41)With this, the time-integrated FPE reads: v ∂P ( x ) ∂r = D ∂ P ( x ) ∂x (42)with the boundary conditions P ( r ) = 0 at r = L and J ( r ) = R ∞ j ( r, t ) dt = R ∞ vP ( r, t ) − D∂ r P ( r, t ) dt = 1 at r = 0 . There, J ( r ) is the time-integrated probability current. Theconstraint J (0) = 1 corresponds to the injection of an unit probability current at t = 0 ,this means a start of the system in the origin. The solution of the differential equation isgiven by P ( r ) = 1 v (cid:0) − e v ( x − L ) /D (cid:1) (43)Then the MET is simply h t i = Z L P ( r ) dr = Lv − Dv (1 − e − vL/D ) , (44)what leads to h t i = N ( − . Nψ +0 . e . Nψ − . ) (0 . Nψ +0 . , s = 0 . . N , s = 0 e − . Nψ N ( e . Nψ (0 . Nψ − . . ) (0 . − . Nψ ) , s = 1 . (45)or in general h t i = 72 e − N ( s − ψ − N (cid:16) e N ( s − ψ + (5 N ( s − ψ − e (cid:17) (36 − N ( s − ψ ) . (46)The proportional dependence of the MET on N and the other parameters is very good,as one can see in Fig. 4. This approach predicts the MET to be proportional to N and independent of the selection strength ω for s = 1 , while for s < it forecasts theMET to be substantially proportional to N e κ ψN , and for s > in good approximationproportional to c N + c e κ ψN . Evolutionary games with three cyclically dominating strategies can stabilize the coexis-tence of strategies superior compared to games with only two strategies which includea dominating strategy. In finite populations, it is still likely that one of these strategiesgoes extinct, though on average it takes a time which is considerably longer than in twostrategy games. We have analyzed the mean times to the extinction of one of the three11trategies numerically in the RPS-game for three different evolutionary processes, the lo-cal update, the frequency-dependent Moran process and the Fermi process. These meanextinction times (MET) have the same fundamental dependencies of the number of agents N , the parameter s and the selection strength ω .For the zero-sum RPS game ( s = 1 ), the MET is proportional to N and independentof ω for all three processes; even the proportionality constants are the same. For s < , wefind a stabilization of biodiversity, as the MET grows exponentially with N . In contrast tothat, for s > , there is a destabilization of biodiversity, and the MET now grows smallerthan in the zero-sum RPS (where it grows ∝ N ), namely, a part of the proportionalityfactor of the zero-sum RPS now decays exponentially with N . In both non zero-sumcases, the coefficient that is multiplied with the N in the exponent, grows with both ω and ( s − .Solving such two-dimensional Markov chains is not possible, but we have developedtwo analytical approaches for the approximation of the dependencies of the MET of N and the parameters s and ω . The first approach can be used for s ≥ only, but servesfor all three processes. In contrast, the second approach gives a good approximation forall parameter values, but it can be used completely analytical only for the local updateprocess, as for the other two processes not all integrals can be done analytical.In summary, the fixation time – as a long-time property of the process – shows acrossover from exponential to polynomial scaling with the population size, being fullyconsistent with the critical population size deived from the “drift reversal” picture whichis based on the short-time dependence of the average drift in phase space - repellingversus attracting towards the fixed point. Thus both approaches to describe stochasticstability in finite populations here lead to the same conclusions: coexistence is preserved(stabilized) for a positive-sum cyclic game, whereas negative-sum games as well as smallpopulations destabilize coexistence. References [1] Antal, T., and Scheuring, I.,
Fixation of Strategies for an Evolutionary Game inFinite Populations , Bulletin of Mathematical Biology 68,1923-1944 (2006)[2] Axelrod R., and May, R.M.,
Infectious Diseases of Humans , (Oxford UniversityPress, Oxford, 1991)[3] Blume, L.E.,
The Statistical Mechanics of Best-Response Strategy Revision , GamesEcon. Behav. , 111-145, (1995)[4] Buss, L.W., Competitive intransitivity and size-frequency distributions of interact-ing populations , Proceedings of the National Academy of Science, USA, , 5355(1980)[5] Claussen, J.C., and Traulsen, A., Cyclic dominance and biodiversity in well mixedpopulations , Phys. Rev. Letters , 058104, (2008)[6] Claussen, J.C.,
Drift reversal in asymmetric coevolutionary conflicts: Influence ofmicroscopic processes and population size,
The European Physical Journal B ,391-399(2007)[7] Dawkins, R., The Selfish Gene , Oxford University Press, (New York, 1976)128] Durrett R., and Levin, S.,
Allelopathy in spatially distributed populations , Journalof Theoretical Biology , 165-174 (1997)[9] Durrett, R., and Levin, S.,
Spatial aspects of interspecific competition , TheoreticalPopulation Biology , 30-43, (1998)[10] Gilg, O., Hanski, I., and Sittler, B., Cyclic dynamics in a simple vertebrate predator-prey community , Science , 866 (2003)[11] Hauert, C., de Monte, S., Hofbauer, J., and Sigmund, K.,
Volunteering as red queenmechanism for cooperation in public goods games , Science , 1129 (2002)[12] Hofbauer J., and Sigmund, K.,
Evolution and the Theory of Games , CambridgeUniversity Press, Cambridge, England, (1998)[13] Kagel J.H., and Roth, A.E., (eds.),
The Handbook of Experimental Economics ,Princeton University Press, Princeton, (1995)[14] Kerr, B., Riley, M.A., Feldman M.W., and Bohannan, B.J.M.,
Local dispersal pro-motes biodiversity in a real-life game of rock-paper-scissors , Nature , 171-174(2002)[15] Kirkup B.C., and Riley, M.A.,
Antibiotic-mediated antagonism leads to a bacterialgame of rock-paper-scissors in vivo , Nature , 412-414 (2004)[16] Moran, P.A.P.,
The statistical processes of evolutionary theory , Clarendon Press(1962)[17] Nowak, M.A.,
Evolutionary Dynamics , Harvard (2007)[18] Nowak, M.A., Sasaki, A., Taylor, C., and Fudenberg, D.,
Emergence of cooperationand evolutionary stability in finite populations , Nature (London) , 646 (2004)[19] Redner, S.,
A Guide to First-Passage Processes , Springer, (Berlin, 1983)[20] Reichenbach, T., Mobilia, M., and Frey, E.,
Coexistence versus extinction in thestochastic cyclic Lotka-Volterra model , Phys Rev E , 051907 (2006)[21] Reichenbach, T., Mobilia, M., and Frey, E., Mobility promotes and jeopardizes bio-diversity in rock-paper-scissors games , Nature , 1046-49 (2007)[22] Sinervo, B, and Lively, C.,
The rock-paper-scissors game and the evolution of al-ternative male strategies , Nature , 240, (1996)[23] Smith, J.M.,
Evolution and the Theory of Games , (Cambridge University Press,Cambridge, England, 1982)[24] Smith, J.M.,
The games lizards play , Nature , 198 (1996)[25] Szab´o, G., and F´ath, G.,
Evolutionary games on graphs , Physics Reports, , 97-216 (2007)[26] Traulsen, A., Claussen, J.C., and Hauert, C.,
Coevolutionary dynamics: from finiteto infinite populations , Physical Review Letters , 238701 (2005)[27] D. Helbing, Interrelations between stochastic equations for systems with pair inter-actions , Physica A , 29 (1992)[28] D. Helbing,
Stochastic and Boltzmann-like models for behavioral changes, and theirrelation to game theory , Physica A , 241 (1993)1329] Chalub F.A., Souza M.O.,
From discrete to continuous evolution models: A unifyingapproach to drift-diffusion and replicator dynamics , Theor. Pop. Biol 76, 268 (2009)[30] Traulsen, A., Pacheco, J. M., and Nowak, M.A.,
Pairwise comparison and selectiontemperature in evolutionary game dynamics , J. Theor. Biol. 246, 522 (2007)[31] McKane, A. J., and Newman, T. J.,
Predator-Prey Cycles from Resonant Amplifica-tion of Demographic Stochasticity , Phys. Rev. Lett. 94, 218102 (2005)[32] Butler, T., and Goldenfeld, N.,
Robust ecological pattern formation induced by de-mographic noise , Phys. Rev. E 80, 030902 (2009)[33] Poncela, J., G´omez-Garde˜nes, J., Traulsen, A., and Moreno, Y.,
Evolutionary gamedynamics in a growing structured population , New J. Phys. 11 083031 (2009)[34] Traulsen, A. Nowak, M.A., and Pacheco, J.M.,
Stochastic dynamics of invasion andfixation , Physical Review E , 011909 (2006)[35] van Kampen, N., Stochastic Processes in Physics and Chemistry , (North Holland,Amsterdam, 1981)[36] Zamudio, K.R., and Sinervo, B.,
Polygyny, mate-guarding, and posthumous fertil-ization as alternative male mating strategies , Proceedings of the National Academyof Science, USA,97