Large Fluctuations in Anti-Coordination Games on Scale-Free Graphs

TL;DR

This paper investigates how the complex topology of scale-free networks influences the dynamics of anti-coordination games, revealing that different update rules lead to distinct effects on metastability and fixation times.

Contribution

It demonstrates that the topology significantly affects game dynamics under voter updates but not under link dynamics, highlighting the role of network structure and update rules.

Findings

Voter update rule causes topology-dependent behavior.

Link dynamics results in behavior similar to complete graphs.

Scale-free topology renormalizes population size under voter updates.

Abstract

We study the influence of the complex topology of scale-free graphs on the dynamics of anti-coordination games (e.g. snowdrift games). These reference models are characterized by the coexistence (evolutionary stable mixed strategy) of two competing species, say "cooperators" and "defectors", and, in finite systems, by metastability and large-fluctuation-driven fixation. In this work, we use extensive computer simulations and an effective diffusion approximation (in the weak selection limit) to determine under which circumstances, depending on the individual-based update rules, the topology drastically affects the long-time behavior of anti-coordination games. In particular, we compute the variance of the number of cooperators in the metastable state and the mean fixation time when the dynamics is implemented according to the voter model (death-first/birth-second process) and the link dynamics (birth/death or death/birth at random). For the voter update rule, we show that the scale-free topology effectively renormalizes the population size and as a result the statistics of observables depend on the network's degree distribution. In contrast, such a renormalization does not occur with the link dynamics update rule and we recover the same behavior as on complete graphs.

Equations

$${\mu }_m\equiv \sum _k{k}^m{n}_k=\sum _i{k}_i^m/N.$$

$$\rho \equiv \sum _i{\eta }_i/N,{\rho }_k\equiv \sum _i^{\prime }{\eta }_i/N_k,$$

$$\omega \equiv \frac{1}{N{\mu }_1}\sum _i{k}_i{\eta }_i=\sum _k\frac{k}{{\mu }_1}{n}_k{\rho }_k.$$

$$\frac{d}{dt}\rho (t)=\rho (t)(1-\rho (t))[{\Pi }^C(\rho (t))-{\Pi }^D(\rho (t))]=(a+d-b-c)\rho (t)(1-\rho (t))(\rho (t)-{\rho }^{\ast }).$$

$$f_i^C=1+s[{\Pi }_i^C-{\bar{\Pi }}_i]=1+s(b-d)(1-\rho ),f_i^D=1+s[{\Pi }_i^D-{\bar{\Pi }}_i]=1+s(c-a)\rho .$$

$$1/f^C=1-s(b-d)(1-\rho )+\mathcal{O}(s^2),1/f^D=1-s(c-a)\rho +\mathcal{O}(s^2).$$

$$f^C=1+s(b-d)(1-\rho ),f^D=1+s(c-a)\rho .$$

$${\Psi }_{ij}^{\mathrm{VM}}=(1-{\eta }_i){\eta }_j/f^D\text{and}{\Psi }_{ji}^{\mathrm{VM}}=(1-{\eta }_j){\eta }_i/f^C,$$

$${\Psi }_{ij}^{\mathrm{LD}}=(1-{\eta }_i){\eta }_j{f}^C\text{and}{\Psi }_{ji}^{\mathrm{LD}}=(1-{\eta }_j){\eta }_i{f}^D,$$

$$T^{+}({\rho }_k)=\sum _i^{\prime }\sum _j\frac{A_{ij}}{N\mathcal{Q}}{\Psi }_{ij}^{\mathrm{VM}/\mathrm{LD}}\text{and}{T}^{-}({\rho }_k)=\sum _i^{\prime }\sum _j\frac{A_{ij}}{N\mathcal{Q}}{\Psi }_{ji}^{\mathrm{VM}/\mathrm{LD}},$$

$$\displaystyle{\cal Q}=\left\{\begin{array}[]{ll}k_{i}&\mbox{for the VM},\\ \mu_{1}&\mbox{for the LD}.\end{array}\right.$$

$$T_k^{+}=n_k\omega (1-{\rho }_k)\left[1-s\rho (c-a)\right],T_k^{-}=n_k(1-\omega ){\rho }_k\left[1-s(1-\rho )(b-d)\right],$$

$$T_k^{+}=n_k\frac{k}{{\mu }_1}\omega (1-{\rho }_k)\left[1+s(1-\rho )(b-d)\right],T_k^{-}=n_k\frac{k}{{\mu }_1}(1-\omega ){\rho }_k\left[1+s\rho (c-a)\right].$$

$${\mathcal{G}}_{\mathrm{f}}(\{{\rho }_k\})$$

$${\mathcal{G}}_{\mathrm{b}}(\{{\rho }_k\})$$

$${\tilde{\mathcal{G}}}_{\mathrm{f}}(\omega )$$

$${\tilde{\mathcal{G}}}_{\mathrm{b}}(\omega )$$

$$\displaystyle N_{{\rm eff}}=N\frac{(\mu_{1})^{2}}{\mu_{2}}\sim\left\{\begin{array}[]{ll}N,&\nu>3,\\ N/\ln{N},&\nu=3,\\ N^{\alpha},&2<\nu<3,\end{array}\right.$$

$$\alpha =\frac{2(\nu -2)}{\nu -1}.$$

$${\tilde{\mathcal{G}}}_{\mathrm{f}}(\rho )=\frac{\partial }{\partial \rho }\left[\tilde{s}\rho (1-\rho )(\rho -{\rho }^{\ast })+\frac{1}{N}\frac{\partial }{\partial \rho }\rho (1-\rho )\right],{\tilde{\mathcal{G}}}_{\mathrm{b}}(\rho )=\rho (1-\rho )\left[\tilde{s}({\rho }^{\ast }-\rho )\frac{\partial }{\partial \rho }+\frac{1}{N}\frac{{\partial }^2}{\partial {\rho }^2}\right].$$

$${\tilde{\mathcal{G}}}_{\mathrm{f}}(\xi )={\rho }^{\ast }(1-{\rho }^{\ast })\left[-\tilde{s}\frac{\partial }{\partial \xi }\xi +\frac{1}{N_{\mathrm{eff}}}\frac{{\partial }^2}{\partial {\xi }^2}\right].$$

$$⟨{\xi }^2(t)⟩=\frac{1}{\tilde{s}{N}_{\mathrm{eff}}}\left[1-e^{-2\tilde{s}{\rho }^{\ast }(1-{\rho }^{\ast })t}\right]\to ⟨{\xi }^2⟩=\frac{1}{\tilde{s}{N}_{\mathrm{eff}}}$$

$$\displaystyle{\rm var}(N\rho)=N^{2}(\langle\rho^{2}\rangle-\langle\rho\rangle^{2})\sim\left\{\begin{array}[]{ll}N,&\nu>3,\\ N\ln{N},&\nu=3,\\ N^{2-\alpha}=N^{2/(\nu-1)},&2<\nu<3.\end{array}\right.$$

$$\displaystyle\sigma=\left\{\begin{array}[]{ll}\tilde{s}N\frac{\mu_{1}^{2}}{\mu_{2}},&\text{(for the VM)},\\ \tilde{s}N,&\text{(for the LD)}.\end{array}\right.$$

$$\displaystyle T_{\rm fix}(\rho)\sim\left\{\begin{array}[]{ll}(1-\phi^{C}(\rho))e^{(\rho^{*})^{2}\,\sigma},&\text{when $\rho>\rho^{*}$},\\ \phi^{C}(\rho)e^{(1-\rho^{*})^{2}\,\sigma},&\text{otherwise.}\end{array}\right.$$

$$\ln T_{\mathrm{fix}}(\rho )\simeq ({\rho }^{\ast }{)}^2\sigma .$$

Full text

I Introduction

Evolutionary game theory (EGT) is a suitable framework to model the dynamics of populations in which the success of one type depends on the actions of the others. In EGT, selection varies with the species densities and is thus “frequency dependent” EGT1 ; EGT2 ; EGT3 ; Nowak . This means that in the realm of EGT the population composition and each species’ fitness change continuously in time. The ensuing evolutionary dynamics is commonly modeled deterministically in terms of the celebrated “replicator equations”, which are nonlinear differential equations EGT1 ; EGT2 ; EGT3 ; Nowak ; freqdepsel1 ; freqdepsel2 well-suited to describe very large populations. The size of a real population is however always finite and is more realistically described by stochastic models whose properties, such as demographic fluctuations, are known to often greatly influence the evolutionary dynamics weaksel1 ; weaksel2 ; weaksel3 . In particular, due to randomness, individuals of one species can take over and fixate the entire population. This stochastic phenomenon, referred to as “fixation”, is characterized by the fixation probability – the probability that a “mutant type” takes over Nowak ; Kimura ; Ewens – and the mean fixation time (MFT) which is the mean time for such an event to occur. Population dynamics is also known to depend on the individuals’ spatial arrangement: while a large body of work has focused on EGT on spatially-homogeneous (well-mixed) populations EGT1 ; EGT2 ; EGT3 , it is known that the outcome of the dynamics may be very different in spatial settings, see e.g. RPS1 ; RPS2 ; RPS3 ; RPS4 ; RPS5 ; RPS6 . For instance, spatial degrees of freedom have been found to promote cooperation in the prisoner’s dilemma game but to hinder coexistence in snowdrift games Spatialcoop1 ; Spatialcoop2 ; Spatialcoop3 ; Xia15a ; EGT1 ; EGT2 ; EGT3 .

EGT was first introduced to study ecological dynamics EGT1 and is also particularly well suited to model the evolution of social behavior of interacting agents EGT2 ; EGT3 . In this context, evolutionary games on graphs NetRef ; Nets1 ; Nets2 ; Nets3 provides a general unifying framework that is able to capture the dynamics of spatially structured populations: it models how individuals interact with their neighbors, to reproduce or die, as prescribed by the underlying game  EGTnets1 ; EGTnets2 ; EGTnets3 ; EGTgraphs ; EGTgraphs-bis . Of particular interest is the question of understanding how the network’s topology affects the fixation properties of a given process, see e.g. Refs. EGTgraphs ; EGTgraphs-bis ; Allen2016 . While it is difficult to give a general answer to this question Ibsen15 , significant progress can be made when the selection pressure is weak, see e.g. weaksel1 ; weaksel2 ; weaksel3 . In fact, most studies have focused on the biologically relevant and tractable limit of weak selection in which exact results have been obtained on regular graphs, see e.g. EGTgraphs-bis ; Taylor07 ; Chen13 . However, much less is known about the fixation properties of evolutionary games on degree-heterogeneous graphs such as scale-free networks NetRef ; Nets1 ; Nets2 ; Nets3 . Most investigations on this important class of graphs have been carried out by means of computer simulations Pacheco1 ; Pacheco2 ; hubs ; Macie14 ; Xia15b , various approximation schemes Tarnita09 ; Konno11 ; Allen2016 , and by considering special graphs Broom11 ; Macie14 . The prisoner’s dilemma game has been extensively studied on scale-free networks and it has been found that the existence of nodes with high connectivity (hubs) can promote cooperation even under adverse conditions Pacheco1 ; Pacheco2 ; EGTgraphs ; hubs ; Xia15b , whereas the presence at the hubs of so-called facilitators facil1 ; facil1bis , which are agents that cooperate with cooperators and defect with defectors, tame the cooperation dilemma facil2 . Recently, tools of statistical mechanics have been used to study simple evolutionary processes in which two types, say cooperators and defectors, interact on degree-heterogeneous graphs under constant fitness and fixed (weak) selection pressure VotNets31 ; VotNets32 ; VotNets11 ; VotNets12 ; VotNets13 ; VotNets21 ; VotNets22 ; VotNets23 . While these models shed light on intriguing properties of evolution on complex graphs, such as the fact that it depends on the microscopic details of the update rules, they cannot describe evolutionary processes characterized by a metastable species coexistence prior to fixation like in the paradigmatic anti-coordination games (ACGs) EGT1 ; EGT2 ; EGT3 ; Gore09 , see e.g. AM1 ; AM2 ; AM3 ; AM4 ; MA10 ; AM10 ; AM16 . In spite of the importance of systems like ACGs, their properties have been investigated mostly on regular lattices Killingback96 ; Doebeli05 ; Spatialcoop2 and on small-world networks Shang06 ; Qiu10 . In fact, there are very few results, mostly based on computer simulations Pacheco1 , on how the scale-free topology affects the metastability and fixation properties of ACGs. Here, the analysis attempted in Refs. PRL12 ; AM13 is critically revisited and generalized.

In this work we study the joint effect of degree-heterogeneous topology and frequency-dependent selection freqdepsel1 ; freqdepsel2 on the dynamics of ACGs on scale-free networks under weak selection. These are well-known EGT models, of particular importance in biology and ecology EGT1 ; EGT2 ; EGT3 ; Gore09 , characterized by a long-lived coexistence state (metastability), and in which fixation is driven by large fluctuations MA10 ; AM10 . In particular, we investigate the influence of the microscopic update rule on the evolutionary dynamics. For the sake of concreteness, here the individual-based dynamics is implemented according to two common update rules, namely with the voter model (VM, death-first/birth-second process) and the link dynamics (LD, birth/death or death/birth at random) VotNets31 ; VotNets32 ; VotNets11 ; VotNets12 ; VotNets13 ; VotNets21 ; VotNets22 ; VotNets23 . By combining analytical and simulation means, we consider systems of large but finite size and determine the effect of the complex topology on the evolutionary dynamics. In particular, with the VM we show that the typical fluctuations in the number of cooperators in the metastable state is anomalous (their variance grows superlinearly in the population size), and that the MFT has a stretched exponential dependence on the population size. We also show that with the LD these quantities coincide in the leading order with the results on complete graphs. It is worth noting that our approach is not limited to scale-free graphs, and is valid for general degree-heterogeneous networks, see, e.g. Refs. VotNets11 ; VotNets12 ; VotNets13 ; VotNets21 ; VotNets22 ; VotNets23 .

This paper is organized as follows: The class of ACGs on networks that we consider is introduced in the next section. Section III is dedicated to a description of the implementation of our computer simulations, while the analytical approach in terms of a multivariate diffusion theory is presented in Sec. IV. The various timescales that characterize the dynamics are presented in Sec. V, where timescale separation is used to derive an effective single-variate diffusion theory. Such an approach, corroborated by extensive simulations of the individual-based system, allows us to characterize the typical fluctuations in the metastable state by determining the variance in the number of cooperators in Sec. VI, and to obtain the MFT in Sec. VII. Finally, we summarize our findings and present our conclusions. Technical details about our computational methods are provided in Appendix A and B.

II The Models: ACGs on complex networks

As in Ref. PRL12 , we consider a scale-free network consisting of $N$ nodes on which population dynamics takes place between two types of agents (see Appendix A), here referred to as cooperators (C’s) and defectors (D’s). Each node is associated with a binary random variable: $\eta_{i}=1$ if the node $i$ is occupied by a C, whereas $\eta_{i}=0$ if it is occupied by a D. The topology of the network is defined by its adjacency matrix ${\bm{A}}=[A_{ij}]$ NetRef ; Nets1 ; Nets2 ; Nets3 , whose elements are $1$ if the nodes $ij$ are connected and [math] otherwise, and its state (or population composition) is described by $\{\eta_{i}\}^{N}=\{\eta_{1},\dots,\eta_{N}\}$. The underlying complex graph is characterized by a degree distribution $n_{k}=N_{k}/N$, where $N_{k}$ denotes the number of nodes of degree $k$, whose $m^{th}$ moment is defined by

[TABLE]

Here $k_{i}$ is the degree of the node $i$, $\mu_{1}$ is the graph’s mean degree and $N\mu_{1}/2$ is the average number of links. There are standard methods to generate random networks whose nodes are distributed according to a prescribed degree distribution, see e.g. Ref. Nets1 ; Nets2 ; Nets3 . Here we use the method outlined in Appendix A to generate the scale-free graphs that we will consider in this work. To study how the population composition changes in time, we introduce the density $\rho$ of cooperators in the network, and the subgraph density $\rho_{k}$ of C’s on nodes of degree $k$. These quantities are defined by

[TABLE]

where $\sum_{i}^{\prime}$ denotes a summation restricted to the nodes $i$ of fixed degree $k$. We note that $\rho=\sum_{k}n_{k}\rho_{k}$. It is also useful to introduce the degree-weighted density of cooperators VotNets11 ; VotNets12 ; VotNets13 :

[TABLE]

Here, we are interested in the class of anti-coordination games (ACGs) which are symmetric two-player two-strategy games. According to the tenets of EGT EGT1 ; EGT2 ; EGT3 , connected cooperators and defectors compete (or “play”) pairwise according to the payoff matrix

[TABLE]

This specifies that two interacting cooperators get a payoff $a$, whereas defectors playing against each other both get a payoff $d$. Moreover, when a C plays against a D the former’s payoff is $b$ and the latter gets a payoff $c$. It is well known that various scenarios, including “cooperation dilemma”, emerge depending on the various values of the entries of (7) EGT1 ; EGT2 ; EGT3 ; EGTnets1 ; EGTnets2 ; EGTnets3 . In this work, we focus on the important class of ACGs for which $c>a$ and $b>d$. The class of ACGs includes the snowdrift game, for which $c>a>b>d$, that is particularly relevant to biological applications, see e.g. Ref. Gore09 .

II.1 Well-mixed setting

In the usual setting of EGT where the population is well-mixed, the expected payoffs (per individual) to cooperators and defectors are respectively $\Pi^{C}(\rho)=a\rho+b(1-\rho)$ and $\Pi^{D}(\rho)=c\rho+d(1-\rho)$, while the population average payoff is $\bar{\Pi}(\rho)=\rho\Pi^{C}(\rho)+(1-\rho)\Pi^{D}(\rho)$. In such a setting, one of the main features of ACGs is a coexistence state in which a fraction $\rho^{*}$ of C’s coexist with a density $1-\rho^{*}$ of D’s over a long period of time. In the game-theoretic language the strategy $(\rho^{*},1-\rho^{*})$ corresponds to playing cooperation C with a frequency $\rho^{*}$ and defection D with frequency $1-\rho^{*}$. This mixed strategy is known to be evolutionary stable (but it is not a strict Nash equilibrium) EGT1 ; EGT2 ; EGT3 . In order to find $\rho^{*}$, we write down the mean-field replicator equation (RE), which reads EGT1 ; EGT2 ; EGT3 ; Nowak :

[TABLE]

This equation is characterized by a stable interior fixed point $\rho^{*}=(b-d)/(b+c-a-d)$ and unstable absorbing states $\rho=0$ (all-D) and $\rho=1$ (all-C). That is, in the deterministic picture, starting from any $0<\rho<1$ the system settles into the stable point $\rho^{*}$ and stays there forever. This picture, however, is altered when the population size is finite ($N<\infty$), due to demographic fluctuations which ultimately drive the system into one of its two absorbing states. As a result, the stable fixed point in the language of the RE, $\rho^{*}$, becomes metastable, and the probability to be in its vicinity slowly decays while the probability to be absorbed in either absorbing states slowly grows. It is well known that the mean decay time of the metastable state, which approximately equals the MFT, grows exponentially with $N$ on complete graphs (well-mixed population, $A_{ij}=1,\forall ij$) MA10 ; AM10 .

II.2 Spatially-structured setting

When the population is spatially-structured and occupies the vertices of a graph, the interactions are among nearest-neighbor agents. The corresponding expected payoffs are thus defined locally: If a node $j$ is occupied by a D individual, its neighbor $i$ receives a payoff $\Pi_{i}^{C}=b$ if the node $i$ is occupied by a C and the payoff of the agent node $i$ is $\Pi_{i}^{D}=d$ if it is a D individual. The local reproductive potential, or fitness, of the agent at node $i$ whose neighbor $j$ is a D individual is proportional to the difference of their expected payoff relative to the population average payoff $\bar{\Pi}_{i}(t)$ as perceived by the agent at node $i$. For the latter, we make the choice to consider $\bar{\Pi}=\rho(t)\Pi_{i}^{C}+(1-\rho(t))\Pi_{i}^{D}$ PRL12 . This mean-field-like form reflects in a simple manner the fact that agents compare their payoffs with those of all others, leading to metastability via a natural and analytically amenable mechanism. As customary in EGT, we also introduce in the definition of the fitness a selection strength $s>0$, accounting for the interplay between demographic fluctuations and selection, as well as a baseline contribution, accounting for the chance contribution to the reproduction, which we set to $1$ Nowak ; weaksel1 ; weaksel2 ; weaksel3 ; Kimura ; Ewens . In this setting, the fitnesses of a C (D) player at node $i$ against a D (C) player at a neighboring node $j$ are 111Naturally, other update rules, such as the Moran-like, the local-update or Fermi-like rules, may also be considered. However, for $s\ll 1$, all these update rules approximately coincide, see, e.g., Refs. weaksel2 ; weaksel3 ; facil1 .

[TABLE]

We now proceed to specify the update rules by which, at an individual-based level, the population evolves. Various types of update rules are possible and, in the case of models without fitness-dependent selection, it has been found that the dynamics may depend crucially on the details of the underlying rules VotNets11 ; VotNets12 ; VotNets13 ; VotNets21 ; VotNets22 ; VotNets23 ; VotNets31 ; VotNets32 . Here, we consider two important choices: the so called (i) “voter model” (VM) rule VotNets11 ; VotNets12 ; VotNets13 ; (ii) and the “link dynamics” (LD) VotNets31 ; VotNets32 ; PRL12 ; AM13 :

(i) In the VM dynamics, a focal agent $i$ is chosen at random with probability $1/N$, and then one of its neighbors is picked with a probability $1/k_{i}$. In this death-first/birth-second process, the focal agent dies and is replaced by the picked neighbor with a probability proportional to the fitness of the latter. Or, equivalently when $0<s\ll 1$, and as implemented in Ref. VotNets11 ; VotNets12 ; VotNets13 for the biased VM and here in our simulations (see below), the focal agent dies and is replaced by an offspring of the picked neighbor with a probability proportional to the inverse of the focal agent’s fitness, see Sec. III. In practice, this means that with the VM at each time increment the population composition changes only when a neighboring CD or DC pair interacts. The reactions $\textsf{C}\textsf{D}\to\textsf{D}\textsf{D}$ (death and replacement of C at $i$ by a D) and $\textsf{D}\textsf{C}\to\textsf{C}\textsf{C}$ (death and replacement of a D at $i$ by a C) thus occur with rates given by the inverse of the individual C’s and D’s fitness respectively, i.e. with rates

[TABLE]

(ii) In the LD, a link is randomly selected at each time step and if it connects a CD pair, one of the neighbors is randomly selected for reproduction with a rate proportional to its fitness, while the other is replaced by the newly produced offspring. In practice, this means that with the LD at each time increment the population composition changes only when a CD or DC pair at neighboring nodes interact. Thus, the reactions $\textsf{D}\textsf{C}\to\textsf{C}\textsf{C}$ and $\textsf{C}\textsf{D}\to\textsf{D}\textsf{D}$ occur with rates given by the fitness $f^{C/D}$ of the agent C and D respectively, i.e. with rates

[TABLE]

In what follows, we consider the evolution of the ACGs under both VM and LD update rules and we show that markedly different behavior emerges.

III Computer simulations of the evolutionary dynamics

Before theoretically analyzing the evolution of ACGs with the VM and LD, in this section we describe how the evolutionary dynamics with the VM and LD have been implemented in our computational individual-based simulations.

We begin by outlining the simulation of the evolution with the LD which has been performed using the Gillespie algorithm Gillespie . Our starting point is a scale-free network, see Appendix A, where each node is in one of two states - C (cooperator) or D (defector). The following steps are repeated until fixation of either C or D occur (or until a prescribed maximal number of time-steps have been performed):

Compute the density $\rho$ of C and the fitnesses of C and D according to Eqs. (9). 2. 2.

Draw random numbers $R_{1}$ and $R_{2}$ uniformly distributed between [math] and $1$. 3. 3.

Pick a random node $i$ and a random neighbor $j$ of $i$. If the states of $i$ and $j$ are different, then if $R_{1}<\frac{f^{C}}{f^{C}+f^{D}}$, both nodes become C. Otherwise, they become D. 4. 4.

Increment time by $\Delta t=-\frac{1}{N}\frac{\ln R_{2}}{f^{C}+f^{D}}$.

The implementation of the VM case goes along the same lines as that of the LD model, except for the update of the randomly chosen node $i$ and its neighbor $j$ in step $3$. For the VM, if node $i$ is a C and node $j$ is a D, then $i$ becomes a D with probability $1/f^{C}$. Otherwise, if node $i$ is a D and node $j$ is a C, then $i$ becomes a C with probability $1/f^{D}$.

IV Diffusion theory

We are particularly interested in the biologically-relevant regime of weak selection strength. In such a limit where $0<s\ll 1$ weaksel1 ; weaksel2 ; weaksel3 ; Kimura ; Ewens , the evolutionary dynamics is generally well described in terms of the so-called diffusion theory using the of appropriate forward and backward Fokker-Planck equations (FPEs) Gardiner . Below, the multivariate FPEs for the subgraph densities are derived. We introduce the following quantities for the evolution with the VM

[TABLE]

while for the LD we use the quantities

[TABLE]

where $(1-\eta_{i})\eta_{j}$ and $(1-\eta_{j})\eta_{i}$ are non-zero only when the nodes $ij$ are occupied by a pair DC and CD, respectively. In an infinitesimal time increment $\delta t=1/N$ the subgraph density $\rho_{k}$ changes by $\pm\delta\rho_{k}=\pm 1/N_{k}$ according to a birth-death process Gardiner defined respectively by the transition rates

[TABLE]

where, from the definition of the VM and LD, we have VotNets11 ; VotNets12 ; VotNets13

[TABLE]

Given an agent at node $i$, the probability to pick one of its neighbors $j$ for an update is $A_{ij}/(N{\cal Q})$, and the transition $\eta_{i}\to 1-\eta_{i}$ occurs with probability $\sum_{j}\frac{A_{ij}}{N{\cal Q}}\left[\Psi_{ij}^{{\rm VM}/{\rm LD}}+\Psi_{ji}^{{\rm VM}/{\rm LD}}\right]$ VotNets11 ; VotNets12 ; VotNets13 ; PRL12 .

To make analytical progress, we assume that the degrees of the nodes of the underlying networks are uncorrelated, i.e. we consider degree-uncorrelated heterogeneous graphs. As explained in Appendix A, such an assumption, which is exact for Molloy-Reed networks Molloy , is here valid because degree-correlations are essentially negligible for the scale-free graphs that we consider. In the realm of this often called “heterogeneous mean-field” or “degree-based mean-field” approximation, see e.g. Refs epidemicsRMP ; Vespignani12 and references therein, we therefore write $A_{ij}=k_{i}k_{j}/(N{\cal Q})$.

We now substitute (11), (12), (13) and (17) into (14), use the above heterogeneous mean-field approximation and the identity $\sum_{i}^{\prime}N^{-1}\eta_{i}=Nn_{k}\rho_{k}$. In the limit of $s\ll 1$, the transition rates, $T^{+}(\rho_{k})\equiv T^{+}_{k}$ and $T^{-}(\rho_{k})\equiv T^{-}_{k}$, thus become

[TABLE]

for the VM, while for the LD the transition rates are

[TABLE]

To remind the reader, $\omega$ is the degree-weighted density of cooperators, and is given by Eq. (3). In the limit of weak selection intensity ($0<s\ll 1$), the birth-and-death process defined by transition rates (18) or (19) MA10 ; AM10 is well described in terms of the multivariate forward and backward FPEs whose generators are respectively Gardiner ; Kimura ; Ewens

[TABLE]

V Timescale separation and effective diffusion approximation

Solving the multivariate FPEs associated with the generators (20) and (21) is a formidable task. Fortunately, the analysis greatly simplifies in the weak selection limit ($0<s\ll 1$) on which we focus our attention. This simplification stems from a separation of timescales which allows us to reduce the dynamics to that of an effective single-variate process VotNets11 ; VotNets12 ; VotNets13 ; VotNets21 ; VotNets22 ; VotNets23 ; PRL12 . Interestingly, timescale separation has also been used to simplify the analysis of the spread of epidemics on degree-heterogeneous networks, see e.g. Parra16 . Here, we first examine the case of the VM update, and then critically revisit the case of the LD discussed in Refs. PRL12 ; AM13 . ‘

V.1 Timescale separation & effective diffusion theory for the VM update rule

For the dynamics with the VM update, the quantity $\omega$ is conserved on a timescale $t\ll s^{-1}$ VotNets11 ; VotNets12 ; VotNets13 . This result, illustrated in Fig. 1, can be understood by computing the rate of change of $\delta\bar{\omega}=\sum_{i}k_{i}(1-2\eta_{i})/(N\mu_{1})$VotNets11 ; VotNets12 ; VotNets13 : $\delta\bar{\omega}/\delta t=\dot{\bar{\omega}}=\sum_{ij}A_{ij}[\bar{\psi}_{ij}-\bar{\psi}_{ji}]/(N\mu_{1})$, where the bar denotes the ensemble average. By using the mean-field approximation $A_{ij}=k_{i}k_{j}/(N{\cal Q})$, we find that $\dot{\bar{\omega}}=(a+d-b-c)s\bar{\omega}(1-\bar{\omega})(\bar{\rho}-\rho^{*})$. This means that on a timescale $t\ll s^{-1}$ the quantity $\omega$ is approximately constant.

Furthermore, on such a timescale, the mean-field rate equation for the subgraph density $\bar{\rho}_{k}$ satisfies $\dot{\bar{\rho}}_{k}=(T_{k}^{+}-T_{k}^{-})/n_{k}=\bar{\omega}-\bar{\rho}_{k}+{\cal O}(s)$. Hence, after a timescale $t={\cal O}(1)$, all $\bar{\rho}_{k}$’s converge to $\omega$ and we have $\bar{\rho}_{k}\approx\bar{\omega}\approx$ constant. Thus, the subgraph densities become independent of $k$. Clearly this implies that $\rho=\sum_{k}\rho_{k}n_{k}\approx\rho_{k}\approx\omega$. From the above, we therefore infer that when $t\gtrsim{\cal O}(1)$ we simply have $\rho_{k}\approx\rho\approx\omega$ and these quantities evolve together towards their common value $\rho^{*}$.

As confirmed in Fig. 1, at timescales $t\gtrsim{\cal O}(1)$ all the densities $\rho$ and $\rho_{k}$ converge to $\omega$, and reach the vicinity of $\rho^{*}$ after a timescale of ${\cal O}(s^{-1})$. This scenario lasts until a chance fluctuation eventually causes the fixation of either C or D on a much longer time scale ($t\gg s^{-1}$). In summary, as illustrated by Fig. 1, with the VM update rule we distinguish three timescales: (i) at $t\gtrsim{\cal O}(1)$, $\omega$ is approximately constant and $\rho$ and $\rho_{k}$ converge to it; (ii) at $t\gtrsim{\cal O}(s^{-1})$, $\rho_{k}\approx\rho\approx\omega$ converge to $\rho^{*}$, and (iii) at $t\gg{\cal O}(s^{-1})$ fixation occurs. The comparison of the top and bottom panels of Fig. 1, illustrates that the separation of the three timescales breaks down when $s$ becomes ${\cal O}(1)$. Overall, the effective diffusion approximation presented below is valid only in the limit of weak selection, i.e. for $0<s\ll 1$.

Since we are interested in metastability and fixation, which occur when $\rho_{k}\approx\rho\approx\omega$, we can legitimately approximate $\rho_{k}$ and $\rho$ by $\omega$ VotNets11 ; VotNets12 ; VotNets13 . We therefore substitute $\rho_{k}$ and $\rho$ by $\omega$ in the transition rates (18) and replace $\partial_{\rho_{k}}$ by $(kn_{k}/\mu_{1})\partial_{\omega}$ in the generators (20) and (21). Under selection of weak intensity, $0<s\ll 1$, this yields the single-variate effective forward and backward generators:

[TABLE]

where $\tilde{s}=(b+c-a-d)s$. Equations (22) and (23) are among the main results of this paper. It has to be noted that in the diffusion term we have neglected the subleading contributions on the order of ${\cal O}(s)$ since we are working in the weak selection limit.

The effective diffusion theory therefore predicts that the main influence of the scale-free topology under the VM update rule is to renormalize the population size $N$ into $N_{{\rm eff}}=N(\mu_{1})^{2}/\mu_{2}$ VotNets11 ; VotNets12 ; VotNets13 ; VotNets21 ; VotNets22 ; VotNets23 . Below we show that numerical simulations fully support this prediction, and we discuss how the effective population size $N_{{\rm eff}}$ affects the metastability and fixation properties of the ACGs. A relevant result for our analysis is the fact that for scale-free networks the degree distribution $n_{k}\sim k^{-\nu}$ with $\nu>2$, such that $\mu_{1}$ is finite, holds up to the maximum degree estimated to scale as $k_{{\rm max}}\sim N^{1/(\nu-1)}$222In principle, the algorithm that we have used generates scale-free graphs and allows for the existence of nodes of degree $k>k_{{\rm max}}$. Yet, as their number is negligible, these nodes play no statistical role, see Appendix A. It is also worth noting that we focus on scale-free graphs with exponent $\nu>2$ since these have a finite mean degree $\mu_{1}$. (size of the largest hub) KS02 . Therefore, the effective population size scales as VotNets11 ; VotNets12 ; VotNets13

[TABLE]

where we have introduced the exponent

[TABLE]

V.2 Timescale separation & effective diffusion theory for the LD update rule

A similar reasoning holds also for the LD update rule. Yet, here it is crucial to realize that the reference observable is the density of cooperators $\rho$. Indeed, under the LD $\rho$ is conserved on the time scale $t\ll s^{-1}$ VotNets11 ; VotNets12 ; VotNets13 ; VotNets31 ; VotNets32 as one can check by noting that $\delta\bar{\rho}=\sum_{i}(1-2\eta_{i})/(N\mu_{1})$ which, proceeding as above, yields $\dot{\bar{\rho}}=\sum_{ij}A_{ij}[\bar{\psi}_{ij}-\bar{\psi}_{ji}]/(N\mu_{1})=(a+d-b-c)s\bar{\rho}(1-\bar{\rho})(\bar{\rho}-\rho^{*})$. This indicates that $\rho$ remains constant when $t\ll s^{-1}$ and then relaxes to its metastable value $\rho^{*}$ on a timescale $t={\cal O}(s^{-1})\gg 1$; this is corroborated by Fig. 2.

Furthermore, at the mean-field level one obtains $\dot{\bar{\rho}}_{k}=(T^{+}_{k}(\bar{\rho}_{k})-T^{-}_{k}(\bar{\rho}_{k}))/n_{k}\sim(\bar{\omega}-\bar{\rho}_{k})k/\mu_{1}$, valid at times $t\gtrsim{\cal O}(1)$. Thus, as illustrated in Fig. 2, this means that on a time scale of $t\gtrsim{\cal O}(1)$, $\rho_{k}$ and $\omega$ approach $\rho$. Then, all these quantities evolve together towards their common metastable value $\rho^{*}$, which is reached at a timescale of $t={\cal O}(s^{-1})\gg 1$, and fluctuate around it before fixation occurs at a much later stage, see below.

To study metastability and fixation in ACGs, we can therefore restrict our attention in the regimes where $\rho_{k}\approx\omega\approx\rho$, and approximate $\rho_{k}$ and $\omega$ by $\rho$ VotNets11 ; VotNets12 ; VotNets13 . A one-body description of the dynamics is thus obtained by substituting $\rho_{k}$ and $\omega$ by $\rho$ in the transition rates (19) and by replacing $\partial_{\rho_{k}}$ by $n_{k}\partial_{\rho}$ in the generators (20) and (21). This yields the single-variate effective forward and backward generators, respectively:

[TABLE]

It has to be noted that these generators are the same as those obtained in a well-mixed population (on a complete graph) of size $N$; here, as opposed to Eqs. (22) and (23), the population size is not renormalized by the graph’s topology. This is similar to what was found for models without selection and with frequency-independent selection evolving with the LD VotNets11 ; VotNets12 ; VotNets13 ; VotNets21 ; VotNets22 ; VotNets23 ; VotNets31 ; VotNets32 . As discussed in more detail below, here we show that the metastability and fixation properties of ACGs evolving with the LD are the same in the leading order as those on a complete graph of size $N$. Therefore, the fixation probability and MFT of the ACGs with the LD are not affected by the scale-free topology, contrary to what was reported in Ref. PRL12 ; AM13 . Yet, a completely different scenario emerges with the VM update rule, where the scale-free topology strongly affects the long-time behavior.

VI Fluctuations in the metastable state

We can use the effective diffusion approximation to study the typical fluctuations in the long-lived metastable state at times $t\gtrsim s^{-1}\gg 1$. In this regime, the one-body forward FPE generators (22) and (29) allow us to describe the dynamics in the metastable state in terms of an Ornstein-Uhlenbeck process (OUP) obtained by linearizing the drift term about the metastable state $\rho^{*}$ in (22) and (29) and by evaluating the diffusion term at $\rho^{*}$.

VI.1 Fluctuations in the metastable state with the VM update

For ACGs evolving with the VM, in the realm of the effective diffusion approximation, the (forward) generator of the OUP for $\xi=\omega-\rho^{*}\approx\rho-\rho^{*}$ is therefore

[TABLE]

Here, we readily recognize the generator of an OUP with drift and diffusion parameters $\tilde{s}\rho^{*}(1-\rho^{*})$ and $2\rho^{*}(1-\rho^{*})/N_{\rm eff}$, respectively. Using the well-known properties of the OUP Gardiner , we readily verify that on average the deviations from the metastable state vanish: $\langle\xi(t)\rangle=\xi(0)e^{-\tilde{s}\rho^{*}(1-\rho^{*})t}\to 0$, where the square bracket denotes the ensemble average over the probability density $p(\xi,t)$ of the forward FPE $(\partial_{t}-\widetilde{{\cal G}}_{\rm f}(\xi))p(\xi,t)=0$ [with zero-flux boundary conditions]. This means that, on average, the deviations from $\rho^{*}$ vanish very quickly in the metastable state. Since $\langle\xi(t)\rangle=0$, the mean number of cooperators in the metasable state is $\langle N\rho\rangle=N\rho^{*}$ and therefore grows linearly with the number of nodes $N$. This result is confirmed in Fig. 3 (obtained from the analysis of stochastic simulations data, see Appendix B), showing that the mean number of cooperators in the metastable state increases linearly in $N$ with a slope $\rho^{*}$.

It is quite interesting to consider the variance of the variable $\xi$. Since initially, ${\rm var}(\xi(0))=\langle\xi^{2}(0)\rangle=0$, we find at times $t\gtrsim s^{-1}\gg 1$

[TABLE]

where $N_{\rm eff}$ is given by Eq. (27). Using (31), one can find the variance of the number of cooperators in the metastable state: ${\rm var}(N\rho)=N^{2}(\langle\rho^{2}\rangle-\langle\rho\rangle^{2})\approx N^{2}(\langle\omega^{2}\rangle-(\rho^{*})^{2})\approx N^{2}\langle\xi^{2}\rangle\sim N^{2}/(\tilde{s}N_{\rm eff})$. Hence, with the definition of $\alpha$ given by (28), we find

[TABLE]

In particular, we notice that for scale-free graphs with $2<\nu<3$ that are characterized by nodes of high degree NetRef , the variance increases faster than linearly in the system size since $2/(\nu-1)>1$. In other words, the effective diffusion theory predicts that $\log{({\rm var}(N\rho))}\sim 2/(\nu-1)\log N$ when $2<\nu<3$. This result is corroborated by extensive computer simulations outlined in Appendix B and illustrated in Fig. 4. Here, the theoretical prediction (35) is in very good agreement (within error bars) with the results of stochastic simulations over a broad range of values $2<\nu<3$. These results mean that the typical fluctuations in the metastable state on scale-free graphs with exponent $2<\nu<3$ are much stronger than on complete graphs where the cooperator variance scales as $N$, see below. In Section VII, we explore how the strong fluctuations arising on scale-free networks with exponent $2<\nu<3$ also dramatically affect the MFT.

VI.2 Fluctuations in the metastable state with the LD

For ACGs evolving with the LD, in the realm of the effective diffusion approximation, the (forward) generator of the OUP for $\xi=\rho-\rho^{*}$ is again given by (30) but with $N$ instead of $N_{\rm eff}$ on the right hand side. This crucial difference means that for the LD the diffusion parameter of (30) reads $2\rho^{*}(1-\rho^{*})/N$ and is therefore independent of the graph’s structure. The OUP drift parameter is still $\tilde{s}\rho^{*}(1-\rho^{*})$ and therefore coincides with that of the VM update rule. Proceeding as above, we readily find that on average the deviations from the metastable state still vanish exponentially in time, as $\langle\xi(t)\rangle=\xi(0)e^{-\tilde{s}\rho^{*}(1-\rho^{*})t}$. We can also compute the second moment of $\xi$: $\langle\xi^{2}(t)\rangle=(1-e^{-2\tilde{s}\rho^{*}(1-\rho^{*})t})/(\tilde{s}N)\sim 1/N$. Thus, the variance of the number of cooperators in the metastable state, ${\rm var}(N\rho)=N^{2}(\langle\rho^{2}(t)\rangle-\langle\rho(t)\rangle^{2})\approx N^{2}\langle\xi^{2}(t)\rangle\sim N$, scales similarly as on a complete graph. In stark contrast with the VM case on scale-free graphs with highly-connected nodes (when $2<\nu<3$), we thus find that the typical fluctuations in the LD case are not affected by the graph’s structure.

VII Fixation properties

The Moran-like evolutionary processes that we are considering are absorbing Markov chains and their fate is to reach one of the states corresponding to the entire network being populated only by cooperators or by defectors. Two important quantities to characterise the underlying evolutionary dynamics are therefore the (unconditional) mean fixation time (MFT), $T_{\rm fix}(\rho)$, which is the mean time to reach either of the absorbing states from an initial density $\rho$ of cooperators, and the fixation probability $\phi^{C}(\rho)$ that cooperation prevails (the final state is all-C, with the extinction of all D’s) starting from a fraction $\rho$ of cooperators. It is well known that on complete graphs, for ACGs these quantities scale exponentially with the population size: $\ln{T_{\rm fix}(\rho)}\sim sN$ and $\ln{\phi^{C}(\rho)}\sim-sN$ or $\ln{(1-\phi^{C}(\rho))}\sim-sN$, see e.g. Antal06 ; MA10 ; AM10 . Here, we use the effective diffusion theory to compute $T_{\rm fix}(\rho)$ and $\phi^{C}(\rho)$ on scale-free graphs and, together with large-scale computer simulations, to uncover under which circumstances, the complex topology affects the MFT.

While we are mainly interested in the MFT, it is useful to start our discussion by outlining the calculation of the fixation probability. In the realm of the effective diffusion approximation, by assuming that fixation occurs after lingering for a long time in the metastable state (see Figs. 1 and 2), we use the backward operators (23) and (29) to compute the fixation probability which satisfies $\widetilde{{\cal G}}_{\rm b}(\rho)\phi^{C}(\rho)=0$ with boundary conditions $\phi^{C}(0)=1-\phi^{C}(1)=0$ weaksel1 ; weaksel2 ; weaksel3 ; Gardiner , and by now approximating $\omega\approx\rho$ in (23). This yields $\phi^{C}(\rho)=\frac{{\rm erfi}\left[\rho^{*}\sqrt{\sigma}\,\right]-{\rm erfi}\left[(\rho^{*}-\rho)\sqrt{\sigma}\,\right]}{{\rm erfi}\left[\rho^{*}\sqrt{\sigma}\,\right]+{\rm erfi}\left[(1-\rho^{*})\sqrt{\sigma}\,\right]}$ PRL12 , where ${\rm erfi}(z)\equiv\frac{2}{\sqrt{\pi}}\int_{0}^{z}e^{u^{2}}du$ and

[TABLE]

This clearly indicates that the VM/LD lead to very different fixation properties on scale-free networks with nodes of high degree: For the VM, the topology yields an effective population size $N_{\rm eff}$ (27) leading to $\sigma\ll\tilde{s}N$ when $2<\nu<3$. The dependence of $\phi^{C}$ on the system size when $2<\nu<3$ is thus a stretched exponential with exponent $-N^{\alpha}$ and $\alpha<1$, whereas the fixation probability with the LD (and for the VM with $\nu>3$) is the same as on a complete graph of size $N$.

In addition to the fixation probability, we can similarly compute the MFT, $T_{\rm fix}(\rho)$. This is done by solving $\widetilde{{\cal G}}_{\rm b}(\rho)T_{\rm fix}(\rho)=-1$ with boundary conditions $T_{\rm fix}(0)=T_{\rm fix}(1)=0$ Gardiner . Using standard methods Kimura ; Ewens ; Gardiner , the solution to this inhomogeneous backward FPE is given in Ref. PRL12 and to leading order we find:

[TABLE]

Hence, with the expression of $\phi^{C}$, when $\rho^{*}<1/2$ and $\rho>\rho^{*}$ this gives

[TABLE]

This shows that when initially the fraction of cooperators is not too low, metastability occurs prior to fixation and the leading contribution to the MFT (42) is independent of the initial condition AM1 ; AM2 ; AM3 ; AM4 ; MA10 ; AM10 ; AM16 . It stems from this result that fixation occurs much more rapidly for the VM on scale-free networks with $2<\nu<3$ than with the LD. In the VM case, the MFT grows as the stretched exponential $\ln{T_{\rm fix}}\sim N^{\alpha}\ll N$ with the exponent $\alpha$ given by (28). This theoretical prediction is confirmed by our stochastic simulations (within error bars) from which we have computed the exponent $\alpha$ reported in Fig. 5 (see Appendix B for technical details). The results reported in Fig. 6 corroborate, within error bars, the theoretical prediction (42), $\ln\left(d\ln\left(T_{{\rm fix}}\right)/ds\right)\sim[2(\nu-2)/(\nu-1)]\ln{N}$ when $2<\nu<3$. A similar analysis with LD, reported in Fig. 7, see Appendix B, confirms our theoretical result (42) that the MFT of ACGs evolving with the LD grows purely exponentially with $N$ as on complete graphs Antal06 ; MA10 ; AM10 .

We therefore emphasize that, contrary to what stated in Refs. PRL12 ; AM13 , $\phi^{C}$ and the MFT of ACGs evolving with the LD are not affected by the scale-free topology. Therefore, in ACGs with the LD, $T_{\rm fix}$ always grows exponentially with $\tilde{s}N$ (within error bars) as illustrated in Fig. 7. Our analysis of the LD in Refs. PRL12 ; AM13 was based on relatively small systems (up to $N=4,000$) which apparently were governed by finite-size effects leading to a nonlinear dependence on $N$ when $2<\nu<3$. Two of us tried to explain these effects, that turn out to disappear at large-enough $N$, in terms of an effective diffusion theory derived by assuming that $\omega$ was approximately conserved 333In Ref. AM13 the fixation probability of coordination games subject to the LD was considered similarly.. Even though the difference between $\omega$ and $\rho$ is negligible for $s\ll 1$ (see Fig. 2), considering $\omega$ to be constant instead of $\rho$ leads to a large error in the MFT. Here, by carrying out extensive simulations on large systems (up to two orders of magnitude larger than in PRL12 ; AM13 ), and by correcting our theoretical analysis, we have confirmed that the fixation properties of ACGs evolving with the LD exhibit a linear dependence on $N$ which, as seen above, is consistent with the (approximate) conservation of $\rho$ by the LD under weak selection (see Fig. 2).

VIII Summary & Conclusion

In this work we have studied the interplay between scale-free topology, demographic fluctuations and frequency-dependent selection on the dynamics of anti-coordination games. This class of paradigmatic games, particularly relevant in theoretical biology, between two competing types (cooperators and defectors) is characterized by their long-lived coexistence, and eventually by the fixation of one type. These features are here investigated by combining analytical methods and extensive stochastic simulations. Since evolutionary dynamics on heterogeneous graphs is known to depend on the underlying microscopic update rule, we have investigated two types of individual-based dynamics: according to the death-first/birth-second voter model (VM), or according to the links dynamics (LD) in which birth/death or death/birth events occur randomly between connected agents. While we here specifically focus on scale-free graphs, it is noteworthy that our approach is valid for any degree-heterogeneous networks.

Our analytical approach, valid in the weak selection limit, is based on a diffusion approximation for the density of cooperators at nodes of a prescribed degree. The resulting multi-variate Fokker-Planck equations are greatly simplified by exploiting a timescale separation: After a transient, when the selection pressure is weak, the multi-body dynamics towards fixation can be described in terms of the (approximately conserved) cooperator degree-weighted density (for the VM) and the cooperator density (for the LD). As a result, the typical fluctuations in the number of cooperators in the coexistence metastable state, and the fixation properties can be determined by means of effective single-variate forward and backward Fokker-Planck equations.

In particular, with the VM update rule the complex scale-free topology is responsible for an effective reduction of the population size $N$ when $2<\nu<3$, into $N_{{\rm eff}}=N^{\alpha}$ with an exponent $0<\alpha<1$ that depends on the first two moments of the graph’s degree distribution $n_{k}\sim k^{-\nu}$. This results in larger typical fluctuations and a superlinear dependence on $N$ of the variance of the number of cooperators in the metastable state. As an outcome, the MFT is exponentially decreased and displays a non-trivial stretched-exponential dependence on the population size. In the absence of “hubs” ($\nu>3$), in the leading exponential order, the dynamics is independent on the scale-free topology and the MFT scales linearly with $N$ as in the well-mixed dynamics.

When the system evolves according to the LD, the scale-free topology does not effectively renormalize the population size and we find $N_{{\rm eff}}=N$ (for $\nu>2$). With the LD, the variance of the number of cooperators in the metastable state grows linearly with $N$ and the MFT displays a pure exponential dependence on $N$ as shown in Fig. 7, similarly to the well-mixed case. These results are in contrast with what was reported in PRL12 ; AM13 . The correction here was made possible due to a subtle amendment that we have made in the theoretical description accompanied by an efficient numerical code that allowed us to attain very large system sizes and hence to carefully analyze the dependence of the fixation properties on $N$.

IX Acknowledgments

DS and MA were supported by Grant No. 300/14 of the Israel Science Foundation. MM is grateful for the hospitality of the LSFrey at the Arnold Sommerfeld Center (University of Munich) where part of this work was done, as well as for the financial support of the Alexander von Humboldt Foundation by Grant No. GBR/1119205 STP.