Cooperation dilemma in finite populations under fluctuating ...
Cooperation dilemma in finite populations under fluctuating environments Michael Assaf,1 Mauro Mobilia,2 and Elijah Roberts3
arXiv:1305.6580v1 [q-bio.PE] 28 May 2013
2
1 Racah Institute of Physics, Hebrew University of Jerusalem, Jerusalem 91904, Israel Department of Applied Mathematics, School of Mathematics, University of Leeds, Leeds LS2 9JT, U.K. 3 Department of Biophysics, Johns Hopkins University, Baltimore, MD 21218, U.S.A.
We present a novel approach allowing the accurate study of rare events in a broad class of evolutionary processes under fluctuating environments. Our approach consists of mapping these processes, not necessarily displaying metastability, onto auxiliary models characterized by metastability that can be analyzed semiclassically. This enables us to study the interplay between extrinsic (environmental) noise and intrinsic (demographic) fluctuations on the statistics of interest. We illustrate our theory for the paradigmatic prisoner’s dilemma game whose evolution is described by the probability that cooperators take over and fixate the entire population. We demonstrate both analytically and by using Monte Carlo simulations that extrinsic noise may drastically enhance such a fixation probability. In particular, we show that under strong, short-correlated extrinsic noise, the cooperation fixation probability decays algebraically with the population size, rather than exponentially. Finally, we show how our approach generalizes earlier works, notably in population genetics. PACS numbers: 05.40.-a, 02.50.Ey, 87.23.Kg, 02.50.Le
Understanding the origin of cooperative behavior and how it is influenced by the population’s intrinsic properties and by environmental factors are major scientific puzzles [1] that are suitably investigated in the framework of evolutionary game theory (EGT) [2]. In EGT, successful species spread at the expense of others, and the optimization of the reproductive potential (fitness) at an individual level can result in the reduction of the population’s overall fitness, a phenomenon metaphorically captured by the prisoner’s dilemma (PD) game [1, 2]. While in EGT the dynamics is traditionally studied in terms of differential equations, demographic fluctuations – intrinsic noise (IN) – are known to affect the evolution in finite populations. In this case, the dynamics is often described by a Markov chain and characterized by the fixation probability of a given trait, which is the probability that the trait invades the entire population [3]. For the classic PD (with IN only), the cooperation fixation probability (CFP) vanishes exponentially with the population size, see e.g. [4], and defection prevails leading to a cooperation dilemma. This prediction, at odds with many observations, has motivated the investigation of various mechanisms that can promote cooperation [5]. In addition to IN, there often is also environmental (extrinsic) noise (EN) which results in random fluctuations of the interaction parameters. In theoretical population genetics [3, 6–8] and ecology [9–11], it has been shown that the combined effect of IN and EN can significantly affect the extinction properties of a population. While many earlier works focused on systems exhibiting metastability and/or short-correlated EN [3, 6–10], here we present a novel approach allowing one to analyze the combined influence of IN and EN, with arbitrary correlation time and magnitude, on a broad class of systems not necessarily displaying metastability. This is done by mapping the system at hand onto one possessing a long-
lived metastable state and treating the latter semiclassically. We illustrate our approach with the paradigmatic PD game, and find that EN drastically enhances the CFP and may even alter its functional dependence on the population size from an exponential to a power law. The paradigm of social dilemma is provided by the classic PD, whose main features are captured by assuming that the pairwise interaction between cooperators and defectors is described in terms of the benefit b and cost c of cooperation, with b > c > 0 [2]. Here, mutual cooperation leads to a payoff b−c > 0 and mutual defection gives a payoff 0, while when one player defects and the other cooperates, the former gets a payoff b and the latter gets −c. The quantity r ≡ c/b is the cost-to-benefit ratio [2] and the dilemma arises from the fact that each individual is better off defecting, even though r < 1 and mutual cooperation enhances the population overall payoff. Here, we consider a finite population comprising N individuals on a complete graph with n cooperators and N − n defectors. At mean field level (N → ∞), defection prevails and the fraction x ≡ n/N of cooperators in the population evolves to extinction according to the replicator rate equation (RE) (d/dt)x ≡ x˙ ∝ x(1−x) [ΠC (x) − ΠD (x)]. Here, ΠC = bx−c and ΠD = bx are the cooperator and defector average payoffs, respectively [2], and we henceforth assume that b, c = O(1). When the population size is finite, fluctuations alter the mean field predictions and the dynamics is usually described in terms of the Moran model (MM) [2, 4, 12, 13], according to a continuous-time Markov chain with the birth/death transition rates Tn+ = [fC (n)/f¯(n)]n(N − n)/N 2 and Tn− = [fD (n)/f¯(n)]n(N − n)/N 2 . Here, the cooperators/defectors fitnesses are fC (n) = 1 + s [b(n/N ) − c] and fD (n) = 1 + sbn/N, (1) where the selection intensity s > 0 has been intro-
2
∂t P (x, t) = −∂x [A(x)P (x, t)]+1/(2N )∂x2 [B(x)P (x, t)], (2) with A(x) = T + (x)−T − (x), B(x) = T + (x)+T − (x), and T ± (x) = Tn± . Note, that for N → ∞, Eq. (2) approaches the RE x˙ = A(x) = sx(1−x)[ΠC (x)−ΠD (x)]/f¯(x) [2, 17]. An important notion to characterize evolutionary dynamics is the fixation probability φC (x0 ) that cooperation prevails starting from an initial fraction x0 of cooperators. In the absence of environmental fluctuations, φC can be calculated exactly [4, 15]. Yet, it is useful to proceed with an asymptotic calculation (N 1) of φC (x0 ) using an auxiliary problem. For this, we consider the modified model obtained by supplementing the original PD model with a reflecting boundary at n0 = N x0 by im− = 0. Hence, the only absorbing state of the posing Tn=n 0 modified model is the state n = N and as a result, since x˙ < 0, a quasi-stationary distribution (QSD) peaked at x0 (for any value of x0 ) forms after a relaxation time of order O(s−1 ). This QSD slowly decays due to a slow leakage of probability into the absorbing state at x = 1, with a rate given by the inverse of the cooperation mean fixation time (MFT). Therefore, after a long time growing exponentially with N , cooperators eventually take over. Importantly, we now demonstrate that φC (x0 ) in the original PD problem coincides to leading order with the inverse of the MFT in the modified problem. The MFT in the modified problem can be computed by employing the metastable ansatz P (x, t) ' π(x)e−t/τ , where π(x) is the QSD and τ is the MFT. Then, Substituting the WKB ansatz π(x) ∼ e−N S(x) in Eq. (2) yields a stationary Hamilton-Jacobi equation (HJE) H(x, px ) = 0, with the Hamiltonian H(x, px ) = px A(x)+(p2x /2)B(x). In this formalism S(x) is called the action function while p(x) ≡ S 0 (x) is called the momentum [18, 19]. Fixation occurs along the zero-energy trajectory (“optimal path”) px (x) = −2A(x)/B(x), where px (x) ∼ O(s) R x 10 in 0the weak-selection limit. This gives S(x) = px (x )dx = (c/b) ln(2 − cs + 2bsx), which yields, for N s 1 [18, 20] ln τ ' N [S(1) − S(x0 )] ' N sc(1 − x0 ).
(3)
Therefore, in the original PD model (no reflecting boundary), the CFP is given by the inverse of Eq. (3), φC (x0 ) ∼ e−N sc(1−x0 ) , which coincides with the exact result [4, 17]. We now add EN into the dynamics in the form of a fluctuating selection intensity in (1), s → s(t) = s0 + ξ(t),
10
φC
duced and 1 accounts for a baseline fitness contribution [2, 13, 14], while the population average fitness is f¯ = (n/N )fC + [1 − (n/N )]fD = 1 + s(b − c)(n/N ). In a finite population, fluctuations always drive the system to either the absorbing states n = 0 or n = N , and the noisy dynamics is described by the corresponding master equation [15]. For weak selection such that N s2 1, the latter can be approximated [14, 16] by a Fokker-Planck equation (FPE) for the probability P (x, t) of having cooperator density x at time t [3, 17]:
-2
10
-4
10
-6
0.1
0.2
0.3
0.4
x
0.5
0.6
0.7
FIG. 1: (Color online) Comparison between φC and ατ −1 for intermediate EN (× and 5 respectively) with σ/s0 = 0.3 and strong EN ( and 4 respectively) with σ/s0 = 1. The line is an eye guide showing the analytical results for φC ∼ τ −1 for IN only. Here, s0 = 0.01, b = 1.25, c = 1, N = 1500, and τc = 20. α was 35.3 for intermediate and 55.1 for strong EN.
where ξ is an Ornstein-Uhlenbeck (OU) process with 0 mean zero, variance hξ(t)ξ(t0 )i = σ 2 e−|t−t |/τc and correlation time τc > 0 [15]. The EN magnitude is assumed to satisfy σ . s0 1 so that s(t) can become negative when σ = O(s0 ). The OU’s Langevin equation (LE) reads p (4) ξ˙ = −ξ/τc + 2σ 2 /τc η1 (t), where η1 (t) is a white Gaussian noise with hη1 (t)i = 0 and hη1 (t)η1 (t0 )i = δ(t − t0 ) [21]. To account for the joint effects of IN and EN, Eq. (4) is coupled with the LE p (5) x˙ = A(x, ξ) + B(x, ξ)/N η2 (t), where η2 (t) is also white Gaussian noise, and A(x, ξ) = T + (x, ξ) − T − (x, ξ) and B(x, ξ) = T + (x, ξ) + T − (x, ξ) read for s 1: A(x, ξ) ' −x(1 − x)c(s0 + ξ)[1 − (b − c)(s0 + ξ)x], B(x, ξ) ' 2x(1 − x)[1 + c(s0 + ξ)(x − 1/2)].
(6)
We now proceed as in the absence of EN and compute the MFT τ of the modified PD model supplemented with a reflecting boundary at x0 . As with only IN, for N s0 1, φC (x0 ) and τ −1 in the original and modified models exhibit the same asymptotic dependence also in the presence of EN. Indeed, see Fig. 1, numerical simulations [22] with various EN strengths indicate that φC (x0 ) = ατ −1 , where α varies slowly with the model parameters. The description in terms of the coupled LEs (4,5) is equivalent to the following bivariate FPE for the probability P (x, ξ, t) to find cooperator density x and selection strength s = s0 + ξ at time t: ∂t P (x, ξ, t) = [−∂x A + ∂ξ ξ]P (x, ξ, t) + (2N )−1 [∂x2 B + (2V /τc )]P (x, ξ, t), (7) where we have defined V ≡ N σ 2 1. For N 1, we can use the semi-classical ansatz for the QSD
3
where we have defined the momenta px = ∂x S and pξ = ∂ξ S. The optimal path to fixation satisfies the HJE whose solutions obey the Hamilton equations (HE)
-2
10
-4
10
10
-6
x˙ = ∂px H = A + px B
10 10
(9)
10
-2
10
σ/s0
-1
10
-2
-4
-6
0
10
-3
10
-2
10
-1
10
σ/s0
0
FIG. 2: (Color online) φC versus relative EN strength σ/s0 in the short-correlated EN regime. The solid line is from Eq. (10) and the symbols are numerical simulations. Here, s0 = 0.01, b = 1.25, c = 1 and N = 2000, τc = 25, x0 = 0.25 in the left panel, and N = 1750, τc = 20, x0 = 0.1 in the right panel. The agreement improves from left to right as the inequalities N s20 1 and τc s0 1 are better satisfied. 10
0
10 10
φC
where the last line has been obtained by combining the HE for ξ˙ and p˙ξ into a single equation and by keeping terms up to O(px ) = O(s), see below. The solution to the HJE for generic EN (with arbitrary τc ) is found by solving numerically (9), R R yielding the action function S(x, ξ) = px (x, ξ)dx + pξ (x, ξ)dξ [23, 24]. Analytical progress can be done in two important limits: short-correlated (white) and long-correlated (adiabatic) EN. In the limit of short-correlated EN characterized by ¨ τc s−1 0 , ξ is negligible in the third of Eqs. (9) [25] yielding the effective noise strength ξ ' ξeff ' −2cV τc px x(1− x) [10, 26, 27]. Below we find px > 0, implying that ξeff < 0. As a result, the EN is exploited to enhance the CFP by decreasing the selection strength. We now assume the noise strength satisfies τc V s0 , in which case the EN markedly affects the dynamics, see below. Substituting ξeff into the first of Eqs. (9) one finds x˙ = −x(1 − x) cs0 − 2px (1 + c2 τc V x(1 − x) + O(s0 )) . The corresponding effective “white-noise” Hamiltonian is H(x, px ) ' −x(1 − x)px cs0 − px {1 + c2 τc V x(1 − x)} . Thus, the CFP of the original PD model is given by the accumulated action along px = cs0 /(1 + c2 τc V x(1 − x)): Z 1 cs0 du = (10) ln φC (x0 ) ' −N 2 x0 1 + c V τc u(1 − u) γ + 1 − 2x0 N s0 − ln 1 + c2 τc V (1 + γ)/2 , cτc V γ γ − 1 + 2x0 p where γ = 1 + 4/(c2 τc V ). In Fig. 2 we compare (10) with stochastic simulations [22] as a function of the relative EN strength σ/s0 and find a very good agreement. For weak EN, V τc 1, Eq. (10) becomes ln φC (x0 ) ' −N s0 c(1 − x0 ) 1 − (1/6)c2 V τc (1 − x0 )(2x0 + 1) , which coincides with the IN-only result to leading order. The most striking effect of the fluctuating environment appears in the limit of strong EN where γ → 1 [28]. In such a limit, the dependence of φC on N is no longer exponential but decays as a pure power-law: when τc V 1 and x0 > 0 in (10), one indeed finds −(s0 /σ2 )/(cτc ) φC ∼ N (σc)2 τc (1 − x0 )/x0 . This result is confirmed by stochastic simulations, see Fig. 3. The behavior of (10) in the limit x0 1 (small initial density of C’s) is also particularly relevant and yields ln φC ' −(2N s0 )/(cτc V γ) ln 1 + c2 τc V (1 + γ)/2 .(11)
-3
10
0
-6
-3
10
φC
p˙x = −∂x H = −px [∂x A + (px /2)∂x B] ξ¨ = ξ/τc2 − 2(V /τc )px ∂ξ A(x, ξ),
10
φC
H = px A(x, ξ) − pξ ξ/τc + (p2x /2)B(x, ξ) + (V /τc )p2ξ , (8)
10
φC
π(x, ξ) ∼ e−N S(x,ξ) in Eq. (7), which yields the HJE H(x, ξ, px , pξ ) = 0 with the Hamiltonian
-12
10 10
10
2
10
-6
-9
0
4
N
10
6
intrinsic noise theory σ/s0=0 extrinsic noise theory σ/s0=1,τC=30
2000
4000
N
6000
8000
10000
FIG. 3: (Color online) φC versus N under short-correlated EN: lines are theoretical results while ×’s/◦’s are simulation results with/without EN, see legend. Parameters are s0 = 0.01, b = 1.25, c = 1, and x0 = 0.25. IN-only results display exponential dependence on N , whereas for strong EN, V τc 1 (see text), φC exhibits a power-law dependence on N . Inset: φC versus N on log-log scale. Results of the theory and simulations are compared to N −3.4 (dotted line). The power-law φC ∼ N −10/3 predicted by (10) is approached when N → ∞.
In the strong EN limit (τc V 1, γ ' 1), Eq. (11) leads −2(s0 /σ2 )/(cτc ) to the power-law φC ' N (σc)2 τc . In the opposite long-correlated (adiabatic) EN limit, ξ(t) varies slowly and we can integrate over the fixation probability φC (ξ0 ) ∼ e−cN (s0 +ξ0 )(1−x0 ) with the Gaussian 2 2 weight e−ξ0 /(2σ ) for the noise ξ0 . Within a saddle-point approximation, we find an optimal noise magnitude at ξ0 = ξopt = −cV (1 − x0 ) so that φC (x0 ) becomes ln φC (x0 ) ' −N cs0 (1 − x0 ) [1 − (c/s0 )V (1 − x0 )] . (12) This result is valid when | ln φC | 1 which requires V < s0 , i.e. not too strong EN. Eq. (12) shows that adiabatic
4 EN can also significantly increase φC , but the dependence on N remains exponential. A different scenario arises when there is strong adiabatic EN. In this case the IN can be neglected and the fixation of C’s is solely governed by EN and thus the OU process (4) [27]. For completeness, we now briefly illustrate how the results obtained for fluctuating s can be directly used to study the case of fluctuating cost-to-benefit ratio r = c/b due to EN, according to r → r(t) = r0 +ξ(t). Here r0 < 1 with r0 ∼ O(1), and ξ is the same OU process as in (4) 0 with hξ(t)ξ(t0 )i = σr2 e−|t−t |/τc , where Vr ≡ N σr2 . In addition we assume σr r0 to guarantee 0 < r(t) < 1, and that b is fixed so that c(t) = br(t) fluctuates. The dynamics is thus described by Eqs. (4) and (5) with A(x, ξ) = −x(1 − x)bs(r0 + ξ){1 − bs[1 − (r0 + ξ)]x} and B(x, ξ) = 2x(1 − x)[1 + (r0 + ξ)bs(x − 1/2)]. As above, in the limit of weak selection N s2 1 and shortcorrelated EN τc s−1 , one finds the effective value for the EN from (9): ξeff ' −2bsVr τc px x(1 − x). Substituting ξ ' ξeff in the first HE (9) leads to an effective Hamiltonian whose accumulated action gives the CFP Z 1 sbr0 du . (13) ln φC (x0 ) ' −N 2 x0 1 + (sb) Vr τc u(1 − u) As before, we have assumed that the EN is sufficiently strong so that Vr τc s−1 [29]. Eq. (13) can be treated in the same manner as (10) and when EN is strong enough it also predicts a power-law dependence of φC on N [24]. Our approach generalizes earlier works in population genetics where the combined role of IN and EN was investigated by considering a fluctuating selection intensity, see e.g. [3, 6–8]. In these studies the dynamics was implemented with the Wright-Fisher model (WFM) with discrete time and non-overlapping generations [3]. In such a setting, a diffusion theory was devised in the weak selection limit to account for IN and time-uncorrelated (white) EN by averaging separately on the two sources of noise [3, 6–8]. When N σ 2 . N s20 1, the results of this approach coincide with (10) where τc = 1 and N → N/2 [30]. Yet, our approach [Eqs. (9)] is more general and not limited to white EN, since it allows to study EN with arbitrary correlation time [24]. We have studied the effect of environmental noise (EN) in a broad class of evolutionary processes not necessarily displaying metastability. Our approach relies on a suitable mapping onto a system exhibiting metastability and on its semiclassical treatment. Using the paradigmatic prisoner’s dilemma as an example, we have demonstrated how EN can drastically enhance the cooperators fixation probability (CFP). In particular, we have shown that for strong short-correlated EN the CFP decays algebraically (instead of exponentially) with the population size, as confirmed by stochastic simulations. Our results may provide further insight for resolving the contradiction between the observed examples of cooperation and predictions of many theoretical models.
[1] E. Pennisi, Science 309, 90 (2005). [2] J. Maynard Smith, Evolution and the Theory of Games (Cambridge University Press, Cambridge, 1982); M. A. Nowak, Evolutionary Dynamics (Belknap Press, 2006); G. Szabo and G. Fath, Phys. Rep. 446, 97 (2007); R. Axelrod, The Evolution of Cooperation (Basic Books, New York, 1984). [3] J. F. Crow and M. Kimura, An Introduction to Population Genetics Theory (Blackburn Press, New Jersey, 2009); W. J. Ewens, Mathematical Population Genetics (Springer, New York, 2004). [4] T. Antal and I. Scheuring, Bull. Math. Biol. 68, 1923 (2006). [5] W. D. Hamilton, J. Theor. Biol., 7, 1 (1964); R. L. Trivers, Quarterly Review of Biology 46, 35 (1971); R. L. Trivers, Quarterly Review of Biology 46, 35 (1971); M. A. Nowak and R. M. May, Nature 359, 826 (1992); M. A. Nowak and K. Sigmund, Nature 364, 56 (1993); R. Ferri`ere, Nature 393, 517 (1998); A. Traulsen and M. A. Nowak, Proc. Natl. Acad. Sci. USA 103, 10952 (2006); Z. Wang, A. Szolnoki, and M. Perc, Sci. Rep. 2, 369 (2012); M. Mobilia, Phys. Rev. E 86, 011134 (2012). [6] L. Jensen, Gen. Res., Camb. 21, 215 (1972). [7] L. Jensen and E. Pollak, J. Appl. Prob. 6, 19 (1969). [8] S. Karlin and B. Levikson, T. Pop. Biol. 74, 383 (1974). [9] E. G. Leigh, J. Theo. Biol. 90, 213 (1981); R. Lande, Am. Nat. 142, 911 (1993); K. Johst and C. Wissel, Theo. Pop. Biol. 52, 91 (1997). [10] A. Kamenev, B. Meerson, and B. Shklovskii, Phys. Rev. Lett. 101, 268103 (2008). [11] U. Dobramysl and U. C. T¨ auber, Phys. Rev. Lett. 110, 048105 (2013). [12] P. A. P. Moran, The statistical processes of evolutionary theory (Clarendon, Oxford, 1962). [13] M. A. Nowak, A. Sasaki, C. Taylor, and D. Fudenberg, Nature 428, 646 (2004); A. Traulsen and C. Hauert, in Reviews of Nonlinear Dynamics and Complexity, edited by H.-G. Shuster, Vol. 2 (Wiley-VCH, New York, 2010). [14] see, e.g., M. Mobilia and M. Assaf, EPL 91, 10002 (2010); M. Assaf and M. Mobilia, J. Stat. Mech., P09009 (2010); M. Assaf and M. Mobilia, J. Theor. Biol. 275, 93 (2011). [15] C. W. Gardiner, Handbook of Stochastic Methods, (Springer, New York, 2002). [16] M. Assaf and B. Meerson, Phys. Rev. Lett. 97, 200602 (2006); Phys. Rev. E 75, 031122 (2007). [17] A. Traulsen, J. C. Claussen, and C. Hauert, Phys. Rev. Lett. 95, 238701 (2005). [18] C. Escudero and A. Kamenev, Phys. Rev. E. 79, 041149 (2009); M. Assaf and B. Meerson, Phys. Rev. E. 81, 021116 (2010). [19] A. D. Wentzell and M. I. Freidlin, Russ. Math. Surveys 25, 1 (1970); M. I. Dykman, E. Mori, J. Ross, and P. M. Hunt, J. Chem. Phys. 100, 5735 (1994). [20] Here, the full calculation including subleading-order corrections yields τ ' [N cs2 x0 (1 − x0 )]−1 eN sc(1−x0 ) . [21] η1 (t) is the time-continuous limit (dt → 0) of a temporally uncorrelated normal random variable with mean 0 and variance 1/dt [15]. [22] To numerically study fixation in the population we use a kinetic Monte Carlo method based on the next-reaction variant of the Gillespie algorithm [31]. During a sin-
5 gle trajectory, the current number of cooperators n is stochastically updated using the birth/death transition propensities (rates) described in the text. EN is added by permitting the selection intensity parameter s to fluctuate. A pseudo-reaction fires at intervals much less than the EN correlation time and s is updated as if it had been following an OU process satisfying Eq. (4) using the method of [32]. At each update both the birth and death transition propensities are recalculated using the new value of s and the next reaction time is adjusted to account for the change. To calculate φC , many trajectories starting at n0 are run to fixation and the fraction in which cooperation results is calculated directly. For the MFT calculation, a reflecting boundary is placed at n0 and the time until the population fixates in the cooperation state is tracked. The mean time for many (≥1000) trajectories to reach fixation is then calculated. [23] D. M. Roma et. al., Phys. Rev. E. 71, 011902 (2005). [24] M. Assaf, M. Mobilia, and E. Roberts, in preparation. [25] Using the expression of ξeff together with Eqs. (6) and
(9), one a posteriori finds that ξ¨ = O(V τc s30 ) V s0 /τc which is indeed negligible in the third equation of (9). [26] E. Y. Levine and B. Meerson, Arxiv 1210.6436 (2012). [27] M. Assaf, E. Roberts, Z. Luthey-Schulten and N. Goldenfeld, Arxiv 1302.2724 (2013). [28] The expression (10) is well-defined also when V → 0 (γ → 1) and x0 → 0, as nh can be seen i by rewriting o ln φC (x0 ) ' −N cs0 2
1−γ 2 4γ
ln
1+
2 γ−1
γ+1−2x0 γ−1+2x0
.
[29] Since Vr s and V scale in the same manner in s, in the case of fluctuating r the EN is considered strong when Vr τc s−1 and weak when Vr τc . s−1 . [30] As time is discrete in the WFM and the EN between two successive generations is uncorrelated, the correspondence requires to set τc = 1 into (10). Furthermore, a population of N individuals in the WFM maps onto a population of size N/2 with the MM [3, 14]. [31] D. T. Gillespie, J. Comput. Phys. 22, 403 (1976); M. A. Gibson and J. Bruck, J. Phys. Chem. 104, 1876 (2000). [32] D. T. Gillespie, Phys. Rev. E 54, 2084 (1996).
arXiv:1305.6580v1 [q-bio.PE] 28 May 2013
2
1 Racah Institute of Physics, Hebrew University of Jerusalem, Jerusalem 91904, Israel Department of Applied Mathematics, School of Mathematics, University of Leeds, Leeds LS2 9JT, U.K. 3 Department of Biophysics, Johns Hopkins University, Baltimore, MD 21218, U.S.A.
We present a novel approach allowing the accurate study of rare events in a broad class of evolutionary processes under fluctuating environments. Our approach consists of mapping these processes, not necessarily displaying metastability, onto auxiliary models characterized by metastability that can be analyzed semiclassically. This enables us to study the interplay between extrinsic (environmental) noise and intrinsic (demographic) fluctuations on the statistics of interest. We illustrate our theory for the paradigmatic prisoner’s dilemma game whose evolution is described by the probability that cooperators take over and fixate the entire population. We demonstrate both analytically and by using Monte Carlo simulations that extrinsic noise may drastically enhance such a fixation probability. In particular, we show that under strong, short-correlated extrinsic noise, the cooperation fixation probability decays algebraically with the population size, rather than exponentially. Finally, we show how our approach generalizes earlier works, notably in population genetics. PACS numbers: 05.40.-a, 02.50.Ey, 87.23.Kg, 02.50.Le
Understanding the origin of cooperative behavior and how it is influenced by the population’s intrinsic properties and by environmental factors are major scientific puzzles [1] that are suitably investigated in the framework of evolutionary game theory (EGT) [2]. In EGT, successful species spread at the expense of others, and the optimization of the reproductive potential (fitness) at an individual level can result in the reduction of the population’s overall fitness, a phenomenon metaphorically captured by the prisoner’s dilemma (PD) game [1, 2]. While in EGT the dynamics is traditionally studied in terms of differential equations, demographic fluctuations – intrinsic noise (IN) – are known to affect the evolution in finite populations. In this case, the dynamics is often described by a Markov chain and characterized by the fixation probability of a given trait, which is the probability that the trait invades the entire population [3]. For the classic PD (with IN only), the cooperation fixation probability (CFP) vanishes exponentially with the population size, see e.g. [4], and defection prevails leading to a cooperation dilemma. This prediction, at odds with many observations, has motivated the investigation of various mechanisms that can promote cooperation [5]. In addition to IN, there often is also environmental (extrinsic) noise (EN) which results in random fluctuations of the interaction parameters. In theoretical population genetics [3, 6–8] and ecology [9–11], it has been shown that the combined effect of IN and EN can significantly affect the extinction properties of a population. While many earlier works focused on systems exhibiting metastability and/or short-correlated EN [3, 6–10], here we present a novel approach allowing one to analyze the combined influence of IN and EN, with arbitrary correlation time and magnitude, on a broad class of systems not necessarily displaying metastability. This is done by mapping the system at hand onto one possessing a long-
lived metastable state and treating the latter semiclassically. We illustrate our approach with the paradigmatic PD game, and find that EN drastically enhances the CFP and may even alter its functional dependence on the population size from an exponential to a power law. The paradigm of social dilemma is provided by the classic PD, whose main features are captured by assuming that the pairwise interaction between cooperators and defectors is described in terms of the benefit b and cost c of cooperation, with b > c > 0 [2]. Here, mutual cooperation leads to a payoff b−c > 0 and mutual defection gives a payoff 0, while when one player defects and the other cooperates, the former gets a payoff b and the latter gets −c. The quantity r ≡ c/b is the cost-to-benefit ratio [2] and the dilemma arises from the fact that each individual is better off defecting, even though r < 1 and mutual cooperation enhances the population overall payoff. Here, we consider a finite population comprising N individuals on a complete graph with n cooperators and N − n defectors. At mean field level (N → ∞), defection prevails and the fraction x ≡ n/N of cooperators in the population evolves to extinction according to the replicator rate equation (RE) (d/dt)x ≡ x˙ ∝ x(1−x) [ΠC (x) − ΠD (x)]. Here, ΠC = bx−c and ΠD = bx are the cooperator and defector average payoffs, respectively [2], and we henceforth assume that b, c = O(1). When the population size is finite, fluctuations alter the mean field predictions and the dynamics is usually described in terms of the Moran model (MM) [2, 4, 12, 13], according to a continuous-time Markov chain with the birth/death transition rates Tn+ = [fC (n)/f¯(n)]n(N − n)/N 2 and Tn− = [fD (n)/f¯(n)]n(N − n)/N 2 . Here, the cooperators/defectors fitnesses are fC (n) = 1 + s [b(n/N ) − c] and fD (n) = 1 + sbn/N, (1) where the selection intensity s > 0 has been intro-
2
∂t P (x, t) = −∂x [A(x)P (x, t)]+1/(2N )∂x2 [B(x)P (x, t)], (2) with A(x) = T + (x)−T − (x), B(x) = T + (x)+T − (x), and T ± (x) = Tn± . Note, that for N → ∞, Eq. (2) approaches the RE x˙ = A(x) = sx(1−x)[ΠC (x)−ΠD (x)]/f¯(x) [2, 17]. An important notion to characterize evolutionary dynamics is the fixation probability φC (x0 ) that cooperation prevails starting from an initial fraction x0 of cooperators. In the absence of environmental fluctuations, φC can be calculated exactly [4, 15]. Yet, it is useful to proceed with an asymptotic calculation (N 1) of φC (x0 ) using an auxiliary problem. For this, we consider the modified model obtained by supplementing the original PD model with a reflecting boundary at n0 = N x0 by im− = 0. Hence, the only absorbing state of the posing Tn=n 0 modified model is the state n = N and as a result, since x˙ < 0, a quasi-stationary distribution (QSD) peaked at x0 (for any value of x0 ) forms after a relaxation time of order O(s−1 ). This QSD slowly decays due to a slow leakage of probability into the absorbing state at x = 1, with a rate given by the inverse of the cooperation mean fixation time (MFT). Therefore, after a long time growing exponentially with N , cooperators eventually take over. Importantly, we now demonstrate that φC (x0 ) in the original PD problem coincides to leading order with the inverse of the MFT in the modified problem. The MFT in the modified problem can be computed by employing the metastable ansatz P (x, t) ' π(x)e−t/τ , where π(x) is the QSD and τ is the MFT. Then, Substituting the WKB ansatz π(x) ∼ e−N S(x) in Eq. (2) yields a stationary Hamilton-Jacobi equation (HJE) H(x, px ) = 0, with the Hamiltonian H(x, px ) = px A(x)+(p2x /2)B(x). In this formalism S(x) is called the action function while p(x) ≡ S 0 (x) is called the momentum [18, 19]. Fixation occurs along the zero-energy trajectory (“optimal path”) px (x) = −2A(x)/B(x), where px (x) ∼ O(s) R x 10 in 0the weak-selection limit. This gives S(x) = px (x )dx = (c/b) ln(2 − cs + 2bsx), which yields, for N s 1 [18, 20] ln τ ' N [S(1) − S(x0 )] ' N sc(1 − x0 ).
(3)
Therefore, in the original PD model (no reflecting boundary), the CFP is given by the inverse of Eq. (3), φC (x0 ) ∼ e−N sc(1−x0 ) , which coincides with the exact result [4, 17]. We now add EN into the dynamics in the form of a fluctuating selection intensity in (1), s → s(t) = s0 + ξ(t),
10
φC
duced and 1 accounts for a baseline fitness contribution [2, 13, 14], while the population average fitness is f¯ = (n/N )fC + [1 − (n/N )]fD = 1 + s(b − c)(n/N ). In a finite population, fluctuations always drive the system to either the absorbing states n = 0 or n = N , and the noisy dynamics is described by the corresponding master equation [15]. For weak selection such that N s2 1, the latter can be approximated [14, 16] by a Fokker-Planck equation (FPE) for the probability P (x, t) of having cooperator density x at time t [3, 17]:
-2
10
-4
10
-6
0.1
0.2
0.3
0.4
x
0.5
0.6
0.7
FIG. 1: (Color online) Comparison between φC and ατ −1 for intermediate EN (× and 5 respectively) with σ/s0 = 0.3 and strong EN ( and 4 respectively) with σ/s0 = 1. The line is an eye guide showing the analytical results for φC ∼ τ −1 for IN only. Here, s0 = 0.01, b = 1.25, c = 1, N = 1500, and τc = 20. α was 35.3 for intermediate and 55.1 for strong EN.
where ξ is an Ornstein-Uhlenbeck (OU) process with 0 mean zero, variance hξ(t)ξ(t0 )i = σ 2 e−|t−t |/τc and correlation time τc > 0 [15]. The EN magnitude is assumed to satisfy σ . s0 1 so that s(t) can become negative when σ = O(s0 ). The OU’s Langevin equation (LE) reads p (4) ξ˙ = −ξ/τc + 2σ 2 /τc η1 (t), where η1 (t) is a white Gaussian noise with hη1 (t)i = 0 and hη1 (t)η1 (t0 )i = δ(t − t0 ) [21]. To account for the joint effects of IN and EN, Eq. (4) is coupled with the LE p (5) x˙ = A(x, ξ) + B(x, ξ)/N η2 (t), where η2 (t) is also white Gaussian noise, and A(x, ξ) = T + (x, ξ) − T − (x, ξ) and B(x, ξ) = T + (x, ξ) + T − (x, ξ) read for s 1: A(x, ξ) ' −x(1 − x)c(s0 + ξ)[1 − (b − c)(s0 + ξ)x], B(x, ξ) ' 2x(1 − x)[1 + c(s0 + ξ)(x − 1/2)].
(6)
We now proceed as in the absence of EN and compute the MFT τ of the modified PD model supplemented with a reflecting boundary at x0 . As with only IN, for N s0 1, φC (x0 ) and τ −1 in the original and modified models exhibit the same asymptotic dependence also in the presence of EN. Indeed, see Fig. 1, numerical simulations [22] with various EN strengths indicate that φC (x0 ) = ατ −1 , where α varies slowly with the model parameters. The description in terms of the coupled LEs (4,5) is equivalent to the following bivariate FPE for the probability P (x, ξ, t) to find cooperator density x and selection strength s = s0 + ξ at time t: ∂t P (x, ξ, t) = [−∂x A + ∂ξ ξ]P (x, ξ, t) + (2N )−1 [∂x2 B + (2V /τc )]P (x, ξ, t), (7) where we have defined V ≡ N σ 2 1. For N 1, we can use the semi-classical ansatz for the QSD
3
where we have defined the momenta px = ∂x S and pξ = ∂ξ S. The optimal path to fixation satisfies the HJE whose solutions obey the Hamilton equations (HE)
-2
10
-4
10
10
-6
x˙ = ∂px H = A + px B
10 10
(9)
10
-2
10
σ/s0
-1
10
-2
-4
-6
0
10
-3
10
-2
10
-1
10
σ/s0
0
FIG. 2: (Color online) φC versus relative EN strength σ/s0 in the short-correlated EN regime. The solid line is from Eq. (10) and the symbols are numerical simulations. Here, s0 = 0.01, b = 1.25, c = 1 and N = 2000, τc = 25, x0 = 0.25 in the left panel, and N = 1750, τc = 20, x0 = 0.1 in the right panel. The agreement improves from left to right as the inequalities N s20 1 and τc s0 1 are better satisfied. 10
0
10 10
φC
where the last line has been obtained by combining the HE for ξ˙ and p˙ξ into a single equation and by keeping terms up to O(px ) = O(s), see below. The solution to the HJE for generic EN (with arbitrary τc ) is found by solving numerically (9), R R yielding the action function S(x, ξ) = px (x, ξ)dx + pξ (x, ξ)dξ [23, 24]. Analytical progress can be done in two important limits: short-correlated (white) and long-correlated (adiabatic) EN. In the limit of short-correlated EN characterized by ¨ τc s−1 0 , ξ is negligible in the third of Eqs. (9) [25] yielding the effective noise strength ξ ' ξeff ' −2cV τc px x(1− x) [10, 26, 27]. Below we find px > 0, implying that ξeff < 0. As a result, the EN is exploited to enhance the CFP by decreasing the selection strength. We now assume the noise strength satisfies τc V s0 , in which case the EN markedly affects the dynamics, see below. Substituting ξeff into the first of Eqs. (9) one finds x˙ = −x(1 − x) cs0 − 2px (1 + c2 τc V x(1 − x) + O(s0 )) . The corresponding effective “white-noise” Hamiltonian is H(x, px ) ' −x(1 − x)px cs0 − px {1 + c2 τc V x(1 − x)} . Thus, the CFP of the original PD model is given by the accumulated action along px = cs0 /(1 + c2 τc V x(1 − x)): Z 1 cs0 du = (10) ln φC (x0 ) ' −N 2 x0 1 + c V τc u(1 − u) γ + 1 − 2x0 N s0 − ln 1 + c2 τc V (1 + γ)/2 , cτc V γ γ − 1 + 2x0 p where γ = 1 + 4/(c2 τc V ). In Fig. 2 we compare (10) with stochastic simulations [22] as a function of the relative EN strength σ/s0 and find a very good agreement. For weak EN, V τc 1, Eq. (10) becomes ln φC (x0 ) ' −N s0 c(1 − x0 ) 1 − (1/6)c2 V τc (1 − x0 )(2x0 + 1) , which coincides with the IN-only result to leading order. The most striking effect of the fluctuating environment appears in the limit of strong EN where γ → 1 [28]. In such a limit, the dependence of φC on N is no longer exponential but decays as a pure power-law: when τc V 1 and x0 > 0 in (10), one indeed finds −(s0 /σ2 )/(cτc ) φC ∼ N (σc)2 τc (1 − x0 )/x0 . This result is confirmed by stochastic simulations, see Fig. 3. The behavior of (10) in the limit x0 1 (small initial density of C’s) is also particularly relevant and yields ln φC ' −(2N s0 )/(cτc V γ) ln 1 + c2 τc V (1 + γ)/2 .(11)
-3
10
0
-6
-3
10
φC
p˙x = −∂x H = −px [∂x A + (px /2)∂x B] ξ¨ = ξ/τc2 − 2(V /τc )px ∂ξ A(x, ξ),
10
φC
H = px A(x, ξ) − pξ ξ/τc + (p2x /2)B(x, ξ) + (V /τc )p2ξ , (8)
10
φC
π(x, ξ) ∼ e−N S(x,ξ) in Eq. (7), which yields the HJE H(x, ξ, px , pξ ) = 0 with the Hamiltonian
-12
10 10
10
2
10
-6
-9
0
4
N
10
6
intrinsic noise theory σ/s0=0 extrinsic noise theory σ/s0=1,τC=30
2000
4000
N
6000
8000
10000
FIG. 3: (Color online) φC versus N under short-correlated EN: lines are theoretical results while ×’s/◦’s are simulation results with/without EN, see legend. Parameters are s0 = 0.01, b = 1.25, c = 1, and x0 = 0.25. IN-only results display exponential dependence on N , whereas for strong EN, V τc 1 (see text), φC exhibits a power-law dependence on N . Inset: φC versus N on log-log scale. Results of the theory and simulations are compared to N −3.4 (dotted line). The power-law φC ∼ N −10/3 predicted by (10) is approached when N → ∞.
In the strong EN limit (τc V 1, γ ' 1), Eq. (11) leads −2(s0 /σ2 )/(cτc ) to the power-law φC ' N (σc)2 τc . In the opposite long-correlated (adiabatic) EN limit, ξ(t) varies slowly and we can integrate over the fixation probability φC (ξ0 ) ∼ e−cN (s0 +ξ0 )(1−x0 ) with the Gaussian 2 2 weight e−ξ0 /(2σ ) for the noise ξ0 . Within a saddle-point approximation, we find an optimal noise magnitude at ξ0 = ξopt = −cV (1 − x0 ) so that φC (x0 ) becomes ln φC (x0 ) ' −N cs0 (1 − x0 ) [1 − (c/s0 )V (1 − x0 )] . (12) This result is valid when | ln φC | 1 which requires V < s0 , i.e. not too strong EN. Eq. (12) shows that adiabatic
4 EN can also significantly increase φC , but the dependence on N remains exponential. A different scenario arises when there is strong adiabatic EN. In this case the IN can be neglected and the fixation of C’s is solely governed by EN and thus the OU process (4) [27]. For completeness, we now briefly illustrate how the results obtained for fluctuating s can be directly used to study the case of fluctuating cost-to-benefit ratio r = c/b due to EN, according to r → r(t) = r0 +ξ(t). Here r0 < 1 with r0 ∼ O(1), and ξ is the same OU process as in (4) 0 with hξ(t)ξ(t0 )i = σr2 e−|t−t |/τc , where Vr ≡ N σr2 . In addition we assume σr r0 to guarantee 0 < r(t) < 1, and that b is fixed so that c(t) = br(t) fluctuates. The dynamics is thus described by Eqs. (4) and (5) with A(x, ξ) = −x(1 − x)bs(r0 + ξ){1 − bs[1 − (r0 + ξ)]x} and B(x, ξ) = 2x(1 − x)[1 + (r0 + ξ)bs(x − 1/2)]. As above, in the limit of weak selection N s2 1 and shortcorrelated EN τc s−1 , one finds the effective value for the EN from (9): ξeff ' −2bsVr τc px x(1 − x). Substituting ξ ' ξeff in the first HE (9) leads to an effective Hamiltonian whose accumulated action gives the CFP Z 1 sbr0 du . (13) ln φC (x0 ) ' −N 2 x0 1 + (sb) Vr τc u(1 − u) As before, we have assumed that the EN is sufficiently strong so that Vr τc s−1 [29]. Eq. (13) can be treated in the same manner as (10) and when EN is strong enough it also predicts a power-law dependence of φC on N [24]. Our approach generalizes earlier works in population genetics where the combined role of IN and EN was investigated by considering a fluctuating selection intensity, see e.g. [3, 6–8]. In these studies the dynamics was implemented with the Wright-Fisher model (WFM) with discrete time and non-overlapping generations [3]. In such a setting, a diffusion theory was devised in the weak selection limit to account for IN and time-uncorrelated (white) EN by averaging separately on the two sources of noise [3, 6–8]. When N σ 2 . N s20 1, the results of this approach coincide with (10) where τc = 1 and N → N/2 [30]. Yet, our approach [Eqs. (9)] is more general and not limited to white EN, since it allows to study EN with arbitrary correlation time [24]. We have studied the effect of environmental noise (EN) in a broad class of evolutionary processes not necessarily displaying metastability. Our approach relies on a suitable mapping onto a system exhibiting metastability and on its semiclassical treatment. Using the paradigmatic prisoner’s dilemma as an example, we have demonstrated how EN can drastically enhance the cooperators fixation probability (CFP). In particular, we have shown that for strong short-correlated EN the CFP decays algebraically (instead of exponentially) with the population size, as confirmed by stochastic simulations. Our results may provide further insight for resolving the contradiction between the observed examples of cooperation and predictions of many theoretical models.
[1] E. Pennisi, Science 309, 90 (2005). [2] J. Maynard Smith, Evolution and the Theory of Games (Cambridge University Press, Cambridge, 1982); M. A. Nowak, Evolutionary Dynamics (Belknap Press, 2006); G. Szabo and G. Fath, Phys. Rep. 446, 97 (2007); R. Axelrod, The Evolution of Cooperation (Basic Books, New York, 1984). [3] J. F. Crow and M. Kimura, An Introduction to Population Genetics Theory (Blackburn Press, New Jersey, 2009); W. J. Ewens, Mathematical Population Genetics (Springer, New York, 2004). [4] T. Antal and I. Scheuring, Bull. Math. Biol. 68, 1923 (2006). [5] W. D. Hamilton, J. Theor. Biol., 7, 1 (1964); R. L. Trivers, Quarterly Review of Biology 46, 35 (1971); R. L. Trivers, Quarterly Review of Biology 46, 35 (1971); M. A. Nowak and R. M. May, Nature 359, 826 (1992); M. A. Nowak and K. Sigmund, Nature 364, 56 (1993); R. Ferri`ere, Nature 393, 517 (1998); A. Traulsen and M. A. Nowak, Proc. Natl. Acad. Sci. USA 103, 10952 (2006); Z. Wang, A. Szolnoki, and M. Perc, Sci. Rep. 2, 369 (2012); M. Mobilia, Phys. Rev. E 86, 011134 (2012). [6] L. Jensen, Gen. Res., Camb. 21, 215 (1972). [7] L. Jensen and E. Pollak, J. Appl. Prob. 6, 19 (1969). [8] S. Karlin and B. Levikson, T. Pop. Biol. 74, 383 (1974). [9] E. G. Leigh, J. Theo. Biol. 90, 213 (1981); R. Lande, Am. Nat. 142, 911 (1993); K. Johst and C. Wissel, Theo. Pop. Biol. 52, 91 (1997). [10] A. Kamenev, B. Meerson, and B. Shklovskii, Phys. Rev. Lett. 101, 268103 (2008). [11] U. Dobramysl and U. C. T¨ auber, Phys. Rev. Lett. 110, 048105 (2013). [12] P. A. P. Moran, The statistical processes of evolutionary theory (Clarendon, Oxford, 1962). [13] M. A. Nowak, A. Sasaki, C. Taylor, and D. Fudenberg, Nature 428, 646 (2004); A. Traulsen and C. Hauert, in Reviews of Nonlinear Dynamics and Complexity, edited by H.-G. Shuster, Vol. 2 (Wiley-VCH, New York, 2010). [14] see, e.g., M. Mobilia and M. Assaf, EPL 91, 10002 (2010); M. Assaf and M. Mobilia, J. Stat. Mech., P09009 (2010); M. Assaf and M. Mobilia, J. Theor. Biol. 275, 93 (2011). [15] C. W. Gardiner, Handbook of Stochastic Methods, (Springer, New York, 2002). [16] M. Assaf and B. Meerson, Phys. Rev. Lett. 97, 200602 (2006); Phys. Rev. E 75, 031122 (2007). [17] A. Traulsen, J. C. Claussen, and C. Hauert, Phys. Rev. Lett. 95, 238701 (2005). [18] C. Escudero and A. Kamenev, Phys. Rev. E. 79, 041149 (2009); M. Assaf and B. Meerson, Phys. Rev. E. 81, 021116 (2010). [19] A. D. Wentzell and M. I. Freidlin, Russ. Math. Surveys 25, 1 (1970); M. I. Dykman, E. Mori, J. Ross, and P. M. Hunt, J. Chem. Phys. 100, 5735 (1994). [20] Here, the full calculation including subleading-order corrections yields τ ' [N cs2 x0 (1 − x0 )]−1 eN sc(1−x0 ) . [21] η1 (t) is the time-continuous limit (dt → 0) of a temporally uncorrelated normal random variable with mean 0 and variance 1/dt [15]. [22] To numerically study fixation in the population we use a kinetic Monte Carlo method based on the next-reaction variant of the Gillespie algorithm [31]. During a sin-
5 gle trajectory, the current number of cooperators n is stochastically updated using the birth/death transition propensities (rates) described in the text. EN is added by permitting the selection intensity parameter s to fluctuate. A pseudo-reaction fires at intervals much less than the EN correlation time and s is updated as if it had been following an OU process satisfying Eq. (4) using the method of [32]. At each update both the birth and death transition propensities are recalculated using the new value of s and the next reaction time is adjusted to account for the change. To calculate φC , many trajectories starting at n0 are run to fixation and the fraction in which cooperation results is calculated directly. For the MFT calculation, a reflecting boundary is placed at n0 and the time until the population fixates in the cooperation state is tracked. The mean time for many (≥1000) trajectories to reach fixation is then calculated. [23] D. M. Roma et. al., Phys. Rev. E. 71, 011902 (2005). [24] M. Assaf, M. Mobilia, and E. Roberts, in preparation. [25] Using the expression of ξeff together with Eqs. (6) and
(9), one a posteriori finds that ξ¨ = O(V τc s30 ) V s0 /τc which is indeed negligible in the third equation of (9). [26] E. Y. Levine and B. Meerson, Arxiv 1210.6436 (2012). [27] M. Assaf, E. Roberts, Z. Luthey-Schulten and N. Goldenfeld, Arxiv 1302.2724 (2013). [28] The expression (10) is well-defined also when V → 0 (γ → 1) and x0 → 0, as nh can be seen i by rewriting o ln φC (x0 ) ' −N cs0 2
1−γ 2 4γ
ln
1+
2 γ−1
γ+1−2x0 γ−1+2x0
.
[29] Since Vr s and V scale in the same manner in s, in the case of fluctuating r the EN is considered strong when Vr τc s−1 and weak when Vr τc . s−1 . [30] As time is discrete in the WFM and the EN between two successive generations is uncorrelated, the correspondence requires to set τc = 1 into (10). Furthermore, a population of N individuals in the WFM maps onto a population of size N/2 with the MM [3, 14]. [31] D. T. Gillespie, J. Comput. Phys. 22, 403 (1976); M. A. Gibson and J. Bruck, J. Phys. Chem. 104, 1876 (2000). [32] D. T. Gillespie, Phys. Rev. E 54, 2084 (1996).
