A unifying framework reveals key properties of multilevel selection

Journal of Theoretical Biology 341 (2014) 41–52

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

A unifying framework reveals key properties of multilevel selection Shishi Luo n Department of Mathematics, Duke University, Durham, NC 27708, USA

H I G H L I G H T S








A ball-and-urn system is presented as a conceptual tool for understanding multilevel selection. General properties of multilevel selection with minimal dependence on parameter values are derived. Fast group-level events relative to individual-level events enhances group-level selection. Mutation enhances individual-level selection. Smaller group sizes or a larger number of groups enhance group-level selection.

art ic l e i nf o

a b s t r a c t

Article history: Received 8 April 2013 Received in revised form 2 September 2013 Accepted 16 September 2013 Available online 1 October 2013

Natural selection can act at multiple biological levels, often in opposing directions. Viral evolution is an important example, with selection occurring both within infected hosts and between hosts via transmission. A fast-replicating virus may outcompete a slower strain within the same host, however, if rapid viral replication incapacitates the host, this fast-replicating virus may not be transmitted as frequently as its slower counterpart. Such examples of antagonistic multilevel selection arise across biological taxa and scales, from microbial public goods production to male mating strategies. A general formalism for describing and analyzing these diverse systems can identify their common underlying properties. Here I introduce such a unifying framework, which can be intuitively visualized as a stochastic ball-and-urn process. This ball-and-urn process illustrates the dynamics of antagonistic selective forces and allows the systematic derivation of properties with little or no dependence on model parameterization. These properties are consistent with previous studies, both theoretical and empirical, of multilevel selection. In particular I show that selection at the group level is favored when group-level events occur frequently relative to individual-level events, when there is little or no mutation, and when there are many groups relative to the number of individuals in each group. This approach demonstrates how multilevel selection can be understood as a general biological phenomenon, and identifies recurring characteristics that may be independent of specific biological contexts. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Stochastic processes Group selection Moran process Public goods game

1. Introduction Multilevel selection theory is a growing field in evolutionary biology which aims to understand biological systems in the framework of selection acting, possibly antagonistically, at multiple hierarchical levels (Okasha, 2006; Wilson and Wilson, 2007; Kerr, 2009). Examples of systems which exhibit antagonistic multilevel selection are diverse (Table 1). In addition to rapid replication in viruses, the production of a beneficial public good (Turner and Chao, 1999; Rumbaugh et al., 2009), such as an extracellular protein that impairs the host's immune response (Bonhoeffer and Nowak, n Corresponding author. Current address: Center for Nonlinear Studies (CNLS) and Theoretical Biology and Biophysics (T-6), Los Alamos National Laboratory MS B258, Los Alamos, NM 87545, USA. Tel.: þ 1 505 606 0211. E-mail address: [email protected]

0022-5193/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jtbi.2013.09.024

1994), is another pathogen trait that undergoes antagonistic selection. The production of the protein facilitates a successful infection, increasing the chance of transmission to a new host. However, within each host, producer pathogens are at a competitive disadvantage against non-producer pathogens because the latter benefit from the protein without paying the cost of producing it. A similar situation arises in the production of public goods in bacteria (Griffin et al., 2004; Chuang et al., 2009, 2010) and yeast (Gore et al., 2009), where producer individuals help increase overall population growth, but can be locally outcompeted by free-loading nonproducers. Antagonistic multilevel selection is also posited to act on plant leaf size, where plants with larger leaves acquire more resources, but to the detriment of plants in the same stand (Stevens et al., 1995; Weinig et al., 2007). Similarly, water strider males that practice an aggressive mating strategy increase their number of successful matings relative to a more restrained strategy, but

42

S. Luo / Journal of Theoretical Biology 341 (2014) 41–52

Table 1 Empirical examples of antagonistic multilevel systems. Trait

Pathogen replication rate

↑

Impact on fitness at Lower level

Higher level

Increased intrahost competitive advantage

Fewer transmissions due to host mortality/ Myxo-rabbit (Levin and Pimentel, 1981), bacteria–plasmid morbidity (Paulsson, 2002), bacteria–phage (Kerr et al., 2006)

Public goods production

✓ Outcompeted by nonproducers

Multicellularity

✓ Vulnerable to exploitation by ‘cheaters’ ↑ Greater individual resource acquisition ↑ Lower number of successful matings ↑ Less intrahost asexual replication

Plant leaf size Male mating restraint Transmission life stage production

System

Fewer transmissions due to early detection Human papillomavirus (Orlando et al., 2012) and clearance by immune system Higher overall population growth Bacteria (Chuang et al., 2009, 2010; Griffin et al., 2004), yeast (Gore et al., 2009) More transmissions due to successful Bacteria–phage (Turner and Chao, 1999), bacteria–mouse (Rumbaugh infection (host pathogen systems) et al., 2009), various pathogens (Bonhoeffer and Nowak, 1994) Increased access to resources/overall Biofilms (Rainey and Rainey, 2003), fruiting bodies (Velicer et al., 2000), survival cancer (Merlo et al., 2006) Less resource for plants in same stand

Jewelweed (Stevens et al., 1995), Arabdopsis (Weinig et al., 2007)

Higher overall reproductive fitness of females in group Higher transmission success

Amazon molly (Kokko and Heubel, 2011), water strider (Eldakar et al., 2009, 2010) Malaria (Pollitt et al., 2011)

↑ denotes an increase of a continuous trait. ✓ denotes the presence of a discrete trait.

are disadvantageous in a larger context because fewer females reside in communities with more aggressive males (Eldakar et al., 2009, 2010). Mathematical models have been developed for many specific examples of multilevel selection (Levin and Pimentel, 1981; Levin and Bull, 1994; Bonhoeffer and Nowak, 1994; Paulsson, 2002; Day and Proulx, 2004). However, the common underlying feature of these examples – a conflict between what is advantageous at the local scale and what is advantageous at a larger scale – motivates a single unifying framework to capture fundamental properties of these systems. A historically popular tool used by evolutionary biologists to address this question is the Price equation (Sober and Wilson, 1998; Rice, 2004; Okasha, 2006; Kerr, 2009). Price (1972) demonstrated that for a large class of evolutionary systems, the change over one generation of the population frequency of a gene undergoing multilevel selection can be decomposed into two components. One component represents selection at the group level and the other, selection at the individual level. Price also suggested that the equation can be rewritten in terms of regression coefficients and that the intensity of selection at the two levels could therefore be empirically measured. A major advantage of the Price equation is that it holds for almost every system that undergoes evolutionary dynamics. However, with this universal applicability comes several shortcomings. First, the Price equation is not in itself a mathematical model; beyond establishing a decomposition it cannot offer more insight without further assumptions. In particular, additional assumptions are required to extrapolate evolutionary dynamics beyond a single generation. Second, the mathematical ambiguity of the notation used in Price's original formulation can lead to confusion and improper mathematical use of the equation (van Veelen, 2005; van Veelen et al., 2012). Finally, it is unclear whether the regression coefficient form of the Price equation can be applied empirically to make predictions, even in a synthetic bacterial system (Chuang et al., 2010). Overall, these shortcomings make it difficult to use the Price equation, in its most general form, to systematically answer questions about evolutionary dynamics beyond one generation. These limitations of the Price equation suggest that complementary mathematical formulations of multilevel selection, where assumptions about population structure and reproduction are more explicit, may be more illuminating and better suited to

investigate multilevel evolutionary dynamics. Indeed Kimura (1983) used a partial differential equation to model antagonistic intra- and inter-group selection under mutation, selection, and migration. The mathematical steps Kimura took to obtain the equilibrium distribution of the group-advantageous type are, however, not completely transparent and therefore it is not clear whether the results he derives from the equilibrium distribution apply to the original differential equation. More recently, Traulsen and Nowak (2006) introduced a stochastic nested Moran model formulation for a public goods game, analytically deriving a simple condition for when group-level selection dominates individual-level selection. To obtain this condition, the authors assume that group splitting events occur rarely, which is not generally the case in the examples presented in Table 1. It would be useful, therefore, to have such a condition where this assumption is relaxed. Simon (2010) introduced a deterministic model of multilevel selection which, along with relaxing the assumption in Traulsen and Nowak (2006), allows for group events such as fission and extinction in addition to selection and migration. This model serves as a formal framework under which concrete definitions of individual and group selection may be established. In its full generality, however, it is difficult to analytically study its properties; computational approaches are often required to obtain results. In this paper, I introduce a way of viewing multilevel selection that builds on these previous models and overcomes several of their limitations. The idea behind this approach is to focus on the essential property of antagonistic multilevel selection: individuallevel selection favors one type, group-level selection favors another type, and use a simple mathematical formulation to systematically explore the consequences of this property. A major advantage of this formulation is that it can be intuitively visualized as a stochastic ball-and-urn process, and therefore has value as a conceptual and visual tool as well as a mathematical one. With this formulation, I deduce explicit relationships between model parameters and the balance of selection at the two scales. I also derive an integro-partial differential equation which represents the change of the population distribution of group-advantageous types over time. This differential equation further reveals the nature of the competing forces of selection as being a flux at the individual level and a forcing term at the group level, with a single parameter controlling the balance of the forces. The majority of

S. Luo / Journal of Theoretical Biology 341 (2014) 41–52

results I obtain are independent of model parameterization. Where approximations or assumptions about parameters are made, their effects are investigated and understood.

2. A simple conceptual model I use individual to refer to the entity at the lowest level of selection (e.g. virus, cell, plant) and group to refer to a collection of individuals (e.g. virions infecting a host, a multicellular organism, plants in a stand). Individuals can be one of the two types: those, such as a fast-replicating/low-transmission virus, which are selectively advantageous at the individual level (type I for individual) and those, such as a slow-replicating/high-transmission virus, which are selectively advantageous at the group level (type G for group). In this formulation there is a fixed number n of individuals in each group and a fixed number m of groups in the population (Fig. 1a). This fixed population size assumption appears frequently in theoretical population genetics, for example in the Wright– Fisher and Moran models (Durrett, 2008), and allows the opposing forces of selection to be studied without population dynamics confounding the conclusions. Replication and selection occur concurrently at the individual and group level according to the Moran process (Durrett, 2008) as follows (further mathematical details are in Appendix A). Type I individuals replicate at rate 1 þ s and type G individuals at rate 1. When an individual gives birth, another individual in the same group is selected uniformly at random to die, maintaining the number of individuals in the group at n (Fig. 1b). To reflect the antagonism at the higher level of selection, groups replicate at a rate which increases with the number of type G individuals they contain. For simplicity, and because it corresponds to a public goods game (Appendix B), this rate is wð1 þ rk=nÞ, where k=n is the fraction of individuals in the group that are type G, r is the selection coefficient at the group level, and w is the ratio of the rate of group-level events to the rate of individual-level events.

43

As with the individual level, the population of groups is maintained at m by selecting a group uniformly at random to die whenever a group replicates (Fig. 1c). The offspring of groups are assumed to be identical to their parent. These above assumptions may be relaxed or interpreted more broadly, a point I will return to in the discussion. This multilevel system is equivalent to a ball-and-urn process (Fig. 1d). Instead of considering each group as a subpopulation of individuals, a ball is placed in urn k for each group containing k type G individuals. This conveniently translates events which act at two different scales, individual and group, into two ways of moving a ball at the same scale. Animations of the dynamics of this stochastic ball-and-urn process (Movie 1) depict individual balls moving between urns in a seemingly random manner with the overall collection of balls displaying a slight tendency to the left or right ends of the row of urns. To separate the underlying behavior of this system from the stochastic noise, I derive a differential equation which applies when the number of groups m and number of individuals per group n are very large (Champagnat et al., 2006). The corresponding behavior in the actual (finite m and n) ball-and-urn process is verified using stochastic simulations (Appendix A.2). This formulation of multilevel selection has several advantages. First, in contrast to models which assume discrete time steps (Rice, 2004), this continuous time stochastic process handles both individual- and group-level events in parallel, where the generation times of the two levels need not be synchronized. w¼a means that, on an average, a group replicates a times for every time all of its n constituent individuals replicate. As an example, if w¼0.1, then all individuals in a group undergo on an average ten generations of replication and selection before the group as whole replicates. Importantly, this means that in the time between when a group is newly formed and when it replicates, it undergoes a period of maturation where individual-level selection acts. Thus all groups constantly experience individual-level selection that decreases the representation of group-advantageous types. The parameter w

Fig. 1. Schematic of the ball-and-urn system and its approximations. (a) A population of m¼ 3 groups, each with n¼ 3 individuals of either type G (blue) or type I (yellow). (b) A type I (yellow) individual replicates in group 3 and a type G (blue) individual is chosen uniformly at random from group 3 to die. (c) Group 1 replicates and produces group 2′. Group 2 is chosen uniformly at random to die. (d) The states in (a)–(c) mapped onto the ball-and-urn framework. Left: Group 2 has no type G (blue) individuals, represented by ball 2 in urn 0. Similarly, group 3 is represented by ball 3 in urn 2 and group 1 by ball 1 in urn 3. Middle: The number of blue individuals in group 3 decreases from two to one, therefore ball 3 moves to urn 1. Right: A group with zero blue individuals dies, while a group with three blue individuals is born. Therefore ball 2 leaves urn 0 and appears in urn 3 as ball 2′. (e) Left: μm;n is the Markov process that describes the microscopic ball-and-urn process. Middle: As the number of balls (m) increases to t infinity, μm;n -μnt , where μnt is described by a system of ordinary differential equations, each equation describing the number of balls in each urn. Right: As the number of urns t (n) increases to infinity, and with a rescaling of the state space to ½0; 1, μnt -μt where μt is described by Eq. (1). Details in Appendix A.1.

44

S. Luo / Journal of Theoretical Biology 341 (2014) 41–52

Table 2 Empirical examples where parameters can be measured. Parameter

Empirical measurement

System

s

Relative growth rates of two types competing in a mixed culture

Microbial systems (Kerr et al., 2006; Turner and Chao, 1999; Chuang et al., 2009; Gore et al., 2009; Rainey and Rainey, 2003; Velicer et al., 2000) Host–pathogen (Rumbaugh et al., 2009) Water strider (Eldakar et al., 2009, 2010) Microbial systems (Bouma and Lenski, 1988; Chuang et al., 2009; Kerr et al., 2006; Turner and Chao, 1999; Velicer et al., 2000) Water strider (Eldakar et al., 2009, 2010) Microbial systems (Griffin et al., 2004; Kerr et al., 2006; Gore et al., 2009; Brockhurst et al., 2007; Ratcliff et al., 2012 ) e.g. Water strider (Eldakar et al., 2009, 2010)

r

w

Relative growth rates of two types during coinfection in vivo Relative number of successful matings per male Relative growth rates of two types growing separately in pure culture (the fitness of a pure type G group is 1þ r) Number of females in a subpopulation Duration of time between dilutions of cultures divided by duration of microbe generation Equals one when a subpopulation ‘replicates’ after each individual in the subpopulation has produced an offspring

encodes how frequently group-level events, which favor groups with more group-advantageous types, occur to counter this process. A second advantage is that parameters w, s, and r correspond to quantities readily measured in empirical multilevel systems (Table 2). These correspondences follow directly from the definitions of the parameters in the stochastic model. Thus, the formulation here could also provide intuition for specific empirical systems where the parameter values are known or can be estimated. The parameters s and r are also readily reformulated to correspond to a public goods game, where individual-level fitness is frequency-dependent. In this case, the selection coefficient at the individual level is a function of the fraction ðk=nÞ of type G (producer) individuals in a group (Appendix B). While the assumption that groups replicate and that offspring of groups are identical to the parent may seem restrictive, it is a natural and convenient approximation which holds for many of the above examples. In host–pathogen systems, if an offspring infection is seeded by a random sample of pathogens in the parent host (e.g. through transmission or cell division), then the initial pathogen composition in the offspring host will, on average, be the same as in the parent host. This similarity increases with sample size by the Central Limit Theorem (Durrett, 2010). The same applies to systems, such as stands of plants, in which new stands are created from propagules of existing stands (Stevens et al., 1995; Weinig et al., 2007). Group-level replication can also account for higher level selection in multilevel systems where groups do not replicate per se. In microbial systems where populations are not serially diluted, subpopulations with a higher proportion of producers of public goods grow faster than those with less, even though non-producers always outcompete producers in direct competition (Turner and Chao, 1999; Chuang et al., 2009). This is phenomenologically captured in the framework, where groups with a higher fraction of producers (type G) have a faster replication rate, thereby increasing the representation of that producer/non-producer ratio in the overall population. Together with the group maturation process described above, this characterization of group replication captures the key group evolutionary dynamics one expects in a multilevel system: group offspring inherit their parent group's composition, but are also subject to individual-level selection during their existence.

3. Results 3.1. A single parameter in the large system limit In the limit m-1, the ball-and-urn process simplifies to a histogram changing in time as described by a system of ordinary differential equations (ODEs; Fig. 1e; Appendix A.4). Taking the limit n-1 (in addition to m-1) results in the following partial

differential equation (Fig. 1e; Appendix A.5): " # Z 1 ∂ ∂ μ ðxÞ ¼ ½xð1 xÞμt ðxÞ þ λμt ðxÞ x  yμt ðyÞ dy ∂t t ∂x 0

ð1Þ

where

λ¼

wr s

ð2Þ

For a fixed time t and x A ½0; 1, μt ðxÞ is the density of groups whose frequency of type G individuals is x at time t. Eq. (1) succinctly illustrates the dynamics of multilevel selection which have previously been described verbally (Kokko and Heubel, 2011) and observed empirically (Chuang et al., 2009). The first term on the right of Eq. (1) represents the continuous flow under individual-level selection for individual-advantageous types. The second term represents the counteracting force of group-level selection that allows groups with more type G individuals to replicate faster. The balance of these two forces depends on only one parameter: λ. In the case where groups are initially uniformly distributed across all type G frequencies, pure type I groups will dominate when λ o 2, whereas a coexistence of type I and G individuals occurs for λ 4 2 (Appendix A.5). The form of this equation, with individual-level selection represented by flux and group-level selection by forcing, is also seen in a more general formulation of multilevel selection (Simon, 2010). This suggests that though the formulation here contains specific assumptions, its depiction of the nature of the competing selection forces is consistent with a more general model. 3.2. Relative replication rates matter The form of λ in (2) reveals that w, the relative rate of group events to individual events, has the same effect on evolutionary dynamics as r, the group-level selection coefficient. Returning to the virus example, suppose there are two host–virus systems which are exactly the same except that the time to transmission to a new host is shorter in the first system than the second system. Then, under the above framework, the virus is expected to evolve to be more transmissible in the first system, i.e. where transmission events occur more frequently, than in the second system. The role of w indicated by the formula for λ is also consistent with a study of multicellularity in yeast (Ratcliff et al., 2012). In this study, yeast were cultured in tubes and yeast clusters that had settled to the bottom of tubes were transferred to seed the next generation. When researchers shortened the time they waited for the yeast to settle, i.e. increased w, they observed an increase in the trait favored by group-level selection: cluster size achieved before propagule detachment. This relationship between w and grouplevel selection also holds for the stochastic ball-and-urn process, where the time to fixation of type G individuals decreases with w (Fig. 2a). Interestingly, a higher value of w does not increase the fixation probability of type G individuals (Fig. S1).

S. Luo / Journal of Theoretical Biology 341 (2014) 41–52

3.3. Mutation benefits the individual-advantageous type A similar large-system derivation yields a partial differential equation for a case with mutation, in which each individual mutates to the other type at rate v Eq. (A.16). In this system, type I individuals dominate in parameter regimes for which, in the absence of mutation, type G individuals would otherwise dominate (Figs. 2b, S2a and S3). This mutation effect is a feature of selection acting at two levels. In a pure type I group, any type G mutant that arises will be quickly driven to extinction by individual-level selection. The group-level advantage gained through this mutation is minor and unlikely to counteract this loss. On the other hand, in a pure type G group, any type I mutant that arises will be favored by individual-level selection to dominate the group. Meanwhile, the decrease in the group-level advantage from this mutation is minor and unlikely to eliminate the group. These results show that mutation should be low to promote selection at higher levels, which has previously been suggested for the evolution of unselfish pathogens (Bonhoeffer and Nowak, 1994). Note that this result assumes the existence of the two types of individuals, I and G. The question of how mutation generates this initial diversity in addition to the role it plays in subsequent antagonistic multilevel selection is more complex. In particular, for a situation where mutation–selection balance generates a continuum of phenotypic types that subsequently evolve to be dominated by group-advantageous types (Wilson and Dugatkin, 1997), mutation may have a more nuanced role, one where a higher mutation rate is advantageous for type I individuals as above, but is possibly also disadvantageous in that it generates the diversity of phenotypes that leads to selection for group-advantageous types. 3.4. Population size matters How do these findings change if the number of groups m and the number of individuals per group n are finite? When n is finite, type G individuals dominate over a wider range of parameters compared to when n is infinite (Figs. 2c and S4; Fig. S5 illustrates a similar result for a public goods game). This is because n determines the amount of individual-level stochasticity (i.e. genetic drift) in the multilevel system. With a small number of individuals per group, type G individuals can dominate a group after a few chance replacements of type I individuals with type G individuals.

45

In the extreme case where the number of individuals per group n is two, type G individuals almost always win (Appendix A.5). Individual-level selection is therefore weaker – group-level selection stronger – when there are fewer individuals per group. An analogous argument applies to the number of groups m: grouplevel selection is weaker – individual-level selection stronger – when there are fewer groups (Fig. 2d). These findings generalize a result obtained from a theoretical model of bacteria–plasmid associations, where intercellular (group-level) selection was found to dominate in large bacterial populations (large m) while intracellular (individual-level) selection was found to dominate in small bacterial populations (small m) (Paulsson, 2002). They also extend a result obtained for a public goods game which found that cooperation will evolve if m is large or n is small (Traulsen and Nowak, 2006).

4. Discussion and conclusion The conflict between traits that are beneficial at a local individual scale and traits that are beneficial at a larger group scale is pervasive throughout biology. I have presented here a unifying mathematical formulation for this phenomenon that enhances our understanding of antagonistic multilevel selection beyond what would be achieved from a system-specific model. This formulation allows the systematic derivation of general ‘rules of thumb’ which have been either difficult to obtain analytically or were subject to stricter assumptions in previous general approaches. I find that a faster rate of group-level events enhances group-level selection, a higher mutation rate enhances individual-level selection, and a larger number of groups enhances group-level selection while a larger number of individuals per group enhances individual-level selection. In addition, I show how multilevel selection can be conveniently visualized as a ball-and-urn process and derive a partial differential equation approximation to this process which concisely illustrates the nature of the counteracting selection forces. Importantly, these results depend minimally on model parameterizations and are therefore stronger than simulation-based or computational-based approaches. The specific assumptions underlying these results capture many properties that are common to examples of multilevel selection. First, group maturation between group-level events is captured through w, the relative rate of group-level and

Fig. 2. General properties of the ball-and-urn system. (a) The time to fixation (T fix ) of type G individuals (red circles) decreases as the relative rate of group replication to individual replication (w) is increased. (b) Upper panel: In the absence of mutation, a population that is initially uniformly distributed over ½0; 1 in a system with λ ¼ 3 will evolve to be dominated by type G individuals. Lower panel: With a moderate level of mutation (v=s ¼ 0:01), the same system will evolve to have a majority of type I individuals. Solid white lines denote the overall fraction of type G individuals in the population over time. (c) Long-term fraction of type G individuals in total population with n finite (n¼ 50) over range of parameter values. White dashed line is the threshold for dominance by type G individuals when n is infinite (λ 42). Heat map generated using numerical solutions of approximating system of ODEs run to time t¼ 50. (d) The fixation probability of type G individuals (red circles) increases with the number of groups (m), asymptoting to the corresponding value for the system of ODEs (i.e. m ¼ 1). In (a) and (d), each point is the mean of 1000 replicate simulations starting from the uniform distribution on f0; 1=n; …; 1g with parameter values (unless otherwise stated) n¼49, m¼ 50, r¼ 0.2, s ¼0.2, w¼ 1, v ¼0.

46

S. Luo / Journal of Theoretical Biology 341 (2014) 41–52

individual-level replication. Second, the assumptions of constant group size n and constant number of groups m allow the model to focus on the changes in frequency of the two types in the population without having to explicitly keep track of group size or population size dynamics. These assumptions are also less restrictive when interpreted as effective population sizes rather than physical population sizes (Durrett, 2008). This broader interpretation allows systems where group sizes or population sizes change to be modeled with constant effective group size and effective number of groups. Third, the assumption that group fitness is linear in the fraction of type G individuals is a minimalist way of capturing antagonistic forces and, as mentioned earlier, corresponds to a public goods game, which constitutes a large class of multilevel systems. Last, the fact that group fitness depends only on group composition is reasonable for most of the examples in Table 1 where the individuals, typically lowerorder organisms, cannot directly influence how they are replicated at the group level. With the results established here as a basis of comparison and as a guide to intuition, one may then explore the behavior of more complex systems by relaxing the assumptions in the model. Some examples of how the above assumptions can be relaxed are as follows. Nonconstant group size can be explicitly modeled by doubly indexed urns: one index for the number of type G individuals and one for the number of type I individuals. Group fitness functions that depend more generally on the fraction of type G individuals are also readily incorporated by replacing ð1 þ rj=nÞ in the formulation with f ðj=nÞ. Similarly, as with the public goods production case, other individual fitness functions can be substituted for ð1 þ sÞ. For systems where group fitness cannot be expressed as a function of the proportion of type G individuals, one can augment the ball-andurn formulation to include additional parameters or variables that determine the rules for group replication and selection. It should be noted that these relaxed assumptions could lead to less tractable formulations and numerical- and simulation-based approaches may be need to obtain results. To close, it should be acknowledged that along with multilevel selection theory, inclusive fitness theory, the idea that an individual's fitness includes both direct reproductive fitness and the indirect fitness derived from benefitting related individuals (West et al., 2007), has also been applied to examples in Table 1. Whether evolution proceeds via multilevel selection or via selection on inclusive fitness continues to be debated (West et al., 2007; Wilson and Wilson, 2007; Kerr, 2009; Okasha, 2010; Nowak et al., 2010; Abbot et al., 2011). In the broader context of this debate, the ball-and-urn process presented here has value as a standard formulation for multilevel selection. Much like the Wright–Fisher and Moran models which are standard models for evolution at one biological scale, its properties are well-characterized and mathematically concrete. By highlighting the fundamental and common properties that may be expected in multilevel selection, it can serve as a basis of comparison for empirical observations, for models that are specific to particular biological examples, or for models based on an inclusive fitness approach. With its transparent mathematical construction and with parameters that correspond to empirically measurable quantities, the ball-and-urn formulation can be an illuminating and natural approach for understanding multilevel selection.

Acknowledgments Thanks to Katia Koelle, Rick Durrett, and Anand Pai for feedback on the manuscript, and Jonathan Mattingly, Mike Reed, and the Koelle Lab for helpful discussions. This work was supported by grant NSF-EF-08-27416 and the Duke Stern Fellowship.

Appendix A. Basic mathematical framework A.1. Mathematical description The basic mathematical framework is a nested Moran process (Durrett, 2008) with two phenotypes for the individual-level Moran process and n þ 1 phenotypes for the group-level Moran process. Consider a population of m groups, each containing n individuals:


each type I individual replicates at rate 1 þ s, s Z 0, and each type G individual at rate 1


the number of individuals in each group is kept constant at n by





selecting an individual from the group uniformly at random to die whenever an individual replicates in that group groups replicate at a rate linear in the proportion ðk=nÞ of individuals they contain that are type G:   k w 1þr n where w 4 0 controls the rate of group replication relative to individual replication and r Z 0 the total number of groups is kept constant at m by selecting a group uniformly at random to die whenever a group replicates

The above can be represented as the vector-valued stochastic j process X t ¼ ðX 1t ; …; X m t Þ, where Xt denotes the fraction of type G j individuals in group j at time t, X t A f0; 1=n; …; 1g. For example, in Fig. 1a, group 1 has three out of three type G individuals, group 2 has zero out of three, and group 3 has two out of three, which corresponds to X t ¼ ð1; 0; 23 Þ. The reformulation of this stochastic process as a ball and urn system (Fig. 1d) is obtained by considering the empirical distribution function:

μm;n t ðÞ ¼

1 m ∑ δ j ðÞ m j ¼ 1 Xt

where δx ðyÞ ¼ 1 if x¼y and 0 otherwise. In other words, μm;n t counts the number of groups that contain a given frequency of type G individuals and divides it by m. E.g. μm;n t ð1=nÞ is the fraction of groups which have 1=n individuals that are type G. I will refer to a group as being class i=n if it contains i=n type G individuals. Continuing the earlier example, X t ¼ ð1; 0; 23Þ corresponds to m;n 2 m;n m;n 1 1 1 μm;n ð0Þ ¼ 13, μm;n is a stochast ð3Þ ¼ 0, μt ð3Þ ¼ 3, and μt ð1Þ ¼ 3. μt tic process and, because its realizations are probability measures, n;m ∑m k ¼ 1 μt ðk=nÞ ¼ 1, it is a probability measure-valued process. To lead up to the mathematical description of the dynamics of this process, fix i A f0; …; ng. μn;m t ði=nÞ will decrease by 1=m if one of the following events occur:

1. In a group of class i=n, a type I individual replicates and a type G individual is chosen to die (a ball in urn i moves one urn to the left). There are mμm;n t ði=nÞ such groups. In each of these, n  i type I individuals replicate at rate 1 þ s. A type G individual will be chosen to die with probability i=n. 2. In a group of class i=n, a type G individual replicates and a type I individual is chosen to die (a ball in urn i moves one urn to the right). There are mμm;n t ði=nÞ such groups. In each of these, i type G individuals replicate at rate 1. A type I individual will be chosen to die with probability 1  i=n. 3. A group of class j=n replicates and a group of class i=n (i aj) is chosen to die (a ball in urn i moves to urn j). There are mμm;n t ðj=nÞ groups of class j=n each replicating at rate wð1 þrj=nÞ. A class i=n group will be chosen to die with probability μm;n t ði=nÞ.

S. Luo / Journal of Theoretical Biology 341 (2014) 41–52

Similarly, μn;m t ði=nÞ will increase by 1=m if one of the following events occur: 4. In a group of class ði þ 1Þ=n, a type I individual replicates and a type G individual is chosen to die. 5. In a group of class ði 1Þ=n, a type G individual replicates and a type I individual is chosen to die. 6. A group of class i=n replicates and a group of class j=n ði a jÞ is chosen to die. The transition rates of the process μm;n are given by R ¼ R1 þ wR2 , t where   8 i i > > ðn  iÞð1 þ sÞ if j ¼ i 1; i o n > mu > n n   > <    1 i i R1 u; u þ ðej  ei Þ ¼ > mu i 1 if j ¼ i þ1; i 4 0 m > > n n > > : 0 otherwise describes individual-level events, and      8 i j j   < u 1þr mu 1 n n n R2 u; u þ ðej  ei Þ ¼ : m 0

Gi is the rate that a ball, initially in urn I1, goes to urn i. I2 is the urn that the ball drawn from I1 will be going to 5. Update u and T:     I1 I1 1 ¼u  u m n n     I2 I2 1 ¼u þ u m n n T ¼ T þτ

Note that this implementation is more efficient than a standard Gillespie (2007) algorithm. In the latter algorithm, the number of rate calculations required is on the order of n2 (a ball moves from urn i to urn j with 0 r i; j r n) for each event. In comparison, for the implementation here, the number required is only on the order of 2n (steps 1 and 4) for each event. Repeated iterations of this algorithm were used to generate Figs. 2a, d and Figs. S1, S3 and S4.

A.3. Taking the limit m-1 iaj otherwise

describes group-level events. ei is the vector of length n þ 1 with a one in the ith position and zeros elsewhere and the term ej  ei corresponds to a ball moving from urn i to urn j. A.2. Stochastic implementation These rules by themselves can be used to stochastically simulate the process using an exact method (i.e. one that mimics the mathematical description precisely). Note that without mutation and with w 4 0, the process has two absorbing states: all type I individuals ðμm;n ð0Þ ¼ 1Þ or all type G individuals ðμm;n ð1Þ ¼ 1Þ. I programmed this process in Cþ þ using the following exact stochastic simulation algorithm:

(initial distribution of balls)
Input: m, n, r, s, w, u ¼ μm;n 0
Initialize time: T¼ 0.
While uð0Þ o 1 or uð1Þ o 1 (while fixation has not occurred): 1. Calculate rate that balls leave each urn: 8        i i i j j > > ðn  iÞ ð2 þ sÞþ mu ∑ mu 1þr mu > > < n n n jai n n    Li ¼   i j j > > > ∑ mu 1þr > : n jai n n

0 o io n i ¼ 0; n

2. Draw τ  PoissonðλÞ where λ ¼ ∑ni¼ 0 LðiÞ to determine the time to the next event. 3. Draw a number, I1, randomly from f0; 1; …; ng where the probability of drawing i is Li =∑nj¼ 0 Lj . I1 is the urn that a ball will be leaving from 4. Draw a number, I2, from f0; 1; …; ng where the probability of drawing i is Gi =∑nj¼ 0 Gj , given by 8    I1 I1  1 I1  1 > > ðn  I 1 Þð1þ sÞ þ u 1þr > > > n n n > >     > <  I1 I1 þ 1 I1 þ 1 þu 1þr Gi ¼ I 1 1  n n n > >  >   > > i i > > 1þ r >u : n n

47

if i ¼ I1  1; I 1 o n if i ¼ I1 þ 1; I 1 4 0 if ji  I 1 j4 1

To derive the differential equations that approximate this system, I compute the infinitesimal mean and variance for this process in a manner akin to diffusion approximations for standard Moran processes (Karlin and Taylor, 1981; Durrett, 2008). To calculate the infinitesimal mean of μm;n for 0 o io n, note t that      i m;n i  μ E μm;n t t þ Δt n n     1 i 1 ¼  P μm;n decreases by t m n m     1 i 1 þ P μm;n increases by t m n m     1 i i ¼ mμm;n i 1  ð2 þ sÞ t m n n     ) j m;n i m;n j ∑μ 1þr þwmμt Δt n jai t n n      1 iþ1 iþ 1 mμm;n ðn  ði þ 1ÞÞð1 þ sÞ þ t m n n     i 1 i1 ði  1Þ 1  þmμm;n t n n      i i i þwmμm;n 1þr 1  μm;n Δt þ oðΔtÞ ðA:1Þ t t n n n where each term in (A.1) corresponds (in order) to the events 1–6 described in Section A.1 (the first two events have been combined in the first term). I have used the fact that the probability of an event occurring in a small time interval, Δt, is equal to the rate that the event occurs multiplied by Δt. The oðΔtÞ term represents compound events (more than one event occurring in the interval Δt), which have a lower probability of occurring and can be ignored in the limit as Δt-0. Thus, after some rearranging of (A.1),      1 i m;n i  E μm;n μ lim t t þ Δt n n Δt-0 Δt           1 i i i i i þ m;n i ¼ D2 μm;n 1  þ sD 1  μ t t 1 n n n n n n n #  "   n i i j j  ∑ μm;n þwr μm;n ðA:2Þ t n n j¼0 t n n

48

S. Luo / Journal of Theoretical Biology 341 (2014) 41–52

where



   f j ¼ D1þ f n

   jþ1 j f n n 1=n

þ sD1þ jon

is a first-order difference quotient and       jþ 1 j j1    f  2f þf j n n n D2 f ¼ n 1=n2

0 o jo n

1

E½μm;n ð0Þ  μm;n t ð0Þ t þ Δt      n 1 1 j m;n j 1  ð1 þ sÞ  wr μm;n ¼ μm;n t t ð0Þ ∑ μt n n n n j¼0

lim

Δt

ðA:3Þ

and for i¼n, lim

Δt-0

1

Δt

E½μm;n ð1Þ  μm;n t ð1Þ t þ Δt

¼ μm;n t



n1 n



  # n 1 j m;n m;n j 1 þ wr μt ð1Þ 1  ∑ μt n n n j¼0 

"

ðA:4Þ

Next, the infinitesimal variance ð0 o io nÞ:     2 i m;n i  μ E μm;n t t þ Δt n n      1 1 m;n i decreases by ¼ þ 2 P μt n m m     i 1 increases by þ P μm;n t n m      1 i i i 1  ð2 þ sÞ ¼ 2 mμm;n t n n n m      j m;n i m;n j ∑ μt 1þr þ wmμt n jai n n     i þ 1 iþ 1 i þ 1 1 ð1 þ sÞ þ mμm;n t n n n     i1 i 1 i1 1 þ mμm;n t n n n      i i m;n i 1 þr 1  μm;n þ wmμt Δt þ oðΔtÞ t n n n

A.4. ODE approximation variance

is

zero,

the

ðA:6cÞ

where (A.6b) applies to 0 o io n. This system has ðn þ 1Þ and depends on two parameters, s and wr. I solve this system numerically using ode45 in Matlab to generate Fig. 2c. I also analytically study the linearized system for n¼ 2. The coefficient matrix in this case is 2 3  1 wr 12μ1 þ μ2  wr μ0 2½ð1 þ sÞ  wr μ0   6 7 6 7 0  12ð2 þ sÞ þ wr 12  μ2  μ1  wr μ1 4 5  1 1 0 wr 1  2μ1  2μ2 2ð1  wr μ2 Þ

Next, I take the limit n-1. Under this limit, ð1=nÞD2 -0 and D1þ -∂=∂x (where x A ½0; 1), and μt ðxÞ≔limn-1 μnt ði=nÞ. Eq. (A.6) then becomes " # Z 1 ∂ ∂ μ ðxÞ ¼ s ðμt ðxÞxð1  xÞÞ þ wrμt ðxÞ x  μt ðyÞy dy ðA:7Þ ∂t t ∂x 0

Remark 2. A mathematically rigorous proof of the weak convergence of this process in the space of probability measure-valued processes is currently in preparation with Jonathan Mattingly.

infinitesimal

" # n d 1 j μn ðtÞ ¼ μn ðtÞðn  1Þ þ wrμn ðtÞ 1  ∑ μj ðtÞ dt n n j¼0

ðA:6bÞ

A.5. Integro-partial differential equation approximation

Remark 1. Unlike standard diffusion approximations for the Moran and Wright–Fisher models (Durrett, 2008; Karlin and Taylor, 1981), I did not rescale the selection coefficients s and r.

the

      d 1 i i i i þ sD1þ μi ðt Þ 1  μi ðtÞ ¼ D2 μi ðt Þ 1  dt n n n n n " # n i j þ wr μi ðtÞ  ∑ μj ðtÞ n j¼0 n

For the state where all groups contain only type I individuals ðμ ¼ ð1; 0; 0ÞÞ, the eigenvalues are 0, wr, and 12ðwr ð2 þsÞÞ. Therefore, as long as wr 4 0 (there is selection at the group level) this state is locally unstable. For the state where all groups contain only type F individuals ðμ ¼ ð0; 0; 1ÞÞ, the eigenvalues of this matrix are  wr (repeated) and  12ð2 þ s þwrÞ which are both negative if wr 4 0. Therefore, this state is locally stable.

Notice that there is a factor of 1=m in front of each term after dividing through by m2 (this is also true for i ¼ 0; n), thus     2 1 i m;n i  E μm;n μ ¼ 0 for i ¼ 0; …; n ðA:5Þ lim t t þ Δt n n m-1;Δt-0 Δt

Since

     "   # n i i i i j n i n j 1 þwr μt  ∑ μ μ n n n n n j¼0 t n n n t

To see that this gives a system of differential equations, it is convenient to set μi ðtÞ ¼ μnt ði=nÞ. The limit as m-1 of Eqs. (A.2)– (A.4) in system is then given by " #   n d 1 j μ0 ðtÞ ¼ μ1 ðtÞ 1  ð1 þ sÞ þ wr μ0 ðtÞ  ∑ μj ðtÞ ðA:6aÞ dt n n j¼0

is the second-order difference quotient. Similarly, for i¼0, Δt-0



process

μm;n behaves like its mean, E½μm;n t t  when m-1. Informally, μnt ≔limm-1 μm;n ¼ limm-1 E½μm;n t t , where ≔ denotes a definition. Using this, take the limit m-1 of both sides of (A.2) to obtain        d n i 1 i i i ¼ D2 μnt 1 μt dt n n n n n

where x has replaced i=n and I have assumed that μt is a density so R1 that ∑nj¼ 0 μnt ðj=nÞj=n- 0 μt ðyÞy dy as n-1. This system can be nondimensionalized by dividing through by s and rescaling time. Rewrite μt ðxÞ ¼ μðt; xÞ to conform with conventional PDE notation to obtain " # Z 1 ∂ ∂ μðt; xÞ ¼ ðμðt; xÞxð1 xÞÞ þ λμðt; xÞ x  μðt; yÞy dy ðA:8Þ ∂t ∂x 0 where λ ¼ wr=s. This is the form of Eq. (1) in the main text. It can be solved numerically using a finite difference method. Analytically, I can further obtain a closed-form solution for an initial condition that is uniform. I use this solution to derive conditions for when group-level selection dominates. As shown below, if λ o 2, there is a delta mass at x¼ 0 (individual-level selection dominates), and if λ Z 2, there will be no delta mass at x ¼0, with the density at x¼ 1 increasing linearly with λ(group-level selection dominates). R1 To solve this equation, set hðtÞ ¼ 0 μðt; yÞy dy and use the method of characteristics (see, for example, Pinchover and Rubinstein, 2005) to obtain a solution in terms of h(t):   Rt xet tλ hðzÞ dz 0 e μðt; xÞ ¼ μ0 ½1 þ xðet  1Þðλ  2Þ 1 þ xðet  1Þ

S. Luo / Journal of Theoretical Biology 341 (2014) 41–52

where μ0 ðxÞ ¼ μð0; xÞ is the initial state of the system. Assuming that the system is initially at uniform distribution (μ0  1 on ½0; 1), then the solution simplifies to Rt μðt; xÞ ¼ et  λ 0 hðzÞ dz ½1 þ xðet  1Þðλ  2Þ ðA:9Þ To find h(t), use the fact that μ is a probability density (it has total mass of 1): Z 1 Rt t λ hðzÞ dz 0 ½1 þ xðet  1Þðλ  2Þ dx ðA:10Þ 1¼e 0

There are then two cases to consider, λ a1 and λ ¼ 1. Case λ a 1. The integral on the right hand side of (A.10) can be calculated explicitly, yielding Rt tðλ  1Þ 1 t λ hðzÞ dz e 0 1¼e ðet 1Þðλ  1Þ Rt ðet  1Þðλ  1Þ t λ hðzÞ dz 0 ðA:11Þ ¼ tðλ  1Þ )e e 1 Plugging (A.11) into (A.9), ðet  1Þðλ  1Þ ½1 þ xðet  1Þðλ  2Þ etðλ  1Þ  1

μðt; xÞ ¼

ðA:12Þ

I use this closed-form expression for μðt; xÞ to plot Fig. 2b (upper panel) and Fig. S6 (upper panel). To investigate the dependence of long-term behavior on λ, consider the limits as t-1 at the endpoints. The density at x¼ 0 has steady-state: 8 if λ 4 2 >0 ðet  1Þðλ 1Þ < - 1 μðt; 0Þ ¼ tðλ  1Þ if λ ¼ 2 > e 1 : 1 if λ o 2

¼

ðet  1Þðλ 1Þ tðλ  2Þ e etðλ  1Þ 1

ðetðλ  1Þ  etðλ  2Þ Þðλ  1Þ etðλ  1Þ  1

(

λ 1 if λ 4 1 if λ o 1

0

Note: By L'Hopital's rule, μðt; xÞ-0 as t-1 if λ o1 for any x a 0. Case λ ¼ 1. Rt t t hðzÞ;dz 0 1¼e et  1 Rt et  1 t hðzÞ dz 0 ðA:13Þ ¼ )e t Plugging (A.13) into (A.9):

μðt; xÞ ¼

et  1 ½1 þ xðet  1Þ  1 t

A.6. Incorporating mutation To incorporate mutation, I proceed as above but assume that each individual mutates at a constant rate v Z 0 to the opposite type. Analytical results can no longer be obtained, but the stochastic process can still be simulated and the approximating differential equations studied numerically. With mutation, in addition to the previous individual-level moves, μm;n t ði=nÞ can also change via an individual mutating from type I to type G (at rate iv per group in class i=n), and from type G to type I (at rate ðn iÞv per group in class i=n). Three addition terms are therefore appended to the infinitesimal mean in (A.2):      1 i i  μm;n E μm;n lim t Δ t t þ n n Δt-0 Δt           1 i i i i i i 1 þ sD1þ μm;n 1 ¼ D2 μm;n t t n n n n n n n  "   # n i i j j  ∑ μm;n þwr μm;n t n n j¼0 t n n     i þ1 i1 ðiþ 1Þv þ μm;n ðn  ði 1ÞÞv þ μm;n t t n n   i ðiv þ ðn  iÞvÞ  μm;n t n           1 i i i i i i 1 þ sD1þ μm;n 1 ¼ D2 μm;n t t n n n n n n n  "   # n i i j j  ∑ μm;n þwr μm;n t n n j¼0 t n n          i þ1 iþ 1 i  m;n i  D 1  μ þv D1 μm;n t t 1 n n n n where, analogous to the definition of D1þ ðf ðj=nÞÞ,     j j1    f f j n n j40 ¼ D1 f 1=n n

and the density at x¼ 1 has steady-state:

μðt; 1Þ ¼

49

ðA:14Þ

Doing the same analysis at the endpoints x ¼0 and x ¼1, I obtain et 1 -1 t 1et -0 μðt; 1Þ ¼ t

μðt; 0Þ ¼

which is consistent with the results for λ a1. These results indicate that when λ o 2, there is a positive Dirac delta mass at x ¼0, i.e. a (possibly large) positive fraction of pure type I groups (such as in Fig. S6). When λ Z2, the steady-state distribution consists of groups which comprise a mixture of types G and I (such as in Fig. 2b (upper panel)), with the density of pure G groups linear in λ. This dependence of long-term population composition on λ appears to hold under different smooth initial conditions (Fig. S7).

The infinitesimal variance will again be 0 as m-1. The resulting system of ODEs is thus " #   n d 1 j μ0 ðtÞ ¼ μ1 ðtÞ 1  ð1 þ sÞ þ wrμ0 ðtÞ  ∑ μj ðtÞ dt n n j¼0 þ vðμ1 ðtÞ  nμ0 ðtÞÞ 





ðA:15aÞ 

  i i μi ðtÞ 1  n n

d 1 i i þ sDjþ μ ðtÞ ¼ Dii μi ðtÞ 1  dt i n n n " # n i j þ wrμi ðtÞ  ∑ μj ðtÞ n j¼0 n       iþ1 i  D1 μi ðtÞ 1  þ v D1 μi þ 1 ðtÞ n n

ðA:15bÞ

" # n d 1 j μn ðtÞ ¼ μn ðtÞðn  1Þ þ wrμn ðtÞ 1  ∑ μj ðtÞ þ vðμn  1 ðtÞ  nμn ðtÞÞ dt n n j¼0 ðA:15cÞ where (A.15b) applies to 0 oi o n and again μi ðtÞ ¼ μ limm-1 μm;n t . The limiting nondimensional integro-partial differential equation is " # Z 1 ∂ ∂ ∂ μ ¼ ðμxð1  xÞÞ þ η ðμð2x  1ÞÞ þ λμ x  μy dy ðA:16Þ ∂t ∂x ∂x 0 n t ði=nÞ≔

where η ¼ v=s, λ ¼ wr=s as before, and the second term comes from the fact that         iþ1 iþ1 i i D1 f 1 D1 f n n n n

50

S. Luo / Journal of Theoretical Biology 341 (2014) 41–52

           iþ 1 iþ 1 i i i þ D1 f  D1 f ¼ D1 f n n n n n ∂ ∂ -2 ðf ðxÞxÞ  ðf ðxÞÞ ∂x ∂x

as n-1

This integro-partial differential equation, solved numerically using a finite difference method, is used to generate Fig. 2b (lower panel) and Fig. S6 (lower panel). Mutation does not seem to change the ultimate outcome in a system that is already favorable to type I individuals (Fig. S6). Indeed, this asymmetric effect of mutation on the long-term fraction of type G individuals is observed for numerical solutions over a range of mutation ðηÞ selection force ðλÞ values (Fig. S2a).

Appendix B. Public goods game framework B.1. Mathematical description In a public goods game, ‘cooperators’ manufacture a public (shared) good which conveys overall benefit b to the population but costs c to produce. ‘Defectors’ do not produce public goods and therefore receive benefit without incurring a cost. In a group of size n with i cooperators, the average payoff to a cooperator is P C ðiÞ ¼

bi c n

bi n

The reproductive fitness of an individual is assumed to be related to this payoff. There are multiple functions for relating fitness to payoff. One is to assume that only the absolute difference in payoff matters for relative reproductive success. In this case, cooperators are at a constant disadvantage of c to defectors: s p c. The basic model (Section A.1) with fitnesses 1 and 1 þ s fits this scenario. Another function, which is used by Traulsen and Nowak (2006), assumes that reproductive success is linear in payoff:   bi SC ðiÞ ¼ 1 þ α c n SD ðiÞ ¼ 1 þ α

i n

 Q1

    8 i i i > > ðn iÞ 1 þ α > mv n n n >  > <      1 i i c i v; v þ ðej ei Þ ¼ > mv i 1þα  1 m > > n n b n > > : 0

if j ¼ i 1; i o n if j ¼ i þ1; i 4 0 otherwise

describes individual-level events, and      8

i j c j   < u 1þβ 1 iaj mv 1 n n b n Q 2 v; v þ ðej  ei Þ ¼ : m 0 otherwise describes group-level events. B.2. Stochastic implementation

bi n

where S stands for reproductive success and α is some constant which is small enough so that SC is nonnegative and therefore makes sense as reproductive fitness. Without loss of generality, and to reduce the number of parameters, I factor out the b and absorb it into the α term, giving   i c ðB:1Þ SC ðiÞ ¼ 1 þ α  n b SD ðiÞ ¼ 1 þ α

1. In a group of class i=n, a type I individual replicates and a type G individual is chosen to die (a ball in urn i moves to the left). There are mνm;n t ði=nÞ such groups. In each of these, n  i type I individuals replicate at rate SD ðiÞ ¼ 1 þ αi=n. A type G group will be chosen to die with probability ni . 2 In a group of class i=n, a type G individual replicates and a type I individual is chosen to die (a ball in urn i moves to the right). i There are mνm;n t ðnÞ such groups. In each of these, i type G individuals replicate at rate SC ðiÞ ¼ 1 þ αði=n  c=bÞ. A type I individual will be chosen to die with probability 1  i=n. 3. A group of class j=n replicates and a group of class i=n (i aj) is chosen to die (a ball in urn i moves to urn j). There are mνm;n t ðj=nÞ groups of class j=n each replicating at rate wð1 þ β ð1  c=bÞj=nÞ. A class i=n group will be chosen to die with probability νm;n t ði=nÞ. The events for a decrease of 1=m are analogous. Again, these verbal descriptions can be summarized mathematically with the transition rates of this process, Q ¼ Q 1 þ wQ 2 , where

and the average payoff to a defector is P D ðiÞ ¼

To keep notation consistent with the previous section, I will henceforth refer to cooperators as type G individuals and to defectors as type I individuals. Let νm;n denote the public goods t game version of the ball-and-urn process. w remains as defined in Section A.1 as the relative rate of group-level events to individuallevel events. The same type of events as in Section A.1 will change νm;n t . The only difference now is that the rates for individual level events are frequency-dependent and the parameter names are different (differences from previous model are boxed). νm;n t ði=nÞ will decrease by 1=m if one of the following events occur:

ðB:2Þ

Then 0 o SC o 1 þ αð1 c=bÞ and 1 o SD o 1 þ α. To reflect a public goods game, the basic framework (Section A.1) is modified by assuming that replication is frequency-dependent and occurs at rate SC(i) for cooperators and SD(i) for defectors. At the group level, I assume that groups replicate at a rate proportional to the average payoff per individual: 1 þ β ðb  cÞi=n, where β is some constant. Again, I absorb the factor b into the constant β to get 1 þ β ð1  c=bÞi=n. This is the same as the basic framework in Section A.1, but with r ¼ βð1 c=bÞ. Note that α and β measure the extent to which payoffs affect reproduction rates.

This process can be stochastically simulated in exactly the same way as with the basic framework (Section A.2) but with the rates modified to reflect frequency-dependent individual-level fitness and with r ¼ β ð1  c=bÞ. B.3. Taking m-1 In the public goods game, the infinitesimal mean for 0 oi o n is      i m;n i E νm;n  ν t t þ Δt n n      1 i i c m;n i mνt i 1 2 þ 2α  α ¼ m n n n b )    

j i j c m;n m;n ∑ν 1þβ 1 Δt þwmνt n jai t n b n       1 iþ1 i iþ1 mνm;n ðn  ðiþ 1ÞÞ 1 þ α þ t m n n n       i 1 c i 1 m;n i 1 ði  1Þ 1 þ α  1 þmνt n n b n     

c i m;n i m;n i 1þβ 1 1  νt þ wmνt Δt þ oðΔtÞ n b n n

S. Luo / Journal of Theoretical Biology 341 (2014) 41–52

stronger group- to individual-level selection (β=α large) all increase λ~ and therefore enhance group-level selection. On the other hand, a relatively higher cost of production, slower grouplevel events, and stronger individual- to group-level selection all decrease λ~ and therefore enhance individual-level selection.

Thus

     1 i i  νm;n E νm;n lim t Δ t t þ n n Δt-0 Δt       1 i i i i 1  1 þ ¼ D2 νm;n α t n n n n n      c  i i i þ α D1 νm;n 1 t b n n n  "   #

 n c m;n i i j m;n j  ∑ ν þ w β 1  νt b n n j¼0 t n n

B.6. Incorporating mutation ðB:3Þ

For i ¼0, 1

E½νm;n ð0Þ  νm;n t ð0Þ Δt  t þΔt    

n 1 1 1 c j m;n j 1 1 þ α  wβ 1  νm;n ¼ νm;n t t ð0Þ ∑ νt

lim

Δt-0

n

n

n

b

j¼0

n n ðB:4Þ

and for i¼n, 1

E½νm;n ð1Þ  νm;n ð1Þ Δt  t þ Δt  t   n1 1 n1 c 1  1 þ  α α ¼ νm;n t

lim

Δt-0

n

n n b "   # n c  m;n j m;n j þ wβ 1  νt ð1Þ 1  ∑ νt b n n j¼0

ðB:5Þ

The infinitesimal variance is calculated analogously and again converges to 0 in the limit m-1. B.4. ODE approximation The infinitesimal mean calculations give the following ODE approximation (based on Eqs. (B3)–(B5)):   

n d 1 1 c j 1 þ α  wβ 1  ν0 ðtÞ ∑ νj ðtÞ ðB:6aÞ ν0 ðtÞ ¼ ν1 ðtÞ 1  dt n n b n j¼0         d 1 i i i c i i 1 1þα þ α D1 νi ðtÞ 1  νi ðtÞ ¼ D2 νi ðtÞ dt n n n n b n n " #

n c i j ðB:6bÞ þwβ 1  νi ðtÞ  ∑ νj ðtÞ b n j¼0 n    d 1 n1 c 1þα α νn ðtÞ ¼ νn  1 ðtÞ 1  dt n n b " #

 n c j þ wβ 1  ν1 ðtÞ 1  ∑ νj ðtÞ b n j¼0

I incorporate mutation in the public goods game in an analogous manner to Section A.6. The ODEs are   

n d 1 1 c j 1 þ α  wβ 1  μt ð0Þ ∑ νj ðtÞ ν0 ðtÞ ¼ ν1 ðtÞ 1  dt n n b n j¼0 þ vðν1 ðtÞ  nν0 ðtÞÞ

ðB:8aÞ



       d 1 i i i c i i 1 1þα þ α D1 νi ðtÞ 1  νi ðtÞ ¼ D2 νi ðtÞ dt n n n n b n n " #   

 n c i j iþ1 þ wβ 1  νi ðtÞ  ∑ νj ðtÞ þ v D1 νi þ 1 ðtÞ b n j¼0 n n   i  D1 ðνi ðtÞ 1  ðB:8bÞ n    d 1 n1 c 1þ α νn ðtÞ ¼ νn  1 ðtÞ 1  dt n n b " #

n c j þ w 1  νn ðtÞ 1  ∑ νj ðtÞ b n j¼0 þ vðνn  1 ðtÞ  nνn ðtÞÞ The integro-partial differential equation is " # Z 1 ∂ ∂ ∂ ~ ν ¼ ðνxð1 xÞÞ þ η~ ðνð2x  1ÞÞ þ λν x  νy dy ∂t ∂x ∂x 0

ðB:8cÞ

ðB:9Þ

where η~ ¼ v=ðαc=bÞ and λ~ is as before. The impact of mutation ðη~ Þ on the long-term fraction of type G individuals in the population for different cost-benefit ratios ðc=bÞ is illustrated in Fig. S2b.

Appendix C. Supplementary data Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.jtbi.2013.09.024. References

ðB:6cÞ

This system of differential equations was numerically solved using ode45 in Matlab and used to generate Fig. S5. B.5. Integro-partial differential equation approximation As n-1, the solution of (B.6a) approaches the solution, νðt; xÞ, of the PDE " # Z 1

∂ c∂ c ν ¼ α ðνxð1  xÞÞ þwβ 1  ν x  νy dy ∂t b∂x b 0 which, non-dimensionalized, is " # Z 1 ∂ ∂ ~ x ν ¼ ðνxð1  xÞÞ þ λν νy dy ∂t ∂x 0

51

ðB:7Þ

where λ~ ¼ ðwβ ð1  c=bÞÞ=ðαc=bÞ. This is exactly the same integropartial differential equation in Section A.5, but with λ~ instead of λ. Applying the earlier result here, a low cost of production relative to benefit (c=b small), faster group-level events (w large), and

Abbot, P., Abe, J., Alcock, J., Alizon, S., Alpedrinha, J.A.C., et al., 2011. Inclusive fitness theory and eusociality. Nature 471, E1–E4. Bonhoeffer, S., Nowak, M.A., 1994. Intra-host versus inter-host selection: viral strategies of immune function impairment. Proceedings of the National Academy of Sciences 91, 8062–8066. Bouma, J.E., Lenski, R.E., 1988. Evolution of a bacteria/plasmid association. Nature 335, 351–352. Brockhurst, M.A., Buckling, A., Gardner, A., 2007. Cooperation peaks at intermediate disturbance. Current Biology 17, 761–765. Champagnat, N., Ferrière, R., Méléard, S., 2006. Unifying evolutionary dynamics: From individual stochastic processes to macroscopic models. Theoretical Population Biology 69, 297–321. Chuang, J.S., Rivoire, O., Leibler, S., 2009. Simpson's paradox in a synthetic microbial system. Science 323, 272–275. Chuang, J.S., Rivoire, O., Leibler, S., 2010. Cooperation and Hamilton's rule in a simple synthetic microbial system. Molecular Systems Biology 6. Day, T., Proulx, S.R., 2004. A general theory for the evolutionary dynamics of virulence. American Naturalist 163, E40–E63. Durrett, R., 2008. Probability Models for DNA Sequence Evolution, 2nd ed. Springer. Durrett, R., 2010. Probability: Theory and Examples, 4th ed. Cambridge University Press 124 pp. Eldakar, O.T., Dlugos, M.J., Pepper, J.W., Wilson, D.S., 2009. Population structure mediates sexual conflict in water striders. Science 326, 816. Eldakar, O.T., Wilson, D.S., Dlugos, M.J., Pepper, J.W., 2010. The role of multilevel selection in the evolution of sexual conflict in the water strider Aquarius remigis. Evolution 64, 3183–3189.

52

S. Luo / Journal of Theoretical Biology 341 (2014) 41–52

Gillespie, D.T., 2007. Stochastic simulation of chemical kinetics. Annual Review of Physical Chemistry 58, 35–55. Gore, J., Youk, H., van Oudenaarden, A., 2009. Snowdrift game dynamics and facultative cheating in yeast. Nature 459, 253–256. Griffin, A.S., West, S.A., Buckling, A., 2004. Cooperation and competition in pathogenic bacteria. Nature 430, 1024–1027. Karlin, S., Taylor, H.M., 1981. A Second Course in Stochastic Processes, 1st ed. Academic Press. Kerr, B., 2009. Theoretical and experimental approaches to the evolution of altruism and the levels of selection. In: Garland, T., Rose, M.R. (Eds.), Experimental Evolution: Concepts, Methods, and Application of Selection Experiments. University of California Press, pp. 585–630. Kerr, B., Neuhauser, C., Bohannan, B.J.M., Dean, A.M., 2006. Local migration promotes competitive restraint in a host–pathogen ‘Tragedy of the commons’. Nature 442, 75–78. Kimura, M., 1983. Diffusion model of intergroup selection, with special reference to evolution of an altruistic character. Proceedings of the National Academy of Science 80, 6317–6321. Kokko, H., Heubel, K.U., 2011. Prudent males, group adaptation, and the tragedy of the commons. Oikos 120, 641–656. Levin, B.R., Bull, J.J., 1994. Short-sighted evolution and the virulence of pathogenic microorganisms. Trends in Microbiology 2, 76–81. Levin, S., Pimentel, D., 1981. Selection of intermediate rates of increase in parasitehost systems. American Naturalist 117, 308–315. Merlo, L.M.F., Pepper, J.W., Reid, B.J., Maley, C.C., 2006. Cancer as an evolutionary and ecological process. Nature Reviews Cancer 6, 924–935. Nowak, M.A., Tarnita, C.E., Wilson, E.O., 2010. The evolution of eusociality. Nature 466, 1057–1062. Okasha, S., 2006. Evolution and the Levels of Selection. Oxford University Press, pp. 40–74. Okasha, S., 2010. Altruism researchers must cooperate. Nature 467, 653–655. Orlando, P.A., Gatenby, R.A., Giuliano, A.R., Brown, J.S., 2012. Evolutionary ecology of human papillomavirus: trade-offs, coexistence, and origins of high-risk and low-risk types. Journal of Infectious Diseases 205, 272–279. Paulsson, J., 2002. Multileveled selection on plasmid replication. Genetics 161, 1373–1384. Pinchover, Y., Rubinstein, J., 2005. An Introduction to Partial Differential Equations. Cambridge University Press. Pollitt, L.C., Mideo, N., Drew, D.R., Schneider, P., Colegrave, N., et al., 2011. Competition and the evolution of reproductive restraint in malaria parasites. American Naturalist 177, 358–367.

Price, G., 1972. Extension of covariance selection mathematics. Annals of Human Genetics 35. Rainey, P.B., Rainey, K., 2003. Evolution of cooperation and conflict in experimental bacterial populations. Nature 425, 72–74. Ratcliff, W.C., Denison, R.F., Borrello, M., Travisano, M., 2012. Experimental evolution of multicellularity. Proceedings of the National Academy of Science 109, 1595–1600. Rice, S.H., 2004. Evolutionary Theory, 1st ed. Sinauer Associates, Inc., pp. 297–321. Rumbaugh, K.P., Diggle, S.P., Watters, C.M., Ross-Gillespie, A., Griffin, A.S., et al., 2009. Quorum sensing and the social evolution of bacterial virulence. Current Biology 19, 341–345. Simon, B., 2010. A dynamical model of two-level selection. Evolutionary Ecology Research 12, 555–588. Sober, E., Wilson, D.S., 1998. Unto Others: The Evolution and Psychology of Unselfish Behavior. Harvard University Press. Stevens, L., Goodnight, C.J., Kalisz, S., 1995. Multilevel selection in natural populations of Impatiens capensis. American Naturalist 145, 513–526. Traulsen, A., Nowak, M.A., 2006. Evolution of cooperation by multilevel selection. Proceedings of the National Academy of Science 103, 10952–10955. Turner, P.E., Chao, L., 1999. Prisoner's dilemma in an RNA virus. Nature 398, 441–443. van Veelen, M., 2005. On the use of the price equation. Journal of Theoretical Biology 237, 412–426. van Veelen, M., García, J., Sabelis, M.W., Egas, M., 2012. Group selection and inclusive fitness are not equivalent; the Price equation vs. models and statistics. Journal of Theoretical Biology 299, 64–80. Velicer, G.J., Kroos, L., Lenski, R.E., 2000. Developmental cheating in the social bacterium Myxococcus xanthus. Nature 404, 598–601. Weinig, C., Johnston, J.A., Willis, C.G., Maloof, J.N., 2007. Antagonistic multilevel selection on size and architecture in variable density settings. Evolution 61, 58–67. West, S.A., Griffin, A.S., Gardner, A., 2007. Evolutionary explanations for cooperation. Current Biology 17, R661–R672. Wilson, D.S., Dugatkin, L.A., 1997. Group selection and assortative interactions. American Naturalist 149, 336–351. Wilson, D.S., Wilson, E.O., 2007. Rethinking the theoretical foundation of sociobiology. The Quarterly Review of Biology 82, 327–348.