#### Abstract

The Snowdrift Game, also known as the Hawk-Dove Game, is a social dilemma in which an individual can participate (cooperate) or not (defect) in producing a public good. It is relevant to a number of collective action problems in biology. In a population of individuals playing this game, traditional evolutionary models, in which the dynamics are continuous and deterministic, predict a stable, interior equilibrium frequency of cooperators. Here, we examine how finite population size and multilevel selection affect the evolution of cooperation in this game using a two-level Moran process, which involves discrete, stochastic dynamics. Our analysis has two main results. First, we find that multilevel selection in this model can yield significantly higher levels of cooperation than one finds in traditional models. Second, we identify a threshold effect for the payoff matrix in the Snowdrift Game, such that below (above) a determinate cost-to-benefit ratio, cooperation will almost surely fix (go extinct) in the population. This second result calls into question the explanatory reach of traditional continuous models and suggests a possible alternative explanation for high levels of cooperative behavior in nature.

#### 1. 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:

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 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 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.

As we discuss, our results provide insight into the evolution of cooperation, particularly from a multilevel perspective (Luo 2014 [22]; Okasha 2006 [23]; Simon et al. 2013 [24]; Traulsen and Nowak 2006 [25]), evolving games (Akçay and Roughgarden 2011 [26]; Hashimoto and Kumagai 2003 [27]; Smead 2014 [28]), and the relationship between discrete and continuous models of evolutionary dynamics (Traulsen et al. 2005 [29]).

#### 2. The Model

Here, we describe the within- and between-group Moran processes in our model.

Payoff Matrix 2 provides the payoff matrix for the two-player 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 , 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 . 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 . 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 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, byThe 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 and the second as . 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 or go extinct . 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]).

##### 2.2. 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. 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 reproduction-extinction 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 and refer, respectively, to these two probabilities, we havewhere 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

#### 3. Methods

##### 3.1. 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 (). The source code is available by request.

###### 3.1.1. Simulations of Within-Group Moran Process

In our first class of simulations, we sought to explore how the within-group 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 , where* b *was set to 10 and* c *ranged from 1 to 9 in increments of 1. We simulated a range of group sizes, , 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 birth-death events was plotted against the equilibrium frequency of cooperators predicted by the deterministic within-group model, as given by (7).

###### 3.1.2. 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 (; ; ). For each -value, we ran simulations where . 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 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 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.

#### 4. Results

##### 4.1. Characterization of Within-Group Moran Process

###### 4.1.1. Diffusion Approximation and Fixation Probability

Following Traulsen 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 that solves the following SDE: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, bywhere

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 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 tends to as . Then, for all and , we have the following.*(1)*Deterministic approximation: Moran process versus ODE:**where is the solution to the ODE (17) starting with , “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 , solves (14), and is such that*for any sequence such that .*

*Remark 2. *Convergence (18) says 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 . That is, once we specify the error tolerance and the duration , we can choose large enough so that uniformly in 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 slower than . This roughly says that at a rate of order at least .

*Remark 4. *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 Kurtz (1978 [38], theorems 2.2 and 3.3).

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 . In mathematical notation, we let be the time it takes for the population to become monomorphic, and we let be the probability law of starting at . Then,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 and*

###### 4.1.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 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 frequency of cooperators, can be extremely high even when group size is relatively small (e.g., ) (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.

Theorem 6 (snowdrift game threshold theorem). *There exists a threshold number such that, for all , we have the limiting probabilitywhere is the probability that cooperation goes to fixation, given that the starting point is . Furthermore, is the unique solution of the following equation:on the positive real axis.*

*Remark 7. *The endpoints are given by and for all .

*Proof. *Asymptotic analysis (see, e.g., Murray 1984 [40]) reveals that the large behavior of the probability is governed by the maximum of the function . Recall the explicit formula and the notation in Lemma 5. For the Snowdrift Game, we haveClearly, and the maximum of is attained at the endpoints and . To derive (23), note thatElementary calculus shows that when , when , and when , where is the solution of on the positive real axis, which was found to be . When (i.e., ), we have, by the Laplace method (see, e.g., Murray 1984 [40]),Result (23) then follows.

###### 4.1.3. 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).

**(a)**

**(b)**

**(c)**

**(d)**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 of 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)).

##### 4.2. 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 between-group 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.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

#### 5. Discussion

##### 5.1. 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 medium-term 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 (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.

##### 5.2. 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 size is small (), 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 within-group 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 (), the within-group dynamics, over 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), 6(b), 6(d), 6(e), 6(g), and 6(h)) and gone extinct in others (Figures 6(c), 6(f), and 6(i)) had we allowed the simulations to continue to evolve, owing to the distribution of* r*-values of the groups in those simulations.

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 will be more likely to replicate. When group sizes are small (e.g., ), cooperators will fix in many of these groups (Figures 7(a) and 7(b)). But even as group size becomes larger (), 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 payoff 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 (Akç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 well-mixed 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.

#### 6. Conclusion

This article presented a model of discrete individuals packaged into discrete groups, where the individual- and group-level 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., ).

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.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.

#### Authors’ Contributions

This study was led by Laurence Loewe, Brian McLoone, and Wai-Tong Louis Fan, all three of whom conceptualized the project and developed the model. The simulations were implemented and analyzed in* R* by Brian McLoone. Wai-Tong Louis Fan performed the main mathematical analyses, including the derivation of the value of . Adam Pham and Rory Smead contributed to the mathematical development of the model and to the writing and organization of the manuscript. Rory Smead was responsible for the model’s link to the literature on evolving games. The manuscript was written by Brian McLoone with input from all, particularly Rory Smead, Wai-Tong Louis Fan, and Laurence Loewe.

#### Acknowledgments

The authors would like to thank Caitlin Pepperell and Tracy Smith for discussion and for detailed comments on earlier versions of this paper and Neil Van Lysel and Christina Koch at the Center for High Throughput Computing at UW-Madison for help in conducting their simulations. Shishi Luo, Burton Simon, and William Sandholm provided valuable comments on the manuscript. John Yin and David Schwartz raised important questions on the project in conversation. Brian McLoone acknowledges support from NSF Career Award 1149123 to Laurence Loewe and from an NHGRI training grant to the Genomic Sciences Training Program 5T32HG002760. Wai-Tong Louis Fan is supported by NSF Career Award DMS-1149312 to Sebastien Roch and by the Wisconsin Institute for Discovery at UW-Madison. Laurence Loewe is supported by NSF Career Award 1149123 and by the Wisconsin Institute for Discovery at UW-Madison.