The Fate of Cooperation during Range Expansions

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.


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,39,[41][42][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 [36,44,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 Lc Lt~D where c(t,x) is the population density at time t and spatial location x, D is the effective diffusion constant, and G c (c) is the per capita growth rate.
Since any habitat has a limited carrying capacity, the per capita growth rate, G c (c), 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 c, 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] where K is the carrying capacity, c Ã is the Allee threshold, i.e. the minimal density required for populations to grow, and g c K 2 sets the overall magnitude of the per capita growth rate. One typically distinguishes a strong Allee effect when c Ã w0 and a weak Allee effect when {Kvc Ã v0. Allee effect is absent when c Ã v{K because G c (c) monotonically decreases for cw0. 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 t~0, we assume that the habitat is colonized for all xv0, but it is empty for xw0. 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,47,48,56,57] and is given by

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 Coordinates in the comoving reference frame (t,f) are then defined in term of (t,x) as t~t and f~x{v c t. For c Ã §{K=2, the shape of the wave profile is known exactly [44,56,57]: For c Ã v{K=2, the front has a qualitatively similar shape with c(f) declining from 1 to 0, as f goes from {? to z?.

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: Lf Lt~D provided the alleles do not affect how organisms migrate and disperse. Here, f (t,x)[(0,1) is the frequency (fraction) of defectors, and G f (f ) 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].
where g f §0 is the strength of selection, and f Ã is the preferred or equilibrium frequency of defectors. In spatially homogeneous populations, cooperators and defectors coexist at a stable fixed point f~f Ã , provided f Ã [(0,1). The coexistence of cooperators and defectors has been observed experimentally; see, e.g., [59]. When f Ã w1, defectors outcompete cooperators, while cooperators prevail when f Ã v0. Negative g f and f Ã [(0,1) 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 In the comoving reference frame, the frequency of defectors changes from 0 to f Ã (or from 0 to 1, if f Ã w1), as f goes from {? to z?. The characteristic width of this transition scales as ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi

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 G c and G f depending on both c and f . Here, we consider a natural extension of equations (2) and (6): 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, G c (c,f ) and G f (c,f ) of an arbitrary functional form can be recast in the form of equation (8) by allowing g c to depend on c and g f on f . 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 g c is only a function of f and g f in only a function of c, provided g c (c,f )w0, g f (c,f )w0, and G f is a decreasing function of f for f [(0,f Ã ). In the following, we will assume that G c and G f 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 f Ã (c) are constrained in any particular population, or how an increase in the death rate affects the evolutionary dynamics.
The behavior of g c (f ) and g f (c) 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 c or f 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 K(f ) is a decreasing function [20][21][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 K is a constant in most numerical solutions.
The Allee threshold c Ã is expected to increase with f , and this was indeed observed both in models [23] and recent experiments [61]. In numerical solutions, we use c Ã (f )~c Ã (0)=(1{f ), 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 f Ã on population density is particularly important for the cooperatordefector dynamics during range expansions. Modeling of public goods games revealed that, under certain conditions, f Ã is a increasing function of c [23,24], which has recently been confirmed by experiments with budding yeast [61]. For simplicity, in numerical solutions we assume that f Ã changes from a high to a low value at a critical population density c c; in other words, we use is the Heaviside step function, which equals one for positive arguments and zero for negative arguments. Note that the value of f (0) can be negative when cooperators outcompete defectors at low population densities.
A spatial model with G c and G f 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 c c~c (1{f ) and c d~c f . Upon using equation (8) to calculate the growth rates for c c and c d and adding the diffusion terms, we obtain that the spatial dynamics of mixed populations can be described by the following equation Here, the terms with G c come from the ecological dynamics controlling the population density, and the terms with G f 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 v m 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 f Ã is a decreasing function of c, 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: Note that the additional advection term 2D L ln(c) Lx Lf Lx appears because changes in f (t,x) due to migration depend on c(t,x). 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 c(t,x) 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) finitedifference method (4-point stencil) [62]. Spatial discretization step was 0:1, and temporal discretization step was minf10 {3 =D, 10 {2 =maxfg f maxfDf Ã D,1g,g c K 2 gg. 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 10% of the habitat. For genetic waves, the whole habitat was occupied by cooperators, while defectors were initially present only the left 10% of the habitat. At t~0, populations were always in local equilibrium, i.e. with c~K(f ) and f~f Ã (c). We used no-flux boundary conditions and computed solutions up to the time when the wave had spread into 90% of the habitat. The habitat size L, i.e. the distance between the left and right boundaries, was typically equal to 100.
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 K. 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 m 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 g c c 2 (Kzc Ã )=K 3 =(1zg c c 2 (Kzc Ã )=K 3 zg c (c 3 z cKc Ã )=K 3 ), and the probability of death was given by g c (c 3 zcKc Ã )=K 3 =(1zg c c 2 (Kzc Ã )=K 3 zg c (c 3 zcKc Ã )=K 3 ). 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 A which was chosen to mimic equation (8). Given that a birth event occurred, a cooperator was born with probability (1{f )(A cc (1{f )zA cd f )=(A cc (1{f ) 2 zA cd (1{f )f zA dc f (1{f )zA dd f 2 ), where f is the frequency of defectors, and we used letters c and d as indices. Otherwise, a defector was born. K 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 v c , invasion of cooperators by defectors with velocity v i , and spreading of mixed waves of defectors and cooperators with velocity v m . 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 (g f %g c K 2 ), mixed waves can experience long periods of acceleration.
The behavior of mixed waves depends on the ratio of cooperator and invasion velocities. When v i wv c , 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 v m . When v c wv i , two qualitatively different outcomes are possible. Either cooperators and defectors can spread together in a mixed wave with velocity v m (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 v c is driven solely by cooperators followed by a slower invasion by defectors with velocity v i .
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 g f (c)~g f and f Ã (c)~f Ã , we can immediately see that f (t,x)~f Ã is a stationary solution of equation (10). With this solution for f (t,x), the dynamics of c(t,x) become identical to that in equation (1) with the parameters c Ã , K, and g c evaluated at f~f Ã . Therefore, the velocity and density profile of the mixed wave can be immediately obtained from equations (3) and (4), e.g., the mixed velocity v m is given by Defectors are expected to decrease the rate of spatial expansions [63]. Consistent with this expectation, equation (11) predicts that v f Ã vv c , provided that all or some of the growth parameters (g c , K, and c Ã ) substantially decrease with f . Note that, because f Ã 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: K(f )~K, g c (f )~g c , and c Ã (f )~c Ã . The partial differential equation for c(t,x) is then independent from the dynamics of f (t,x), so the velocity of the population wave and wave profile are given by equations (3) and (4)  Lf Lt~D where the term with v c appears because we changed to the comoving reference frame. From equation (12), one can see that The effective growth term is always positive and peaks at the middle of the population wave front where c~K=2 and f~0 (by definition). The effective advection term is also positive, provided f decreases with f, 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 v i vv c .
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 f (t,f) in equation (12). When a steady state does not exist, population and genetic waves split, and f (t,f) shifts to negative f with velocity v c {v i . Since G f decreases with f , the existence of a steady state requires that operator L, obtained by linearizing the right hand side of equation (12) with respect to f , has a positive eigenvalue l. Indeed, if all eigenvalues of L are negative, then Lf Lt v0 because G f (c,f )vG f (c,0), while a positive eigenvalue ensures that small f (t,f) will increase until G f (c,f ) is sufficiently diminished so that Lf Lt~0 . Thus, we look for a solution of the following equation with where g l (f)~G f (c(f),0), and we used primes to denote derivatives with respect to f. Equation (14) can be transformed to a canonical form by the following change of variables which also insures that y(+?)~0; see Text S1. The result reads We can now use the standard theory of second order differential equations to establish conditions necessary for the existence of a solution for lw0; see [64,65]. These conditions require that there must be a sufficiently large region where the potential V (f) is negative and the values of V (f) in this region must be sufficiently low. To show this, we multiply both sides of equation (16) 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 V (f)v0 because y 2 (f) §0. This region has a finite width because V (+?)w0, which follows from equations (8) and (13). Indeed, , which is positive since v c wv i . In agreement with the earlier discussion, when v c vv i , the potential at {? is negative ensuring the existence of eigenfunction for lw0 and, therefore, the existence of mixed waves. For f?z?, we find that V (z?)~{g f (0)f Ã (0)z½v c {v a (z?) 2 =(4D), which is positive as long as defectors are sufficiently disfavored at low densities, e.g. when f Ã (0)v0. Thus, V (f) is a potential well, and the existence of an eigenfunction for lw0 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 V(f). The first term lowers V (f) for cw c c in the population bulk (f?{?), but it increases V (f) for cv c c at the front (f?z?). The second term is always negative; it deepens the potential well around f~0 and vanishes for f?+?. The third term is always positive, but the reduction in the depth of the potential well due to v c is lessened by v a (at least for some f). 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 c c, the density at which natural selection starts to favor cooperators over defectors, as shown in figures 2 and 3a. Decreasing f Ã (0) or v i =v c 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 c Ã~0 , and f Ã (c)~f Ã for cw c c, but defectors are not viable for cv c c. In other words, we impose an absorbing boundary condition at f f such that c( f f)~ c c. With these simplifications, we find that splitting occurs provided 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 v i approaches v c . 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 g f . In wellmixed populations, it only determines the rate of approach to the equilibrium not the equilibrium itself, while, in expanding populations, g f affects the invasion velocity and the conditions for wave splitting, both of which are expected to be under selection. The critical density c c is another examples. It has little effect in the well-mixed populations, but large c c 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 v c and v i . 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 lw0 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 c c, is varied so that the system behavior changes from non-splitting to splitting. Increasing c c directly increases the abundance of cooperators at the front and increases the velocity of the mixed wave from v f Ã , when c c~0, to v c , 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, G f is barely sufficient for defectors to keep up with cooperators, and, since G f decreases with f , 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 v i 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 c(0,x)w0 and f (0,x)w0) has a nonzero projection on the eigenfunction of the operator L with the largest positive eigenvalue. The growth of this projection ensures the establishment of defector population at the front even when v i vv c as shown in figure 3b. The ability of defectors to catch up, however, relies on the assumption that c(t,x) and f (t,x) 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 v c wv i 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 v i wv c because defectors can always catch up with cooperators. In contrast, when v i vv c , 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. 10 4 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 frequencyindependent 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) [35,67]. On long time scales, mutations will change the dynamics in two ways. First, the coexistence fraction f Ã will be determined by both natural selection and the mutational pressures similar to the classic theory of two genetic variants [35,67]. This shift in f Ã does not change the dynamics qualitatively and could be included in our theory by modifying G f (c,f ) 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 L pc 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 L pc =(v c {v i ) for deterministic splitting. While the latter time scales as the inverse of the product of the mutation rate from cooperators to defectors m d , population size K(0)L pc , and fixation probability g f ½K(0)f Ã ½K(0), i.e. as fm d L pc K(0)g f ½K(0)f Ã ½K(0)g {1 ; see [35,67]. For deterministic splitting, this balance yields , for stochastic splitting, L pc 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 [30,31,68,69]. 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 [26,30,31,69]. 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 ecoevolutionary 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  figure 3c with the exception of L~2000 for Kv100 and L~3000 for K~100. The error bars show the standard deviation of the mean. Similar increase in the average time to splitting was observed for increasing migration because larger m leads to longer wave profiles containing more organisms, and, therefore, smaller number fluctuations. doi:10.1371/journal.pcbi.1002994.g004 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 [47,48,56,57]. 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.

Supporting Information
Text S1 Text S1 contains the derivation of the boundary conditions for equation (15), the derivation of equation (18), and the discussion of the effects of effective advection and growth on the ability of defectors to travel faster during a range expansion with cooperators than during an invasion into a population of cooperators. (PDF) Figure 5. Population waves can accelerate due to evolutionary changes at the front. (a) A mixed wave of cooperators and defectors spreads into empty territories. The shape of this invasion is not linear in the (t,x) coordinates indicating wave acceleration. The acceleration is caused by a slow increase in the number of cooperators (blue) at the front compared to the bulk (magenta). (b) The velocity of the front increases from v f Ã to its asymptotic value v m , which is smaller than v c because population and genetic waves do not split in this case. (c) Cooperator and defector density profiles at different times show both acceleration of the expansion and the increase in the frequency of cooperators at the front. For this figure, we used K~1, g c~5 , g f~0 :1, c Ã (f )~0:24=(1{f ), f Ã~3 q (c{0:5){2:5, and D~0:5. doi:10.1371/journal.pcbi.1002994.g005