The fate of cooperation during range expansions

Kirill S. Korolev

1 Massachusetts Institute of Technology, Department of Physics, Cambridge, MA 02139

E-mail: papers.korolev@gmail.com

## Abstract

Species expand their geographical ranges following an environmental change, long range dispersal, or a new adaptation. Range expansions not only bring an ecological change, but also affect the evolution of the expanding species. Although the dynamics of deleterious, neutral, and beneficial mutations have been extensively studied in expanding populations, the fate of alleles under frequency-dependent selection remains largely unexplored. The dynamics of cooperative alleles are particularly interesting because selection can be both frequency and density dependent resulting in a coupling between population and evolutionary dynamics. This coupling leads to an increase in the frequency of cooperators at the expansion front, and, under certain conditions, the entire front can be taken over by cooperators. Thus, a mixed population wave can split into an expansion wave of only cooperators followed by an invasion wave of defectors. After the splitting, cooperators increase in abundance by expanding into new territories faster than they are invaded by defectors. Our results not only provide an explanation for the maintenance of cooperation, but also elucidate the effect of eco-evolutionary feedback on the maintenance of genetic diversity during range expansions. When cooperators do not split away, we find that defectors can spread much faster with cooperators than they would be able to on their own or by invading cooperators. This enhanced rate of expansion in mixed waves could counterbalance the loss of genetic diversity due to the founder effect for mutations under frequency-dependent selection. Although we focus on cooperator-defector interactions, our analysis could also be relevant for other systems described by reaction-diffusion equations.

## Author Summary

Cooperation is beneficial for the species as a whole, but, at the level of an individual, defection pays off. Natural selection is then expected to favor defectors and eliminate cooperation. This prediction is in stark contrast with the abundance of cooperation at all levels of biological systems: from bacterial biofilms to ecosystems and human societies. Several explanations have been proposed to resolve this paradox, including direct reciprocity and group selection. Our work, however, builds upon an observation that natural selection on cooperators might depend both on their relative frequency in the population and on the population density. We find that this feedback between the population and evolutionary dynamics can substantially increase the frequency of cooperators at the front of an expanding population, and can even lead to a splitting of cooperators from defectors. After splitting, only cooperators colonize new territories, while defectors slowly invade them from behind. Since range expansions are very common in nature, our work provides a new explanation of the maintenance of cooperation. More generally, the phenomena we describe could be of interest in other situations when coexisting entities spread in space, be it species in ecology or diffusing and reacting molecules in chemical kinetics.

## Introduction

Cooperation between organisms has always interested evolutionary biologists [1, 2, 3]. On the one hand, cooperative interactions are widespread in living systems. Microbes cooperate to digest food, scavenge for scarce resources, or build a protective biofilm [4]; animals cooperate to hunt or avoid predation [5]; and individual cells cooperate to enable multicellular life [6, 7]. On the other hand, evolutionary theory has struggled to explain the mere existence of cooperation [1, 2, 3]. Although cooperation is beneficial to the species, it is susceptible to invasion by defectors, individuals who reap the benefits of cooperation without paying the costs. Defectors, having a higher relative fitness, are expected to take over the population leading to the demise of cooperation. The breakdown of cooperation often has devastating consequences such as the tragedy of the commons [8] or cancer [9, 10]. The break down of cooperation can also be desirable, e.g., when it destroys biofilms protecting pathogens from antibiotics and the immune system [4]. Understanding the evolution and maintenance of cooperation is therefore an important problem in economics, medicine, and biology.

Several mechanisms have been proposed that can stabilize cooperation against defection [1, 2, 11], including direct reciprocity [12, 13], group and kin selection [14, 15, 16], and spatial structure [17, 18, 19]. Most of the previous studies have focused exclusively on the changes in the relative frequencies of cooperators and defectors and neglected possible changes in the population size. The naive expectation, however, is that a reduced level of cooperation should lead to a lower average fitness of the population and, therefore, to a lower population size. Indeed, this effect of evolutionary on ecological dynamics has been observed in several experimental populations [20, 21, 22] as well as in ecological public goods games [23, 24, 25]. In the latter studies, the authors also found that, under certain conditions, cooperators are favored by natural selection at low population densities while defectors are favored at high population densities. This dependence of evolutionary dynamics on population density suggests that the low-density edges of population ranges might be conducive to the evolution of cooperation. The edge of an expanding population is of particular interest because the new territories might be colonized mostly by cooperators.

Range expansions and range shifts are common in nature [26, 27, 28, 29, 30, 31, 32, 33]. Examples include the recolonizations of temperate latitudes between glaciations, the invasion of North American forests by the Asian long-horned beetle (Anoplophora glabripennis), and the spread of the western corn rootworm (Diabrotica virgifera) in the Midwestern United States. Apart from ecological and sometimes economic impact, range expansions bring about significant changes in the genetic diversity and the evolutionary dynamics of the expanding species [34, 35, 36, 37, 38, 39, 40]. In particular, genetic diversity is typically lost because of the founder effect, and mutations appearing close to the expansion edge are very likely to reach high frequencies in the population even if they are neutral or slightly deleterious [34, 41, 42, 39, 43]. Although the dynamics of neutral, deleterious, and beneficial mutations arising at the edge of a range expansion have been extensively studied, the fate of mutations subject to frequency-dependent selection, e.g. encoding cooperative traits, has received little attention.

Here, we study how the interplay between ecological and evolutionary processes affects the evolution of cooperation during range expansions. To this end, we formulate a model that combines the effect of genetic composition on population growth and the effect of population density on the fitnesses of different genotypes. We find that cooperators are favored at the edge of an expanding population and, under certain conditions, cooperators can outrun defectors and spread into unoccupied territories, leaving defectors behind. This mechanism of maintaining cooperation might play an important role in populations that experience frequent disturbances (e.g. forest fires) followed by range expansions. More generally, we report the splitting of a mixed two-species (or two-allele) traveling wave into a wave of one species followed by an invasion wave of the other species. This wave-splitting phenomenon might also be relevant to other processes described by traveling reaction-diffusion waves, which are widely used in biology and chemical kinetics [44, 36, 45].

Our works complements two previous studies of the evolutionary dynamics of cooperation in spatial populations [25, 46]. In [25], the intricate spatial patters of coexisting cooperators and defectors were reported for ecological public goods games. In [46], the authors found that cooperation can persist in a spatial prisoner’s dilemma game because of the cyclic turnover of cooperators, defectors, and extinctions. Similar to our work, their results rely on the fact that cooperators grow faster than defectors when spatially isolated from each other.

When splitting does not occur, defectors slow down cooperators, and the rate of spreading is primarily controlled by the evolutionary dynamics at the front. For example, we find that a mixed wave of cooperators and defectors can experience long periods of acceleration, as the genetic composition at the front adjusts to the low-density conditions. We also find that defectors can spread faster in mixed waves than they would be able in isolation or by invading cooperators. This finding could be important for conservation efforts to help ecosystems shift cohesively in response to a rapid climate change or habitat deterioration.

## Models

To understand the fate of cooperation during range expansions, we need a spatial model of a mixed population of cooperators and defectors colonizing new territories. We first discuss range expansions of pure cooperators and then consider defectors invading cooperators. Although these two spreading phenomena have a similar mathematical description, they have rarely been studied together, because one is an ecological, and the other is an evolutionary process. Finally, we introduce a complete model that allows for changes in both population density and allele frequencies and explicitly includes the coupling between these two variables. For simplicity, only one-dimensional expansion waves are considered, which should be a good approximation when number fluctuations and the curvature of the wave front can be neglected. The theoretical studies of range expansions were pioneered in [47, 48], and a good introduction to this topic can be found in [44].

### Ecological dynamics

We begin by considering populations of only cooperators. Range expansions are driven by migration from colonized to new territories and by population growth. When migration (or dispersal) is short-range and isotropic, it can be approximated by a diffusion term leading to the following reaction-diffusion model of range expansions

(1) |

where is the population density at time and spatial location , is the effective diffusion constant, and is the per capita growth rate.

Since any habitat has a limited carrying capacity, the per capita growth rate, , must decline and become negative at high population densities. Populations with cooperative growth may also experience reduced or negative growth rates at low population densities because, for small , the probability of forming cooperating groups or the size of these groups is too small. Such nonmonotonic dependence of the per capita growth rate on population density is called an Allee effect and has been observed in different species ranging from budding yeast to desert bighorn sheep [49, 50, 51, 52, 53, 54, 55]. The most common model of an Allee effect assumes the following growth rate [44, 50]

(2) |

where is the carrying capacity, is the Allee threshold, i.e. the minimal density required for populations to grow, and sets the overall magnitude of the per capita growth rate. One typically distinguishes a strong Allee effect when and a weak Allee effect when . Allee effect is absent when because monotonically decreases for . Thus, equation (2) is sufficiently flexible to describe populations with and without an Allee effect.

Even if the growth rate is negative at small densities, populations can spread into unoccupied territories. At , we assume that the habitat is colonized for all , but it is empty for . After an initial transitory period, expansion waves typically move at a constant velocity and the density profile does not change in the comoving reference frame, i.e. the reference frame moving along with the expansion wave; see figure 1ab. The expansion velocity is known exactly [44, 48, 47, 56, 57] and is given by

(3) |

Coordinates in the comoving reference frame are then defined in term of as and . For , the shape of the wave profile is known exactly [44, 56, 57]:

(4) |

### Evolutionary dynamics

Similar to organisms, alleles encoding cooperator or defector phenotypes can spread in populations. Under the assumption of short-range and isotropic migration, the spreading of genetic changes can also be modeled by a reaction-diffusion equation:

(5) |

provided the alleles do not affect how organisms migrate and disperse. Here, is the frequency (fraction) of defectors, and is the relative growth rate of defectors describing the force of frequency-dependent natural selection. Note that the diffusion constant in equation (5) is the same as in equation (1) because both genetic and population spreading are due to the same migration process.

The following model of frequency-dependent selection is most commonly used because of its simplicity and because it appears naturally as a weak-selection limit of evolutionary game theory [2, 37, 58]

(6) |

where is the strength of selection, and is the preferred or equilibrium frequency of defectors. In spatially homogeneous populations, cooperators and defectors coexist at a stable fixed point , provided . The coexistence of cooperators and defectors has been observed experimentally; see, e.g., [59]. When , defectors outcompete cooperators, while cooperators prevail when . Negative and describe a bistable behavior, e.g. due to a chemical warfare, and is not considered here. In game theory, these four scenarios are known as the snowdrift, prisoner’s dilemma, harmony, and coordination games respectively [58].

To understand the maintenance of cooperation, we need to know how defectors invade cooperators; see figure 1cd. The velocity of this invasion can be calculated using the results of [48] and is given by

(7) |

In the comoving reference frame, the frequency of defectors changes from to (or from to , if ), as goes from to . The characteristic width of this transition scales as .

### Coupling between ecological and evolutionary dynamics

In the preceding discussion of population and genetic waves, we avoided the coupling between ecology and evolution either by assuming a constant genetic composition (no defectors) or neglecting the changes in population properties (e.g. carrying capacity) due to defector invasion. More general situations require a combined model with and depending on both and . Here, we consider a natural extension of equations (2) and (6):

(8) |

where the parameters of population dynamics depend on the genetic composition, and the parameters of evolutionary dynamics depend on population density.

Two comments on the functional form of our model are in order. First, equation (8) allows for the most general dynamics that reduces to the classic models of frequency-dependent selection and cooperative growth with an Allee effect. Second, and of an arbitrary functional form can be recast in the form of equation (8) by allowing to depend on and on . We expect these dependences to be small in many situations because the other terms in equation (8) describe the most important aspects of the population dynamics. The analysis presented below, however, does not depend on the assumption that is only a function of and in only a function of , provided , , and is a decreasing function of for . In the following, we will assume that and satisfy these conditions and will illustrate our results in the context of a simpler model given by equation (8). The advantage of this approach is that many calculations can be carried out explicitly in the simpler model thus allowing us to provide an intuitive interpretation of the results.

We also note that equation (8) is phenomenological in nature and is not derived from a more mechanistic description of species interactions. One the one hand, this approach allows us to present a very general analysis that is valid for a large number of populations. On the other hand, our model cannot answer more mechanistic questions, e.g., how the dynamics of are constrained in any particular population, or how an increase in the death rate affects the evolutionary dynamics.

The behavior of and has not been previously investigated; therefore, we typically assume that these two functions are constants. The dependencies of other model parameters in equation (8) on or are known qualitatively from previous experimental and theoretical studies. Note, however, that most of our results are derived for general functional forms and do not depend on any specific assumptions about the parameters.

Several experiments have confirmed the naive expectation that is a decreasing function [21, 20, 22]. Interestingly, the opposite effect of defectors has also been observed in recent experiments with budding yeast, where it was established that mixed populations of defectors and cooperators have a higher carrying capacity than pure populations of cooperators [60]. For simplicity, we assume that is a constant in most numerical solutions.

The Allee threshold is expected to increase with , and this was indeed observed both in models [23] and recent experiments [61]. In numerical solutions, we use , which is equivalent to requiring a certain density of cooperators to produce a sufficient amount of public goods for population growth.

As we show below, the dependence of the preferred frequency on population density is particularly important for the cooperator-defector dynamics during range expansions. Modeling of public goods games revealed that, under certain conditions, is a increasing function of [23, 24], which has recently been confirmed by experiments with budding yeast [61]. For simplicity, in numerical solutions we assume that changes from a high to a low value at a critical population density ; in other words, we use , where is the Heaviside step function, which equals one for positive arguments and zero for negative arguments. Note that the value of can be negative when cooperators outcompete defectors at low population densities.

A spatial model with and defined in equation (8) cannot be obtained by simply combining equations (1) and (5) because changes in population density and defector frequency due to migration are not independent. Instead, migration terms have to be added to the dynamics of cooperator and defector densities, defined as and . Upon using equation (8) to calculate the growth rates for and and adding the diffusion terms, we obtain that the spatial dynamics of mixed populations can be described by the following equation

(9) |

Here, the terms with come from the ecological dynamics controlling the population density, and the terms with come from the evolutionary dynamics controlling the relative abundances of cooperators and defectors. Note that we used the same diffusion constants for both cooperators and defectors. This is justified as long as the mutations causing defector phenotype do not affect dispersal.

A range expansion of a mixed population is shown in figure 1ef. Similar to the pure cooperator and invasion waves, mixed waves can spread with a constant velocity after an initial transient. Both cooperator and defector density profiles reach a constant shape in the comoving reference frame, but these shapes are different. In particular, the relative frequency of cooperators is much higher at the front than in the population bulk. This behavior is expected because is a decreasing function of , and population density declines at the front.

It is convenient to separate evolutionary and ecological dynamics, so we recast equation (9) in terms of population density and defector frequency:

(10) |

Note that the additional advection term appears because changes in due to migration depend on . In particular, migration increases the relative frequency of organisms in low-density regions if these organisms are more abundant in nearby high-density regions because regions with higher density send out more migrants. This effect of density gradients is particularly important at the wave front, where is changing rapidly, and there is a constant imbalance between migrants coming from the high-density colonized regions and migrants coming from low-density uncolonized regions. A more general discussion of the advection terms in populations with many different species (or alleles) can be found in [36].

### Numerical solutions and simulations

Numerical solutions of equation (9) can be easily obtained by standard methods [62]. Equations (1) and (5) did not require a separate solver because they are special cases of equation (9). We used an explicit forward-time centered-space (FTCS) finite-difference method (4-point stencil) [62]. Spatial discretization step was , and temporal discretization step was . This level of discretization was sufficient to ensure that numerical wave velocities did not differ from the analytical results in equations (3) and (7) by more than a percent. For cooperator and mixed waves, the initial population typically occupied the left % of the habitat. For genetic waves, the whole habitat was occupied by cooperators, while defectors were initially present only the left % of the habitat. At , populations were always in local equilibrium, i.e. with and . We used no-flux boundary conditions and computed solutions up to the time when the wave had spread into % of the habitat. The habitat size , i.e. the distance between the left and right boundaries, was typically equal to .

We also performed individual-based simulations to demonstrate the possibility of stochastic splitting. The simulations were done on a one-dimensional lattice of sites each with a carrying capacity . Every time step consisted of one possible migration event and one reproduction or death event at every site. During a migration event, a randomly chosen organism migrated with probability to one of the two neighboring sites. Migration was isotropic, and we imposed no-flux boundary conditions. During a reproduction or death event, the population at a site could increase by one, decrease by one, or remain unchanged. The probability of birth was given by , and the probability of death was given by . Death events were equally likely to eliminate cooperators and defectors, but birth events created either a new cooperator or a new defector with probabilities determined by the interaction matrix which was chosen to mimic equation (8). Given that a birth event occurred, a cooperator was born with probability , where is the frequency of defectors, and we used letters and as indices. Otherwise, a defector was born. time steps constituted a generation.

## Results

In the previous section, we formulated a reaction-diffusion model that includes the coupling between ecological and evolutionary dynamics of cooperators and defectors. This model, equation (10), and its special cases, equations (1) and (5), describe range expansions of cooperators with velocity , invasion of cooperators by defectors with velocity , and spreading of mixed waves of defectors and cooperators with velocity . For mixed waves, we find that the frequency of cooperators is higher at the wave front because cooperators are favored at low population densities; see figure 1ef. In this section, we show that, under certain conditions, cooperators can take over the entire wave front and split from defectors by colonizing new territories faster than defectors can invade from behind. We also investigate how the speed of mixed waves depends on the parameters of the model and show that, when evolutionary dynamics are much slower than ecological dynamics (), mixed waves can experience long periods of acceleration.

The behavior of mixed waves depends on the ratio of cooperator and invasion velocities. When , defectors invade faster than cooperators can spread into new territories; therefore, any initial condition leads to a mixed population of defectors and cooperators. This mixed population will expand into empty territories with both defectors and cooperators spreading at the same velocity . When , two qualitatively different outcomes are possible. Either cooperators and defectors can spread together in a mixed wave with velocity (figure 2abc), or a mixed wave can split into an ecological expansion wave and an evolutionary invasion wave (figure 2def). In the latter scenario, the population expansion with velocity is driven solely by cooperators followed by a slower invasion by defectors with velocity .

The analysis of equation (10) is complicated both by the coupling between the two differential equations and by the unknown dependences of the parameters on population density and defector frequency. To reduce this complexity, we first consider the effect of evolution on ecology and the effect of ecology on evolution separately and then analyze the general case.

### Effects of evolutionary on ecological dynamics

By neglecting the effect of ecology on evolution and setting and , we can immediately see that is a stationary solution of equation (10). With this solution for , the dynamics of become identical to that in equation (1) with the parameters , , and evaluated at . Therefore, the velocity and density profile of the mixed wave can be immediately obtained from equations (3) and (4), e.g., the mixed velocity is given by

(11) |

Defectors are expected to decrease the rate of spatial expansions [63]. Consistent with this expectation, equation (11) predicts that , provided that all or some of the growth parameters (, , and ) substantially decrease with . Note that, because is an attracting fixed point everywhere in space, equation (10) does not allow wave splitting when evolutionary dynamics are decoupled from population dynamics.

### Effects of ecological on evolutionary dynamics

To understand the effect of ecology on evolution, we assume that the ecological parameters do not depend on the frequency of defectors: , , and . The partial differential equation for is then independent from the dynamics of , so the velocity of the population wave and wave profile are given by equations (3) and (4) respectively. As a result, equation (10) is reduced to a single equation for the defector frequency by substituting the solution for . In the reference frame comoving with the population wave, this substitution results in an explicit dependence of the growth dynamics on through the -dependence of and . Compared to equation (5), the additional term comes from the effect of density gradients on spatial diffusion of alleles. This term is similar to advection in reaction-diffusion equations with the medium moving at a velocity . To make this analogy explicit, we need to take the advection velocity inside the partial derivative, which results in the following equation for

(12) |

where the term with appears because we changed to the comoving reference frame. From equation (12), one can see that results in both effective advection and effective growth. The effective terms are only functions of and can be computed from the known density profile . For , this computation can be carried out explicitly with the following results

(13) |

The effective growth term is always positive and peaks at the middle of the population wave front where and (by definition). The effective advection term is also positive, provided decreases with , which is expected when cooperators are favored at low population densities. As we show below and in Text S1, these two terms in some cases allow defectors to keep pace with cooperators even when .

The existence of a mixed traveling wave, where both cooperators and defectors spread at the same velocity, is equivalent to the existence of a steady state for in equation (12). When a steady state does not exist, population and genetic waves split, and shifts to negative with velocity . Since decreases with , the existence of a steady state requires that operator , obtained by linearizing the right hand side of equation (12) with respect to , has a positive eigenvalue . Indeed, if all eigenvalues of are negative, then because , while a positive eigenvalue ensures that small will increase until is sufficiently diminished so that . Thus, we look for a solution of the following equation with

(14) |

where , and we used primes to denote derivatives with respect to . Equation (14) can be transformed to a canonical form by the following change of variables

(15) |

which also insures that ; see Text S1. The result reads

(16) |

We can now use the standard theory of second order differential equations to establish conditions necessary for the existence of a solution for ; see [64, 65]. These conditions require that there must be a sufficiently large region where the potential is negative and the values of in this region must be sufficiently low. To show this, we multiply both sides of equation (16) by and integrate over , which gives

(17) |

after an integration by parts in the first term. Since the first two terms in equation (17) are positive, the third term must be negative, which in turn requires that there exists a region where because . This region has a finite width because , which follows from equations (8) and (13). Indeed, for , we find that , which is positive since . In agreement with the earlier discussion, when , the potential at is negative ensuring the existence of eigenfunction for and, therefore, the existence of mixed waves. For , we find that , which is positive as long as defectors are sufficiently disfavored at low densities, e.g. when . Thus, is a potential well, and the existence of an eigenfunction for requires that this potential well be sufficiently wide and deep.

We now interpret the effects of the three terms on the right hand side of equation (16) in the context of the depth and width of the potential well . The first term lowers for in the population bulk (), but it increases for at the front (). The second term is always negative; it deepens the potential well around and vanishes for . The third term is always positive, but the reduction in the depth of the potential well due to is lessened by (at least for some ). The transition from non-splitting to splitting behavior can then be achieved by reducing the potential well. In particular, the transition threshold can be crossed by increasing , the density at which natural selection starts to favor cooperators over defectors, as shown in figures 2 and 3a. Decreasing or has a similar effect.

To understand wave splitting better, we solve a special case of equation (14) exactly in Text S1 and find how the threshold between splitting and non-splitting behavior depends on model parameters. An exact solution can be found when , and for , but defectors are not viable for . In other words, we impose an absorbing boundary condition at such that . With these simplifications, we find that splitting occurs provided

(18) |

We draw three important conclusions from this result. First, the exact solution confirms that both splitting and non-splitting behaviors are possible depending on the parameters of the model. Second, the severity of selection against defectors required for splitting increases as approaches . Third, defectors can spread faster in mixed waves than they can invade cooperators when splitting does not occur. These three conclusions do not depend on the simplifying assumptions used to derive equation (18); see the discussion below and figures 2 and 3a. We now discuss the biological significance of wave splitting.

The possibility of wave splitting has important implications for the evolution of cooperation. When splitting is possible, cooperators outrun defectors, leading to a constant increase in the relative abundance of cooperators for as long as there are uncolonized territories. Frequent local extinctions followed by recolonizations could therefore maintain high levels of cooperation in natural populations. Equation (18) also suggests that traits that have little effect on the dynamics in well-mixed populations might be under selection during range expansions. One such trait is . In well-mixed populations, it only determines the rate of approach to the equilibrium not the equilibrium itself, while, in expanding populations, affects the invasion velocity and the conditions for wave splitting, both of which are expected to be under selection. The critical density is another examples. It has little effect in the well-mixed populations, but large allows cooperators to escape from defectors during range expansions.

The possibility of non-splitting highlights the importance of the coupling between ecology and evolution in the maintenance of genetic or species diversity during range expansions. Quite surprisingly, we find that defectors can spread in a mixed wave faster than they can invade cooperators despite the fact that defectors are disfavored or even eliminated by natural selection at the front. In other words, the community of cooperators and defectors is able to move to new territories while preserving the coexistence between the two phenotypes. Our results could, therefore, have important implications for conservation efforts to preserve genetic diversity or ecosystem integrity during potentially rapid range shifts due to climate change or habitat deterioration. For negative frequency-dependence discussed here, we found that diversity is more stable than one would naively predict from just measuring and . Other types of interactions could be less resilient and would require managed interventions to prevent splitting. In Text S1, we show that interventions increasing advection or relative growth rate at the front can achieve that goal.

### General case of eco-evolutionary feedback

In the general case when ecology and evolution affect each other, the nature of the transition from non-splitting to splitting behavior remains the same. In particular, the condition for splitting is still given by the existence of a solution of equation (14) with because, close to the splitting transition, the expansion front is populated almost exclusively by cooperators. To show this, let us consider how the dynamics at the front changes as one of the model parameters, say , is varied so that the system behavior changes from non-splitting to splitting. Increasing directly increases the abundance of cooperators at the front and increases the velocity of the mixed wave from , when , to , when splitting occurs. In addition to that, defectors lag more and more behind the front, as the splitting transition is approached. Indeed, close to the transition, is barely sufficient for defectors to keep up with cooperators, and, since decreases with , the fraction of defectors at the front must approach zero right before splitting occurs. These dynamics are illustrated in figure 3a and further discussed in Text S1. Another effect of the coupling between ecology and evolution is that is no longer given by equation (7) and one has to solve equation (10) to describe defector invasion.

We would now like to discuss the effect of initial conditions on wave splitting. One may naively think that deterministic wave splitting can be achieved simply by creating a region of pure cooperators in front of a mixed population, even if populations do not split when started from a well-mixed state. Initially, the range expansion started from such an initial condition indeed resembles splitting with cooperators transiently outrunning defectors and increasing in relative abundance; see figure 3b. Nevertheless, defectors eventually catch up because any initial condition (with and ) has a nonzero projection on the eigenfunction of the operator with the largest positive eigenvalue. The growth of this projection ensures the establishment of defector population at the front even when as shown in figure 3b. The ability of defectors to catch up, however, relies on the assumption that and vary continuously and can increase from arbitrarily small values. When the discrete nature of organisms is taken into account, one finds that wave fronts have a finite width [66]. Therefore, initial conditions can force splitting provided and the region of pure cooperators is sufficiently large compared to the width of expansion and invasion wave profiles.

Number fluctuations must also be considered for traveling waves of discrete entities. For biological systems, these fluctuations could arise due to random fluctuations of the environment or due to the randomness of births and deaths, which is known as genetic or ecological drift. Number fluctuations are largely irrelevant when because defectors can always catch up with cooperators. In contrast, when , splitting is the ultimate outcome because, given enough time, a fluctuation will create a region of pure cooperators, which is sufficiently large to prevent defectors from catching up. This region of pure cooperators will then grow with time making it exceedingly unlikely for another fluctuation to destroy it. Although splitting is inevitable, it may take a very long time (e.g. generations in figure 3c) because the deterministic dynamics discussed above create an effective activation barrier for stochastic splitting; see figure 3c and 4. The magnitude of this barrier goes to zero as the conditions for deterministic splitting are approached. Stochastic splitting is also possible for frequency-independent selection, provided expansion and invasion velocities are different; see for example [42].

In the preceding discussion, we neglected possible transitions between cooperators and defectors due to mutations or other heritable changes. This is justified on short time scales because mutation rates are typically small and even beneficial mutations struggle to survive number fluctuations (genetic drift) [67, 35]. On long time scales, mutations will change the dynamics in two ways. First, the coexistence fraction will be determined by both natural selection and the mutational pressures similar to the classic theory of two genetic variants [67, 35]. This shift in does not change the dynamics qualitatively and could be included in our theory by modifying accordingly. The second and more important difference is that the region of pure cooperators will not expand indefinitely following a splitting event because defectors will appear in the interior of the region of pure cooperation due to a mutation in addition to invading the region of pure cooperators from behind. In the limit of rare mutations, the average length of the region of pure cooperators will be determined by the balance between the time necessary to create this region after splitting and the time to the next successful defector mutation in the region of pure cooperators. The former time scales as for deterministic splitting. While the latter time scales as the inverse of the product of the mutation rate from cooperators to defectors , population size , and fixation probability , i.e. as ; see [67, 35]. For deterministic splitting, this balance yields , while, for stochastic splitting, will also be affected by the average waiting time to stochastic splitting. In the other limit of very frequent mutations, splitting would not be able to occur before additional defectors arise at the front, and the dynamics would be primarily determined by the balance between mutational transitions leading to mixed populations.

Another interesting consequence of the coupling between ecological and evolutionary dynamics is wave acceleration. We found that cooperators are favored at the wave front (see figure 1c). As the frequency of cooperators is increasing at the front, the instantaneous velocity of the expansion wave must also increase because defectors slow down expansions (see figure 2). The evolutionary change could however be very slow compared to the colonization dynamics leading to a gradual acceleration of the range expansion. Indeed, long periods of dramatic wave acceleration are possible in our model and are shown in figure 5. The acceleration of expansion fronts has been observed in a number of species and is typically explained by the evolution of shorter generation times, greater dispersal abilities, or specific adaptations to the environment in the newly colonized regions [68, 31, 69, 30]. Our analysis of cooperator and defector waves suggests that a changing rate of expansion could simply be a consequence of the slow adjustment of the genetic composition at the wave front to a new low-density optimum, which is different from that in the population bulk.

## Discussion

Understanding the link between ecology and evolution remains a major open problem [70]. This coupling is particularly important during range expansions, which are often accompanied by both demographic and evolutionary changes [31, 69, 30, 26]. We formulated a simple two-allele model that is capable of describing a wide range of frequency and density dependencies of selection and growth. Our analysis revealed the necessity of taking eco-evolutionary feedback into account in order to make accurate predictions. In particular, we found that the genetic composition of the population bulk may be a poor predictor for the genetic composition of the population front. Similarly, the measurements at the expansion front may be a poor predictor of the properties of the new population once it is fully established. These differences between the population bulk and expansion front make it also harder to predict the rates of expansions. Indeed, we found that the rates of expansions could be accelerating as the result of slow evolutionary changes at the expanding population edge and that species can spread faster in mixed waves than in isolation or by invading already established populations.

Our main result is that colonization of new territories can proceed in two qualitatively different ways. The first possibility is a single mixed wave, where both alleles (or species) move with the same velocity and their relative abundances reach a steady state in the reference frame comoving with the expansion. The second possibility is that one of the alleles (or species) outruns the other. The faster allele is solely responsible for the colonization, while the slower allele invades the faster one from behind with a smaller velocity. As a result, there is a growing region occupied exclusively by the faster allele. The effect of wave splitting could be especially dramatic for species that have markedly smaller migration rates in the population bulk compared to population front [71, 72] because the invasion by the slower allele would be significantly slower than the range expansion. The existence of secondary genetic waves could lead to unexpected population and genetic dynamics, which cannot be described by the classic models of range expansions [48, 47, 57, 56]. It would be interesting to know when the commonly observed changes behind the expansion front [69] are caused by de novo adaptation to the new environment and when they are caused by the secondary genetic waves, which could reach and alter populations at the newly colonized territories many generations after the arrival of the population wave.

In the context of cooperator-defector interactions, wave splitting and the increase in cooperation at the front could stabilize cooperation against defectors in populations that experience frequent disturbances followed by recolonizations. This mechanism of maintaining cooperation does not rely on reciprocity or multi-level selection, but is instead grounded in the density dependence of the evolutionary dynamics. More importantly, we have demonstrated that the coupling between ecological and evolutionary dynamics can have a profound effect on the fate of cooperation and should, therefore, be considered in both theoretical and experimental studies.

## Acknowledgments

We would like to acknowledge useful discussions with Jeff Gore, Ivana Cvijovic, and Manoshi Datta.

## Text S1

In this Supporting Text, we first expand our discussion of the boundary conditions for equation (16). We then outline the derivation of the exact solution, given by equation (18), for the transition from non-splitting to splitting behavior in a special case of equation (12). Finally, we show that both effective advection and effective growth can allow defectors to expand faster in mixed waves than they can do by invading cooperators.

Boundary conditions for equation (16)

Here, we show that are the proper boundary conditions for equation (16). Since the boundary conditions for are and , it is sufficient to show that and ; see equation (15). The former condition follows from the fact that and . The latter condition follows from the fact that , which we prove below.

The density profile at large decays exponentially as ; moreover, for a weak Allee effect [73]. We then substitute this scaling form in the equation (9) for , change to the comoving reference frame, and require that the wave profile does not depend on . This calculation yields

(S1) |

This expression for can now be compared to the advection velocity . For a strong Allee effect, is negative and, therefore, . For a weak Allee effect, , where the last inequality follows from ; see [73].

Derivation of equation (18)

Here, we derive equation (18). Since, for decreasing with , wave splitting is affected only by the linearization of for small , we are free to replace with any other function that has identical behavior for . It is particularly useful to approximate with a piece-wise linear function because this allows for an exact solution for everywhere in space, which captures the qualitative dependence of the solution on the parameters of the model. In contrast, eigenfunctions of only describe the behavior of close to the front and only up to a normalization factor. To this purpose, we redefine as

(S2) |

Note that this redefinition does not have a fixed point at , which is appropriate as long as we are not interested in cooperators invading defectors.

We now turn to equation (12), which can be further simplified by the change of variables from to using equation (4). The result reads

(S3) |

where primes now denote derivatives with respect to . We further assume that and obtain that

(S4) |

where we also set because vanishes at the splitting transition. We now use equation (S2) and the assumption that for and for , i.e. defectors are not viable at very low densities. The right hand side of equation (S4) is then proportional to for and to for ; therefore, further simplifications occur upon integrating both sides of equation (S4) over and defining a new variable for and for . For example, for the region where and , we obtain

(S5) |

The general solution of this equation is given by

(S6) |

Similarly, for , we find

(S7) |

Then, the solution for is given by from equation (S7) for and by from equation (S6) for and . The four constants , , , and and the matching point are determined by the boundary conditions and and the matching conditions at , where and . Upon using four of the five conditions, one can easily express the four constants in terms of because the conditions are linear with respect to , , , and . The remaining condition determines and can also be solved exactly:

(S8) |

We then obtain the condition for splitting, stated in equation (18), by requiring that is real, and . Note that, as the splitting state is approached, say, by increasing . Since corresponds to , we see that the wave of defectors lags more and more behind the wave of cooperators. As we show in the main text, this is a general result, and the front is populated mostly by the cooperators close to the splitting transition. As a consequence, the expansion velocity at the splitting transition is given by equation (3).

Increased spreading through effective advection and growth

Here, we demonstrate the effects of effective advection and growth on the ability of defectors to keep up with cooperators. When ecological dynamics are not affected by the evolutionary dynamics, the ability of defectors to follow cooperators in a mixed wave is mathematically equivalent to the ability of a Fisher wave to follow a moving boundary between good and bad conditions. In particular, the assumption that for makes this moving boundary absorbing. One can easily show that classic Fisher waves can only follow an absorbing boundary that moves with a velocity smaller than the Fisher velocity. Here, we show that even boundaries moving with higher velocities can be followed provided there is an effective advection or growth. This analysis illustrates the role of the term in equation (10) and suggests that both effective advection and growth could be useful management strategies to maintain genetic diversity during range expansions or to help species shift in response to a rapid climate change.

To avoid confusion with the earlier discussion, we introduce new notation. Let be the velocity of the absorbing barrier, and is the per capita growth rate. The Fisher velocity is then , where is the effective diffusion constant. When , the population cannot keep up with the barrier without external interventions. One possible intervention is to increase by a factor , say by supplying extra food, in a region of length behind the barrier, which is analogous to the effective growth discussed earlier. Another possibility, analogous to effective advection, is to move the population with velocity in a region of length behind the front, say by capturing and transporting individuals within the population. The analysis of both of these scenarios follows the same procedure as the derivation of (18), but is much simpler because the resulting linearized equations have constant coefficients; therefore they can be easily solved by the standard methods. Below we just summarize the results.

For effective growth, we find that barriers moving with any velocity can be followed, but the region of effective growth should be sufficiently large:

(S9) |

Note that the minimal size is finite when , and that it decreases only as for large . This suggests that there is an optimal size that minimizes , which could be interpreted as the total cost of the intervention.

For effective advection, we also find that barriers moving with any velocity can be followed provided and the region of advection is sufficiently large:

(S10) |

Similar to the previous case, minimal is finite when . More importantly, the minimal size is the smallest when and are about equal. Indeed, when , advection has a small effect on the dynamics, while, when , the point where the region of advection ends becomes an effectively absorbing boundary because organisms entering advection region are quickly moved towards the front of the wave.

### References

- Gardner A, Foster KR (2008) The evolution and ecology of cooperation–history and concepts. In: Korb J, Heinze J, editors, Ecology of Social Evolution, Springer, Berlin. pp. 1–36.
- Nowak M (2006) Evolutionary dynamics: exploring the equations of life. Belknap Press.
- Maynard Smith J (1982) Evolution and the Theory of Games. Cambridge University Press, Cambridge UK.
- West SA, Diggle SP, Buckling A, Gardner A, Griffin AS (2007) The social lives of microbes. Annu Rev Ecol Evol Syst 38: 53–77.
- Hölldobler B, Wilson E (1990) The ants. Belknap Press.
- Michod RE, Roze D (2001) Cooperation and conflict in the evolution of multicellularity. Heredity 86: 1–7.
- Sachs JL (2008) Resolving the first steps to multicellularity. Trends Evol Evol 23: 245–248.
- Hardin G (1968) The tragedy of the commons. Science 162: 1243–1248.
- Crespi B, Summers K (2005) Evolutionary biology of cancer. Trends Ecol Evol 20: 545–552.
- Merlo LMF, Pepper JW, Reid BJ, Maley CC (2006) Cancer as an evolutionary and ecological process. Nat Rev Cancer 6: 924–935.
- Nowak MA (2006) Five rules for the evolution of cooperation. Science 314: 1560–1563.
- Axelrod R, Hamilton WD (1981) The evolution of cooperation. Science 211: 1390–1396.
- Nowak MA, Sigmund K (1992) Tit for tat in heterogeneous populations. Nature 355: 250–253.
- Wilson D (1975) A theory of group selection. Proc Natl Acad Sci U S A 72: 143–146.
- Hamilton WD (1964) The genetical evolution of social behaviour. i. J Theor Biol 7: 1–16.
- Hamilton WD (1964) The genetical evolution of social behaviour. ii. J Theor Biol 7: 17–52.
- Nowak MA, May RM (1992) Evolutionary games and spatial chaos. Nature 359: 826–829.
- Roca C, Cuesta J, Sanchez A (2009) Effect of spatial structure on the evolution of cooperation. Phys Rev E Stat Nonlin Soft Matter Phys 80: 046106.
- Nadell CD, Foster KR, Xavier J (2010) Emergence of spatial structure in cell groups and the evolution of cooperation. PLoS Comput Biol 6: e1000716.
- Turner PE, Chao L (1999) Prisoner’s dilemma in an RNA virus. Nature 398: 441–443.
- Rainey PB, Rainey K (2003) Evolution of cooperation and conflict in experimental bacterial populations. Nature 425: 72–74.
- Griffin AS, West SA, Buckling A (2004) Cooperation and competition in pathogenic bacteria. Nature 430: 1024–1027.
- Hauert C, Holmes M, Doebeli M (2006) Evolutionary games and population dynamics: maintenance of cooperation in public goods games. Proc Biol Sci 273: 2565–2571.
- Hauert C, Wakano JY, Doebeli M (2008) Ecological public goods games: Cooperation and bifurcation. Theor Popul Biol 73: 257–263.
- Wakano JY, Nowak MA, Hauert C (2009) Spatial dynamics of ecological public goods. Proc Natl Acad Sci U S A 106: 7910–7914.
- Shine R, Brown GP, Phillips BL (2011) An evolutionary process that assembles phenotypes through space rather than through time. Proc Natl Acad Sci U S A 108: 5708–5711.
- Templeton A (2002) Out of Africa again and again. Nature 416: 45–51.
- Parmesan C, Ryrholm N, Stefanescu C, Hill JK, Thomas CD, et al. (1999) Poleward shifts in geographical ranges of butterfly species associated with regional warming. Nature 399: 579–583.
- Pateman RM, Hill JK, Roy DB, Fox R, Thomas CD (2012) Temperature-Dependent Alterations in Host Use Drive Rapid Range Expansion in a Butterfly. Science 336: 1028-1030.
- Gray ME, Sappington TW, Miller NJ, Moeser J, Bohn MO (2009) Adaptation and Invasiveness of Western Corn Rootworm: Intensifying Research on a Worsening Pest. Annu Rev Entomol 54: 303–321.
- Phillips BL, Brown GP, Webb JK, Shine R (2006) Invasion and the evolution of speed in toads. Nature 439: 803–803.
- Hitch AT, Leberg PL (2007) Breeding Distributions of North American Bird Species Moving North as a Result of Climate Change. Conserv Biol 21: 534–539.
- Bronnenhuber JE, Dufour BA, Higgs DM, Heath DD (2011) Dispersal strategies, secondary range expansion and invasion genetics of the nonindigenous round goby, Neogobius melanostomus, in Great Lakes tributaries. Mol Ecol 20: 1845–1859.
- Excoffier L, Foll M, Petit RJ (2009) Genetic consequences of range expansions. Annu Rev Ecol Evol Syst 40: 481–501.
- Korolev KS, Avlund M, Hallatschek O, Nelson DR (2010) Genetic demixing and evolution in linear stepping stone models. Rev Mod Phys 82: 1691–1718.
- Vlad MO, Cavalli-Sforza LL, Ross J (2004) Enhanced (hydrodynamic) transport induced by population growth in reaction-diffusion systems with application to population genetics. Proc Natl Acad Sci U S A 101: 10249–10253.
- Korolev KS, Nelson D (2011) Competition and Cooperation in One-Dimensional Stepping-Stone Models. Phys Rev Lett 107: 0881031–0881035.
- McInerny GJ, Turner JRG, Wong HY, Travis JMJ, Benton TG (2009) How range shifts induced by climate change affect neutral evolution. Proc Biol Sci 276: 1527–1534.
- Roques L, Garnier J, Hamel F, Klein EK (2012) Allee effect promotes diversity in traveling waves of colonization. Proc Natl Acad Sci U S A 109: 8828–8833.
- Arenas M, Ray N, Currat M, Excoffier L (2011) Consequences of Range Contractions and Range Shifts on Molecular Diversity. Mol Biol Evol 29: 207–218.
- Hallatschek O, Nelson D (2008) Gene surfing in expanding populations. Theor Popul Biol 73: 158–170.
- Hallatschek O, Nelson DR (2010) Life at the front of an expanding population. Evolution 64: 193–206.
- Travis JMJ, Munkemuller T, Burton OJ, Best A, Dytham C, et al. (2007) Deleterious Mutations Can Surf to High Densities on the Wave Front of an Expanding Population. Mol Biol Evol 24: 2334–2343.
- Murray JD (2003) Mathematical Biology. Springer.
- Douglas JF, Efimenko K, Fischer DA, Phelan FR, Genzer J (2007) Propagating waves of self-assembly in organosilane monolayers. Proc Natl Acad Sci U S A 104: 10324–10329.
- Sella G, Lachmann M (2000) On the Dynamic Persistence of Cooperation: How Lower Individual Fitness Induces Higher Survivability. J Theor Biol 206: 465–485.
- Fisher RA (1937) The wave of advance of advantageous genes. Ann Eugen 7: 355–369.
- Kolmogorov AN, Petrovsky N, Piscounov NS (1937) A study of the equation of diffusion with increase in the quantity of matter, and its application to a biological problem. Moscow University Bulletin of Mathematics 1: 1.
- Allee W, Bowen E (1932) Studies in animal aggregations: mass protection against colloidal silver among goldfishes. J Exp Zool 61: 185–207.
- Courchamp F, Clutton-Brock T, Grenfell B (1999) Inverse density dependence and the Allee effect. Trends Ecol Evol 14: 405–410.
- Berec L, Angulo E, Courchamp F (2007) Multiple Allee effects and population management. Trends Ecol Evol 22: 185–191.
- Dai L, Vorselen D, Korolev KS, Gore J (2012) Generic indicators for loss of resilience before a tipping point leading to population collapse. Science 336: 1175–1177.
- Courchamp F, Macdonald DW (2001) Crucial importance of pack size in the African wild dog Lycaon pictus. Animal Conservation 4: 169–174.
- Clutton Brock TH, Gaynor D, McIlrath GM, Maccoll A, Kansky R, et al. (1999) Predation, group size and mortality in a cooperative mongoose, Suricata suricatta. J Anim Ecol 68: 672–683.
- Mooring MS, Fitzpatrick TA, Nishihira TT, Reisig DD (2004) Vigilance, predation risk, and the Allee effect in desert bighorn sheep. J Wildl Manage 68: 519–532.
- Fife PC, McLeod JB (1977) The approach of solutions of nonlinear diffusion equations to travelling front solutions. Arch Ration Mech Anal 65: 335–361.
- Aronson DG, Weinberger HG (1975) Nonlinear diffusion in population genetics, combustion and nerve propagation Lectures Notes Math, volume 446. Springer, New York, 5â49 pp.
- Frey E (2010) Evolutionary game theory: Theoretical concepts and applications to microbial communities. Physica A 389: 4265–4298.
- Gore J, Youk H, van Oudenaarden A (2009) Snowdrift game dynamics and facultative cheating in yeast. Nature 459: 253–256.
- MacLean RC, Fuentes-Hernandez A, Greig D, Hurst LD, Gudelj I (2010) A Mixture of “Cheats” and “Co-Operators” Can Enable Maximal Group Benefit. PLoS Biol 8: e1000486.
- Sanchez A, Gore J (2013) Feedback between population and evolutionary dynamics determines the fate of social microbial populations Available: http://arxiv.org/abs/1301.2791. Accessed 14 February 2013.
- Press WH, Teukolsky SA, Vetterling WT, Flannery BP (2007) Numerical Recipes: The Art of Scientific Computing. Cambridge University Press, Cambridge.
- Hilker FM, Lewis MA, Seno H, Langlais M, Malchow H (2005) Pathogens can slow down or reverse invasion fronts of their hosts. Biol Invasions 7: 817–832.
- Titchmarsh E (1946) Eigenfunction expansions associated with second-order differential equations. Clarendon Press, Oxford.
- Merzbacher E (1997) Quantum Mechanics. Wiley.
- Hallatschek O, Korolev KS (2009) Fisher Waves in the Strong Noise Limit. Phys Rev Lett 103: 108103–108106.
- Gillespie J (2004) Population genetics: a concise guide. Johns Hopkins University Press.
- Thomas CD, Bodsworth EJ, Wilson RJ, Simmons AD, Davies ZG, et al. (2001) Ecological and evolutionary processes at expanding range margins. Nature 411: 577–581.
- Phillips BL, Brown GP, Greenlees M, Webb JK, Shine R (2007) Rapid expansion of the cane toad (Bufo marinus) invasion front in tropical Australia. Austral Ecology 32: 169–176.
- Kokko H, López-Sepulcre A (2007) The ecogenetic link between demography and evolution: can we bridge the gap between theory and data? Ecol Lett 10: 773–782.
- DeGiorgio M, Jakobsson M, Rosenberg N (2009) Explaining worldwide patterns of human genetic variation using a coalescent-based serial founder model of migration outward from africa. Proc Natl Acad Sci U S A 106: 16057–16062.
- Korolev KS, Müller MJI, Karahan N, Murray O Andrew W Hallatschek, Nelson DR (2012) Selective sweeps in growing microbial colonies. Phys Biol 9: 026008–026023.
- Van Saarloos W (2003) Front propagation into unstable states. Phys Rep 386: 29–222.