Stochasticity , Selection , and the Evolution of Cooperation in a Two-Level Moran Model of the Snowdrift Game

1Wisconsin Institute for Discovery, University of Wisconsin-Madison, 330 North Orchard Street, Madison, WI 53715, USA 2Laboratory of Genetics, University of Wisconsin-Madison, 425-G Henry Mall, Madison, WI 53706, USA 3School of Philosophy, National Research University Higher School of Economics, Staraya Basmannaya st., No. 21/4, Moscow 105066, Russia 4Department of Mathematics, University of Wisconsin-Madison, 408 Lincoln Drive, Madison, WI 53706, USA 5Department of Philosophy, University of Wisconsin-Madison, 600 North Park Street, Madison, WI 53706, USA 6Department of Philosophy and Religion, Northeastern University, 360 Huntington Avenue, 371 Holmes Hall, Boston, MA 02115, USA


Introduction
Evolutionary game theory (Hofbauer and Sigmund 1998 [1]; Smith 1982 [2]; Smith and Price 1973 [3]) allows one to analyze the evolutionary dynamics of a population of individuals in which individual fitness is frequency dependent.A model in evolutionary game theory is based on a payoff matrix, which describes the payoff an individual will receive given its own behavior and the behavior of its partner(s).Payoff Matrix 1 represents a symmetric, two-player game.The entries in the matrix-, , , and -represent the payoff the row player will receive when it plays strategy A or B with its partner, the column player, who can also play either strategy A or B. Payoff Matrix 1.A two-player, symmetric game is as follows: . (1) The Prisoner's Dilemma, in which  >  >  > , is perhaps the most widely studied form of social interaction, having been used to model systems ranging from microorganisms (Conlin et al. 2014; [4], Frey 2010 [5]; West et al. 2006 [6]) to human societies (Axelrod and Hamilton [7]; Bowles and Gintis [8]).However, The Prisoner's Dilemma does not 2 Complexity accurately represent social dilemmas in which there is no strictly dominant strategy.
The Snowdrift Game (Sugden 1986 [9]), also known as the Hawk-Dove Game (Smith 1982 [2]), is a social dilemma in which an individual can participate ("cooperate") or not ("defect") in producing a public good.In this game, two individuals, each in her own car, want to get home, but their cars are blocked by a large pile of snow.For either to get home, at least one person must get out of her car to shovel the snow.But there is a cost to doing so, creating a strategic dilemma in which  >  >  > .The game is relevant to a number of phenomena in biology, such as collective defense and resource extraction among microorganisms (Gore et al. 2009 [10]; Conlin et al. 2014 [4]), behavioral contests (Smith 1982 [2]), distributive justice in humans (Sugden 1986 [9]), behavioral diversification (Doebeli et al. 2004 [11]), and branching events (Wakano and Lehmann 2014 [12]).
In what we might think of as a "standard" model in evolutionary game theory, there is an infinitely large, well-mixed population within which individuals either cooperate or defect in one-time, pairwise games with each other.One generally assumes that individuals reproduce asexually and that the number of offspring each individual has is proportional to its fitness.This makes it easy to track how the frequencies of cooperators and defectors change in the population over time.The rate of change of a strategy is given by the replicator dynamics (Hofbauer and Sigmund 1998 [1]; Hofbauer et al. 1979 [13]; Taylor and Jonker 1978 [14]), an ordinary differential equation.
If a population of individuals is playing this game, then traditional evolutionary models, in which the dynamics are continuous and deterministic, predict a stable, interior equilibrium frequency of cooperators (Doebeli and Hauert 2005 [15]; Doebeli et al. 2004 [11]; Hauert and Doebeli 2004 [16]) (see (7)).Thus, the standard, deterministic model of the Snowdrift Game describes a scenario in which, even in the absence of facultative trait expression or heterozygote superiority, a stable polymorphism of behaviors can emerge.
However, all biological communities are finite, and many are small and organized into groups that together form a metapopulation (Gilpin 2012 [17]; Hanski 1999 [18]).This was certainly the case for ancestral hominins, which has plausibly influenced the evolution of human cooperation (Bowles and Gintis 2011 [8]), and it also appears to characterize many other taxa, including bacteria (Lieberman et al. 2016 [19]).It is therefore important to understand how a model with finite population size and metapopulation structure can change the predictions generated from models that assume a single population whose size tends to infinity.
Here we explore a discretization of the standard model of the Snowdrift Game.We consider a metapopulation composed of a finite set of discrete, nonintermixing groups, which are themselves composed of a finite set of discrete individuals who either cooperate or defect in the Snowdrift Game.The evolutionary dynamics both between and within groups are governed by a discrete time Moran process (Moran 1958 [20], 1962 [21]), a special case of a discrete time Markov chain.
Our analysis has two main results.First, we show that the combination of within-group stochasticity and group selection can promote the evolution of cooperation in the metapopulation and can even result in cooperation's fixation.This would be an impossible result were the evolutionary dynamics deterministic.
Second, we describe a phase transition for the fixation and extinction probabilities of cooperation in any finite group of a constant size whose members play the Snowdrift Game.Letting r stand for the cost-to-benefit ratio in a given Snowdrift Game, this threshold quantity, which we call  * , is approximately equal to 0.796.If  <  * , the probability cooperators will fix tends to 1 as group sizes go to infinity.If  >  * , the probability cooperators will go extinct tends to 1 as group size goes to infinity.This is true so long as the starting frequency of cooperators is strictly between 0 and 1. (As we detail, more complicated dynamics occur when  =  * .)The existence of this threshold quantity allows us to state a sufficient condition for cooperators to fix in a metapopulation.Moreover, while this threshold result comes about from taking the limit of group size, we in fact show that even when group size is fairly small (e.g., 100), the value of  effectively determines whether cooperators will fix or go extinct within a group.This threshold has no analogue in a deterministic model of the Snowdrift Game, and it makes the important differences between discrete and continuous models of evolution salient.

The Model
Here, we describe the within-and between-group Moran processes in our model.
Payoff Matrix 2 provides the payoff matrix for the twoplayer Snowdrift Game we will assume throughout this paper.Let b refer to the benefit of getting home and let c refer to the cost of shoveling snow, where  > .Each driver can either get out of her car to shovel ( for "cooperate") or stay in her car ( for "defect").If both drivers cooperate, then each receives a net payoff of  − /2, since both receive the benefit of going home, while the cost of shoveling is divided in half.If one driver cooperates and the other defects, then both get to go home, but the cooperator must pay the full cost of shoveling snow, receiving a net payoff of −, while the defector pays no cost, receiving a net payoff of .If each driver stays in her car, neither pays the cost of shoveling, but neither gets to go home, so neither receives any reward.Following Zheng et al. (2007 [30]), we let  stand for the ratio of the cost of cooperating in the Snowdrift Game to the benefit of doing so ( = /).We can thereby speak of the "-value" of an instantiation of the game.
Payoff Matrix 2. The Snowdrift Game is as follows ( = /): Note that "cooperation" in this context refers to the shoveling snow behavior-that is, strategy C in Payoff Matrix 2. This strategy is not a form of altruism since cooperation can be in the interest of the actor, depending on the behavior of the other player.Strategy C in Payoff Matrix 2 coincides with a technical definition of cooperation (West et al. 2007 [31]): cooperative behaviors (i) increase the payoff to others and (ii) carry a benefit (or cost) to the actor contingent on the behavior of others.
2.1.Within-Group Moran Process.Suppose there is a metapopulation composed of a finite number  of discrete groups, indexed by , so that  ∈ {1, 2, . . ., }.The size of a given group  (  ) is finite and constant, and each group in the metapopulation is of the same size.We assume the individuals in each group are hard-wired to either cooperate or defect in the Snowdrift Game.The evolutionary dynamics within each group are governed by a discrete time Moran process (Moran 1958 [20], 1962 [21]).At each time step, an individual in a group j is chosen to reproduce and have one offspring.The probability that a given individual is chosen to reproduce is proportional to its fitness relative to the average fitness of the individuals in its group.The behavior of an offspring is always identical to the behavior of its parent-that is, cooperators always beget cooperators, while defectors always beget defectors.At the moment it is born, an offspring replaces uniformly at random some individual in j, perhaps its parent, but not itself.The probability an individual will be replaced is unaffected by that individual's phenotype or fitness and is always 1/  .For our purposes here, we assume there is no migration between groups and no mutations during reproduction.
Since we are interested in modeling interactions that generate public goods that can be used by all, we will assume individuals play the Snowdrift Game with all of the members of their respective groups, including themselves, simultaneously.This is equivalent to using the expected fitness of random pairwise interactions among members of the group allowing for self-interaction.
Formally, we can represent the fitness of cooperators and defectors in a group j as follows.Let   and   index the parameters b and c for a given group j, let   stand for the frequency of cooperators in j, and let 1 −   stand for the frequency of defectors in j.The fitness of a cooperator and a defector in group j for a given value of   are given, respectively, by The average individual fitness of the members of a group j is given by The composition of a group can change in one of two ways: a cooperator can replace a defector, or a defector can replace a cooperator.Letting   stand for the number of cooperators in a group j, we can represent the first transition as (  →   + 1) and the second as (  →   − 1).As described elsewhere (Fudenberg et al. 2006 [32]; Nowak 2006 [33]; Taylor et al. 2004 [34]), we can use the fitness given in (3)-(4) to calculate the probability, , of each of these two transitions: Since no other changes within a group are possible, the probability that the state of the group will not change is given by, Were each group well-mixed and infinitely large, each group j would converge to a stable, internal equilibrium frequency of cooperators ( *  ), which is given by the following (Hauert and Doebeli 2004 [16]): However, because group size is finite in our model and there are no mutations, the only truly stable states of a group are the two absorbing states, in which cooperators fix (  = 1) or go extinct (  = 0).Nevertheless, the frequency of cooperators in a group in our model will often temporarily oscillate around what its internal equilibrium value would be were group size to be infinite; for a finite population, this is sometimes called its "quasi-equilibrium" (Shpak et al. 2013 [35]).

Between-Group
Moran Process.Our model also involves a discrete time Moran process that occurs between groups.Within the metapopulation, a "parent" group is chosen to replicate, thereby producing a "daughter" group, which replaces uniformly at random some group in the metapopulation, perhaps its parent, but not itself.A daughter group of some group j that is created at time t has the same frequency of cooperators as j at t and has the same payoff matrix as j.(See Figure 1.) The probability that a given group is chosen to reproduce is proportional to the average individual fitness of its members relative to the average individual fitness of the members of all the groups in the metapopulation.Taking the first and second derivatives of (4), we see that the average individual fitness of a group increases monotonically but nonlinearly as its frequency of cooperators increases: In our model, a group's fitness is equal to the probability it will give rise to a daughter group at the next time step.
Figure 1: Two-Level Moran Process.This is a condensed representation of the two-level, discrete time Moran process in our model.Each isolated cluster of red and blue ovals is a "group, " while the ovals themselves are "individuals," which come in two varieties, cooperators (blue) and defectors (red).The set of six groups constitutes a metapopulation.Within a group, an individual is chosen to reproduce and replace some individual in the same group: the arrow represents an individual having an identical offspring, which replaces another individual in the group.(The replaced individual has a dashed border around it.)Across groups, a group is chosen to reproduce and replace some group in the metapopulation.The group in the bottom-left corner reproduces and replaces the group in the top-right corner.(The replaced group has a dashed border around it.) Hence, a group's fitness increases along with its frequency of cooperators.(One can of course consider other configurations of the relationship between individual fitness and group fitness-e.g., where a group's fitness decreases as its frequency of cooperators increases-but we do not pursue such extensions here.) In our model, we assume individual-level birth-death events occur more frequently than group-level reproductionextinction events.Again, the abstract structure of the model itself does not require this.Presumably, the relative rates of the individual-and group-level events in a model should be governed by the target system one is modeling.
There is a straightforward way to calculate the probability that a given group j will replicate or die at the next time step.Letting   (+1) and   (−1) refer, respectively, to these two probabilities, we have where   is the average fitness of the individuals in group , given by equation ( 4), refers to the average fitness of all groups other than group , and   is the average of the average individual fitness of each group in the metapopulation, which is given by The probability that there will be no change in the metapopulation is given by

Simulations of Moran Process.
We ran two classes of simulations.In the first class, we were concerned exclusively with the Moran process within a group.In the second class, we were concerned with our "complete" model, in which Moran processes occur both within and between groups.Both classes of simulations were implemented in R version 3.2.3(R Core Team 2016 [36]).High Performance Computing (HPC) was necessary to perform some of the metapopulation simulations.This is because our code "counts" each individual in the metapopulation at each individual-level birth-death event and at each group-level replication-extinction event.
This creates a computational bottleneck that can be nontrivial when group size is large (  ⩾ 1000).The source code is available by request.

Simulations of Within-Group Moran
Process.In our first class of simulations, we sought to explore how the withingroup evolutionary dynamics were affected by the size of a group and the group's -value.We therefore simulated the within-group evolutionary dynamics of nine groups, where each group was associated with a different r-value.The range of r-values that we simulated was  ∈ {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9}, where b was set to 10 and c ranged from 1 to 9 in increments of 1.We simulated a range of group sizes,   ∈ {10, 100, 1000, 10,000}, and we allowed each group to evolve for 100 generations, where 1 generation for a group of size   is equal to  birth-death events.The observed frequency of cooperators within each group after a given number of birthdeath events was plotted against the equilibrium frequency of cooperators predicted by the deterministic within-group model, as given by (7).

Simulations of Two-Level Moran Process.
In our second class of simulations, we simulated the complete version of our model, in which there are discrete time Moran processes that occur both within groups and between groups.
In our first set of simulations within this class, we considered a scenario in which the groups in the metapopulation do not vary in their r-values.We seeded the metapopulation with 25 groups, where the frequency of cooperators in each was initially set to 0.10.We ran three versions of these simulations, which differed in the r-value we assigned to each group ( = .25; = .67; = .90).For each -value, we ran simulations where   ∈ {10, 100, 1000}.We ran 500 independent simulations for each combination of -value and group size.Within a given simulation, there were approximately 100 within-group generations and 20 between-group generations, where one generation for a metapopulation of size  is equal to  group replication events.
In our second set of simulations within this class, we considered a scenario in which the groups in the metapopulation initially varied in their r-values.To explore the sensitivity of the metapopulation dynamics to the r-values of each group, we ran two versions of this type of simulation.In the first version, we seeded the metapopulation with 25 groups and associated each with a unique r-value so that the average r-value was relatively low: the benefit b of cooperating was set to 10, while c ranged from 0.25 to 6.25, in increments of 0.25, resulting in an initial average r-value of 0.33.In the second version, we performed the same analysis, except that we associated each group with a unique r-value so that the average r-value was relatively high: b was set to 10, while c ranged from 7.79 to 9.94, in increments of approximately 0.09 on average, resulting in an initial average -value in the metapopulation of 0.89.These simulations were "truncated" in that they could stop before the metapopulation became monomorphic.We conducted truncated simulations because we were interested in documenting the evolutionary dynamics in the short and medium term, not only as one takes the limit of time.As our results show, the groups in the metapopulation, and the metapopulation itself, often spend a good deal of time at a polymorphic state.This is because the fixation times can be long, which underscores the need to not only consider the end state of the system.Eventually, this polymorphism will disappear, and so the simulations are an inaccurate representation of the dynamics of the system for an arbitrarily long sequence of events.
Comparison with Null Models.For each set of simulations within this class, we plotted the observed frequency of cooperators in the simulations against what the long-term, expected frequency of cooperators in the metapopulation (E[  ]) would be given two alternative, null models.In the first null model, there is a single, well-mixed population of individuals playing the Snowdrift Game, within which the dynamics are deterministic.In the second, there is a metapopulation of groups whose members play the Snowdrift Game, where the within-group dynamics are deterministic, there is group reproduction, and there is no group selection (i.e., groups do not vary in their probability of replicating).Conveniently, the value for (E[  ]) for each of these null models is given by the same equation: where   represents the number of groups with a given value in the metapopulation,  *  stands for the equilibrium frequency of cooperators in a group of type-r, and the notation  ∈  indicates that we perform this operation for each r-value that is in the set of r-values that are represented in the metapopulation.In the case where there is a single, infinitely large population, we can think of this population as a "metapopulation" with a single, infinitely large group. ) ≥0 that solves the following SDE:

Characterization of Within-Group Moran Process
where  is the size of the group (note that we have left out the  subscripts from  for ease of reading),   is the Wiener process (standard Brownian motion), and  and  are the "drift" and the "volatility" functions given, respectively, by where Equation ( 14) describes how a change in the frequency of cooperators is dependent both on selection and on the stochasticity of the system that emerges because the population is finite.Observe that when the group's size tends to infinity, the SDE reduces to an ODE: The precise connections between the Moran process and the corresponding SDE and the Moran process and the Complexity corresponding ODE are given by the following lemma.A rigorous proof follows from standard stochastic analysis (Durrett 1996 [37]).

Lemma 1.
Let  be the total number of individuals and let   () be the number of cooperators at the th step in the Moran process.Suppose the initial fraction of cooperative players   (0)/ tends to  0 ∈ [0, 1] as  → ∞.Then, for all  > 0 and  < ∞, we have the following.
(1) Deterministic approximation: Moran process versus ODE: where () is the solution to the ODE (17) starting with (0) =  0 , "sup" is the supremum, and [] is the largest integer smaller than or equal to . (2) Diffusion approximation: Moran process versus SDE: there exists a diffusion process  ()  defined on the same probability space as that of the Moran process, which starts at  () 0 =  0 , solves (14), and is such that for any sequence   such that lim →∞ (log /  ) = 0.
In evolutionary biology, one is often interested in the probabilities of fixation and extinction.Let ℎ() be the probability that the diffusion process  () hits 1 before 0 (i.e., cooperators fix), given that the process starts at  () 0 = .In mathematical notation, we let  fl inf{ ≥ 0 :  ()  = 0 or 1} be the time it takes for the population to become monomorphic, and we let P  be the probability law of  When  = 0.3, which is below  * , the probability that cooperators will fix gets very large, regardless of the starting frequency of cooperators, as group size increases.
Note that ℎ() depends on  and the parameters in the game matrix.This probability can be found in standard textbooks on the subject (e.g., Otto and Day 2007 [39], Chapter 15).We record the exact formula for our case in the following lemma.Lemma 5. Let ℎ() be the probability that cooperation fixes, defined in (20).Then, where () = −2 ∫

Asymptotic Analysis Reveals a Threshold Cost-to-Benefit
Ratio.Through an application of Laplace's method to the explicit formula for the fixation probability obtained in Lemma 5, we discover a threshold quantity  * ≈ 0.796 such that when  <  * the probability cooperators will fix tends to 1 as group size goes to infinity (Figure 2).Likewise, when  >  * , the probability cooperators will go extinct tends to 1 as group size goes to infinity (Figure 3).More complicated dynamics occur when  is roughly equal to  * (Figure 4).Remarkably, this holds irrespective of the initial frequency of cooperators in the group.Moreover, the probability cooperation will fix or go extinct, given an arbitrary starting When  = 0.9, which is above  * , the probability that cooperators will fix is negligible, regardless of the starting frequency of cooperators, as group size increases.frequency of cooperators, can be extremely high even when group size is relatively small (e.g.,   = 100) (Figures 2  and 3).To the best of our knowledge, this threshold quantity is unreported in the extant literature.We now describe its derivation in the "Snowdrift Game Threshold Theorem."Note that this theorem applies only to the form of the Snowdrift Game presented in Payoff Matrix 2, not to the more general form of the game presented in Lemma 5.

Simulations of Within-Group Moran
Process.The results of our purely within-group simulations are presented in Figure 5.The jagged lines are the observed frequency of cooperators within a group with a particular r-value after a given number of birth-death events.The horizontal lines of the corresponding color are the equilibrium frequency of cooperators one would expect in that group were the group to be infinitely large, as given by (7).
In accord with analytic predictions, over a given number of generations, relatively smaller groups are more susceptible to larger fluctuations in group composition.Moreover, given a sufficient amount of time, and large enough groups, the within-group dynamics tend to follow a given pattern: in the short term, the frequency of cooperators in a group migrates toward the equilibrium frequency of cooperators given by (7), and in the medium term the frequency cooperators in a group tends to oscillate around this equilibrium frequency, creating a quasi-equilibrium.This pattern becomes more evident as one moves from the simulations with the smallest group size (Figure 5(a)) to the simulations with the largest group size (Figure 5(d)).

Characterization of Metapopulation Dynamics.
Theorem 6 allows us to state a sufficient condition for whether cooperators will fix (go extinct) in our two-level Moran model of the Snowdrift Game.As one takes the limit of group size, cooperation will fix (go extinct) in the metapopulation if each group in the metapopulation has an r-value below (above)  * .Note that, in the limit of group size, the betweengroup selection process is irrelevant to the long-term dynamics of cooperation in the metapopulation.While this sufficient condition is a limit result, we stress that a group's size need not be very large for its -value, relative to  * , to substantially influence the probability cooperators will fix or go extinct.
When groups' -values straddle  * , when group sizes are small, or when one is concerned with evolutionary dynamics over shorter intervals of time, simulations are informative, owing to the challenge of formal analysis of these cases.Even though the two-level Moran process is a finite state Markov chain that counts the number of cooperators in each group, and so the probability of fixation can be computed explicitly, it is not obvious how to write down a closed-form expression in terms of (i) the number of groups, (ii) the number of individuals in each group, and (iii) the initial frequencies of cooperators in the groups when the number of groups and/or the size of the groups are as large as we are concerned with here.
The results of our metapopulation simulations are presented in Figures 6 and 7.

Stochasticity and the Snowdrift Game.
There is a good deal of work that explores how stochasticity affects the evolutionary dynamics of a population of individuals playing a game in general (Fudenberg et al. 2006 [32]; Nowak et al. 2004 [41]; Sandholm 2011 [42]; Traulsen and Hauert 2009 [43]).In our model, we explore a particular form of "demographic stochasticity" (Rice 2008 [44]; Shaffer 1981 [45]; Shpak et al. [35]).Demographic stochasticity occurs when there is a difference between an individual's (or type's) expected fitness and its realized fitness.In our model, this occurs when the type ("cooperator" or "defector") that is chosen to reproduce within a group is not the type with the higher fitness, and when the group that is chosen to replicate is not the group with the highest group fitness.Other forms of demographic stochasticity can emerge if an individual plays a random sample of the members of a finite population (Shpak et al. 2013 [35]), or if there are "shocks" to the payoffs of a game (Fudenberg and Harris 1992 [46]).
There is, in fact, a considerable amount of work that explores how stochasticity affects the evolution of cooperation in a single, unstructured population of individuals playing the Snowdrift Game (Ficici and Pollack 2000 [47]; Fogel et al. 1997 [48]; Fogel et al. 1998 [49]; Fudenberg et al. 2006 [32]; Shpak et al. 2013 [35]; Taylor et al. 2004 [34]; Voelkl 2010 [50]).While the results of these analyses are consistent with our mathematical and simulation results, to the best of our knowledge we are the first to report a threshold quantity for the fixation and extinction probabilities of a finite group of individuals playing the version of the Snowdrift Game presented in Payoff Matrix 2.
It is worth commenting on the general approach we have used to derive the value of this threshold quantity.In many traditional analyses in evolutionary game theory, one takes the limit of population size ( → ∞) and then time ( → ∞) to arrive at the long-term behavior of an evolutionary system.To derive the value of  * , we have switched the order in which these limits are taken.We take ( → ∞) to derive the formula for the fixation probability of cooperation.Then we take ( → ∞) by applying Laplace's method (Murray 1984 [40]) to this formula.Determining the "correct" order in which limits are taken presumably depends on the nature of the system one has in mind.(For a discussion of this issue as it relates to mutations and population size, see Sandholm 2010 [51].) Within the context of the Snowdrift Game, if the population is large and one is concerned with a short-or mediumterm time horizon, then the evolutionary trajectory of the system may be adequately represented by a "deterministic approximation" of the microscopic dynamics (Sandholm 2011 [42]).Indeed, our simulation results show a marked correspondence to the deterministic dynamics after 100 generations when  ⩾ 1000 (Figures 5(c) and 5(d)).If, in contrast, one is concerned with smaller populations or the longer-term behavior of a larger population, the stochasticity of the system becomes increasingly important-see Figures 5(a) and 5(b) and Theorem 6, respectively.

Multilevel Selection and the Snowdrift Game. Our results
provide insight into when within-group stochasticity and group selection do (and do not) result in evolutionary dynamics that are different from those in our two null models.
For instance, our results show that when the members of each group in the metapopulation play precisely the same parameterization of the Snowdrift Game, and when group  The jagged lines represent the frequency of cooperators in a group.As r-values get larger, the equilibrium frequency of cooperators gets smaller; thus,  = 0.1 is black, while  = 0.9 is teal.Each horizontal line of a corresponding color represents the equilibrium frequency of cooperators for that r-value given by (7).
size is small ( = 10), the combination of within-group stochasticity and group selection can promote the evolution of cooperation in the metapopulation above what its value would be given either of the two null models we consider (Figures 6(a)-6(c)).Why does this occur?Within each group, the frequency of cooperators is oscillating.Because the probability that a group will replicate at the next group replication event is proportional to its frequency of cooperators at the time of the group replication event, those groups with a relatively higher frequency of cooperators during the group replication event are more likely to replicate.When this occurs, the frequency of cooperators in the metapopulation will "jump" upward, as in a jump diffusion process.Indeed, within-group stochasticity can result in cooperators fixing in a group, which results in that group having the highest possible group fitness, making it possible that this group will then overtake the entire metapopulation.It is for this reason that we say withingroup stochasticity in combination with group selection can promote the evolution of cooperation in the metapopulation, even its fixation, for a range of parameter values.(Note, as the number of groups in the metapopulation becomes increasingly large, the between-group selection dynamics become increasingly deterministic; thus, that there are a finite, rather small number of groups in our metapopulation simulations in fact dampens the efficacy of group selection; this modeling choice was intentional.We wanted to, if anything, load the die against the evolution of cooperation.)However, our results also show that as the size of each group increases ( ⩾ 100), the within-group dynamics, over  (a, d, g) show the results for  = .25;panels (b, e, h) show  = .67;panels (c, f, i) show  = .9.Each panel shows results for 100 independent simulations with populations of 10 (a, b, c), 100 (d, e, f), and 1000 (g, h, i).In each simulation, there were 25 groups with an initial cooperator frequency of .10,500 group replication events, and 100 within-group generations.The red, black, and light blue (cyan) dots are the first, second, and third quartiles of the observed frequencies of cooperators in the metapopulation, while the thin black line is the mean frequency.The blue lines are the maximum and minimum observed frequency across 100 simulations.The dashed, horizontal line is the long-term, expected frequency of cooperation in the metapopulation assuming a null model given by (13).Note that the observed frequency of cooperators is substantially above its value predicted by the null models for smaller populations (a)-(c) and for one larger population with a high quasi-equilibrium value (d).
the time intervals we have simulated here, become increasingly deterministic.Over these time intervals, within-group stochasticity is less capable of resulting in jumps in the frequency of cooperators in the metapopulation after a group replication event, or resulting in cooperators fixing within a group (Figures 6(d)-6(h)).To be sure, because each metapopulation we simulated is composed of finite groups, the deterministic approximation of the dynamics will prove to be inaccurate in the long term.Our  * results tell us that it is extremely likely cooperators would have overtaken the entire metapopulation in some of our simulations (Figures 6(a When we allow the groups in the metapopulation to initially vary in their r-values, group selection can be quite efficacious, even for large group sizes (Figures 7(c)-7(f)).Regardless of group size, those groups with quasi-equilibria frequencies that are relatively closer to  = 1 will be more likely to replicate.When group sizes are small (e.g.,  = 10), cooperators will fix in many of these groups (Figures 7(a) and 7(b)).But even as group size becomes larger ( ⩾ 100), and the within-group dynamics become increasingly deterministic, those groups with higher quasi-equilibria are more likely to replicate.For this reason, we again see the jump diffusion process of cooperators in the metapopulation that we saw in the version of our model in which each group initially had the same -value.Within-group stochasticity and group selection are pushing the observed frequency of cooperators above what its expected frequency would be, given our null models.
Why have we considered a case in which we allow groups to vary in the parameter values of the Snowdrift Game their members play?First, in a biological context, if a metapopulation is composed of groups whose members are engaged in collective action that is qualitatively characterized by the Snowdrift Game, it would in fact be somewhat surprising if there were no quantitative differences between the payoff matrices associated with each group.(To be sure, there would also likely be variation in payoffs within groups, an interesting possibility we did not explore here.)For instance, these differences could emerge for reasons intrinsic to the group if, as an example, the members of some groups can more efficiently expend energy, so that the cost of participating in collective action is mitigated.Differences could also emerge because of ecological facts about the location of a group within the metapopulation, such as its proximity to energy sources or to external threats.To model the latter case, which we have not done, one would need to allow the payoff matrix to vary as a function of a group's location, which is an interesting avenue for future research.
Second, by allowing group payoff matrices to vary in our model, we learn something about the evolution of games themselves.Note that a central change brought about by selection and replacement of groups is a change to the payof f matrix of the underlying game.That is, as groups reproduce and go extinct, the payoffs of the game evolve.The possibility of an evolving game has been explored in other contexts, especially that of the Prisoners Dilemma, where it tends to promote cooperation in some contexts (Akc ¸ay and Roughgarden 2011 [26]; Worden and Levin 2007 [52]) and inhibit cooperation in others (Stewart and Plotkin 2014 [53]).Additionally, Hashimoto and Kumagai (2003 [27]) and Smead (2014 [28]) both present frameworks for the evolution of payoffs across a broad range of games.With respect to evolving games, our model is novel in two ways.First, we specifically target the Snowdrift Game and find that evolution of the payoffs promotes cooperation.Second, our model provides a mechanism for the evolution of the payoffs based on multilevel selection.This framework could be extended to a number of other games or settings for the purposes of modeling evolving payoffs.
Our model and results connect in illuminating ways to other work on multilevel selection in general, as well as to applications of multilevel selection to the Snowdrift Game in particular.Others have used a multilevel Moran model, or close variants, to explore the evolution of social traits (Hauert and Imhof 2012 [54]; Luo 2014 [22]; Simon et al. 2013 [24]; Traulsen and Nowak 2006 [25]).Notably, Traulsen and Nowak consider a model in which within-group dynamics are governed by the payoffs of playing a game, while in the metapopulation a group fissions (with some probability) when a group reaches a maximum size, replacing some other group in the metapopulation.Traulsen and Nowak provide conditions for when cooperation will be "favored" in terms of a payoff matrix, the maximum size of a fissioning group, and the number of groups in the metapopulation.In their work, selection favors cooperation when a mutant cooperator has a higher probability of fixing than does a mutant defector.Traulsen and Nowak's results extend to cases in which the within-group dynamics are governed by a Snowdrift Game (SI Traulsen and Nowak 2006 [25]).(See Simon and Pilosov 2016 [55] for discussion of limitations of Traulsen and Nowak's model.) Ours is a different model than Traulsen and Nowak's in that we consider a "strict" two-level Moran process, wherein both group size and metapopulation size are held constant.While two-level Moran models have been discussed before, and close variants of such a model have been described within the context of the Snowdrift Game (SI Traulsen and Nowak 2006 [25]), ours appears to be the most in-depth application of a two-level Moran model to the Snowdrift Game in the literature.
In fact, the dynamics of a strict two-level Moran model can differ in nontrivial ways from more general two-level Markov chain models.For instance, both Traulsen and Nowak (2006 [25]) and Simon et al. (2013 [24]) consider a model in which, when a group fissions, both the size and composition of the two resultant daughter groups are random variables.Because of this, a daughter group may end up with all cooperators (or none), even though its parent group was polymorphic.All else being equal, compared to a model, like ours, in which an entire group replicates, the stochasticity that emerges from the random seeding of daughter groups would appear to increase the probability that cooperation will fix in some groups and, hence, in the metapopulation as a whole.While such a model of group fissioning is biologically justifiable, in many cases more so than our own, our results show that one does not need the size or composition of daughter groups to be random variables in order for cooperation to evolve.
Finally, there has been other work that explores how group structure can affect the evolution of cooperation in the Snowdrift Game, where it has been shown to sometimes promote the evolution of cooperation, though not always (Hauert and Doebeli 2005 [15]; Killingback and Doebeli 1996 [56]).Hauert and Doebeli's work is particularly relevant to our own because they contrast the equilibrium of frequency of cooperators that would be found in a well-mixed population, given particular parameter values for a Snowdrift Game, with the frequency of cooperators observed in a population with individuals arranged on a lattice, given those same parameter values.The authors found that, for small  * values, the equilibrium frequency of cooperators in a population with a lattice structure is larger than what is expected in a wellmixed population.However, for intermediate and large  * values, the equilibrium frequency of cooperators found in the population is smaller than what is expected in a well-mixed population.They conclude, "unexpectedly, spatial structure reduces the proportion of cooperators for a wide range of parameters," adding, "our results caution against the common belief that spatial structure is necessarily beneficial for cooperative behavior" (p.643).
In contrast, our analysis shows that quite a different kind of population structure-namely, in which the groups are discrete units that are spatially isolated within a metapopulation and replicate as a function of their internal composition-can, in certain cases, facilitate the evolution of cooperation, rather than hinder it.Whether a lattice structure or the form of population structure we consider above is a better representation of a given biological system of course depends on the details of that system.

Conclusion
This article presented a model of discrete individuals packaged into discrete groups, where the individual-and grouplevel dynamics were stochastic.We showed that within-group stochasticity and group selection promote the evolution of cooperation in the metapopulation when group size is small, and also when group size is large and groups are allowed to vary in the payoff matrices their members play.We further showed that the long-term fate of cooperation as one takes the limit of group size is determined by the cost-to-benefit ratio of cooperating in the Snowdrift Game and that this ratio effectively determines the long-term fate of cooperation even when group size is fairly small (e.g.,  = 100).
All biological populations are finite (though they can be quite large), and they are often nested within a larger metapopulation.When these factors are incorporated into a model that would otherwise assume deterministic dynamics and no group structure (or group structure but no group selection), the differences in the predicted levels of cooperation can be substantial.

4. 1 . 1 .
Diffusion Approximation and Fixation Probability.FollowingTraulsen et al. (2005 [29]), we can obtain the diffusion approximation of the Moran process for a group whose members are playing a two-player game of the form given in Payoff Matrix 1.It is a diffusion process  () = (()

Figure 2 :
Figure2: The plot of ℎ().When  = 0.3, which is below  * , the probability that cooperators will fix gets very large, regardless of the starting frequency of cooperators, as group size increases.

Figure 3 :
Figure3: The plot of ℎ().When  = 0.9, which is above  * , the probability that cooperators will fix is negligible, regardless of the starting frequency of cooperators, as group size increases.

Figure 4 :
Figure4: The plot of ℎ().When  is almost exactly equal to  * , the probability that cooperators will fix is approximately  * /4 ≈ 0.2, regardless of the starting frequency of cooperators, as group size increases.

Figure 5 :
Figure5: Simulation results of within-group Moran process.Each panel shows the evolution of cooperation in each of nine different groups over the course of 100 generations within each group.Within a panel, groups are all the same size but vary in their r-values.The range of r-values that were simulated was  ∈ {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9}, where b was set to 10 and c ranged from 1 to 9 in increments of 1.The jagged lines represent the frequency of cooperators in a group.As r-values get larger, the equilibrium frequency of cooperators gets smaller; thus,  = 0.1 is black, while  = 0.9 is teal.Each horizontal line of a corresponding color represents the equilibrium frequency of cooperators for that r-value given by(7).

Figure 6 :
Figure6: Two-level Moran model with no variation in r-values.Each panel represents simulation results for our metapopulation model.Panels (a, d, g) show the results for  = .25;panels (b, e, h) show  = .67;panels (c, f, i) show  = .9.Each panel shows results for 100 independent simulations with populations of 10 (a, b, c), 100 (d, e, f), and 1000 (g, h, i).In each simulation, there were 25 groups with an initial cooperator frequency of .10,500 group replication events, and 100 within-group generations.The red, black, and light blue (cyan) dots are the first, second, and third quartiles of the observed frequencies of cooperators in the metapopulation, while the thin black line is the mean frequency.The blue lines are the maximum and minimum observed frequency across 100 simulations.The dashed, horizontal line is the long-term, expected frequency of cooperation in the metapopulation assuming a null model given by(13).Note that the observed frequency of cooperators is substantially above its value predicted by the null models for smaller populations (a)-(c) and for one larger population with a high quasi-equilibrium value (d).
This roughly says that sup ∈[0,] |  ([])/ −  ()  | → 0 at a rate of order at least (log /).Lemma 1 remains true if we replace the discrete time Moran process by the continuous time version of the Moran process with rate .This follows directly from (19) that when the population size  is large, most sample paths of the process   ([])/ stay within distance  from the solution () of the ODE for any time  ∈ [0, ].That is, once we specify the error tolerance  and the duration , we can choose  large enough so that |  ([])/ − ()| <  uniformly in  ∈ [0, ] with a probability close to one.Remark 3. Convergence(19)says that  () is a better approximation than () and precisely quantifies the extent to which it is a better approximation.Specifically,(19)holds as long as   → 0 slower than log /.