Evolution of cooperation in deme-structured populations on graphs

Understanding how cooperation can evolve in populations despite its cost to individual cooperators is an important challenge. Models of spatially structured populations with one individual per node of a graph have shown that cooperation, modeled via the prisoner’s dilemma, can be favored by natural selection. These results depend on microscopic update rules, which determine how birth, death and migration on the graph are coupled. Recently, we developed coarse-grained models of spatially structured populations on graphs, where each node comprises a well-mixed deme, and where migration is independent from division and death, thus bypassing the need for update rules. Here, we study the evolution of cooperation in these models in the rare migration regime, within the prisoner’s dilemma. We find that cooperation is not favored by natural selection in these coarsegrained models on graphs where overall deme fitness does not directly impact migration from a deme. This is due to a separation of scales, whereby cooperation occurs at a local level within demes, while spatial structure matters between demes.


I. INTRODUCTION
Cooperation between individuals is commonly observed in living systems, at the cell scale [1][2][3][4], between animals [5][6][7][8], and in human societies [9][10][11].For instance, bacteria often produce diffusible resources (e.g.enzymes, toxins or signaling molecules) that can benefit neighboring bacteria, including those that do not produce this resources (cheaters) [12,13].Because producing these resources is generally costly, natural selection a priori works against the producers of public goods, and favors cheaters.Therefore, understanding how cooperation can evolve and spread is an important question [13,14], which has led to intense debates.Kin selection can lead to the evolution of cooperation under some conditions [15,16]: if the individuals that benefit from cooperation are related to cooperators, then cooperation can be selected even if it carries a cost.Spatial structure can also favor cooperation in some models [13,[17][18][19][20][21].Since spatial structure generally entails that closely located individuals are genetically related, its ability to promote cooperation is related to kin selection [22].In evolutionary game theory, individuals play games such as the prisoner's dilemma [23,24] or the snowdrift game [23,25] using fixed strategies.Their fitnesses depend on the composition of the surrounding population, via the payoff they receive from the games played with their neighbors [26][27][28][29].This allows to model cooperation in a straightforward way and to study its evolution.In these models, each individual is located on a node of a graph [30], allowing to investigate the impact of complex spatial structures on the evolution of cooperation [17,18].It was found that in the prisoner's dilemma game, a cooperator mutant can be favored by natural selection in some graphs under specific update rules that prescribe how individuals are chosen to divide, die, and replace each other [17,18].
Recently, we developed more general models of spatially structured populations on graphs, where each node comprises a well-mixed deme instead of a single individual, and where migration between demes is independent from division and death [31,32].These models therefore do not require an update rule, contrary to other models with demes on graphs [33][34][35].We showed that for mutants with constant (non-frequency-dependent) fitness, some spatial structures with specific migration asymmetries can amplify natural selection in the regime of rare migration [31].We further found that suppression of selection is pervasive with more frequent migrations, in the branching process regime [32].Importantly, these results do not depend on update rules, contrary to models with one individual per node where suppression or amplification of natural selection for a given graph strongly depends on the update rule choice [36][37][38][39].In the rare migration regime, this dependence on update rules was replaced by a dependence on migration asymmetry, which can be experimentally tuned and measured, while update rules model assumptions about the microscopic dynamics at the cell scale, which are generally not strictly followed by natural or experimental populations.
Can spatial structure favor the evolution of cooperation in these deme-structured models, as in models with one individual per node of the graph?To address these questions, we consider the prisoner's dilemma in spatially structured populations within the frameworks of [31,32].In these models, the fate of mutants is impacted by the choice of graph structure and specific migration rates between demes, but the average fitness of a deme does not directly impact the intensity of its outgoing (or incoming) migrations.In other words, selection is essentially soft [40], as we aimed for a minimal incorporation of spatial structure.Note however that cooperation can be favored via a coupling between composition and population size [13,[19][20][21].
We first lay out our models of deme-structured populations on graphs with cooperation in the case of the prisoner's dilemma.Next, we consider a single cooperator (mutant), placed in a population entirely composed of defectors (wild-types), and we ask what the probability that it takes over (i.e.fixes) is.We address this question in different graphs, for rare migrations.We compare the fixation probability of a single cooperator to that of a single defector, and to that of a neutral mutant, to determine whether cooperation is favored or not.We also evaluate the impact of structure by comparing the fixation probability of a cooperator mutant in a structure to that in a well-mixed population with the same total size.We find that cooperators are not favored by natural selection in our models.The impact of spatial structure on their fixation probability reduces to that observed for mutants with fixed fitness disadvantage.This is due to a separation of scales, with cooperation occurring within demes, while spatial structure matters between demes.We make complete comparisons with models with one individual per node of the graph, including formal calculations in the Appendix, where we show the points that differ between these two classes of models.

A. Spatially structured populations with demes on graphs
We model spatially structured populations on directed graphs.Each of the D nodes of the graph contains a well-mixed subpopulation or deme.The specific graphs considered here are shown in Fig. 1.Two types of individuals make up the population: mutants (M ) and wildtypes (W ).We compare two different models, which differ by how deme size is regulated.In the first one, each individual can divide, die and migrate with given rates, deme sizes are limited by a carrying capacity, and we consider the regime where they fluctuate around a steady-state value [31].In the second one, inspired by serial transfer experiments, the population undergoes an alternation of growth and bottleneck phases.Starting from a bottleneck state, each deme undergoes exponential growth, and then some individuals are sampled and migrate to form a new bottleneck state [32].Formally, the first model is close to a Moran model, but deme size is not strictly fixed [31].Meanwhile, the second one is close to a Wright-Fisher model [32].We mainly study the first one, which is more directly comparable to models with one individual per node [17,18,30].The second model allows us to extend our results to more frequent migrations.
Within our models of spatial structure, we consider an evolutionary game theory setup that provides a simple description of cooperation, following [17,18,28].Our two types of individuals then correspond to cooperators and defectors, interacting within a prisoner's dilemma with payoff matrix where b > c > 0. If a cooperator interacts with another cooperator, then both receive a benefit b but also lose a cost c; this yields b − c in the upper left corner of P .If a cooperator interacts with a defector, then the cooperator does not receive anything, but it loses c because of the interaction (upper right corner), while the defector receives b (bottom left corner).Finally, if two defectors interact, both get 0. Note that here, being a cooperator or a defector in our model is an intrinsic (genetic) characteristic of an individual.Individual fitnesses are expressed as f = 1 − w + wπ where w is the intensity of selection and π the total payoff received by the individual upon interactions with all its partners.We will usually consider cooperators as mutants (M ) and defectors as wild-types (W ), and ask how likely it is that cooperator mutants take over.We consider that interactions between individuals occur inside demes, which are well-mixed.Thus, each individual interacts with all other individuals in the same deme.Denoting by k the number of mutants in a deme of size N , the payoff π M (k, N ) for a mutant individual and π W (k, N ) for a wild-type individual in this deme read [28]: In these expressions, (k−1)/(N −1) is the fraction of mutants among the individuals that a mutant of focus can interact with (i.e.all but itself), while (N − k)/(N − 1) is the fraction of wild-types among them.Similarly, k/(N − 1) is the fraction of mutants among the individuals that a wild-type of focus can interact with.In each of our two models, fitness can be expressed from these payoffs, giving rise to frequency-dependent fitnesses (see below).

B. Model with carrying capacities
We first consider the model introduced in [31].In this model, each deme has maximum carrying capacity K. Besides, each individual has fitness f M or f W , and death rate g M = g W = g, as we focus on selection on division.Here, fitness represents the maximal division rate of microorganisms, reached in exponential growth.The division rate of individuals of type A = M, W in deme i is given by the logistic function f A (1 − N i /K), where N i is the number of individuals in deme i.As in [31], we focus on the rare migration regime, where fixation or extinction within a deme occurs on much shorter time scales than migrations (see Fig. 2).In this regime, we describe the evolution of the microbial population as a Markov process where each event is a migration and where each deme is either fully mutant or fully wild-type [31,41].Note that we will present an extension to more frequent migrations using our second model (see below).We focus on the regime where deme sizes N i fluctuate weakly around their deterministic steady-state value K(1 − g/f W ) (where we have neglected differences between f M and f W , which is acceptable for weak selection, and where we have neglected migrations, which is acceptable if they are rare).Note that the deterministic steady-state size is equal to the maximum carrying capacity K if g = 0.In what follows, we will refer to K as the deme carrying capacity, for simplicity.We approximate the fixation process of a type within a deme by a Moran process [42,43] for a population of size K(1 − g/f W ), thus neglecting size fluctuations.Here, we implement logistic regulation in the division rate, which decreases if N i grows.Appendix C briefly discusses a variant where it is instead implemented in the death rate.Our conclusions are not affected.
Denoting by k the number of mutants in a deme of size N as above, the frequency-dependent fitness f M (k, N ) for a mutant individual and f W (k, N ) for a wild-type individual in this deme read: where payoffs π M and π W are given by Eq. (2).Under this convention, where fitnesses depend linearly on payoffs, we will focus on the weak selection regime defined by w ≪ 1 and wN ≪ 1, where both f M and f W are close to the base fitness f = 1 [17,18].We will then briefly discuss an alternative convention where fitnesses depend exponentially on payoffs [44], which allows an extension of our results beyond the weak selection regime.
FIG. 2. Fixation of a mutant in a clique.We start with a single mutant individual (yellow) in one deme, while all the other individuals are wild-types (red).In the rare migration regime, mutants first fix in the initial deme before one migration event to another deme occurs.This process is repeated until mutants take over the complete structure.
We perform analytical calculations and stochastic simulations.For the latter, we use a Gillespie algorithm [31,45,46], described in detail in Appendix D. We will refer to this model in figures as the "carrying-capacity" (C.C.) one.

C. Model with serial dilutions
Recently, in Ref. [32], we introduced a new model of spatially structured populations on graphs, directly inspired by the batch culture setups with serial transfers that are used in many evolution experiments [47][48][49][50][51][52][53][54] and are important in ecology [55].This serial dilution model is a variant of the Wright-Fisher model, and features demes placed on the nodes of a graph with possible migrations along its edges, like the carrying-capacity model defined above.Each deme holds a well-mixed subpopulation, and within-deme interactions implement cooperation between individuals following the payoff matrix (1).In this model, all demes have the same bottleneck size B, and repeatedly undergo a two-step process.To clarify notations, quantities common to the two models will be denoted with a prime for the serial dilution model.1.First, subpopulations grow exponentially in each deme during a fixed growth time t.Time within this growth phase is denoted by τ ∈ [0, t].Starting at τ = 0 from k ′ (0) mutants and l ′ (0) wild-types in a deme at bottleneck, such that k ′ (0) + l ′ (0) = N ′ (0) = B, these numbers grow as: with rates f M (τ ) and f W (τ ) reading Here π M (τ ) is a shorthand for the payoff π M (k ′ (τ ), N ′ (τ )) given in Eq. ( 2), with population size N ′ (τ ) = k ′ (τ ) + l ′ (τ ) [and similarly for π W (τ )].The growth rates in Eq. ( 5) are chosen to have a similar form to the fitnesses given in Eq. ( 3).
Here too, the payoffs obtained from the prisoner's dilemma within demes directly impact the ability of cells to divide.
2. Then, a dilution and migration step occurs (see details in Section B of the Appendix), bringing the system to a new bottleneck state.For each deme, we perform a multinomial sampling to pick exactly B individuals that form the new bottleneck state of this deme.These incoming individuals are sampled from all demes, according to migration probabilities (all incoming migration probabilities to a given deme sum to 1).
In addition to bridging with experiments, this model has the advantage of facilitating the treatment of frequent migrations.Indeed, in our carrying capacity model under frequent asymmetric migrations, some demes can have a steady-state size that exceeds K, which prevents division there.Choosing a division rate independent from N and a death rate proportional to N/K would instead result in substantially faster turnover in these demes.These situations lack realism and strongly depend on the exact form of deme size regulation.The serial dilution model does not suffer from this drawback.We perform stochastic simulations under this serial dilution model [32], and we refer to it in figures as the "dilution" one.Considering the serial dilution model allows us to test the robustness of our conclusions from the carrying-capacity model in the rare migration regime, and to extend our study beyond this regime.
To quantitatively compare results from the two models, we choose the bottleneck size B in the dilution model to be equal to the steady-state size of the demes K(1 − g/f W ) in the carrying-capacity model.Furthermore, in the dilution model, we introduce the effective intensity of selection 2w ′ t.Indeed, it plays the same role as the intensity of selection w in the carrying-capacity model, see Section B 3 of the Appendix (in particular, the factor of 2 comes from the difference between a Moran and a Wright-Fisher description in the diffusion limit).

A. Fixation probability in a deme and in a structured population with rare migrations
In the rare migration regime, fixation of a mutant in the whole population starts by fixation in the deme where the mutant was introduced (see Methods and Fig. 2).Thus, we can write the fixation probability ρ struct M of a mutant in the structure as where ρ M is the fixation probability of one mutant in a well-mixed deme where all other individuals are wildtype, while Φ 1 is the fixation probability of the mutant type in the whole population starting from one fully mutant deme and D − 1 wild-type demes.
Let us first focus on ρ M , and let us compute it in the carrying-capacity model.Neglecting deme size fluctuations in the fixation process within a deme (see Methods) and denoting by N the number of individuals in the deme, we can follow the calculation of [28], which yields: Note that N corresponds to the deterministic steadystate deme size N = K(1 − g), where we have neglected the impact of fitness differences on N and used f W ≈ f M ≈ 1, as we focus on the weak-selection regime (this amounts to neglecting terms of order wg/(f − g), which are explicitly written in [31] -this is acceptable if g ≪ f , which is realistic).Similarly, the fixation probability ρ W of one wild-type in a well-mixed deme where all other individuals are mutant reads Inserting the fitness expressions from Eq. (3) into Eqs.( 7) and ( 8), and performing an expansion in w in the weak selection regime w ≪ 1 and wN ≪ 1 yields: We note that ρ M < 1/N < ρ W : we recover the result that in a well-mixed population, natural selection favours defectors over cooperators.
To calculate the fixation probability ρ struct M of a mutant in the whole structured population, we need to express the probability Φ 1 that the mutant type fixes, starting from one fully mutant deme, see Eq. ( 6).Importantly, Φ 1 depends on the specific graph that is considered for the structured population.Below, we focus on the simple graphs shown in Fig. 1, and investigate how each of them impacts the evolution of cooperation.To compute Φ 1 , we notice that in the rare migration regime, the state of the population can be fully described by specifying which demes are fully mutant and which ones are fully wildtype.Furthermore, this state evolves due to migration events, which may change it if they result into fixation in the destination deme.The evolution of the state of the population is a Markov chain [31,41].

B. Evolution of cooperation in the clique and cycle
Let us first consider the clique and the cycle, where all demes are equivalent.To compute Φ 1 , and more generally the fixation probability Φ i when starting from i fully mutant demes (consecutive in the case of the cycle) and D − i fully wild-type demes, we follow the method we used in [31], which is based on [28].Briefly, we write down a recurrence relation on Φ i by discriminating over all possible outcomes of the first migration event: either i increases by one if a mutant migrates to a wild-type deme and fixes there, or it decreases by one in the opposite case, or it stays constant in all other cases (see Appendix A 1 and [31] for details).Solving this recurrence relation yields with The reasoning made with frequency-independent fitness in [31] thus extends to the present case with cooperation, and accordingly, the formal expression of Φ 1 is the same in both cases.However, the expression of γ differs because the fixation probabilities within a deme are impacted by the cooperation model.In the present case, γ can be expressed using Eqs.( 9) and (10).Note that Eq. ( 11) is independent from migration rates.In particular, in the case of the cycle, Eq. ( 11) does not depend on the ratio m A /m C of the anticlockwise migration rate m A to the clockwise one m C , see Fig. 1.
In [31], we further showed that all circulation graphs (such that the sum of incoming migration rates to any given deme is equal to the sum of outgoing ones from that deme), including the cycle, have the same fixation probability, thereby extending the circulation theorem which holds for graphs with one individual per node [30].The circulation theorem further extends to the present case with cooperation.Indeed, the exact same proof as in [31] holds, albeit with a different γ.Therefore, all circulation graphs have the same fixation probability of cooperator mutants, given by Eqs. ( 11) and ( 12).We will henceforth denote it by Φ circ. 1 .Using Eqs. ( 9) and ( 10) to express Eq. ( 11) yields which gives the fixation probability of a mutant in the whole population, using Eq. ( 6) and Eq. ( 9): Importantly, this fixation probability is smaller than the neutral value 1/(N D) for all b > c > 0: cooperation is never favored by natural selection in circulation graphs, including the clique or the cycle.It is also very similar to the fixation probability of a mutant in a well-mixed population of size N D [obtained from Eq. ( 9) by replacing N with N D] -both coincide for N ≫ 1 (see also [31]).In Fig. 3(a), we compare numerical simulation results and analytical predictions for the fixation probability ρ cycle M in the cycle.We find that stochastic simulation results do not depend on migration asymmetry α = m A /m C , consistently with our analytical calculations of A slight discrepancy can however be observed between analytical predictions and simulation results.It arises from the constant-size approximation we made to compute the analytical expressions of ρ M and ρ W , which are both involved in the expression of ρ cycle M .Indeed, using Gillespie simulation results obtained on well-mixed demes for ρ M and ρ W instead of their analytical expressions for fixed size within the expression yields an excellent agreement with simulation results [see Fig. 3(a)].
In Fig. 3(b), we explore more frequent migrations, using the serial dilution model introduced above.Our simulation results suggest that the fixation probability for a cycle remains the same for more frequent migrations.Besides, as expected, matching the effective intensity of selection 2w ′ t in this model to w in our first carryingcapacity model, simulations from both models give close results for rare migrations.The rare migration analytical prediction for the dilution model is obtained with Eq. ( 11), but replacing ρ M and ρ W with their expressions ρ ′ M and ρ ′ W derived in the dilution model (see details in Section B of the Appendix).Note that [32] extended the circulation theorem to the dilution model without cooperation, both for rare migrations, and for frequent migrations within the branching process approximation.
Our results stand in contrast to those of [18], where the present model of cooperation was studied in a cycle graph with one individual per node.There, it was found that cooperators could be favored by natural selection under the death-Birth (dB) and imitation update rules, but not under the Birth-death (Bd) update rule (the uppercase B indicates that selection on fitness happens at the birth step).Specifically, natural selection was found to favor cooperation in the dB case if b/c > 2 + 4/(D − 4) [18].To understand the formal origin of this difference, we revisit the derivation of [18] in Appendix A 1. The key difference is that in our model, cooperative interactions occur within a deme and are not involved in migration events.Conversely, in the model of [18], they impact the spread of mutants on the graph since birth, death and migration events are coupled via update rules.

C. Evolution of cooperation in the star
Let us now consider the star graph, where the center sends out individuals to any leaf with rate m O , and receives individuals from any leaf with rate m I [see Fig. 1(c)].This graph structure has been extensively studied for mutants with constant fitness.In evolutionary graph theory with one individual per node, the star amplifies or suppresses natural selection depending on the update rule [30,[36][37][38].For instance, it is an amplifier of selection using the Bd update rule, but a suppressor with the dB rule.In evolutionary graph theory with nodes containing subpopulations, the evolutionary outcome also depends on the specific update mechanism [35].In our model with rare migrations [31] that does not rely on a specific update rule, the star's behavior depends on migration asymmetry α = m I /m O , as α > 1 amplifies selection while α < 1 suppresses it.This same outcome holds for the rare migration regime in the serial dilution model [32].However, for more frequent migrations, the star becomes a suppressor for all migration asymmetries in the branching process regime [32].How does the star graph impact the fixation of a cooperator, in the rare migration regime and beyond?
To address this question, we extend the calculation of [31] to our cooperation model, and we express the fixation probability Φ star 1 starting from one fully mutant deme, uniformly chosen among all demes, in the star graph.As in [31], we obtain where γ is given by Eq. ( 12), while α = m I /m O is the ratio of the migration rate incoming to the center from a leaf m I , to that outgoing from the center to a leaf m O , see Fig. 1.Eq. (15) shows that migration asymmetry α directly impacts the fixation probability in the star.In [31], we demonstrated that the star amplifies natural selection for α > 1 and suppresses it for α < 1.For α = 1, the fixation probability reduces to Eq. ( 11), as the star is then a circulation graph.
Using Eqs. ( 9) and ( 10) to express Eq. ( 15) yields Studying the ratio between Eq. ( 16) and Eq.(13) shows that the fixation probability of a cooperator mutant is larger in the star than in the clique for α < 1 (i.e. when the star behaves as a suppressor of selection [31]) and smaller for α > 1 (i.e. when the star behaves as an amplifier of selection [31]).This result is confirmed by stochastic simulations in Fig. 4(a).These effects are purely due to the population structure and not to the cooperation model.In other words, what we find here holds as well for a deleterious mutant without frequency-dependent fitness (see [31]).Simulations in Fig. 4(b) show that cooperation does not impact the star's behavior in the dilution model either, compared to the case of mutants with constant fitness [32].Indeed, the rare migration regime recovers results from the carrying-capacity model, but as migrations become frequent, the star becomes a suppressor of natural selection regardless of migration asymmetry [32].

D. Evolution of cooperation in the square lattice
In [18], a general rule for cooperation on graphs was found for the model with one individual per node, where interactions are defined by a specific update rule.The key result of [18] is that under the death-Birth (dB) update rule, natural selection favours cooperators (i.e. mutant individuals here) if b/c > k, where k here is the average degree of the graph (i.e., the average number of neighbours per node).No such effect exists under the Bd update rule [18].In Appendix A 2, we revisit stepby-step the derivation of [18], which uses a pair approximation, and adapt it to our model.We find that the fixation probability Φ 1 starting from one fully mutant deme, given in Eq. (A25), does not involve the average degree k and that natural selection never favors cooperation.Again, the key differences are that in our model, birth, death and migration are all independent, and that cooperation occurs within demes.
The proof of [18] is exact for Bethe lattices, i.e. infinite graphs without cycles where each individual has exactly k neighbors.Graphs possessing this last property are said to be regular.Because Bethe lattices pose some difficulties for simulations, finite regular graphs, including square lattices, were simulated in [18], yielding good overall agreement with the prediction that cooperators are favored if b/c > k.Thus motivated, here, we study the square lattice with D = 9 demes shown in Fig. 1, adding periodic boundary conditions, so that each deme is connected to exactly 4 other demes.This graph is a circulation, since for each deme, the incoming migration rates sum to m N + m S + m E + m W , and the outgoing ones too.Therefore, we predict that the fixation probability of a mutant in the rare migration regime is the same as in other circulation graphs (see Section III B), and does not depend on migration asymmetry.Fig. 5(a) shows that our simulation results are in good agreement with this prediction.In addition, this result seems to extend to frequent migrations in the dilution model, see Fig. 5(b).

E. Beyond weak selection
So far, we assumed that fitnesses depend linearly on payoffs [see Eq. ( 3)], and we remained in the weak selection regime where w ≪ 1.This allowed us to give simple analytical expressions of fixation probabilities in demes, and to conduct a direct comparison with Refs.[17,18].
It is analytically difficult to go beyond weak selection if fitnesses depend linearly on payoffs.However, Ref. [44] showed that calculations beyond the weak se-lection regime simplify if fitnesses depend exponentially on payoffs.This amounts to writing where w is a parameter that determines the intensity of selection (as w in the linear case), instead of Eq. ( 3).Payoffs π M (k, N ) and π W (k, N ) remain expressed by Eq. ( 2).Using the same reasoning as before, we obtain the same expressions for the probabilities of fixation of cooperator mutants in subdivided populations for rare migrations, namely Eq. ( 11) for the clique, the cycle, and all circulations, and Eq. ( 15) for the star.What changes is the expression of the fixation probabilities within demes, which now read [44] and Their ratio, which enters Eq. ( 11) and Eq. ( 15), is: where we neglected the dependence of steady-state deme size on fitness, which is acceptable if w ≪ 1 (see Methods and Ref. [31]).Hence, this exponential convention for fitnesses yields simple analytical expressions of fixation probabilities for rare migrations in our model of spatially structured populations.Here, we assumed w ≪ 1 but not N w ≪ 1: in this sense, these results hold beyond weak selection.In addition, the first-order expansions of Eqs. ( 18) and ( 19) when w ≪ 1 and N w ≪ 1 exactly match Eqs. ( 9) and ( 10), respectively, with w standing in for w.Thus, the exponential and the linear conventions for fitnesses are equivalent in the weak selection regime.Importantly, because only the fixation probabilities within demes are modified by changing from linear to exponential fitnesses, our main conclusions remain valid in the exponential case.This includes our results for circulations and stars which are based respectively on Eq. ( 11) and Eq. ( 15), as well as our reasoning for general graphs based on the pair approximation in Appendix A 2. Our main result that natural selection does not favor cooperation in our model of spatially structured populations therefore extends beyond weak selection, at least to the regime where w ≪ 1 without needing N w ≪ 1.

IV. DISCUSSION
In all structures we considered, the fixation probability of a cooperator in a population of defectors is smaller than vice-versa.Thus, cooperation is never favored in our model.The reason for this is that cooperation only takes place within nodes.Migration events do not directly involve cooperation or any interaction between individuals from different demes, not even competition.Cooperation is therefore a local interaction.In contrast, cooperative interactions happen at the level of the graph in [17].Indeed, each individual (placed on one node) interacts with its nearest neighbors on the graph, which affects its fitness and thus its probability of being picked at the birth step.This separation of scales entails that cooperation cannot be favored by natural selection in our model, including in other graphs, complex networks [56], and structures with multiple network layers [57].
More formally, our model decouples migration, birth, and death events.In the rare migration regime, fixing in one deme and spreading to the whole spatially structured population happen on two different time scales.Because of this, fitnesses enter overall fixation probabilities only through the fixation probabilities in one wellmixed deme, namely ρ M and ρ W . Assuming that fitnesses depend linearly on payoffs, in the weak selection regime, we can compare Eqs. ( 9) and ( 10) to the fixation probabilities that would be obtained in the frequencyindependent case, i.e. f M = f W (1 − ϵ) with a slightly deleterious mutant with fitness disadvantage ϵ satisfying 0 < ϵ ≪ 1, N ϵ ≪ 1.We get a perfect mapping when choosing ϵ = w[c + b/(N − 1)].Moreover, assuming that fitnesses depend exponentially on payoffs as in Section III E, Eqs. ( 18) and ( 19) also match the frequencyindependent fixation probabilities obtained in that convention with f M /f W = exp {− w [c + b/(N − 1)]}, and this holds beyond weak selection, assuming w ≪ 1 without needing N w ≪ 1.Thus, the fixation probabilities ρ M and ρ W have the same expression as with constant fitnesses, as long as we pick an effective fitness disadvantage for the mutant, which depends on the parameters of cooperation.This confirms that at the level of the whole structure, the frequency-dependence of the model does not matter.Because of this, cooperation cannot benefit from the specificities of the structure.Thus, our previous results regarding the impact of spatial structure on mutant fixation for constant fitnesses [31,32] extend to the present model with cooperation.
In our framework, selection is essentially soft, i.e. the contributions of demes are not affected by their average fitnesses [40].Indeed, in our dilution model, the total size of each deme after growth does not impact its contribution to the next bottleneck state.In our carryingcapacity model, the carrying capacity of each deme is fixed and independent of its composition, and equilibrium sizes depend only weakly on it.We implemented soft selection because we aimed to assess the impact of population structure in its most minimal form.A promising direction for future research would be to study such a model with hard selection.Indeed, in models where subpopulations contribute to the next bottleneck proportionally to their sizes, the fraction of cooperators can in-creasing overall despite decreasing within each subpopulation [19]: this is known as the Simpson paradox [58][59][60][61].Such a coupling between composition and population size can favor the evolution of cooperation [13,20,21].It will thus be important to understand how this effect is impacted by graph structure.A first interesting step in this direction was conducted in Ref. [62], which determined conditions under which cooperation can be favored at mutation-selection equilibrium, in a model with constant total population size but variable deme sizes.
In this context, we can summarize our results as follows.The ability of specific graphs and update rules to favor the evolution of cooperation in models with one individual per node of the graph [17,18] does not survive coarse-graining in a model with soft selection.This result does not arise from considering larger deme sizes, but from separating scales.Indeed, in our model, migration is independent from birth and death, and cooperation occurs locally within demes, while spatial structure matters between demes.
In addition to considering hard selection, possible extensions include addressing other ways through which migration may couple with cooperation [63], studying other games such as the snowdrift game [25], and explicitly modeling the dynamics of a diffusible resource [12].
Appendix A: Formal comparison with models with one individual per node 1. Cycle [17] In [17], evolutionary games, and especially the prisoner's dilemma, were studied in the context of the cycle graph with one individual per node and with different update rules (birth-death, death-birth, imitation).Analytical expressions were obtained for the fixation probability of a cooperator in a population of defectors and vice-versa, in the weakselection regime w ≪ 1.The key result of [17] is that selection can favour cooperators over defectors for death-birth and imitation updates, but not birth-death updates.To understand the similarities and differences with our model, let us compare the derivations of the fixation probabilities in the cycle within each of the two models.
In our model with a well-mixed deme on each node of the cycle and migration events independent from birth and death events, and in the rare migration regime, let i denote the number of fully mutant demes.Upon each migration event, the number of mutant demes can either increase by one with probability λ i , due to a migration of one mutant individual to a wild-type deme and its subsequent fixation in it, or decrease by one with probability µ i upon the opposite succession of events, or stay constant in all other cases.A recurrence relationship can then be written on the fixation probability Φ cycle i in a cycle of D demes starting from i consecutive fully mutant demes.For this, one considers all the possibilities at the first migration event that occurs, which allows to relate Φ cycle i to Φ cycle i−1 and Φ cycle i+1 .Solving this recurrence relationship yields where γ j = µ j /λ j = ρ W /ρ M , which is independent of j.Therefore, Eq. (A1) gives A more detailed derivation (which is presented without frequency-dependent fitness, but extends to the present case) is given in [31], and the general method is presented in [28].Note that the expression in Eq. (A2) also holds in the case of the clique, where it can be obtained with an extremely similar derivation [31].
In the model with one individual per node of a cycle with N nodes considered in [17], the same approach was employed to express the fixation probability ρ cycle M of one mutant, yielding which is formally equivalent to Φ cycle 1 in our model (see Eq. (A1)).However, in this model, γ j = µ j /λ j , where λ j (resp.µ j ) represents the probability that the number of mutant nodes increases (resp.decreases) by one upon an update, takes different expressions depending on the update rule [17].Indeed, with the birth-death update rule, for 2 ≤ j ≤ N −2, where f e W (resp. f e M ) is the fitness of a wild-type (resp.mutant) at the edge of the cluster of mutant nodes, i.e. which has a mutant neighbor and a wild-type neighbor.Here we expressed each fitness as 1 − w + wπ where the total payoff π received by an individual on the cycle is computed by summing the payoffs they receive from interactions with their two neighbors, with the payoff matrix in Eq. (1).We also used the fact that during the fixation process, starting from one mutant, there is a cluster of consecutive mutant nodes that cannot break.With the death-birth update rule, the fitnesses involved in γ j depend on whether the mutant that was at an edge of the cluster of mutant nodes (individual 1) or its wild-type neighbor (individual 2) is selected to die.If individual 1 is selected to die, which can lead to a decrease of j by one, the fitness of the wild-type (resp.mutant) neighbor of individual 1 is denoted by f e W (resp. f e M ).If individual 2 is selected to die, which can lead to an increase of j by one, the fitness of the wild-type (resp.mutant) neighbor of individual 1 is denoted by f ′ e W (resp. f ′ e M ).We have for 3 ≤ j ≤ N − 3. The expressions in Eqs.(A4) and (A5) differ for all b > c > 0, which leads to different probabilities of fixation under these two different update rules [17].
This comparison of the derivation of fixation probabilities within our model and within that of [17] shows that the key difference is that in our model, migration events do not depend on birth or death and are thus decoupled from interactions between individuals, while in models with one individual per node, any move on the graph is coupled to birth and death via the update rule.
2. General graphs with pair approximation [18] While [17] focused on the cycle, [18] derived a general rule for cooperation on graphs in the same model with one individual per node.The key result of [18] is that under the death-birth (dB) update rule, natural selection favours cooperators (i.e. mutant individuals here) if b/c > k, where k here is the average degree of the graph (i.e., the average number of neighbours per node).No such effect exists under the Bd update rule [18].To understand the origin of the difference between these results and ours, we here follow the main analytical derivation of [18] and adapt it to our model.Again, the key difference is that our model does not involve any update rule, as birth, death and migration are all independent.We will point out where this matters.
Following [18], let p M (resp.p W ) denote the frequency of fully mutant M (resp.wild-type W ) demes in a given graph comprising D nodes.In our model, each of these nodes contains a well-mixed deme, while it comprises a single individual in [18].Let us also introduce the frequencies p M M , p W W , p M W and p W M of neighbouring M M , W W , M W and W M pairs of demes, and the conditional probability q X|Y to find an X deme at a certain position given that one randomly chosen nearest-neighbour of this X deme (i.e. a deme connected to it by an edge) is a Y deme.We focus on the rare migration regime.As in [18], we will make a pair approximation where only the pair correlation is accounted for, meaning that higher-order correlations are fully determined by these pair correlations.The pair approximation requires Bethe lattices (or Cayley trees), which are regular graphs (i.e., all nodes have the same degree) with no loops.We thus focus on regular graphs, and we further assume k > 2 (the case k = 2 corresponds to the cycle, treated in the previous section).The following identities hold: (A6) Thus, the system is fully described by p M and q M |M .We will focus on these two quantities.
To derive the fixation probability of a mutant starting from a fully mutant deme (and D − 1 fully wild-type demes), we first compute the probabilities that p M changes upon one migration event [18].The probability that the total number of M demes increases by one (meaning that one M individual has migrated from a fully M deme to a fully W deme and has fixed there) is given by: where k is the degree of our regular graph, while k M (resp.k − k M ) is the number of M (resp.W ) neighbours of a deme.Indeed, for this to happen, a W deme needs to be picked as the deme of arrival of the migration (probability p W ), one of its M neighbors needs to be picked as the deme of origin of the migration (probability k M /k), and the mutant migrant has to fix in the wild-type deme (probability ρ M ).In addition, we need to take into account all different possibilities for the numbers of neighbors of each type that the W deme where the migration arrives possesses, leading to a sum with binomial weights.Recognising a binomial sum, we obtain which has a simple interpretation: we need to pick a W M neighbouring pair for our migration event, and then the mutant needs to fix in the wild-type neighboring deme.Importantly, this last simplification cannot be made in the proof of [18], because the replacement probabilities upon an update (which are replaced by fixation probabilities ρ M in our model) depend upon the state of the neighborhood (i.e., on the value of k M ) due to the cooperative interaction.
In the model of [18], this dependence involves the update rule, and vanishes for the Bd update rule while it remains for the dB one.In fact, here, the reasoning we just made starts by picking the arrival deme of the migration (in a dB spirit), but an equivalent one can be made by first picking the deme where the migration originates (in a Bd spirit): In the deterministic limit, we have: where we introduced By computing the probabilities that p M and p M M change in one migration step, we derived two differential equations on our two variables of interest, p M and q M |M , in the deterministic limit.In the weak selection regime, we suppose that the conditional probability q M |M which describes the local interactions reaches an equilibrium much more quickly than the frequency p M (which describes the overall state of the structure) [18].Therefore, the system quickly converges onto the space defined by F 2 p M , q M |M = 0.This gives Under this assumption, the system is fully described by p M .In order to compute probabilities of fixation, we need to go beyond the deterministic limit, and to work in the diffusion limit [18].To write a diffusion equation satisfied by p M , we need the second moment ∆p 2 M of ∆p M , which reads where we have used Eq.(A20).We can now write down a Kolgomorov backward equation on the fixation probability ϕ M (y) of M starting with an initial frequency p M (t = 0) = y: where m(y) and v(y) are given by: Recall that k > 2 and that we focus on the weak selection regime w ≪ 1.This Kolgomorov backward equation can be solved with the two conditions ϕ M (0) = 0 and ϕ M (1) = 1 (which describe the two absorbing states): Starting from one single fully mutant deme, the probability that mutants take over the population is thus: Similarly, starting from one single fully wild-type deme, the probability that wild-types take over the population is: Note that the same equation holds with mutants with constant fitness [32], where wild-types have fitness f W = 1, and mutants have fitness f M = 1 + s, but with s instead of −w ′ c.Thus, here, with the approximations we made, −w ′ c plays the part of a constant fitness advantage.The evolution of the mutant fraction x = k ′ /(k ′ + l ′ ) can be directly deduced from Eq. (B3).Under Kimura's diffusion approximation, the fixation probability ρ ′ (x 0 ) starting from a fraction x 0 of mutants in a population of size N ≫ 1, with w ′ ≪ 1, then reads: This yields the probability ρ M of fixation of a single mutant The probability ρ ′ W of fixation of a single wild-type individual in a population where other individuals are mutants can be obtained similarly and reads: We can use these expressions to calculate the fixation probability in a structure in the rare migration regime, via Eq.( 6).

Mapping between carrying-capacity and dilution model in the diffusion approximation
For simplicity, let us first consider the case where individuals have constant fitness.In our work, evolution within a deme is modeled via the Moran model in the carrying-capacity model (neglecting deme size fluctuations), and via a dilution model close to the Wright-Fisher model, where the e s ′ t − 1 ≈ s ′ t plays the role of the fitness advantage s WF in the Wright-Fisher model [32].In a well-mixed population of constant size N , the mutant fixation probability starting from a fraction x 0 of mutants is well-known in the diffusion limit, both for the Moran and the Wright-Fisher models [64][65][66].Let f W = 1 denote the wild-type fitness, s the mutant's fitness advantage in the Moran model and s ′ in the Wright-Fisher model.If N ≫ 1 and s, s ′ ≪ 1, the diffusion approximation predicts a fixation probability 1 − e N σ , (B7) where σ = s in the Moran model, while σ = 2s WF in the Wright-Fisher model, and σ = 2(e s ′ t − 1) ≈ 2s ′ t in our dilution model.Thus, there is a mapping between s in the carrying-capacity model and 2s ′ t in the dilution model.Note that the factor of 2 differs between the Moran model and the Wright-Fisher model: it arises from their different variances in offspring number [67].
Let us now include cooperation.In our serial dilution model, the probability of fixation of a single cooperator mutant is given by Eq. (B5) in the diffusion limit.Let us write the first-order expansion of the previous equation in w ′ , assuming that w ′ ct ≪ 1 and w ′ ctN ≪ 1: In the regime of parameters we consider here, where N ≫ 1, while b and c are of order 1, this is the same as the probability of fixation ρ M obtained in Eq. ( 9) within the carrying-capacity model, but replacing w by 2w ′ t.Thus, 2w ′ t is the effective selection intensity we use in the dilution model for our quantitative comparisons with the carrying-capacity model.

Appendix C: Variant of the carrying capacity model
In the logistic model used throughout, the division rate of individuals of type A = M, W in a deme is given by the logistic function f A (1 − N/K), where f A is the fitness given in Eq. ( 3) and N is the number of individuals in the deme.Meanwhile, the death rate g is independent from type and deme size (see Model and methods).An alternative way of implementing logistic regulation of population size is to consider division rates that do not depend on deme size but a death rate that does.Let us consider such a variant of the carrying capacity model, by taking a division rate f A for individuals of type A = M, W and a death rate f N/K, where f is the mean fitness in the deme, for all individuals.Since f = [kf M (k, N ) + (N − k)f W (k, N )] /N , the death rate reads [kf M (k, N ) + (N − k)f W (k, N )] /K.In this model, the steady state of the deme is equal to its carrying capacity K, which thus plays the role of K(1 − g/f ) in our usual model.
In the rare migration regime, within-deme events enter the fixation probabilities in graph-structured populations only through the fixation probabilities in demes.Therefore, to compare the two variants of the carrying capacity models in this regime, it is sufficient to compare the fixation probabilities in demes.In Fig. 6, we present simulation results corresponding to these two variants, which are in good agreement.As expected, they are also both close to the analytical prediction obtained within the Moran model.The slight discrepancy arises from deme size fluctuations, which do not exist in the Moran model.Note however that the two models yield different time scales.For instance, in the simple case where f W = f M = 1, at steady state a division event occurs every 1/g time unit on average in the first model, and every time unit in the second model.

1/N
• W i mij −−→ W j : migration of a wild-type from deme i to deme j.
Here, we have denoted by N M,i and N W,i the numbers of mutant and wild-type individuals in deme i.The total rate of possible events is then given by: m ij (N W,i + N M,i ) .(D1) Reactions occur at time intervals picked in an exponential distribution with mean 1/k tot , and the specific reaction that occurs is selected randomly proportionally to the ratio of its rate and k tot .Reactions are performed until one of the two types takes over the entire population.
In the dilution model, we performed our simulations as described in [32], except that multinomial sampling was employed throughout, as described in Appendix B, and that the deterministic growth was performed using a numerical resolution of Eq. (B1) to include cooperation.

FIG. 1 .
FIG. 1. Highly symmetric population structures: (a) clique, (b) cycle, (c) star and (d) 3 × 3 regular lattice.(a), (b), and (c) have 5 demes, while (d) has 9 demes, represented by large white circles.Each deme holds a well-mixed population of microorganisms represented by red circles.Individuals can migrate from one deme to another along the depicted arrows, following the written migration rates.

FIG. 3 .
FIG. 3. (a) Mutant fixation probability ρ cycle M in the cycle versus intensity of selection w for different migration asymmetries α = mA/mC , in the rare migration regime.Solid lines represent analytical predictions; red "+" markers are semianalytical predictions; other markers are simulations.The well-mixed population case (black dashed line), and the neutral mutant case (gray line) are shown for reference.Analytical and semi-analytical predictions are obtained using Eqs.(6,11), with ρM and ρW calculated either using Eqs.(7,8) that hold for fixed deme size, or using simulations on wellmixed demes (10 7 replicates) to avoid errors from assuming fixed deme size.Pure simulation results are obtained over 10 6 replicates with (mA, mC ) ×10 5 = (4,1), (1,1) and (1,4).(b) Mutant fixation probability ρ cycle M in the cycle versus effective intensity of selection 2w ′ t in the dilution model, for different migration asymmetries and intensities.The solid red line represents the analytical prediction in the rare migration regime in the dilution model, from Eqs. (6,11), using Eqs.(B5, B6) for ρM and ρW .The well-mixed population case (black dashed line) in the dilution model and simulations for the carrying-capacity (C.C.) model in the rare migration regime (versus w) are shown for comparison.Dilution simulation results are obtained over 10 6 realizations for (mA, mC ) ×10 and ×10 3 = (4,1), (1,1) and (1,4).In this figure, we consider D = 5 demes with carrying capacity K = 20, death rate g = 0.1 yielding steady-state size N = 18, and benefit and cost b = 2 and c = 1.The wellmixed population holds DN = 90 individuals.

FIG. 5 .
FIG. 5. (a) Mutant fixation probability ρ lattice M in the square lattice with periodic boundary conditions versus intensity of selection w for different values of migration rates, in the rare migration regime.(b) Mutant fixation probability ρ lattice M in the square lattice with periodic boundary conditions for the dilution model versus effective intensity of selection 2w ′ t, for different migration asymmetries and intensities.Here, we consider D = 9 demes with carrying capacity K = 20 and steady-state size N = 18.Other parameter values and conventions are the same as in Fig. 3.All simulation results are obtained over 10 6 replicates.