# Cooperation dilemma in finite populations under fluctuating environments

###### Abstract

We present a novel approach allowing the study of rare events like fixation under fluctuating environments, modeled as extrinsic noise, in evolutionary processes characterized by the dominance of one species. Our treatment consists of mapping the system onto an auxiliary model, exhibiting metastable species coexistence, that can be analyzed semiclassically. This approach enables us to study the interplay between extrinsic and demographic noise on the statistics of interest. We illustrate our theory by considering the paradigmatic prisoner’s dilemma game whose evolution is described by the probability that cooperators fixate the population and replace all defectors. We analytically and numerically demonstrate that extrinsic noise may drastically enhance the cooperation fixation probability and even change its functional dependence on the population size. These results, which generalize earlier works in population genetics, indicate that extrinsic noise may help sustain and promote a much higher level of cooperation than static settings.

###### pacs:

05.40.-a, 02.50.Ey, 87.23.Kg, 02.50.LeUnderstanding the origin of cooperative behavior and how it is influenced by the population’s intrinsic properties and by environmental factors are major scientific puzzles Pennisi () that are suitably investigated in the framework of evolutionary game theory (EGT) EGT (). In EGT, successful species with a high reproductive potential (fitness) spread and the optimization of the fitness at an individual level can result in the reduction of the population’s overall fitness, a phenomenon suggestively captured by the prisoner’s dilemma (PD) game Pennisi (); EGT (). 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 (or “pure strategy”), which is the probability that the trait invades the entire population Kimura (). For the classic PD (with IN only), the cooperation fixation probability (CFP) vanishes exponentially with the population size, see e.g. Antal06 (), and defection prevails leading to a cooperation dilemma. This prediction, at odds with many experimental observations, has motivated the investigation of various mechanisms that can promote and sustain cooperation various ().

Besides IN, an important source of fluctuations in such systems is extrinsic noise (EN) mostly due to the inherent environmental fluctuations, and from being coupled to other fluctuating systems. Such EN can be aptly modeled in the form of random fluctuations in one or more interaction parameters. In theoretical population genetics Kimura (); Jensen69 (); Jensen72 (); Karlin74 (), ecology Leigh (); KMS08 (); DT13 () and cellular biology AR2013 (), it has been shown that the combined effect of IN and EN can significantly affect the lifetime of the long-lived metastable coexistence state the system dwells in prior to escape. In this work, we go beyond these and other works that focused on systems exhibiting metastability, and present a novel approach that allows us to analyze the combined influence of IN and EN, with arbitrary correlation time, magnitude and statistics, in systems characterized by the dominance of one species instead of metastability. This is done by a suitable mapping onto an auxiliary model possessing a long-lived metastable state and treating the latter semiclassically. We illustrate our approach on the prototypical example of the PD game. We show that EN can drastically enhance the CFP and may even change its functional dependence on the population size. These results may be interpreted as the evolutionary signature of noisy environments on population diversity bethedging ().

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 and cost of cooperation, with EGT (). Here, mutual cooperation leads to a payoff and mutual defection gives a payoff , while when one player defects and the other cooperates, the former gets a payoff and the latter gets . The quantity is the cost-to-benefit ratio EGT () and the dilemma arises from the fact that, while and mutual cooperation enhances the population overall payoff, each individual is better off defecting.

We consider a finite and well-mixed population of size , with cooperators and defectors. At mean field level (), defection always prevails and the fraction of cooperators evolves to extinction, , according to the replicator rate equation . and are the cooperator and defector average payoffs, respectively EGT (), and we assume that .

When the population size is finite, demographic fluctuations always drive the system to either the absorbing states or , and the stochastic dynamics is described by the master equation , where and are the respective birth and death rates. As often, these are given in terms of the Moran model Moran (); EGT (); weaksel (); Antal06 (): and , where the cooperators/defectors fitnesses are

(1) |

and the population average fitness is . In Eqs. (1) the term accounts for a baseline fitness contribution and the selection strength is denoted by EGT (); weaksel (); MA10 (). While our approach applies to arbitrary selection strength, throughout the paper we focus on the biologically relevant limit of weak selection, Kimura (); weaksel (), which ensures that in Eqs. (1).

Furthermore, it is convenient to work in the regime where . In this regime, one can accurately approximate the master equation by a Fokker-Planck equation (FPE) MA10 (); Assaf06 () for the probability of having cooperator density at time Kimura (); TCH ():

(2) |

where , giving a relaxation time , , and .

An important notion to characterize evolutionary dynamics is the CFP – the probability that cooperation fixates starting from an initial fraction of cooperators. In the absence of EN, can be calculated exactly Gardiner (); Antal06 (), and one finds in the leading exponential order . Here, we purposely adopt another route and show how to asymptotically calculate via an auxiliary problem. For this, we consider the modified model obtained by supplementing the original PD system with a reflecting boundary at by imposing . Hence, the only absorbing state of the modified model is the state . Therefore, as for any , a quasi-stationary distribution (QSD) peaked at (for any value of ) forms after an relaxation time. This metastable state, however, slowly decays due to a slow leakage of probability into the absorbing state at , with a rate given by the inverse of the cooperation mean fixation time (MFT).

Employing the metastable ansatz in Eq. (2), where is the QSD, the MFT of the auxiliary model can be computed using the semiclassical ansatz, . Here is called the action function, while is the momentum Dykman1994 (); AM10 (). This yields a stationary Hamilton-Jacobi equation, , with Hamiltonian . Fixation occurs along the zero-energy trajectory , where . This gives , from which the QSD at , , is found. Since , we have AM10 (); subleading ()

(3) |

where this result is valid when , which ensures a long-lived metastable state MA10 (). Importantly, we find that for the MFT of the modified problem (3) coincides to leading order with the inverse of the CFP in the original PD model Antal06 (); TCH (). We now use this finding to study the CFP in the presence of EN.

To this end, we incorporate EN in the form of one or more fluctuating parameters. For concreteness we take a fluctuating selection strength, . By directly affecting the fitness of individuals, this choice is particularly relevant in population genetics Fisher47 (); Kimura (); Jensen72 (); Jensen69 (); Karlin74 (), ecology ecol () and cellular biology bethedging (). Here, is taken as an Ornstein-Uhlenbeck (OU) process with mean zero, variance and correlation time Gardiner (); OU (). We assume that is arbitrary so that can become negative for . The OU process satisfies the following Langevin equation

(4) |

where is white Gaussian noise eta1 ().

We now proceed as in the absence of EN and compute of the modified PD model supplemented with a reflecting boundary at . We have numerically confirmed (see SM SM () for details) that for , and exhibit the same asymptotic behavior in the original and modified models also in the presence of EN, see Fig. 1.

To account for the joint effects of IN and EN, we couple Eq. (4) with FPE (2), and arrive at the following bivariate FPE for the probability to find cooperator density and selection strength at time :

(5) | |||||

Here, and read for :

(6) |

and we have defined . For , we can use the semi-classical ansatz for the QSD in Eq. (5), which yields the Hamilton-Jacobi equation with Hamiltonian

(7) |

where we have defined and . The corresponding Hamilton equations are

(8) |

where the third equation has been obtained by combining the equations for and into a single equation and by keeping terms up to , see below. The solution to the Hamilton-Jacobi equation for generic EN (with arbitrary ) is found by solving numerically Eqs. (Cooperation dilemma in finite populations under fluctuating environments), yielding the action function Roma (). Here, we focus on two important and analytically amenable regimes: short-correlated (white) EN, when , and long-correlated (adiabatic) EN, when .

For short-correlated EN, , we find that is negligible in the third of Eqs. (Cooperation dilemma in finite populations under fluctuating environments) ddotxi () yielding the effective noise strength KMS08 (); Meerson2012 (); AR2013 (). Since , see below, , thus EN is exploited to enhance the CFP by decreasing the selection strength.

Substituting into the first of Eqs. (Cooperation dilemma in finite populations under fluctuating environments) one finds . It appears that EN markedly affects the dynamics when its magnitude satisfies . In this regime the corresponding effective white-noise Hamiltonian is . Solving , we find . This yields the MFT in the modified model, and therefore, the CFP of the original PD model:

(9) | |||

where . In Fig. 2 we compare Eq. (9) with numerical simulations as a function of the relative EN strength and find a very good agreement for both (left panel) and (right panel). One can clearly see that EN, by effectively decreasing the selection strength , enhances the CFP compared to the IN-only case with (see also Fig. 3 and Fig. S2 where we respectively plot the CFP versus and ).

For a given short-correlated EN, , there are two interesting limits to (9): (i) strong and (ii) weak EN. (i) The most striking effect of EN appears in the limit of strong EN, , which yields . Here, for finite values of , the dependence of on becomes a power-law, and Eq. (9) gives way to . This result is confirmed by numerical simulations, see Fig. 3.

(ii) For weak EN, , Eq. (9) can be approximated as , which coincides with the IN-only result to leading order.

The behavior of Eq. (9) for a small initial density of C’s () is particularly relevant in EGT EGT (). In this case, for arbitrary EN strength and , the CFP is

(10) |

Again, for strong EN, , Eq. (10) becomes a power-law other-form ().

Note that, while Eq. (9) has been formally derived in the regime , its predictions also hold when with , as illustrated by the numerical results in Fig. 3. This is because the leading correction to due to EN is independent of when [see the denominator of the integrand of (9)]. Thus, our results due to EN are applicable as long as , and are expected to hold also in the non-FPE regime where SM ().

The case of long-correlated EN, , is investigated in the SM SM (). Here, for weak EN, , we find that . Under strong EN, , the intrinsic fluctuations are negligible AR2013 () and is solely governed by Eq. (4), yielding , see SM SM () for the details.

The various EN parameter regimes for fluctuating are summarized in a diagram, see Fig. S1 in the SM SM ().

For completeness, we have also considered the case of external fluctuations in the cost-to-benefit ratio , with and [where ]. In this case, the dynamics of is given by (4) with , where . In addition we assume to guarantee , and that is fixed so that fluctuates. Performing the calculations along the same lines as for fluctuating , we find for short-correlated EN,

(11) |

Similarly as before, for strong EN , Eq. (11) also predicts that decays algebraically with .

Our approach generalizes earlier works in population genetics where the combined role of IN and EN was investigated by considering a fluctuating selection strength, see e.g. Fisher47 (); Jensen69 (); Jensen72 (); Karlin74 (); Kimura (). In these studies the dynamics was implemented with the Wright-Fisher model with discrete time and non-overlapping generations Kimura (). 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 Jensen69 (); Jensen72 (); Karlin74 (); Kimura (). When , the results of this approach coincide with Eq. (9) for and WFM (). Yet, our approach is more general, since it allows to study EN with arbitrary correlation time and statistics, as well as in the presence of frequency-dependent selection.

In this work, we have analyzed fixation in evolutionary processes characterized by the dominance of one species. Our approach relies on a semi-classical treatment applied to an auxiliary model exhibiting metastability. This allows to study how fixation is affected by the interplay between intrinsic and extrinsic noise (EN). Our theory is general in the sense that it can deal with EN of arbitrary statistics, correlation time and magnitude, with one or multiple fluctuating parameters, and can be also used for systems exhibiting metastable coexistence. Using the prototypical prisoner’s dilemma game we have shown that EN is exploited to effectively reduce the selection strength and thereby, to drastically enhance cooperation, whose fixation probability is otherwise vanishingly small. This indicates that EN may be vital in sustaining a certain level of cooperation and population diversity by effectively opposing single-type dominance, as reported in recent microbial experiments bethedging (). Therefore, EN may contribute to reconcile the theoretical predictions with observed examples of cooperative behaviors.

Supplemental Material for:

Cooperation dilemma in finite populations under fluctuating environments

In this supplemental material, we summarize in a schematic diagram the results obtained in the main text for the cooperation fixation probability (CFP) and discuss the various extrinsic noise (EN) parameter regimes. We also outline the derivation of the CFP under adiabatic EN. Finally, we briefly explain our simulation method.

## I Parameter regime diagram

In this section we map the results for the CFP in the various regimes of parameter space when the selection strength fluctuates, . Here, EGT (); weaksel () is the mean selection intensity and is the Ornstein-Uhlenbeck (OU) process [Eq. (4)] with mean zero, variance and correlation time .
When the selection strength fluctuates, in addition to the benefit and cost of cooperation, there are four essential parameters that control the system’s dynamics: the population size , , and the EN magnitude and correlation time . In order to present the results in a two-dimensional diagram, we fix and the relative EN magnitude , so that , and draw the schematic diagram of versus . As shown in Fig. S1, this diagram is characterized by distinct regimes which are discussed below.

## Ii Discussion of the various EN regimes

Weak short-correlated noise regime (I): Here, , and . In this regime the CFP is given by Eq. (9), and the dependence on is exponential.

Strong short-correlated noise regime (II): Here, , and . In this regime is also given by Eq. (9), but the dependence on becomes algebraic, see Fig. 3.

Short-correlated EN / non-Fokker-Planck regime (III): Here, and . In this regime the Fokker-Planck equation is generally not an accurate approximation of the underlying master equation. However, as argued in the main text, and as corroborated by the numerical simulations of Fig. 3 (see the large- results), the correction to the CFP due to EN predicted by our theory [given by Eq. (9)] also holds well into the non-Fokker-Planck regime of .

Weak long-correlated noise regime (IV): Here, and . In this regime is given by Eq. (S1), and the dependence on is exponential as explained in Section 3 of this supplemental material.

Strong long-correlated noise regime (V): Here, and . In this regime the dynamics is solely governed by the OU process (4) and the CFP satisfies to leading order as explained in Section 3 of this supplemental material (see also Fig. S2 and the main text).

Quasi-neutral regime (VI): Here, Kimura (). In this regime the dynamics is close to neutral, and the CFP scales as EGT (); weaksel ()

Strong-selection regime (VII): Here, (for the sake of illustration, this regime is shown as the shaded region of
in Fig. S1). This strong-selection regime, which is of marginal biological relevance Kimura (), can be treated by considering other expressions [than Eqs. (1)] for the fitnesses EGT (); weaksel (); Antal06 (). In addition, one needs to proceed with a direct analysis of the master equation (instead of the Fokker-Planck equation), as, e.g., in Ref. MA10 ().

## Iii CFP under long-correlated (adiabatic) EN

In this section we calculate the CFP under fluctuating selection strength in the case of long-correlated (adiabatic) EN, . Here, the selection strength fluctuates slowly and can be considered as almost constant while a rare fluctuation leads to the fixation of cooperators. As a result, the CFP can be found by integrating over the fixation probability given noise , , with the noise’s Gaussian weight . A saddle-point approximation gives the optimal noise strength yielding

(S1) |

This result is valid when which requires , i.e. not too strong EN (Regime IV in Fig. S1). In addition, since in the original PD model the fixation time is Antal06 (), the adiabatic regime holds provided that . Eq. (S1) shows that can be exponentially enhanced by adiabatic EN.

A different scenario arises under adiabatic noise of strong intensity (Regime V in Fig. S1): When the intrinsic noise is negligible and the CFP is solely governed by the OU process (4) AR2013 (). As a result, the mean fixation time (MFT) in the auxiliary model is determined by the mean first passage time (MFPT) it takes to reach the value starting from at . This MFPT, denoted by , is governed by the following equation Gardiner ()

where , and we assume absorbing and reflecting boundaries at and , respectively, such that . The solution of this equation is given by

where , is the generalized hypergeometric function, , and in the regime of interest .

We are interested in the MFPT to reach noise magnitude starting from . Once crosses , the selection pressure vanishes and the auxiliary model rapidly fixates (compared to the fixation time when ). As a result, the fixation time is governed by . For strong selection and for we have and . Thus, we find , and as a result, the CFP under strong adiabatic noise in the original problem (see main text) satisfies , which is confirmed by Fig. S2.

## Iv Stochastic simulations with EN

To study fixation we use a kinetic Monte Carlo method based on the next-reaction variant of the Gillespie algorithm GillespieGibsonBruck (). During a single trajectory, the current number of cooperators is stochastically updated using the birth/death transition rates, , described in the main text. EN is added by permitting the selection strength parameter to fluctuate. A pseudo-reaction fires at intervals much less than the EN correlation time and is updated as if it had been following an OU process satisfying Eq. (4) using the method of Gillespie1996ens (). To calculate , the fraction of many trajectories starting at and resulting in the cooperation state is calculated directly. For the MFT calculation, a reflecting boundary is placed at and the mean time for many (1000) trajectories to reach the cooperation state is then calculated.

## References

- (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. Szabó and G. Fáth, 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ère, 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 (1973).
- (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äuber, Phys. Rev. Lett. 110, 048105 (2013).
- (12) M. Assaf, E. Roberts, Z. Luthey-Schulten and N. Goldenfeld, Phys. Rev. Lett. 111, 058102 (2013).
- (13) H. J. E. Beaumont et. al., Nature 462, 90 (2009).
- (14) P. A. P. Moran, The statistical processes of evolutionary theory (Clarendon, Oxford, 1962).
- (15) 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).
- (16) 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).
- (17) M. Assaf and B. Meerson, Phys. Rev. Lett. 97, 200602 (2006); Phys. Rev. E 75, 031122 (2007).
- (18) A. Traulsen, J. C. Claussen, and C. Hauert, Phys. Rev. Lett. 95, 238701 (2005).
- (19) C. W. Gardiner, Handbook of Stochastic Methods, (Springer, New York, 2002).
- (20) C. Escudero and A. Kamenev, Phys. Rev. E. 79, 041149 (2009); M. Assaf and B. Meerson, Phys. Rev. E. 81, 021116 (2010).
- (21) 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).
- (22) Here, the full calculation including subleading-order corrections yields .
- (23) R. A. Fisher and E. B. Ford, Heredity 1, 143 (1947).
- (24) Since the work of Fisher and Ford Fisher47 (), it is conjectured that noisy selection, , can explain the fluctuations other than those produced by demographic noise Kimura (); Jensen72 (); Karlin74 (). Because most fixation properties depend on the product (for ), our findings can also be related to fluctuations in the population size, , with kept fixed, a case particularly relevant for ecology.
- (25) Other non-Gaussian statistics are possible. Yet, without specific knowledge of the EN properties, we chose the OU noise which is arguably the simplest form of EN with Gaussian statistics and arbitrary correlation time.
- (26) is the time-continuous limit () of a temporally uncorrelated normal random variable with mean and variance Gardiner ().
- (27) D. M. Roma et. al., Phys. Rev. E. 71, 011902 (2005).
- (28) Using the expression of together with Eqs. (Cooperation dilemma in finite populations under fluctuating environments) and (Cooperation dilemma in finite populations under fluctuating environments), one finds a posteriori that which is negligible in the third of Eqs. (Cooperation dilemma in finite populations under fluctuating environments).
- (29) E. Y. Levine and B. Meerson, Phys. Rev. E 87, 032127 (2013).
- (30) In the Supplemental Material (SM) available along with this e-print, the different results for the CFP are summarized in a phase diagram, the various EN parameter regimes are discussed, and the CFP under strong adiabatic EN is derived.
- (31) It can be checked that Eq. (9) is well-defined in the double limit of () and .
- (32) As time is discrete in the Wright-Fisher model and the EN between two successive generations is uncorrelated, the correspondence requires to set in Eq. (9). Also, a population of size in the Wright-Fisher model maps onto a population of size in the Moran model Kimura (); MA10 ().
- (33) D. T. Gillespie, J. Comput. Phys. 22, 403 (1976); M. A. Gibson and J. Bruck, J. Phys. Chem. 104, 1876 (2000).
- (34) D. T. Gillespie, Phys. Rev. E 54, 2084 (1996).