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

Authors: Brian McLoone; Wai-Tong Louis Fan; Adam Pham; Rory Smead; Laurence Loewe
Journal: Complexity
 (2018)

---

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


---

## Body

## 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 strategyA orB with its partner, the column player, who can also play either strategyAorB.Payoff Matrix 1. A two-player, symmetric game is as follows:(1)ABAαβBγδ.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. Lettingrstand for the cost-to-benefit ratio in a given Snowdrift Game, this threshold quantity, which we call r∗, is approximately equal to 0.796. If r<r∗, the probability cooperators will fix tends to 1 as group sizes go to infinity. If r>r∗, 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 r=r∗.) 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 r 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. Letbrefer to thebenefit of getting home and letcrefer to thecost of shoveling snow, where b>c. Each driver can either get out of her car to shovel (C for “cooperate”) or stay in her car (D for “defect"). If both drivers cooperate, then each receives a net payoff of b-c/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 b-c, while the defector pays no cost, receiving a net payoff of b. 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 r stand for the ratio of the cost of cooperating in the Snowdrift Game to the benefit of doing so (r=c/b). We can thereby speak of the “r-value" of an instantiation of the game.Payoff Matrix 2. The Snowdrift Game is as follows (r=c/b):(2)CDCb-c2b-cDb0.Note that “cooperation” in this context refers to the shoveling snow behavior—that is, strategyC 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. StrategyC 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 numbern of discrete groups, indexed by j, so that j∈{1,2,…,n}. The size of a given group j (Nj) 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 groupjis 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 inj, 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/Nj. 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 groupj as follows. Let bj and cj index the parametersbandcfor a given groupj, let xj stand for the frequency of cooperators inj, and let 1-xj stand for the frequency of defectors inj. The fitness of a cooperator and a defector in groupj for a given value of xj are given, respectively, by(3)wC,j=bj-cj2xj+bj-cj1-xj,wD,j=bjxj.The average individual fitness of the members of a groupjis given by(4)w¯j=xjwC,j+1-xjwD,j.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. Lettingij stand for the number of cooperators in a groupj, we can represent the first transition as (ij→ij+1) and the second as (ij→ij-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, P, of each of these two transitions:(5)Pij⟶ij+1=wC,jw¯j×ijNj×Nj-ijNj,Pij⟶ij-1=wD,jw¯j×ijNj×Nj-ijNj.Since no other changes within a group are possible, the probability that the state of the group will not change is given by,(6)1-Pij⟶ij+1-Pij⟶ij-1.Were each group well-mixed and infinitely large, each groupj would converge to a stable, internal equilibrium frequency of cooperators (xj∗), which is given by the following (Hauert and Doebeli 2004 [16]):(7)xj∗=bj-cjbj-cj/2.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(xj=1) or go extinct (xj=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]).


### 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 groupj that is created at timet has the same frequency of cooperators asj att and has the same payoff matrix asj. (See Figure 1.)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.)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:(8)dw¯dx=c-2bx+2b-c,d2w¯dx2=c-2b<0.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 groupjwill replicate or die at the next time step. Letting Pj(+1) and Pj(-1) refer, respectively, to these two probabilities, we have(9)Pj+1=w¯jw¯M×1n×n-1n,Pj-1=w¯j′w¯M×1n×n-1n,where w¯j is the average fitness of the individuals in group j, given by equation (4),(10)w¯j′=1n-1∑k:k≠jw¯krefers to the average fitness of all groups other than group j, and w¯M is the average of the average individual fitness of each group in the metapopulation, which is given by(11)w¯M=1n∑j=1nw¯j.The probability that there will be no change in the metapopulation is given by(12)1-∑j=1nPj+1+Pj-1.


## 2.1. Within-Group Moran Process

Suppose there is a metapopulation composed of a finite numbern of discrete groups, indexed by j, so that j∈{1,2,…,n}. The size of a given group j (Nj) 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 groupjis 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 inj, 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/Nj. 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 groupj as follows. Let bj and cj index the parametersbandcfor a given groupj, let xj stand for the frequency of cooperators inj, and let 1-xj stand for the frequency of defectors inj. The fitness of a cooperator and a defector in groupj for a given value of xj are given, respectively, by(3)wC,j=bj-cj2xj+bj-cj1-xj,wD,j=bjxj.The average individual fitness of the members of a groupjis given by(4)w¯j=xjwC,j+1-xjwD,j.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. Lettingij stand for the number of cooperators in a groupj, we can represent the first transition as (ij→ij+1) and the second as (ij→ij-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, P, of each of these two transitions:(5)Pij⟶ij+1=wC,jw¯j×ijNj×Nj-ijNj,Pij⟶ij-1=wD,jw¯j×ijNj×Nj-ijNj.Since no other changes within a group are possible, the probability that the state of the group will not change is given by,(6)1-Pij⟶ij+1-Pij⟶ij-1.Were each group well-mixed and infinitely large, each groupj would converge to a stable, internal equilibrium frequency of cooperators (xj∗), which is given by the following (Hauert and Doebeli 2004 [16]):(7)xj∗=bj-cjbj-cj/2.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(xj=1) or go extinct (xj=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]).


## 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 groupj that is created at timet has the same frequency of cooperators asj att and has the same payoff matrix asj. (See Figure 1.)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.)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:(8)dw¯dx=c-2bx+2b-c,d2w¯dx2=c-2b<0.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 groupjwill replicate or die at the next time step. Letting Pj(+1) and Pj(-1) refer, respectively, to these two probabilities, we have(9)Pj+1=w¯jw¯M×1n×n-1n,Pj-1=w¯j′w¯M×1n×n-1n,where w¯j is the average fitness of the individuals in group j, given by equation (4),(10)w¯j′=1n-1∑k:k≠jw¯krefers to the average fitness of all groups other than group j, and w¯M is the average of the average individual fitness of each group in the metapopulation, which is given by(11)w¯M=1n∑j=1nw¯j.The probability that there will be no change in the metapopulation is given by(12)1-∑j=1nPj+1+Pj-1.


## 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 (Nj⩾1000). 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’sr-value. We therefore simulated the within-group evolutionary dynamics of nine groups, where each group was associated with a differentr-value. The range ofr-values that we simulated was r∈{0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9}, wherebwas set to 10 andcranged from 1 to 9 in increments of 1. We simulated a range of group sizes, Nj∈{10,100,1000,10,000}, and we allowed each group to evolve for 100 generations, where 1 generation for a group of size Nj is equal to N 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 theirr-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 ther-value we assigned to each group (r=.25; r=.67; r=.90). For each r-value, we ran simulations where Nj∈{10,100,1000}. We ran 500 independent simulations for each combination of r-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 n is equal to n 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 theirr-values. To explore the sensitivity of the metapopulation dynamics to ther-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 uniquer-value so that the averager-value was relatively low: the benefitb of cooperating was set to 10, whilec ranged from 0.25 to 6.25, in increments of 0.25, resulting in an initial averager-value of 0.33. In the second version, we performed the same analysis, except that we associated each group with a uniquer-value so that the averager-value was relatively high:b was set to 10, whilec ranged from 7.79 to 9.94, in increments of approximately 0.09 on average, resulting in an initial average r-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[xM]) 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[xM]) for each of these null models is given by the same equation:(13)ExM=1n∑r∈Rnrxr∗,where nr represents the number of groups with a given r-value in the metapopulation, xr∗ stands for the equilibrium frequency of cooperators in a group of type-r, and the notation r∈R indicates that we perform this operation for eachr-value that is in the set ofr-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.


## 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 (Nj⩾1000). 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’sr-value. We therefore simulated the within-group evolutionary dynamics of nine groups, where each group was associated with a differentr-value. The range ofr-values that we simulated was r∈{0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9}, wherebwas set to 10 andcranged from 1 to 9 in increments of 1. We simulated a range of group sizes, Nj∈{10,100,1000,10,000}, and we allowed each group to evolve for 100 generations, where 1 generation for a group of size Nj is equal to N 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 theirr-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 ther-value we assigned to each group (r=.25; r=.67; r=.90). For each r-value, we ran simulations where Nj∈{10,100,1000}. We ran 500 independent simulations for each combination of r-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 n is equal to n 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 theirr-values. To explore the sensitivity of the metapopulation dynamics to ther-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 uniquer-value so that the averager-value was relatively low: the benefitb of cooperating was set to 10, whilec ranged from 0.25 to 6.25, in increments of 0.25, resulting in an initial averager-value of 0.33. In the second version, we performed the same analysis, except that we associated each group with a uniquer-value so that the averager-value was relatively high:b was set to 10, whilec ranged from 7.79 to 9.94, in increments of approximately 0.09 on average, resulting in an initial average r-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[xM]) 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[xM]) for each of these null models is given by the same equation:(13)ExM=1n∑r∈Rnrxr∗,where nr represents the number of groups with a given r-value in the metapopulation, xr∗ stands for the equilibrium frequency of cooperators in a group of type-r, and the notation r∈R indicates that we perform this operation for eachr-value that is in the set ofr-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.


## 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’sr-value. We therefore simulated the within-group evolutionary dynamics of nine groups, where each group was associated with a differentr-value. The range ofr-values that we simulated was r∈{0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9}, wherebwas set to 10 andcranged from 1 to 9 in increments of 1. We simulated a range of group sizes, Nj∈{10,100,1000,10,000}, and we allowed each group to evolve for 100 generations, where 1 generation for a group of size Nj is equal to N 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 theirr-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 ther-value we assigned to each group (r=.25; r=.67; r=.90). For each r-value, we ran simulations where Nj∈{10,100,1000}. We ran 500 independent simulations for each combination of r-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 n is equal to n 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 theirr-values. To explore the sensitivity of the metapopulation dynamics to ther-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 uniquer-value so that the averager-value was relatively low: the benefitb of cooperating was set to 10, whilec ranged from 0.25 to 6.25, in increments of 0.25, resulting in an initial averager-value of 0.33. In the second version, we performed the same analysis, except that we associated each group with a uniquer-value so that the averager-value was relatively high:b was set to 10, whilec ranged from 7.79 to 9.94, in increments of approximately 0.09 on average, resulting in an initial average r-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[xM]) 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[xM]) for each of these null models is given by the same equation:(13)ExM=1n∑r∈Rnrxr∗,where nr represents the number of groups with a given r-value in the metapopulation, xr∗ stands for the equilibrium frequency of cooperators in a group of type-r, and the notation r∈R indicates that we perform this operation for eachr-value that is in the set ofr-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 X(N)=(Xt(N))t≥0 that solves the following SDE:(14)dXtN=μXtNdt+σXtNNdWt,where N is the size of the group (note that we have left out the j subscripts from N for ease of reading), Wt is the Wiener process (standard Brownian motion), and μ and σ are the “drift" and the “volatility” functions given, respectively, by(15)μx=T+x-T-x,σx=T+x+T-x,where(16)T+x=αx+β1-xw¯x1-x,T-x=γx+δ1-xw¯x1-x.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:(17)dxtdt=μxt.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.
LetN be the total number of individuals and let CN(i) be the number of cooperators at the ith step in the Moran process. Suppose the initial fraction of cooperative players CN(0)/N tends to x0∈[0,1] as N→∞. Then, for all ϵ>0 and T<∞, we have the following.(1)
Deterministic approximation: Moran process versus ODE:(18)limN→∞⁡Psupt∈0,T⁡CNtNN-xt<ϵ=1,wherex(t) is the solution to the ODE (17) starting with x(0)=x0, “sup" is the supremum, and [s] is the largest integer smaller than or equal to s.(2)
Diffusion approximation: Moran process versus SDE: there exists a diffusion process Xt(N) defined on the same probability space as that of the Moran process, which starts at X0(N)=x0, solves (14), and is such that(19)limN→∞⁡Psupt∈0,T⁡CNtNN-XtN<aNϵ=1for any sequenceaN such that limN→∞⁡(log⁡N/NaN)=0.Remark 2.
Convergence (18) says that when the population size N is large, most sample paths of the process CN([tN])/N stay within distance ϵ from the solution x(t) of the ODE for any time t∈[0,T]. That is, once we specify the error tolerance ϵ and the duration T, we can choose N large enough so that |CN([tN])/N-x(t)|<ϵ uniformly in t∈[0,T] with a probability close to one.Remark 3.
Convergence (19) says that X(N) is a better approximation than x(t) and precisely quantifies the extent to which it is a better approximation. Specifically, (19) holds as long as aN→0 slower than log⁡N/N. This roughly says that supt∈[0,T]CNtN/N-XtN→0 at a rate of order at least O(log⁡N/N).Remark 4.
Lemma1 remains true if we replace the discrete time Moran process by the continuous time version of the Moran process with rate N. 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. Leth(x) be the probability that the diffusion process X(N) hits 1 before 0 (i.e., cooperators fix), given that the process starts at X0(N)=x. In mathematical notation, we let τ≔inf⁡{t≥0:Xt(N)=0or1} be the time it takes for the population to become monomorphic, and we let Px be the probability law of X(N) starting at x∈[0,1]. Then,(20)hx≔PxXτN=1.Note that h(x) depends on N 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.
Leth(x) be the probability that cooperation fixes, defined in (20). Then, (21)hx=∫0xeNfydy∫01eNfydy,where f(y)=-2∫0y(Cu+D)/(Au+B)du and(22)A=α+γ-β-δ,B=β+δ,C=α-β-γ+δ,D=β-δ.


#### 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 Lemma5, we discover a threshold quantity r∗≈0.796 such that when r<r∗ the probability cooperators will fix tends to 1 as group size goes to infinity (Figure 2). Likewise, when r>r∗, the probability cooperators will go extinct tends to 1 as group size goes to infinity (Figure 3). More complicated dynamics occur when r is roughly equal to r∗ (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., Nj=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.Figure 2
The plot ofh(x). When r=0.3, which isbelowr∗, the probability that cooperators will fix gets very large, regardless of the starting frequency of cooperators, as group size increases.Figure 3
The plot ofh(x). When r=0.9, which isabover∗, the probability that cooperators will fix is negligible, regardless of the starting frequency of cooperators, as group size increases.Figure 4
The plot ofh(x). When r is almost exactly equal to r∗, the probability that cooperators will fix is approximately r∗/4≈0.2, regardless of the starting frequency of cooperators, as group size increases.Theorem 6 (snowdrift game threshold theorem).
There exists a threshold numberr∗∈(0,1) such that, for all x∈(0,1), we have the limiting probability(23)limN→∞⁡hx=0,ifr>r∗,r∗4,ifr=r∗,1,ifr<r∗,where h(x)≔Px(Xτ(N)=1) is the probability that cooperation goes to fixation, given that the starting point is X0(N)=x∈[0,1]. Furthermore, r∗≈0.79559250931 is the unique solution of the following equation:(24)1-r24-21-rln⁡4-r2-2r=0on the positive real axis.Remark 7.
The endpoints are given byh(0)=0 and h(1)=1 for all N.Proof.
Asymptotic analysis (see, e.g., Murray 1984 [40]) reveals that the large N behavior of the probability h(x) is governed by the maximum of the function f. Recall the explicit formula and the notation in Lemma 5. For the Snowdrift Game, we have(25)fy=-2∫0yCu+DAu+Bdu=-2CAy+AD-BCA2ln⁡1+ABy,f′y=-2Cy+DAy+B,f′′y=2AD-BCAy+B2.Clearly, f′′>0 and the maximum of f is attained at the endpoints f(0) and f(1). To derive (23), note that(26)f1-f0=2b2A2θrwhereθr=1-r24-21-rln⁡4-r2-2r.Elementary calculus shows that f(0)<f(1) when r∈0,r∗, f0=f1 when r=r∗, and f0>f1 when r∈r∗,1, where r∗ is the solution of θ(r)=0 on the positive real axis, which was found to be r∗≈0.79559250931. When f0=f1 (i.e., r=r∗), we have, by the Laplace method (see, e.g., Murray 1984 [40]),(27)limN→∞⁡hx=1f′0-1f′1-11f′0=r∗4∀x∈0,1.Result (23) then follows.


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

The results of our purely within-group simulations are presented in Figure5. The jagged lines are the observed frequency of cooperators within a group with a particularr-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).Figure 5
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 theirr-values. The range ofr-values that were simulated was r∈{0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9}, wherebwas set to 10 andcranged from 1 to 9 in increments of 1. The jagged lines represent the frequency of cooperators in a group. Asr-values get larger, the equilibrium frequency of cooperators gets smaller; thus, r=0.1 is black, while r=0.9 is teal. Each horizontal line of a corresponding color represents the equilibrium frequency of cooperators for thatr-value given by (7).
(a)
N
=
10
(b)
N
=
100
(c)
N
=
1000
(d)
N
=
10
,
000In 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

Theorem6 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 anr-value below (above) r∗. 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 r-value, relative to r∗, to substantially influence the probability cooperators will fix or go extinct.When groups’r-values straddle r∗, 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 Figures6 and 7.Figure 6
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 r=.25; panels (b, e, h) show r=.67; panels (c, f, i) show r=.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).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 7
Two-level Moran model with variation in r-values. Each panel represents simulation results for our model when groups in the metapopulation vary in their initial r-values. The panels (a, e, c) show the results for r=.33, while the panels (b, d, f) show r=0.89. Each panel represents summary results for 100 independent simulations with populations of 10 (a, b), 100 (c, d), and 1000 (e, f) individuals. 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, respectively, 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 greater than that predicted by the null models in all cases, sometimes substantially.
(a)
(b)
(c)
(d)
(e)
(f)


## 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 X(N)=(Xt(N))t≥0 that solves the following SDE:(14)dXtN=μXtNdt+σXtNNdWt,where N is the size of the group (note that we have left out the j subscripts from N for ease of reading), Wt is the Wiener process (standard Brownian motion), and μ and σ are the “drift" and the “volatility” functions given, respectively, by(15)μx=T+x-T-x,σx=T+x+T-x,where(16)T+x=αx+β1-xw¯x1-x,T-x=γx+δ1-xw¯x1-x.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:(17)dxtdt=μxt.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.
LetN be the total number of individuals and let CN(i) be the number of cooperators at the ith step in the Moran process. Suppose the initial fraction of cooperative players CN(0)/N tends to x0∈[0,1] as N→∞. Then, for all ϵ>0 and T<∞, we have the following.(1)
Deterministic approximation: Moran process versus ODE:(18)limN→∞⁡Psupt∈0,T⁡CNtNN-xt<ϵ=1,wherex(t) is the solution to the ODE (17) starting with x(0)=x0, “sup" is the supremum, and [s] is the largest integer smaller than or equal to s.(2)
Diffusion approximation: Moran process versus SDE: there exists a diffusion process Xt(N) defined on the same probability space as that of the Moran process, which starts at X0(N)=x0, solves (14), and is such that(19)limN→∞⁡Psupt∈0,T⁡CNtNN-XtN<aNϵ=1for any sequenceaN such that limN→∞⁡(log⁡N/NaN)=0.Remark 2.
Convergence (18) says that when the population size N is large, most sample paths of the process CN([tN])/N stay within distance ϵ from the solution x(t) of the ODE for any time t∈[0,T]. That is, once we specify the error tolerance ϵ and the duration T, we can choose N large enough so that |CN([tN])/N-x(t)|<ϵ uniformly in t∈[0,T] with a probability close to one.Remark 3.
Convergence (19) says that X(N) is a better approximation than x(t) and precisely quantifies the extent to which it is a better approximation. Specifically, (19) holds as long as aN→0 slower than log⁡N/N. This roughly says that supt∈[0,T]CNtN/N-XtN→0 at a rate of order at least O(log⁡N/N).Remark 4.
Lemma1 remains true if we replace the discrete time Moran process by the continuous time version of the Moran process with rate N. 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. Leth(x) be the probability that the diffusion process X(N) hits 1 before 0 (i.e., cooperators fix), given that the process starts at X0(N)=x. In mathematical notation, we let τ≔inf⁡{t≥0:Xt(N)=0or1} be the time it takes for the population to become monomorphic, and we let Px be the probability law of X(N) starting at x∈[0,1]. Then,(20)hx≔PxXτN=1.Note that h(x) depends on N 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.
Leth(x) be the probability that cooperation fixes, defined in (20). Then, (21)hx=∫0xeNfydy∫01eNfydy,where f(y)=-2∫0y(Cu+D)/(Au+B)du and(22)A=α+γ-β-δ,B=β+δ,C=α-β-γ+δ,D=β-δ.


### 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 Lemma5, we discover a threshold quantity r∗≈0.796 such that when r<r∗ the probability cooperators will fix tends to 1 as group size goes to infinity (Figure 2). Likewise, when r>r∗, the probability cooperators will go extinct tends to 1 as group size goes to infinity (Figure 3). More complicated dynamics occur when r is roughly equal to r∗ (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., Nj=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.Figure 2
The plot ofh(x). When r=0.3, which isbelowr∗, the probability that cooperators will fix gets very large, regardless of the starting frequency of cooperators, as group size increases.Figure 3
The plot ofh(x). When r=0.9, which isabover∗, the probability that cooperators will fix is negligible, regardless of the starting frequency of cooperators, as group size increases.Figure 4
The plot ofh(x). When r is almost exactly equal to r∗, the probability that cooperators will fix is approximately r∗/4≈0.2, regardless of the starting frequency of cooperators, as group size increases.Theorem 6 (snowdrift game threshold theorem).
There exists a threshold numberr∗∈(0,1) such that, for all x∈(0,1), we have the limiting probability(23)limN→∞⁡hx=0,ifr>r∗,r∗4,ifr=r∗,1,ifr<r∗,where h(x)≔Px(Xτ(N)=1) is the probability that cooperation goes to fixation, given that the starting point is X0(N)=x∈[0,1]. Furthermore, r∗≈0.79559250931 is the unique solution of the following equation:(24)1-r24-21-rln⁡4-r2-2r=0on the positive real axis.Remark 7.
The endpoints are given byh(0)=0 and h(1)=1 for all N.Proof.
Asymptotic analysis (see, e.g., Murray 1984 [40]) reveals that the large N behavior of the probability h(x) is governed by the maximum of the function f. Recall the explicit formula and the notation in Lemma 5. For the Snowdrift Game, we have(25)fy=-2∫0yCu+DAu+Bdu=-2CAy+AD-BCA2ln⁡1+ABy,f′y=-2Cy+DAy+B,f′′y=2AD-BCAy+B2.Clearly, f′′>0 and the maximum of f is attained at the endpoints f(0) and f(1). To derive (23), note that(26)f1-f0=2b2A2θrwhereθr=1-r24-21-rln⁡4-r2-2r.Elementary calculus shows that f(0)<f(1) when r∈0,r∗, f0=f1 when r=r∗, and f0>f1 when r∈r∗,1, where r∗ is the solution of θ(r)=0 on the positive real axis, which was found to be r∗≈0.79559250931. When f0=f1 (i.e., r=r∗), we have, by the Laplace method (see, e.g., Murray 1984 [40]),(27)limN→∞⁡hx=1f′0-1f′1-11f′0=r∗4∀x∈0,1.Result (23) then follows.


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

The results of our purely within-group simulations are presented in Figure5. The jagged lines are the observed frequency of cooperators within a group with a particularr-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).Figure 5
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 theirr-values. The range ofr-values that were simulated was r∈{0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9}, wherebwas set to 10 andcranged from 1 to 9 in increments of 1. The jagged lines represent the frequency of cooperators in a group. Asr-values get larger, the equilibrium frequency of cooperators gets smaller; thus, r=0.1 is black, while r=0.9 is teal. Each horizontal line of a corresponding color represents the equilibrium frequency of cooperators for thatr-value given by (7).
(a)
N
=
10
(b)
N
=
100
(c)
N
=
1000
(d)
N
=
10
,
000In 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.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 X(N)=(Xt(N))t≥0 that solves the following SDE:(14)dXtN=μXtNdt+σXtNNdWt,where N is the size of the group (note that we have left out the j subscripts from N for ease of reading), Wt is the Wiener process (standard Brownian motion), and μ and σ are the “drift" and the “volatility” functions given, respectively, by(15)μx=T+x-T-x,σx=T+x+T-x,where(16)T+x=αx+β1-xw¯x1-x,T-x=γx+δ1-xw¯x1-x.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:(17)dxtdt=μxt.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.
LetN be the total number of individuals and let CN(i) be the number of cooperators at the ith step in the Moran process. Suppose the initial fraction of cooperative players CN(0)/N tends to x0∈[0,1] as N→∞. Then, for all ϵ>0 and T<∞, we have the following.(1)
Deterministic approximation: Moran process versus ODE:(18)limN→∞⁡Psupt∈0,T⁡CNtNN-xt<ϵ=1,wherex(t) is the solution to the ODE (17) starting with x(0)=x0, “sup" is the supremum, and [s] is the largest integer smaller than or equal to s.(2)
Diffusion approximation: Moran process versus SDE: there exists a diffusion process Xt(N) defined on the same probability space as that of the Moran process, which starts at X0(N)=x0, solves (14), and is such that(19)limN→∞⁡Psupt∈0,T⁡CNtNN-XtN<aNϵ=1for any sequenceaN such that limN→∞⁡(log⁡N/NaN)=0.Remark 2.
Convergence (18) says that when the population size N is large, most sample paths of the process CN([tN])/N stay within distance ϵ from the solution x(t) of the ODE for any time t∈[0,T]. That is, once we specify the error tolerance ϵ and the duration T, we can choose N large enough so that |CN([tN])/N-x(t)|<ϵ uniformly in t∈[0,T] with a probability close to one.Remark 3.
Convergence (19) says that X(N) is a better approximation than x(t) and precisely quantifies the extent to which it is a better approximation. Specifically, (19) holds as long as aN→0 slower than log⁡N/N. This roughly says that supt∈[0,T]CNtN/N-XtN→0 at a rate of order at least O(log⁡N/N).Remark 4.
Lemma1 remains true if we replace the discrete time Moran process by the continuous time version of the Moran process with rate N. 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. Leth(x) be the probability that the diffusion process X(N) hits 1 before 0 (i.e., cooperators fix), given that the process starts at X0(N)=x. In mathematical notation, we let τ≔inf⁡{t≥0:Xt(N)=0or1} be the time it takes for the population to become monomorphic, and we let Px be the probability law of X(N) starting at x∈[0,1]. Then,(20)hx≔PxXτN=1.Note that h(x) depends on N 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.
Leth(x) be the probability that cooperation fixes, defined in (20). Then, (21)hx=∫0xeNfydy∫01eNfydy,where f(y)=-2∫0y(Cu+D)/(Au+B)du and(22)A=α+γ-β-δ,B=β+δ,C=α-β-γ+δ,D=β-δ.


## 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 Lemma5, we discover a threshold quantity r∗≈0.796 such that when r<r∗ the probability cooperators will fix tends to 1 as group size goes to infinity (Figure 2). Likewise, when r>r∗, the probability cooperators will go extinct tends to 1 as group size goes to infinity (Figure 3). More complicated dynamics occur when r is roughly equal to r∗ (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., Nj=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.Figure 2
The plot ofh(x). When r=0.3, which isbelowr∗, the probability that cooperators will fix gets very large, regardless of the starting frequency of cooperators, as group size increases.Figure 3
The plot ofh(x). When r=0.9, which isabover∗, the probability that cooperators will fix is negligible, regardless of the starting frequency of cooperators, as group size increases.Figure 4
The plot ofh(x). When r is almost exactly equal to r∗, the probability that cooperators will fix is approximately r∗/4≈0.2, regardless of the starting frequency of cooperators, as group size increases.Theorem 6 (snowdrift game threshold theorem).
There exists a threshold numberr∗∈(0,1) such that, for all x∈(0,1), we have the limiting probability(23)limN→∞⁡hx=0,ifr>r∗,r∗4,ifr=r∗,1,ifr<r∗,where h(x)≔Px(Xτ(N)=1) is the probability that cooperation goes to fixation, given that the starting point is X0(N)=x∈[0,1]. Furthermore, r∗≈0.79559250931 is the unique solution of the following equation:(24)1-r24-21-rln⁡4-r2-2r=0on the positive real axis.Remark 7.
The endpoints are given byh(0)=0 and h(1)=1 for all N.Proof.
Asymptotic analysis (see, e.g., Murray 1984 [40]) reveals that the large N behavior of the probability h(x) is governed by the maximum of the function f. Recall the explicit formula and the notation in Lemma 5. For the Snowdrift Game, we have(25)fy=-2∫0yCu+DAu+Bdu=-2CAy+AD-BCA2ln⁡1+ABy,f′y=-2Cy+DAy+B,f′′y=2AD-BCAy+B2.Clearly, f′′>0 and the maximum of f is attained at the endpoints f(0) and f(1). To derive (23), note that(26)f1-f0=2b2A2θrwhereθr=1-r24-21-rln⁡4-r2-2r.Elementary calculus shows that f(0)<f(1) when r∈0,r∗, f0=f1 when r=r∗, and f0>f1 when r∈r∗,1, where r∗ is the solution of θ(r)=0 on the positive real axis, which was found to be r∗≈0.79559250931. When f0=f1 (i.e., r=r∗), we have, by the Laplace method (see, e.g., Murray 1984 [40]),(27)limN→∞⁡hx=1f′0-1f′1-11f′0=r∗4∀x∈0,1.Result (23) then follows.


## 4.1.3. Simulations of Within-Group Moran Process

The results of our purely within-group simulations are presented in Figure5. The jagged lines are the observed frequency of cooperators within a group with a particularr-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).Figure 5
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 theirr-values. The range ofr-values that were simulated was r∈{0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9}, wherebwas set to 10 andcranged from 1 to 9 in increments of 1. The jagged lines represent the frequency of cooperators in a group. Asr-values get larger, the equilibrium frequency of cooperators gets smaller; thus, r=0.1 is black, while r=0.9 is teal. Each horizontal line of a corresponding color represents the equilibrium frequency of cooperators for thatr-value given by (7).
(a)
N
=
10
(b)
N
=
100
(c)
N
=
1000
(d)
N
=
10
,
000In 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

Theorem6 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 anr-value below (above) r∗. 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 r-value, relative to r∗, to substantially influence the probability cooperators will fix or go extinct.When groups’r-values straddle r∗, 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 Figures6 and 7.Figure 6
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 r=.25; panels (b, e, h) show r=.67; panels (c, f, i) show r=.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).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 7
Two-level Moran model with variation in r-values. Each panel represents simulation results for our model when groups in the metapopulation vary in their initial r-values. The panels (a, e, c) show the results for r=.33, while the panels (b, d, f) show r=0.89. Each panel represents summary results for 100 independent simulations with populations of 10 (a, b), 100 (c, d), and 1000 (e, f) individuals. 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, respectively, 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 greater than that predicted by the null models in all cases, sometimes substantially.
(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 (N→∞) and then time (t→∞) to arrive at the long-term behavior of an evolutionary system. To derive the value of r∗, we have switched the order in which these limits are taken. We take (t→∞) to derive the formula for the fixation probability of cooperation. Then we take (N→∞) 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 N⩾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.


### 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 (N=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 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 (N⩾100), 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 r∗ 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 ofr-values of the groups in those simulations.When we allow the groups in the metapopulation to initially vary in theirr-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 x=1 will be more likely to replicate. When group sizes are small (e.g., N=10), cooperators will fix in many of these groups (Figures 7(a) and 7(b)). But even as group size becomes larger (N⩾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 r-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 isqualitatively characterized by the Snowdrift Game, it would in fact be somewhat surprising if there were noquantitative 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 x∗ 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 x∗ 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.


## 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 (N→∞) and then time (t→∞) to arrive at the long-term behavior of an evolutionary system. To derive the value of r∗, we have switched the order in which these limits are taken. We take (t→∞) to derive the formula for the fixation probability of cooperation. Then we take (N→∞) 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 N⩾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.


## 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 (N=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 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 (N⩾100), 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 r∗ 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 ofr-values of the groups in those simulations.When we allow the groups in the metapopulation to initially vary in theirr-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 x=1 will be more likely to replicate. When group sizes are small (e.g., N=10), cooperators will fix in many of these groups (Figures 7(a) and 7(b)). But even as group size becomes larger (N⩾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 r-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 isqualitatively characterized by the Snowdrift Game, it would in fact be somewhat surprising if there were noquantitative 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 x∗ 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 x∗ 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.,N=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.



---

*Source: 9836150-2018-02-13.xml*