Stochastic dynamics of the prisoner’s dilemma with cooperation facilitators
Abstract
In the framework of the paradigmatic prisoner’s dilemma game, we investigate the evolutionary dynamics of social dilemmas in the presence of “cooperation facilitators”. In our model, cooperators and defectors interact as in the classic prisoner’s dilemma, where selection favors defection. However, here the presence of a small number of cooperation facilitators enhances the fitness (reproductive potential) of cooperators, while it does not alter that of defectors. In a finite population of size , the dynamics of the prisoner’s dilemma with facilitators is characterized by the probability that cooperation takes over (fixation probability) and by the mean times to reach the absorbing states. These quantities are computed exactly and using FokkerPlanck equations. Our findings, corroborated by stochastic simulations, demonstrate that the influence of facilitators crucially depends on the difference between their density and the game’s costtobenefit ratio . When , the fixation of cooperators is likely in a large population and, under weak selection pressure, invasion and replacement of defection by cooperation is favored by selection if , where is the cooperation payoff benefit. When , the fixation probability of cooperators is exponentially enhanced by the presence of facilitators but defection is the dominating strategy.
pacs:
05.40.a, 02.50.r, 87.23.Kg, 87.23.GeI Introduction
Understanding the origin of cooperative behavior is a central issue in the life and behavioral sciences, and has recently been listed among the major scientific puzzles to be elucidated (1). Evolutionary game theory (EGT) provides the ideal framework to study the competition between species and there is a long tradition of modeling the evolution of cooperation using evolutionary games (2); (3). In recent years, these processes have increasingly been investigated using the methods of statistical physics, see e.g. (2) and references therein. In EGT, successful species spread at the expense of the others, and each individual’s reproductive potential (fitness) varies with the population’s composition that continuously changes in time. The interaction between the species is thus accounted for by a fitnessdependent (or “frequencydependent”) selection pressure (2), as observed in various experiments (4). Quite intriguingly, in such a setting the optimization of the fitness at an individual level can result in the reduction of the population overall fitness (2); (3). An influential example of such a paradoxical behavior, is provided by the celebrated prisoner’s dilemma (PD) game that serves as a metaphor for social dilemmas. In fact, in the classic PD individual interest leads to defection, even though mutual cooperation would be socially more beneficial (3); (2). While the PD is the paradigmatic model for the evolution of cooperation, its main prediction is at odds with the cooperative behavior that is commonly observed in experimental realizations (5); (4). This has motivated an upsurge of research aiming to identify the possible mechanisms capable of promoting cooperation in biological and social systems (6). Notably, it has been proposed that cooperation can be promoted by kin and group selection (7), as well as by conditional behavioral rules leading to direct or indirect reciprocity (8); (9). It has also been found that local interactions may promote cooperation in some social dilemmas (10). Furthermore, it has been shown that cooperation is supported in games with voluntary participation, or with punishment for noncooperation (11).
In this work, we investigate an alternative scenario for the spread of cooperation in social dilemmas: we consider the evolution of the prisoner’s dilemma in a finite population comprising a small number of “cooperation facilitators”. The facilitators participate in the dynamics only by enhancing the reproductive potential of cooperators, while they do not affect the fitness of defectors (see Sec. II below). To study the influence of cooperation facilitators on the prisoner’s dilemma dynamics, the evolution is modeled in terms of a birthdeath process and the fixation properties are studied analytically. In fact, it is well established that the evolutionary dynamics in finite populations is efficiently characterized by the probabilities of reaching the absorbing states, where the extinction of one or more species and the fixation of another occur (2); (12); (15); (14); (13). Here, we are particularly interested in the probability that, from a given initial composition, the population eventually comprises only cooperators and a small fraction of facilitators, but no defectors (“cooperation fixation probability”). The mean times for these events (mean fixation times) are also studied and our results are checked against stochastic simulations. This approach allows us to (i) discuss how demographic fluctuations alter the mean field predictions of the classic replicator equations (2), and (ii) thoroughly analyze the circumstances under which facilitators and selection favor a single cooperator invading and replacing a population of defectors.
This paper is organized as follows: The PD with cooperation facilitators is introduced in the next section, where some of its properties are discussed. In Section III the dynamics with the Fermi process is characterized by the fixation probability (Sec. III.A) and the mean fixation times (Sec. III.B). The dynamics with the Moran process is studied in Section IV, while we summarize our findings and present our conclusions in Section V.
Ii Prisoner’s dilemma with cooperation facilitators: model and dynamics
In evolutionary game theory, twoplayer games can be interpreted as dilemmas of cooperation. In fact, the two possible strategies can be interpreted as “cooperation” (C) and “defection” (D). The paradigm of social dilemma is provided by the classic prisoner’s dilemma (PD), whose main features are captured by the following payoff matrix giving the pairwise interaction between cooperators and defectors (2); (3); (10)(16):
(1) 
where and respectively represent the benefit and the cost of cooperation, with . Here, without loss of generality, we assume that . According to (1), mutual cooperation leads to a payoff and mutual defection gives a payoff ; whereas when one player defects and the other cooperates, the defector receives a payoff and the cooperators gets . In the (classic) PD, the dilemma arises from the fact that each individual is better off not cooperating, even though mutual cooperation enhances the population overall payoff. Hence, while cooperation is socially beneficial, defection is the only (strict) Nash equilibrium in the PD (2); (3).
In this work, we consider a finite population comprising individuals on a complete graph (no spatial structure). The number of cooperators and defectors is respectively denoted by and . In addition to cooperators and defectors, we consider that the population also comprises a fixed (small) number of “cooperation facilitators” (). These facilitators cooperate with players and therefore enhance the reproductive potential (fitness) of cooperators, while they leave the fitness of defectors unaltered, see below. Hence, while the number of cooperators and defectors in the population changes in time ( and vary), the total number of cooperators and defectors is conserved. According to the tenets of EGT, the variation in time of the number of cooperators and defectors depends on their average payoffs, and respectively, obtained from the payoff matrix (1). Here, since facilitators enhance by cooperating with C individuals and have no (direct) influence on , one has
(2) 
where we have excluded selfinteractions from the definition of the payoffs (2). The population average payoff is given by . It is worth noticing that the expression of now comprises a term reflecting the positive contribution of facilitators to the cooperators payoff. In evolutionary dynamics, it is customary to add a baseline constant, here set to , to the payoffs of the spreading species (2); (13), yielding the fitness of species C and D, respectively given by
(3) 
where, we have introduced the costtobenefit ratio (with ) and have used . Similarly, the average fitness of the entire population reads and grows linearly with the density of cooperators.
The size of the population being finite, the evolutionary dynamics is modeled as a continuoustime birthdeath process (2); (17); (18). In this model, only pairs of cooperators and defectors interact (according to (1)) and the stochastic dynamics is implemented as follows: (i) at each time step a pair of individuals is randomly chosen from the entire population; (ii) unless a pair of cooperatordefector is drawn, nothing happens; and (iii) if one picks a cooperatordefector pair, one of these individuals is randomly chosen for reproduction (proportionally to its fitness) and the other is replaced by the newborn offspring. Hence, at each interaction the number of cooperators increases or decreases by one. The time evolution of this birthdeath process can therefore be described by the random variable giving the number of cooperators and by the rates associated with the transitions , respectively. Here, we consider
(4) 
where accounts for the probability of picking a cooperatordefector pair, while are functions of the fitnesses (II) that encode the interactions (selection) according to the chosen “microscopic” update rule (2). We here discuss the cases where correspond to (i) the Fermi process (FP) (19); (20) and (ii) the Moran process (MP) (21); (12); (2); (14); (20) that are commonly used in EGT (2).
Stochastic evolutionary dynamics and the influence of selection are generally characterized by the fixation properties, namely the probability that a given species fixates (takes over) the whole population and by the mean time for such an event to occur (2); (12); (13). In the absence of facilitators, fixation happens when only one species survives and the population composition is uniform. Here, as the number of facilitators remains constant, fixation will be achieved when one of the absorbing states is reached and either all cooperators are replaced by defectors, or vice versa, resulting in a (nonuniform) population comprising facilitators and cooperators or defectors. In this work, we are particularly interested in the probability that, starting with cooperators, all defectors are eventually removed from the population and replaced by cooperators. As discussed in Sec. III.A., the fixation probability is necessary to establish when selection favors cooperation replacing defection (13). In the framework of the above birthdeath process (4), this probability obeys the backward master equation (18); (14); (2)
(5) 
with absorbing boundaries and . The formal solution of Eq. (5) reads (18); (14); (2)
(6) 
Since the above birthdeath process is a onedimensional Markov chain, other quantities like the mean fixation times (MFTs) can, in principle, be obtained exactly, but yield unwieldy expressions (18); (14). When the population size is large, it is often much more useful to describe the fixation properties in terms of the diffusion approximation obtained in the continuum limit () by a secondorder sizeexpansion of the master equation resulting in a FokkerPlanck equation (17); (18); (12); (15). By denoting and the initial density of cooperators and defectors, respectively; and with being the fraction of facilitators in the population, the (backward) FokkerPlanck equation (FPE) associated with (5) reads (17); (18)
(7) 
where and
(8)  
with and, as usual, the density changes by at each cooperatordefector interaction. In the realm of the FokkerPlanck equation, the timescale is such that the timestep is . The formulation in terms of the FPE allows a neat connection with the mean field treatment of the dynamics. In fact, when and all demographic fluctuations are negligible, the time variation of the density of cooperators is given by the drift term of (8) (17), i.e.
(9)  
As for the classic PD, this rate equation admits two absorbing fixed points, (no cooperators) and (no defectors), but possesses no interior fixed point since for the Fermi and Moran processes, see below. As discussed in what follows, the stability of these fixed points depends on the difference between the costtobenefit ratio and the fraction of facilitators.
Iii Dynamics with the Fermi Process
The stochastic dynamics of evolutionary games is often conveniently modeled in terms of the socalled Fermi process (FP), see, for example, (19). In the FP, at each timestep two individuals are randomly drawn from the entire population and one of them reproduces at the expense of the other that is replaced by the newborn offspring. This happens with a probability proportional to the difference between the fitness of the interacting individuals and given by the Fermi function from statistical physics. Since only the CD pairs interact, the dynamics with the FP is described by the birthdeath process defined by (4) and (19). With these expressions of , one checks that (since , see Eqs. (II),(24)), which confirms the absence of an interior fixed point in the mean field (continuum) limit.
iii.1 Fixation probability
With (II), the transition rates (4) for the Fermi process read
(10) 
with
(11) 
This quantity measures the selection pressure (22). Clearly, and . In the continuum limit , the densities , and are treated as continuous quantities, and the absence of selfinteraction is ignored yielding the transition rates (10)
(12) 
The rate equation corresponding to the mean field dynamics with the Fermi process is obtained by using (12) into (9), and is characterized by a single stable (absorbing) fixed point corresponding to a stationary density
(13) 
of cooperators. This means that cooperation prevails (, no defectors) in an infinitely large population comprising a fraction of facilitators higher than the costtobenefit ratio , i.e. when . However, as in the traditional PD, defection wins (, no cooperators) if is less than (). In other words, for cooperation to prevail () at mean field level, it is necessary that the density of facilitators compensates the cost of cooperation relative to its benefit.
When the population size is finite, demographic fluctuations are important and the evolution is thus no longer aptly described by the mean field dynamics (9). In particular, the mean field results (13) do not account for the nonzero probability that a single cooperator can invade and replace a (finite) population of defectors. Here, to investigate how the above mean field picture (9),(13) is altered by fluctuations arising in a finite population, we compute the probability that defection is eventually replaced by cooperation in a population comprising initially cooperators and defectors. Since , this probability can be obtained explicitly using (6) and, when (22), one finds
(14) 
while fixation probability of defectors is simply given by . As shown in Fig. 1, where results of stochastic simulations obtained using the Gillespie algorithm (23) are reported, the predictions of (14) are in excellent agreement with numerical simulations. The expression (14) implies that
(15) 
when . In particular, the cooperation fixation probability starting with a single cooperator reads
(16)  
The findings (14)(16), summarized in Fig. 1, illustrate how a small fraction of cooperation facilitators affects the fixation probabilities in a large, yet finite, population with an initial density of cooperators comparable to, or larger than, the density of defectors: When (), the fixation probability of cooperators is much higher than that of defectors, , and the spread of cooperation is thus efficiently promoted by facilitators. Yet, it is worth noticing that defectors still have finite probability to fixate even when (and ), contrary to the mean field prediction (13). The opposite situation arises when (), as shown in Fig. 1.
The results (15) and (16) can also be used to assess the influence of selection on the evolutionary dynamics (2): Following the seminal work of Ref. (13), we can establish when selection favors cooperation (C) invading and replacing defection (D). Selection is said to favor the replacement of D by C if the fixation probability of a single cooperator in a population of defectors is greater than in absence of selection pressure () when (22). With (16), this yields the condition that is generally satisfied in large populations under nonvanishing selection pressure. An interesting result arises when the selection intensity is weak and the population size is large, i.e. and . In such a limit, when (see Fig. 2 where ) and selection favors cooperation replacing defection provided that . Moreover, selection favors C invading D when (13). With (II), this yields the condition . Therefore, under weak selection ( and ) selection favors the invasion and replacement of D by C if . Since , one has and cooperation invading and replacing defection is favored by selection provided that
(17) 
One can also use the results (16) to determine the circumstances under which defection is evolutionary stable. In fact according to (13), and as natural extension of the concept of evolutionary stability for infinitely large populations and deterministic evolutionary dynamics (2), D is evolutionary stable in a finite population if (i) selection opposes C invading D, implying , i.e. ; and if (ii) selection opposes C replacing D, i.e. . The condition (ii) is clearly always satisfied when is finite, and in this case defection is evolutionary stable if . In the weak selection limit where (with ), the condition (ii) yields . Hence, defection is evolutionary stable under weak selection in a large population if . Since , this clearly means that defection is evolutionary stable and is the dominating strategy when . It is worth noticing that in the limit of an infinite population, , one recovers the mean field results (13): cooperation prevails only if , according to (17), and defection dominates otherwise.
The meaning of the results (15)(17) is illustrated in Fig. 2 where has been computed in populations comprising a small initial number of cooperators () and an excellent agreement with (15) and (16) has been found. In Fig. 2, and we notice that increases linearly in , with a slope steeper than when and selection favors cooperation replacing defection. The slope is less than when and the fixation of cooperation is opposed by selection. To further appreciate the implications of (14)(17), it is useful to compare these results with those obtained in the absence of facilitators. Putting in (14)(17), one recovers the results for the classic PD when cooperation fixation probability vanishes exponentially with the population size : and (2); (14) (see Fig. 1).
Our findings therefore demonstrate that facilitators greatly influence the probability that cooperation prevails and are summarized in Fig. 1. As illustrated in that figure, the influence of facilitators crucially depends on the difference between their density and the costtobenefit ratio :

When , the fixation of cooperators is likely (but not certain) even when they are initially in minority, i.e. even when initially .

When and , the fixation probability of a single cooperator is generally higher than in the absence of selection pressure (). In this case, with , selection favors cooperation invading and replacing defection, see (16) and Figs. 1 and 2. Furthermore, under weak selection pressure and in a large population ( and ), the fixation probability of a single cooperator is independent of , . In this case invasion and replacement of defection by cooperation is favored by selection if (17) is satisfied.

When , selection always opposes cooperation replacing defection. In this case, while defection is evolutionary stable and is the dominating strategy when , the cooperation fixation probability is exponentially enhanced by a small fraction of facilitators. Yet, cooperation is likely to fixate only if defectors are initially outnumbered by cooperators, i.e. if , as illustrated in Fig. 1.
iii.2 Mean fixation times
Another quantity of great interest to unveil the influence of facilitators in the evolutionary dynamics of the PD is the (unconditional) mean fixation time. This quantity gives the average time necessary to reach one of the absorbing boundaries, i.e. a population composition with either or cooperators. The unconditional mean fixation time (MFT), , for a system comprising initially cooperators obeys the following backward master equation (18); (14); (2) (where the timestep is )
(18) 
with boundary conditions . In principle, this equation can be solved exactly but the final result is cumbersome and not enlightening. Here, in the continuum limit , we work with the continuous quantities and , and adopt the approach of the diffusion theory (12); (17). The diffusion approximation is known to be particularly suited to analyze the dynamics under weak selection, which here corresponds to the regime where (13); (2); (15). Exact methods (when available) or other approximations (14), e.g. the WKB approach (20), are particularly useful to deal with the case of strong selection intensity and/or with phenomena like metastability. In the realm of the diffusion theory, the transition rates of the FP are given by (12) and the fixation probability of cooperation is obtained by solving (8) which yields
while for defection the probability is .
Similarly, the unconditional MFT is obtained by solving the backward FPE (17); (15), i.e.
(19) 
with the absorbing boundary conditions . When the drift and diffusive terms are of the same order, i.e. when , it follows from Eq. (19) that the MFT scales linearly with , i.e.
(20) 
The scaling function can be obtained explicitly by solving Eq. (19) using standard methods, see e.g. (17). For instance, when the initial density of cooperators and defectors is the same, , and , one finds
(21)  
where , denotes the usual exponential integral and is EulerMascheroni constant. While the expression of is usually cumbersome, some useful properties can be directly inferred from (19). In fact, as Eq. (19) is invariant under the transformation , one has when is kept fixed. The unconditional MFT in the Fermi process is therefore characterized by the symmetry
(22) 
where, on the righthandside is replaced by and transformed into , with kept fixed. Furthermore, when is kept fixed but varies, (19) is invariant under the transformation and , while the boundary conditions become and . In the weak selection regime , the second boundary condition can be approximated by , which allows a mapping onto (19) that yields:
(23) 
with fixed. The comparison between the solution of (19) and the results of stochastic simulations (for the FP with rates (10)) reported in Fig. 3 shows that the diffusion approximation aptly captures the functional dependence of , even though some deviations (of about ) can be noticed. These deviations stem from the selfinteraction terms that are excluded from (10) but not in the continuum limit (12) [e.g. in Fig. 3 one has and when , and and for ]. More importantly, the scaling (20) and the relationship (23) are confirmed by the numerical simulations of Fig. 3. In fact, in Fig. 3 we notice that is a humped function with a maximum well separated from the absorbing boundaries and located at when and, while scales linearly with , the presence of facilitators increases the unconditional MFT and its maximum value at the hump.
In addition to the unconditional MFT, it is also relevant to consider the mean time to specifically reach one of the absorbing boundaries. Hence, the conditional mean fixation times and respectively give the average time to reach the absorbing boundaries and (14); (20). As for the unconditional MFT, these quantities can be obtained from a backward FPE in the realm of the diffusion approximation. In fact, obeys , with the absorbing boundaries (12). Since and, from (14), , the conditional MFTs in the regime (weak selection) are related by the relationship where is kept fixed, as illustrated in Fig. 3. Furthermore, one has when and (), while when and (). This implies that As shown in Fig. 3, decreases while increases monotonically with .
The influence of facilitators on the unconditional and conditional MFTs is summarized in Fig. 3. We have found that in the PD with cooperation facilitators all conditional and unconditional MFTs scale linearly with the population size when (weak selection). While a similar scaling is also obtained in the absence of facilitators, the MFTs at a fixed value are found to be significantly increased by the presence of facilitators. Hence, the presence of cooperation facilitators has the quantitative effect to prolong the coexistence and the competition between cooperators and defectors before an absorbing state is reached, see Fig. 3.
Iv Dynamics with the fitnessdependent Moran process
The stochastic dynamics of evolutionary games is often implemented in terms of the Moran process, see e.g. (13); (2), that was originally introduced in population genetics (21); (12). In its essence, the Moran model is a birthdeath process where one randomly picked individual produces an offspring proportionally to its fitness relative to the population average fitness. The resulting offspring then replaces another individual that is randomly picked to be removed from the population whose size is therefore conserved. Here, as the interactions are between cooperators and defectors, the Moran process is implemented with and in (4). Since when [see (24) and (22)], one verifies that implying the absence of an interior fixed point in the mean field (continuum) limit. The Moran process is usually investigated when the selection intensity is weak, both for technical convenience (the mathematical treatment simplifies greatly) and for the biological relevance of such a limit (2); (12); (13). In this section, the stochastic dynamics with the Moran process is investigated in the weak selection limit, where , using the diffusion approximation.
iv.1 Fixation probability
In the continuum limit, the fitnesses (II) become
(24) 
with . The transition rates for the Moran process thus read
(25) 
With (24) and (25), the mean field dynamics is described by the rate equation (9) whose properties are similar to those discussed for the Fermi process. In particular, the rate equation (9) for the Moran process is also characterized by a single stable (absorbing) fixed point (no defectors) if and (no cooperators) if [see (13)].
To understand how the combined effect of nonlinear selection and demographic noise alters the mean field description, we now compute the cooperation fixation probability in the realm of the diffusion approximation. In such a setting, the fixation probability is given by the FPE (7,8) with the boundary conditions and . The solution of (7) is given by (17)
(26) 
(27)  
Introducing (27) into (26) and performing the integrals, one obtains
(28) 
As shown in Fig. 4, this result is in excellent agreement with numerical simulations and exhibits the same qualitative features obtained for the Fermi process (compare with Fig. 1). The finding (28) implies that in the weak selection limit where and , one has
(29) 
In particular, the probability that cooperation fixates starting with a single cooperator, when is given by . We therefore recover the result derived from (16) for the Fermi process. Clearly, this implies that under weak selection the fixation of a single cooperator is favored by selection if the nontrivial condition (17) is satisfied. Again, it is instructive to compare (28), (29) with the result obtained in the absence of facilitators, when decays to zero exponentially with . The influence of the facilitators on the fixation probabilities for the Moran process is summarized in Fig. 4, where the same features as in Fig. 1 are recognized and summarized as follows:

The fixation of cooperators is likely (but not certain) when the density of facilitators is higher than the costtobenefit ratio ().

When and , selection favors cooperation invading and replacing defection if (17) is satisfied. In particular, the fixation probability of a single cooperator is .

When , selection opposes cooperation replacing defection but the fixation probability of cooperators is exponentially enhanced by the presence of facilitators.
iv.2 Mean fixation times
In the realm of the diffusion approximation, the unconditional mean fixation time obeys the backward FPE , with the absorbing boundary conditions . In the weak selection regime and continuum limit, with (25), one has
With these expression, the backward FPE for the unconditional MFT reads
(30) 
with . When , Eq. (30) coincides with the FPE (19) for the Fermi process with an effective population size . The solution to (30) can therefore readily be obtained from (19) and (20). In particular, we infer from (20) that the MFT scales linearly with when , yielding
(31) 
where is the scaling function (20) obtained for the Fermi process. This function still satisfies the symmetry yielding , when is kept fixed, as in the Fermi process.
In the same manner, from (31) and (23), we infer
(32) 
when is kept fixed and is transformed into . The solution of (30), as well as the relationships (31) and (32), are in excellent agreement with the results of stochastic simulations reported in Fig. 5. As for the FP, we can also consider the conditional mean fixation times and it follows from (31) and (23) that for the Moran process the conditional MFTs are related by where is kept fixed, in agreement with the results of Fig. 5.
The influence of facilitators on the MFTs with the Moran process is summarized in Fig. 5, where the MFTs rescaled by a factor reproduce the same qualitative behavior obtained for the Fermi process (compare with Fig. 3) and is a humped function with a pronounced maximum. Again, all MFTs scale linearly with (in the weak selection limit). Yet, the comparison with the results for reveals that, at a fixed value of , the presence of facilitators increases the MFTs, see Fig. 5. Also, we notice that the monotonic dependence of and on is essentially independent of the sign of (in Fig. 5, and ).
V Summary and conclusion
In this work, we have proposed and investigated an alternative scenario leading to the spread of cooperation in social dilemmas. We have considered the evolutionary dynamics of the prisoner’s dilemma (PD) game in the presence of a small number of cooperation facilitators. These individuals participate in the dynamics only by enhancing the fitness of cooperators. The influence of facilitators on the evolutionary dynamics has been characterized by computing the model’s fixation properties in a finite population of size . Here, fixation occurs either in the state with only defectors (as in the classic PD), or in the state where the entire population is comprised of cooperators and facilitators. The dynamics has been implemented with the Fermi and Moran processes and the same qualitative results have been found, which demonstrates the robustness of our findings. Our analytical approach, corroborated by stochastic simulations, is based on an exact treatment and on the diffusion approximation (FokkerPlanck equation) of the underlying birthdeath process.
Our main results concern the fixation probabilities, whose properties crucially depend on whether the fraction of facilitators is more or less than the game’s costtobenefit ratio . When , we have shown that facilitators are very efficient in promoting the spread of cooperators whose fixation is likely (but not certain, contrary to the mean field predictions) in a large population with comparable initial densities of defectors and cooperators. Furthermore, when the selection intensity is weak and , we have demonstrated that the invasion and replacement of defectors by a single cooperator is favored by facilitators and selection if (where is the cooperation payoff benefit). When , defection is evolutionary stable and is the dominating strategy. In this case, while cooperation is unlikely to fixate, the fixation probability of cooperators is still exponentially enhanced by the presence of facilitators. We have also studied the (unconditional and conditional) mean fixation times in the weak selection limit and found that these quantities grow linearly with the population size. While a similar scaling is also obtained in the absence of facilitators, their presence has the effect of significantly increasing all the mean fixation times and hence to prolong the coexistence of cooperators and defectors.
In conclusion, this work demonstrates that the presence of a small number of cooperation facilitators can effectively enhance the spread of cooperation in a simple model of social dilemmas and prolong the coexistence of competing species. The influence of facilitators is particularly drastic when their abundance exceeds the game’s costtobenefit ratio, in which case cooperation is generally the strategy favored by selection in large populations. These findings pave the way to further investigations of the influence of facilitators in other social dilemmas, e.g. with mixed strategies and/or in spatial settings.
References
 E. Pennisi, Science 309, 90 & 93 (2005).
 J. Maynard Smith, Evolution and the Theory of Games (Cambridge University Press, Cambridge, 1982); J. Hofbauer and K. Sigmund, Evolutionary Games and Population Dynamics (Cambridge University Press, Cambridge, 1998); H. Gintis, Game Theory Evolving (Princeton University Press, Princeton, 2000); M. A. Nowak, Evolutionary Dynamics (Belknap Press, 2006); G. Szabó and G. Fáth, Phys. Rep. 446, 97 (2007); A. Traulsen and C. Hauert in Reviews of Nonlinear Dynamics and Complexity, edited by H.G. Schuster Vol. 2 (WileyVCH, 2010).
 R. Axelrod and W. D. Hamilton, Science 211, 1390 (1981); R. Axelrod, The Evolution of Cooperation (Basic Books, New York, 1984).
 P. E. Turner and L. Chao, Nature (London) 398, 441 (1999); J. Gore, H. Youk, and A. van Oudenaarden, ibid 459, 253 (2009).
 D. Semmann, H. J. Krambeck and M. Milinski, Nature (London) 425, 390 (2003); A. Traulsen, D. Semmann, R. D. Sommerfeld, H. J. Krambeck and M. Milinski, Proc. Natl. Acad. Sci. USA 107, 2962 (2010).
 M. Doebeli and C. Hauert, Ecol. Lett. 8, 748 (2005); M. A. Nowak, Science 314, 1560 (2006).
 J. A. Fletcher and M. Zwick, J. Theor. Biol. 228, 303 (2004); W. D. Hamilton,ibid 7, 1 (1964); W. D. Hamilton, ibid. 7, 17 (1964); D. S. Wilson, Proc. Natl. Acad. Sci. USA 72, 143 (1975); A. Traulsen and M. A. Nowak, Proc. Natl. Acad. Sci. USA 103, 10952 (2006).
 R. L. Trivers, Quar. Rev. Biol. 46, 35 (1971); M. A. Nowak and K. Sigmund, Nature (London) 364, 56 (1993); L. A. Imhof, D. Fudenberg, and M. A. Nowak, Proc. Natl. Acad. Sci. USA 102, 10797 (2005); M. A. Nowak and K. Sigmund, Nature (London) 437, 1291 (2005); L. A. Imhof, D. Fudenberg, and M. A. Nowak, J. Theor. Biol. 247, 574 (2007); J. M. Pacheco, A. Traulsen, H. Ohtsuki and M. A. Nowak, J. Theor. Biol. 250, 723 (2008); A. J. Bladon, T. Galla, and A. J. McKane, Phys. Rev. E 81, 066122 (2010).
 M. A. Nowak and K. Sigmund, Nature (London) 393, 573 (1998); R. Ferrière, ibid. 7 393, 517 (1998); O. Leimar and P. Hammerstein, Proc. R. Soc. Lond. B 268, 745 (2001); K. Panchanathan and R. Boyd, J. Theor. Biol. 224, 115 (2003); H. Ohtsuki and Y. Iwasa, J. Theor. Biol. 231, 107 (2004).
 M. A. Nowak and R. M. May, Nature (London) 359, 826 (1992); C. Hauert and G. Szabó, Complexity 8, 31 (2003); C. Hauert and M. Doebeli, Nature (London) 428, 643 (2004); F. C. Santos and J. M. Pacheco, Phys. Rev. Lett. 95, 098104 (2005); H. Ohtsuki and M. A. Nowak, J. Theor. Biol. 243, 86 (2006), ibid. 251, 698 (2008); C. E. Tarnita, H. Ohtsuki, T. Antal, F. Fu, and M. A. Nowak, J. Theor. Biol. 259, 570 (2009); Z. Wang, A. Szolnoki, and M. Perc, Sci. Rep. 2, 369 (2012); M. Assaf and M. Mobilia, arXiv:1202.3231v1.
 C. Hauert, S. De Monte, J. Hofbauer, and K. Sigmund, Science 296, 1129 (2002); C. Hauert, S. De Monte, J. Hofbauer, and K. Sigmund, J. Theor. Biol. 218, 187 (2002); C. Hauert, A. Traulsen, H. Brandt, M. A. Nowak, and K. Sigmund, Science 316, 1905 (2007); M. Perc and A. Szolnoki, New J. Phys. 14, 043013 (2012).
 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).
 M. A. Nowak, A. Sasaki, C. Taylor, and D. Fudenberg, Nature (London) 428, 646 (2004); C. Taylor, D. Fudenberg, A. Sasaki, and M. A. Nowak, Bull. Math. Biol. 66, 1621 (2004).
 T. Antal and I. Scheuring, Bull. Math. Biol. 68, 1923 (2006).
 A. Traulsen, J. C. Claussen, and C. Hauert, Phys. Rev. Lett. 95, 238701 (2005); A. Traulsen, J. M. Pacheco, and L. A. Imhof, Phys. Rev. E, 74, 021905 (2006); A. Traulsen, J. C. Claussen, and C. Hauert, Phys. Rev. E 74, 011901 (2006); M. Mobilia, EPL 95, 50002 (2011).
 Such a formulation is commonly used since it represents the most convenient way to capture the essential features of the classic PD (10). It is an example of “equal gains from switching” games, see e.g. (2) and references therein.
 C. W. Gardiner, Handbook of Stochastic Methods, (Springer, New York, 2002); H. Risken, The FokkerPlanck equation, (Springer, New York, 1989).
 N. G. van Kampen, Stochastic Processes in Physics and Chemistry (Elsevier, Amsterdam, 1997).
 L. E. Blume, Games Econ. Behav. 5, 387 (1993); G. Szabó and C. Töke, Phys. Rev. E 58, 69 (1998); C. Hauert and G. Szabó, Am. J. Phys. 73, 405 (2005); A. Traulsen, M. A. Nowak, and J. M. Pacheco, Phys. Rev. E 74, 011909 (2006).
 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) and references therein.
 P. A. P. Moran, The statistical processes of evolutionary theory (Clarendon, Oxford, 1962).
 When the density of facilitators is , the selective pressure due to the costtobenefit ratio is balanced and one has . In this special case, the dynamics is neutral (there is effectively no selection) with , and (6) yields the expected result .
 D. T. Gillespie, J. Comput. Phys. 22, 403 (1976).