Genetic drift in range expansions is very sensitive to density feedback in dispersal and growth
Gabriel Birzu Sakib Matin Oskar Hallatschek Kirill S. Korolev

Open and closed magnetic configurations of twisted flux tubes
    [
    December 30, 2023
=============================================================







Abstract: 

Theory predicts rapid genetic drift in expanding populations due to the serial founder effect at the expansion front. Yet, many natural populations maintain high genetic diversity in the newly colonized regions. Here, we investigate whether density-dependent dispersal could provide a resolution of this paradox. We find that genetic drift is dramatically suppressed when dispersal rates increase with the population density because many more migrants from the diverse, high-density regions arrive at the expansion edge. When density-dependence is weak or negative, the effective population size of the front scales only logarithmically with the carrying capacity. The dependence, however, switches to a sublinear power law and then to a linear increase as the density-dependence becomes strongly positive. To understand these results, we introduce a unified framework that predicts how the strength of genetic drift depends on the density-dependence in both dispersal and growth. This theory reveals that the transitions between different regimes of diversity loss are controlled by a single, universal parameter: the ratio of the expansion velocity to the geometric mean of dispersal and growth rates at expansion edge. Importantly, our results suggest that positive density-dependence could dramatically alter evolution in expanding populations even when its contributions to the expansion velocity is small.




§ INTRODUCTION


Range expansions and range shifts have become a central theme in population biology because they offer unique opportunities to study rapid adaptation and eco-evolutionary feedbacks <cit.>. Range expansions are also of great practical interest because they describe the spread of pathogens <cit.>, disease vectors <cit.>, and agricultural pests <cit.>. Predicting and controlling the outcome of range expansions requires a firm understanding of their evolutionary dynamics. Indeed, expansion velocities could increase more than five-fold due to the evolution of faster dispersal <cit.> or decrease due to the rapid accumulation of deleterious mutations at the expansion edge <cit.>. 


The latter scenario is caused by repeated bottlenecks at the expansion frontier—a phenomenon known as the founder effect <cit.>. Genetic drift due to the founder effect controls not only the efficacy of natural selection, but also the amount of genetic diversity on which selection can act. In particular, large reductions in diversity are known to negatively impact the resilience of species to environmental stress <cit.>, resistance to parasites <cit.>, and adaptation to climate change <cit.>. Thus, the founder effect is a key determinant of the evolutionary potential of an expanding population.


The strength of the founder effect is often inferred from the reduction in neutral genetic diversity at the range margin compared to the population core <cit.>. In some case, this reduction is dramatic, and a single clone takes over the expansion front <cit.>. 
Most empirical studies, however, report only a marginal reduction in the genetic diversity of 10-30 % <cit.>. 
These observations are inconsistent with many theoretical studies and the naive expectation of low effective population size due to a small number of propagules at the leading edge <cit.>.
A number of mechanisms have been proposed to resolve this paradox <cit.>. For example, diversity loss could be mitigated by multiple introductions <cit.> or the bottlenecks could be less sever if the expansion rate is limited by the rate at which suitable habitat becomes available <cit.>. The latter scenario could be applicable to range shifts during the climate change or re-expansion of the biosphere following an ice age. While these and other mechanisms could resolve the discrepancy in some specific instances, here we focus on density-dependence as a more general factor that can control the rate of genetic drift at the expansion front.


Density-dependent growth, commonly described as an Allee effect, offers a promising explanation for the high genetic diversity of expanding populations <cit.>.
Previous studies have shown that this influx from the core region reduces the size of the bottleneck and increases the effective population size of the front <cit.>. 
Recently, it was found that the effective population size could be proportional to the carrying capacity or scale sublinearly, for example, logarithmically or as a power law <cit.>. The type of the scaling is very sensitive to the strength of the Allee effect. As a result, effective population size can increase by orders of magnitude in response to a moderate suppression of low-density growth. 
Since growth density-dependence is rarely examined at the expansion front, it is possible that Allee effects are responsible for high genetic diversity in many expanding populations.


Compared to density-dependent growth, density-dependent dispersal received very little attention in the context of diversity loss during expansions. This lack of attention most likely reflects the technical difficulties of studying dispersal both analytically and in the wild, rather than its ecological relevance. 
Indeed, nontrivial density-dependence of dispersal on population density has been documented in many organisms, including plants (Echinops sphaerocephalus, Lythrum alatum, Monarda [istulosa) <cit.>, arachnids(Erigone arctica) <cit.>, birds (Parus major, Sitta europeae, Larus ridibundus) <cit.>, and mammals (Cervus elaphus, Capreolus capreolus, Marmota olympus) <cit.>.
The origins of these strategies are also quite diverse, from cooperative responses to avoid starvation <cit.> to mechanisms for population regulation to reduce overcrowding <cit.>.
Here, we address this open problem and show that positive density-dependent dispersal can preserve genetic diversity. This result is first established via computer simulations with two neutral alleles invading a patchy one-dimensional landscape. We then formulate a continuum description based on the generalized Fisher-Kolmogorov equation and derive an explicit formula for the effective population size of the front. 
Our analytical treatment differs from previous work which relied on spatial discretization of the population density or the separation of the expansion front into the fully-colonized and newly established populations <cit.>. Therefore, our results could be used to
suggest new ways to study population dynamics in continuous space. 


Our theory unifies density-dependence in dispersal and growth and 
shows that both processes can be treated on the same footing, with only a minor increase in complexity.
More importantly, we predict that the scaling of the effective population size with the carrying capacity depends on a single parameter—the ratio of the expansion velocity to the geometric mean of the growth and dispersal rates at low densities.
In other words, most aspects of density-dependence are fully captured just by the expansion velocity itself. This emergent universality is confirmed by numerical simulations of a diverse set of models with density-dependence in growth and dispersal.


The main implication of our work is that even relatively weak positive density-dependence is sufficient to prevent rapid diversity loss during range expansion. For example, the rate of genetic drift can be two orders of magnitude slower in species that expanded only 10% faster than they would in the absence of positive feedback in growth or dispersal.
Hence, weak density-dependence could be an overlooked explanation for the surprisingly high diversity of expanding species.




§ METHODS




 §.§ Model

We simulate range expansions using a variant of the stepping-stone model <cit.>. The model consists of a one-dimensional array of islands (demes) with carrying capacity N. To quantify genetic drift, it is sufficient to consider a model with two alleles.
We denote the density of each allele by n_i(t,x), where i ∈{1, 2}. 
Each generation, the densities are updated in three steps. In the first step, we compute the densities ñ_i after dispersal


    ñ_i (t, x) = n_i(t,x) - m[n(t, x)] n_i(t, x) + m[n (t, x-1)]/2n_i(t, x-1) + m[n(t, x + 1)]/2n_i(t, x + 1),


where n(t,x) =  n_1(t,x) + n_2(t,x) is the total population density, and m(n) is the dispersal probability.
Throughout this paper, we focus on parameter regimes which are well described by a continuum model.
In this regime, the assumption that m(n) only depends on the density in the source deme n(t, x) is not restrictive because alternative models can be cast in the same form as Eq. (<ref>).

In the second step, the expected densities n̂_i after growth are computed according to


    n̂_i(t, x) = [ 1 + r(ñ(t, x)) ] ñ_i(t, x),


where r(n) is the per capita growth rate. 

The third step accounts for stochasticity in dispersal and growth. This step is performed by drawing n_i for the next generation from a trinomial distribution


    {n_1(t+1), n_2(t+1), 1-n(t+1) } = Trinomial(N, n̂_1(t)/N, n̂_2(t)/N),


where Trinomial (M, p_1, p_2) is a trinomial distribution with M tries and success probabilities equal to p_1 and p_2 for alleles 1 and 2, respectively <cit.>. The above procedure is equivalent to two successive draws from a binomial distribution. The first draw accounts for demographic fluctuations and determines total population density n = Bin(N, n̂_1 + n̂_2/N). The second draw accounts for genetic drift and determines the genetic composition: n_1 = Bin(N, n̂_1/n) and n_2 = n - n_1.
This procedure is computationally efficient and enforces a strict carrying capacity in the population bulk.
At the edge of the expansion, where the densities are small, the stochastic dynamics are equivalent to the commonly used Poisson distribution for the number of offspring.


Many studies, both theoretical and experimental, have used linear or power law dependencies to describe density-dependent dispersal <cit.>.
To capture a wide range of possible dependencies, we choose the following functional form for the dispersal probability:


    m(n) = m_0[ 1 + A_b( n/N)^b],


where m_0 is a baseline dispersal probability, and A_b quantifies the strength of density-dependence. 
The sign of A_b controls whether the density-dependence is negative (A_b < 0) or positive (A_b > 0). The former is known to enhance unevenness in the spatial distribution of the population and leads to the formation of spatial patterns <cit.>. In contrast, the latter describes strategies that reduce population clustering and helps reduce competition within the population <cit.>. 
The exponent b controls the shape of m(n). For b≪1, the function increases sharply at low densities, compared to b≫1, when most of the increase occurs when n is close to the carrying capacity (Fig. <ref>a).


Simple logistic growth provides an excellent description of the growth dynamics for many species <cit.>.
However, this simple model neglects cooperative interactions and other mechanisms that generate an Allee effect <cit.>. To capture this additional complexity we chose the following widely used functional functional form for the growth rate <cit.> 



    r(n) = r_0  (1 - n/N)  (1 + B n/N) ,


where r_0 is the low-density growth rate, and B is the strength of the Allee effect.
For B ≤ 0 there is no Allee effect. 
For B > 0 we have a weak Allee effect that becomes more pronounced as B increases. 
We do not consider a strong Allee effect (negative growth at low densities) here since previous work has shown that it lead to the same dynamics as those of Eq. (<ref>) with large B <cit.>. 
In such cases, growth rates at low densities are substantially reduced and multiple introductions take place before the population is established at a new location. As a result, the founder effect is greatly suppressed.


The expansion are performed in a fixed box of size L=300 demes, which moves along with the expansion such that the box is approximately half full (see SI for details). This simulation procedure saves computational resources by excluding demes that are either empty or too far behind the front to affect the genetic composition of the expansion edge. Although this is an approximation, its effect on our results is small if the front width is much less than the box size (see SI). In our case, the typical size of the front width is ten demes, which is an order of magnitude less than L.




 §.§ Measuring the strength of genetic drift


We start with a state of maximum diversity by assigning one allele to each individual with equal probability (Fig. <ref>a). 
A complete description of this process involves the frequency of the first allele f(t,x)=n_1(t,x)/n(t, x) that depends on both space and time. 
The spatial variation, however, can be neglected on the long time scale of diversity loss (see Refs. <cit.> and <cit.> for derivation). This approximation assumes that the spatial relaxation of f(t,x) is much faster than the time to fixation, which is valid for realistically large N. The accuracy of replacing f(t,x) by f(t) is illustrated in Fig. <ref>a, which shows that spatial variation in f is much smaller than the change in allele frequency which occurs over time.
Thus, we can replace f(t,x) by the spatial average defined by


    f(t) = ∑_xn(t, x)f(t, x)/∑_xn(t, x).







To quantify the rate of diversity loss we introduce the “heterozygosity”

    H(t)=⟨ f(t)[1 - f(t)]⟩,


where ⟨·⟩ is the ensemble average over independent simulations. 
For populations of haploid asexuals considered here, the heterozygosity represents the probability of randomly choosing two different alleles from the population
After a short transient, H(t) decays exponentially in time <cit.> 



    H(t) ∼ e^-t/N_e.


We refer to the timescale of this decay as the effective population size N_e by drawing an analogy with the dynamics in a well-mixed population. This analogy is known to be imperfect <cit.> because some aspects of population dynamics could be distinct from the Wright-Fisher model. Nevertheless, N_e serves as a convenient and intuitive measure for the rate of genetic drift at the front.




§ RESULTS




 §.§ Positive density-dependent dispersal suppresses the founder effect





How does density-dependent dispersal affect genetic drift in an expanding population?
To answer this question, we simulated range expansions with m(n) = m_0(1 + An/N) and determined N_e as a function of A (see Methods).
We found that while negative density-dependence has a negligible effect on N_e, the effect of positive density-dependence is quite dramatic.
Indeed, increasing A from 1 to 2.5 resulted in N_e changing by an order of magnitude (Fig. <ref>).
Moreover, the rate of change in the effective population size increases as the (density) feedback becomes stronger. 


These results show that the strength of the founder effect can be influenced just as much by dispersal as by an Allee effect.
What is the connection between these two mechanisms and how dispersal affects diversity loss is what we aim to uncover next.
In order to gain a deeper insight into these processes we extend the theoretical framework used in previous studies <cit.> to the case of density-dependent dispersal.




 §.§ Continuum limit


Because of the discrete nature of our simulations, an analytical treatment of the dynamics of the heterozygosity is difficult. However, we can make progress by taking the continuum limit of our model, which is a valid approximation when r_0 ≪ 1 and √(m_0/r_0)≫ 1.
In the SI, we derive the continuum equations and find


    n_it = x[ D(n)n_i ] + r(n)n_i + √(γ_n(n)n_i) η_i(t, x),


where η_i(t, x) is a Gaussian white noise with unit variance, n_i(t,x) is the density of allele i at time t and position x (which are now continuous variables).


The growth term r(n) is the same as in Eq. (<ref>) in the discrete model apart from an extra factor of Δ t^-1, where Δ t is the generation time. The dispersal and noise terms in Eq. (<ref>) are given by


 

    D(n)    = m(n)Δ x^2/2Δ t, 
    γ_n(n)    = 1/Δ t( 1 - n/N),




where Δ x is the distance between neighboring demes.
The placement of the dispersal function D(n) inside both derivatives emerges naturally from our microscopic model, but it is not an extra assumption in our model.
The alternative formulation, in which the first term in Eq. (<ref>) would be replaced by x[D(n) nx], is equivalent to change in the definition of D(n) (see Ref. <cit.> and SI).



 §.§ Deterministic dynamics: pulled and pushed fronts


Reaction-diffusion fronts of the type given by Eq. (<ref>) have been studied extensively in the deterministic limit, when noise is absent. In this limit, the total population density obeys the following equation:


    nt = x[ D(n)n ] + r(n)n.



Depending on the expansion velocity, these fronts fall into two classes: pulled and pushed waves. The expansion velocity of pulled waves is given by the classic Fisher formula v_F=2√(D(0)r(0)), which can be derived by linearizing Eq. (<ref>) for small n. Since the velocity of the low density edge of the front determined the expansion velocity, these fronts are said to be “pulled" by the leading edge. Conversely, pushed expand with v>v_F, as if they are “pushed" by the nonlinearities in the bulk.
Most investigations into these expansions have focused on D(n)=const, and have examined how feedback in r(n) affects the velocity <cit.>. 
However, it has been suggested that the transition between pulled and pushed fronts could occur due to either growth or dispersal feedback <cit.>. Several recent studies corroborate this conclusion <cit.>. 
In particular, exact analytical solution for v was obtained in Ref. <cit.> for a linear D(n), which showed the existence of a transition between pulled and pushed waves.
A more extensive discussion of the differences between pulled and pushed fronts, with an emphasis on density-dependent dispersal, can be found in the SI.




 §.§ Stochastic dynamics: Effective population size


Most of the previous work on genetic drift in expanding populations has focused on pulled fronts <cit.>. These studies have shown that pulled expansions exhibit large fluctuations in allele frequencies, leading to small effective population sizes, which scale as N_e ∼log^3N <cit.>. 
Recent work extended these results to pushed waves and identified two distinct subclasses of pushed waves, with different scaling of N_e with N <cit.>.
However, these studies only considered the case of constant dispersal rates. Here, we derive a complete theory of genetic drift during range expansions that accounts for arbitrary density-dependence in dispersal and growth.


In the SI, we generalize the calculations from Ref. <cit.> to density-dependent dispersal. The calculation relies on treating the front deterministically up to an appropriately chosen cutoff at low densities. To lowest order in N^-1, the effective population size is given by


    N_e = {    C_1 ln^3N,     ν = 1,
        C_2 N^2√(1 - 1/ν^2)/1 - √(1 - 1/ν^2),      1 < ν < 3/2√(2),
        C_3 N,     ν≥3/2√(2),
    .


where we have defined the normalized expansion velocity


    ν = v/2√(r(0)D(0)).




An intuitive interpretation of ν is the ratio between the actual velocity of the expansion v and the velocity inferred from Fisher's formula. Clearly, pulled waves have ν=1, while pushed waves have ν>1.
The constants C_1, C_2, and C_3 in Eq. (<ref>) are independent of N and are given in the SI.


Two critical values of ν divide expansions into three regimes (Fig. <ref>a). For low values of the feedback, the expansions are pulled and the logarithmic scaling of N_e with N is now well established <cit.>. 
At higher values of the feedback, we find two types of pushed waves: semi-pushed waves, for which N_e scales as a sublinear power law in N, and fully-pushed waves, for which we recover the well-mixed scaling N_e ∼ N (Fig. <ref>b).


The surprising implication of Eq. (<ref>) is that the type of the expansion and the scaling of N_e with N depend on a single parameter, the normalized velocity ν. In particular, it is of no consequence whether high values of ν arise due to nonlinear dispersal or an Allee effect. Furthermore, all complexity of D(n) and r(n) are fully encoded in ν.



 §.§ Testing universality


We performed two tests of the universal behavior of N_e(N,ν).
First, we checked that different dispersal models lead to the same scaling and that the scaling exponent matches the prediction from Eq. (<ref>).
Second, we confirmed that the scaling is unchanged when both the density-dependence in dispersal and the strength of an Allee effect are combined in different proportions to keep ν constant.




To carry out the first test we measure N_e(N) in simulations with m(n) given by Eq. (<ref>) for three different values of the exponent b (see Fig. <ref>a).
For each model, we choose the feedback A_b such that the normalized velocities of the expansions are equal in each of the three regimes. 


For pulled waves, the dynamics are fully captured by the linearized equation so N_e should be independent of the nonlinearities in the dispersal model. Indeed, Fig. <ref>b shows that the results for pulled waves collapse on the same curve, which is well fit by log^3N.
For semi-pushed waves, we find a quantitative agreement between the simulations and the power law scaling predicted by our theory (see SI). For fully-pushed waves, the three models exhibit linear scaling, as expected. Thus, our simulations confirm the predicted functional dependence of N_e on N in all three regimes, independent of the shape of m(n).




To carry out the second test of our theory, we vary both the growth and dispersal parameters, while keeping ν constant, for m(n) given by Eq. (<ref>) and r(n) given by Eq. (<ref>). 
This is conceptually illustrated by Fig. <ref>a.
According to our theory, such changes in parameters should not affect the scaling N_e with N. 
This expectation is confirmed by our simulations in Fig. <ref>b.
Similar to the previous test, all pulled waves collapse to the same curve.
In contrast, for semi-pushed and fully-pushed waves only the scaling is predicted to be universal. Thus the intercepts in Fig. <ref>b are given the nonuniversal constants C_2 and C_3.
Across the isocline, the change in C_2 is small and the curves nearly collapse, while C_3 changes by a factor of 2. Still, this change is much smaller than the change in N_e due to N or ν, which is of several orders of magnitude.




§ DISCUSSION


Density-dependent dispersal has been found in empirical studies across the natural world—from plants and animals <cit.> to insects and microbes <cit.>.
Both theoretical and empirical studies have shown that optimal dispersal strategy can offer significant selective advantage, or and even ensure the ultimate survival of the species <cit.>. 
While the evolution of dispersal strategies has been extensively studied, the reverse problem—the effect of dispersal on evolution—has largely been unexplored. 
Here we have addressed this question in the context of range expansions.


We found that positive density-dependent dispersal dramatically reduces the rate of diversity loss, while negative density-dependence has little effect.
The intuition behind this result is that the flux of immigrants from the population core become comparable to the flux of early colonizers.
To explain this behavior quantitatively, we computed N_e for a general model with arbitrary density-dependence in dispersal and growth.
Using this theory, we showed that the change in the effective population size is caused by two transitions, both marked by changes in the scaling of N_e with the carrying capacity N.


Weak or negative density-dependence results in pulled expansions with ν=1 and N_e∼ln^3N. Expansions become semi-pushed as the density dependence increases, and ν begins to exceed 1. Although semi-pushed waves rely on dispersal and growth behind the expansion edge, their population dynamics are still highly stochastic, and N_e grows as a sublinear power law of N. Further increase in the strength of the density dependence leads to a second transition from semi-pushed to fully-pushed waves (ν>3/2√(2)>1). The fluctuations in fully-pushed waves are suppressed, and N_e∼ N. Thus, the naive expectation that the effective population size is given by the number of organisms at the front is only valid for range expansions with strong positive density dependence. Since the strength of the density-dependence is fully encoded in ν, one can determine the class of the range expansion by comparing the observed expansion velocity to the Fisher velocity expected from the growth and dispersal rates at low population densities <cit.>.


Many of our results parallel those for density-independent dispersal in Ref. <cit.>. However, it is a highly nontrivial result of our work that density-dependent dispersal and growth have similar effects on population dynamics and can be treated in a unified way. In particular, the functional dependence of N_e on N and many other ecological and evolutionary processes in an expanding population are controlled by a single parameter, the normalized velocity ν. To determine ν, one needs to measure only three quantities: the expansion velocity v, the dispersal rate at low densities D(0), and the growth rate at low densities r(0). 
Once ν known it can be used to predict the fixation probabilities of neutral mutations or the amount genetic diversity in the population.
Even more importantly, ν can be used to compare different species or different environments. Such standardization would allow researchers to separate the intrinsic dynamics of the population from their specific biological characteristics such as the genetic architecture or mating system. 


Our results suggest an enticing although yet untested possibility of different genealogies in pulled, semi-pushed, and fully-pushed expansions. In well-mixed populations, distinct scaling of N_e with N is known to be a manifestation of fundamentally different structure of the genealogical trees. The N_e∼ln^3N scaling in pulled waves is identical to that in large populations evolving due to weak, but numerous beneficial mutations. The genealogies in such rapidly adapting populations are described by the Bolthausen-Sznitman coalescent with multiple mergers <cit.>. The standard N_e∼ N scaling is observed, for example, for strictly neutral evolution, which leads to the Kingman coalescent for ancestral lineages <cit.>. These two analogies suggest that the genealogies of expanding populations are described by a family of Λ-coalescents that has Bolthausen-Sznitman and Kingman coalescents as two limiting cases for pulled and fully-pushed waves respectively. If this conjecture is correct, we can gain a better understanding of the eco-evolutionary dynamics of range expansions by analyzing the properties of genealogical trees that go beyond the rate of diversity loss. 


In addition to genetics, other properties of the front could be controlled by density feedback.
and calculate the effective diffusion constant of the expansion front, D_f. Similar to previous results <cit.>, we find that D_f^-1 scales with N in the same way as N_e. Thus, the behavior of both evolutionary (N_e) and ecological (D_f) properties of the expansion falls into one of the three distinct classes, depending on the value of the normalized velocity ν. 


The linear scaling of N_e with the population size we observe in fully-pushed waves is apparently very similar to the main result of Ref. <cit.>. However, in deriving this result, Slatkin et al. relied on approximating the propagation by successive founder events, in which the population at the edge of the expansion reaches carrying capacity before the next deme is colonized. While this approximation is frequently made in the literature <cit.>, our results highlight an important limitation of this approach. The approximation effectively divides the front into fully-colonized and newly established populations, and reduces the evolutionary dynamics to those of a two deme model. The loss of diversity in the population can then either be limited by migration—and hence independent of N—or by the time for diversity loss within one deme—which is proportional to N. The nontrivial scaling we observe for pulled and semi-pushed waves is a consequence of the broad distribution of fixation probabilities across the expansion front (see SI from Ref. <cit.>). Thus, to capture the correct behavior of these waves, it is necessary to consider the full dynamics across the width of the front.


How populations maintain genetic diversity during expansions—also known as the genetic paradox of invasions—has been a major question in evolutionary ecology <cit.>. 
This paradox hinges on the studies of pulled expansions, which indeed have a small effective population size that grows only logarithmically with the number of organisms at the expansion front. For fully-pushed expansions, however, the paradox disappears completely because the effective population size is comparable to the census population size of the front, and the rate of diversity loss is therefore exceedingly slow. In this paper, we have shown that pushed expansions and weak genetic drift are quite generic in the presence of positive feedback in growth and dispersal. The prevalence of pulled, semi-pushed, and fully-pushed expansions in the wild is an important and fascinating topic for future investigations. However,
given the prevalence of density-dependence in both dispersal and growth reported in field studies <cit.>, it is possible that a strong founder effect is more of an exception rather than a rule.





Acknowledgements

This work was supported by Simons Foundation Grants #409704 (to K.S.K.) and #327934 (to O.H.); by the Research Corporation for Science Advancement through Cottrell Scholar Award #24010 (to K.S.K.) and the Scialog grant #26119 (to K.S.K.); by the Moore foundation grant #6790.08 (to K.S.K.); by the National Science Foundation Career Award #1555330 (to O. H.); and by the the National Institute of General Medical Sciences of the National Institutes of Health Award R01GM115851 (to O. H.). Simulations were carried out on the Boston University Shared Computing Cluster.









Supplemental Information
    [
    
========================




This Supplemental Information (SI) contains an expanded description of our computer simulations and data analysis, and detailed derivations of the results from the main text. In addition, exact analytical results are derived for some models of density-dependence.




	


§ DENSITY-DEPENDENCE AND POSITIVE FEEDBACK IN MODELS OF RANGE EXPANSIONS
 

We model population expansions using a costless, short-ranged dispersal in the continuum limit. These three assumptions imply that the change in the population density due to dispersal is given by the divergence of a flux:


    nt = -x[ J(n, nx) ]


The flux J can be written in two equivalent forms:


    J = -D(n)nx,


or


    J = -x[ D̃(n) n ].


Both D(n) and D̃(n) can be interpreted as effective diffusion constants and are related by


    D(n) = D̃(n) + D̃'(n)n.


In our analytical derivation, we use D(n) because this form of dispersal can be more easily converted to a Hermitian operator (see Sec. <ref>). In our simulations, we use the other form because dispersal that only depends on the population density at the source deme is naturally described by D̃(n) = m(n)/2. Note that we used the simpler notation D(n) in the main text for the model given by Eq. (<ref>). However, throughout the SI, we follow the convention given by Eqs. (<ref>) and (<ref>).




§ TRAVELING WAVE SOLUTIONS: PULLED AND PUSHED EXPANSIONS
 


In this section, we explain the difference between pulled and pushed waves emphasizing the role of dispersal.
Prior work on this topic largely assumes density-independent dispersal and considers how the fronts change due to nonlinear growth.
Most of the results, however, can be extended to density-dependent dispersal.
Below, we only provide a minimal, heuristic discussion based on the approach from
Ref. <cit.> (SI, Sec. II).
A more thorough presentation of these ideas can be found in Ref. <cit.>.


The equation we consider in the main text is


    ∂ n/∂ t = x[ D(n) nx] + r(n) n.



Assuming sufficiently localized initial conditions, this equation admits a traveling wave solution of the following form:


    n(t,x) = n(ζ), 
       ζ = x - vt,


where without loss of generality we focus on the right-moving part of the expansion.


Upon substituting Eq. (<ref>) into Eq. (<ref>), we obtain the equation for v and n(ζ):


    [D(n)n']' + vn' + r(n)n = 0,


where primes denote derivatives with respect to ζ. 
The boundary conditions for Eq. (<ref>) are given by


    n(-∞) = N,
       n(+∞) = 0.



The solutions of Eq. (<ref>) clearly have a translational degree of freedom, i.e., if n(ζ) is a solution, then n(ζ + const) is a solution. We can eliminate this degree of freedom by choosing the references frame such that n(0)=N/2. The general solution of Eq. (<ref>) then has only one remaining degree of freedom, which corresponds to the value of n'(0). For a given v, the existence of the solution depends on whether the value of n'(0) can be adjusted to match the boundary conditions specified by Eq. (<ref>).


Depending on the behavior of the solution at ζ→±∞, each boundary condition may or may not provide a constraint on the solution and, therefore, remove either one or zero degrees of freedom. To determine the number of constraints, we linearize Eq. (<ref>) near each of the boundaries.



 §.§ Behavior near the boundaries


For ζ→-∞, we let u=N-n and obtain:


    D(N)u” + vu' + r'(N)Nu = 0,


which has the following solution:


    u = A_b e^k_b ζ + C_b e^q_b ζ,


where


    k_b = √(v^2 - 4D(N)r'(N)N) -v/2D(N),
       q_b = -√(v^2 - 4D(N)r'(N)N) - v/2D(N).


Here, we used subscript b to indicate that we refer to the behavior of the population bulk. Because the carrying capacity is an attractive fixed point, dr/dn is negative at n=N, while D(N) > 0, and therefore k>0 and q<0. We then conclude that the boundary condition at ζ→ -∞ requires that C_b=0 and thus selects a unique value of n'(0). 

Next, we analyze the behavior at the front and linearize Eq. (<ref>) for small n:


    D(0)n” + vn' + r(0)n = 0,


which has the following solution:


    n = A e^-k ζ + C e^-q ζ,


where


    k = v+√(v^2 - 4D(0)r(0))/2D(0),
       q = v-√(v^2 - 4D(0)r(0))/2D(0).



The implications of Eq. (<ref>) depend on the sign of r(0). When r(0) < 0 (strong Allee effect), q<0 and, therefore, one needs to set C=0. This constraint uniquely determines the wave velocity. Such waves are classified as pushed because the linearization of Eq. (<ref>) has only solutions that decay in time. 


For r(0)>0, the analysis is more subtle because, naively, the behavior at the front does not fully constrain the velocity of the wave.
Indeed, one can find a solution of Eq. (<ref>) satisfying Eq. (<ref>) for arbitrary v ≥ 2 √(r(0) D(0)).
For pulled waves, q=k and the minimum velocity is selected by the growth dynamics of the linearization of Eq. (<ref>), which can be proved exactly (see Sec. II in the SI of Ref. <cit.>). In this case, the velocity and front decay rate are given by


    v_f = 2 √(r(0) D(0)).


and


    k_f = √(r(0)/D(0)),


To understand the origin of pushed waves, we need to consider the behavior of the population density near the front for v>v_f. 
Equations (<ref>) and (<ref>) predict that n(ζ) is a sum of two exponentially decaying terms with different decay rates: one term decays with the rate k>k_f, but the other with the rate q<k_f.
Note that we use v_F= and k_F defined above as a convenient parameter combination, even though they do not describe the velocity or the decay rate of the population density of the pushed front.


Rigorous analysis of the linearized equation for the front shows that n(t,x) must decay to zero with a rate at least as large as k_f <cit.>. 
Therefore C=0. This extra requirement makes the problem of satisfying the boundary conditions overdetermined. As a result, there are two alternatives. First, there could be no solutions for any v>v_f. In this case, the expansion must be pulled because it is the only feasible solution. Second, a special value of v exists for which the boundary conditions at ζ→±∞ can be satisfied simultaneously. In this case, the wave is pushed because the pulled expansion at low n is quickly overtaken by the faster expansion from the bulk. 


To summarize, expansions can be either pulled or pushed depending on the relative strength of the growth and dispersal behind vs. at the front. Pushed waves expand with v>v_f, and their profile decays as e^-kζ. In contrast, pulled waves expand with v=v_f, and their profile decays as e^-k_fζ. In general, the class of a wave cannot be determined without solving the full nonlinear problem given by Eq. (<ref>). 







§ DYNAMICS OF NEUTRAL MARKERS AND FIXATION PROBABILITIES
 

Heritable neutral markers provide an easy way to study the changes in genetic diversity that occur during a population expansion. These dynamics can be studied either forward-in-time or backward-in-time. Here, we mostly focus on forward-in-time dynamics. A more complete discussion of the differences between the two approaches can be found in Ref. <cit.> (SI, Sec. III). The goal of this section is to introduce the notation that we use in the following sections and explain how density-dependent dispersal affects the dynamics of neutral alleles.




 §.§ Deterministic dynamics of a neutral marker

We consider an expanding population composed of two subpopulations labeled by a neutral genetic marker inherited from parent to offspring. In the following we study two alleles for simplicity, but the results can easily be generalized to an arbitrary number of alleles <cit.>. Since the growth and dispersal rates are the same for both subpopulations, both n_i obey the equation


    n_it = x[ D(n) n_ix] + r(n)n_i,


where n=∑_i n_i is the total population density. It is convenient to define the relative frequency of one of the markers: f = n_1/n. Taking the time derivative of f and using (<ref>) and (<ref>) gives the equation for the relative frequency:


    ft = x[ D(n)fx] + 2D(n)ln nxfx.


The new advection-like term arrises from the nonlinear change of variables and reflects the larger change in f due to dispersal from locations with higher population density. The net effect of this term is a “flow" of f from the back to the tip of the expansion. It is convenient for further analysis to shift to the comoving reference frame:



    ft = ζ[ D(n)fζ] + [ v + 2D(n)ln nζ] fζ.





 §.§ Stochastic dynamics


The equations derived in the previous subsection apply to deterministic fronts. However, genetic diversity is lost due to genetic drift, an inherently stochastic process. Here, we present the equations for n and f taking into account demographic fluctuations and genetic drift. 
Because the fluctuations depend on the details of the reproductive process, we introduce two new functions of the population density: γ_n(n) and γ_f(n), that describe the strength of fluctuations in n and f. The derivation of the stochastic equations for n(t,ζ) and f(t,ζ) is analogous to the D(n)=const case presented in Ref. <cit.> (SI, Sec. III). Below we give the final result:



    nt = ζ[ D(n) nζ] + r(n)n + √(γ_n(n)n) η(t, ζ),



    ft = ζ[ D(n)fζ] + [ v + 2D(n)ln nζ] fζ  + √(γ_f(n)f(1 - f)/n) η_f(t, ζ).








§ EFFECTIVE POPULATION SIZE OF THE FRONT
 



 §.§ Forward-in-time dynamics of the heterozygosity


We measure the genetic diversity in the expansion front using H(t)—the heterozygosity averaged over the front. To derive an expression for H(t), we follow Refs. <cit.> and introduce the spatially-resolved heterozygosity



    H(t, ζ_1, ζ_2) = ⟨ f(t, ζ_1)[ 1 - f(t, ζ_2) ] + f(t, ζ_2)[ 1 - f(t, ζ_1) ] ⟩,


where ⟨·⟩ denotes the average over independent realizations. As defined, H(t, ζ_1, ζ_2) is the probability of sampling different genotypes from location ζ_1 and ζ_2 in the co-moving reference frame. 
We use the shorter term heterozygosity throughout to refer to H(t,ζ_1,ζ_2) when it is unlikely to cause confusion.


To solve for H, we assume the fluctuations are small enough that we can approximate the front profile by the deterministic solution given by Eq. (<ref>), i.e. ρ(t, ζ) ≈ρ_d(ζ). Taking the time derivative of (<ref>) and substituting ft from (<ref>) we obtain


    Ht = [ℒ_ζ_1 + ℒ_ζ_2]H - γ_f(ρ_d)/Nρ_dδ (ζ_1 - ζ_2) H,


where


    ℒ_ζ = ζ[ D(ρ_d)ζ] + [ v_d + 2D(ρ_d) lnρ_dζ] ζ.



The equation for H is a linear differential equation, which can be solved by eigenvector decomposition. At long times, the solution is described by the eigenvector with the largest eigenvalue, which we identify with -1/N_e, following Eq. (8). In general, it is difficult to find the eigenvalue of the full problem, but in the limit of large N, the last term in Eq. (<ref>) can be treated as a perturbation. Using standard techniques, we obtain the following equation for the eigenvalue:



    1/N_e = 1/N∫ dζ_1 dζ_2 δ (ζ_1 - ζ_2) γ_f[ρ_d(ζ_1)]/ρ_d(ζ_1) L(ζ_1, ζ_2) R(ζ_1, ζ_2) /∫ dζ_1 dζ_2 L(ζ_1, ζ_2)R(ζ_1, ζ_2),


where L(ζ_1, ζ_2) and R(ζ_1, ζ_2) are the left and right eigenvectors of the operator ℒ_ζ_1 + ℒ_ζ_2 with the largest eigenvalues. 




 §.§ Determining the left and right eigenvectors


The expression derived previously in Eq. (<ref>) represents a formal solution for N_e. In order to compute a solution we need to find R(ζ_1,ζ_2) and L(ζ_1,ζ_2). Note that we must use only the negative spectrum of the operator ℒ_ζ_1 + ℒ_ζ_2, since H is bounded. The spectrum, therefore, has an upper bound at zero. It is relatively easy to show by inspection that an eigenvector with zero eigenvalue exists, so we use this as a starting point of our analysis.


We begin with the right eigenvector, which is the solution to the equation


    [ℒ_ζ_1 + ℒ_ζ_2]R(ζ_1,ζ_2) = 0,


which can be solved immediately:


    R(ζ_1, ζ_2) = 1.



From Eq. (<ref>) we can derive L(ζ_1, ζ_2) by standard methods. This derivation can be found in standard textbooks on the Fokker-Planck equation, such as Ref. <cit.> (see also Sec. V in the SI of Ref. <cit.> for an analogous derivation for constant  D). Here, we skip these technical details and give the final result:



    L(ζ_1, ζ_2) = ρ_d^2(ζ_1)ρ_d^2(ζ_2)e^v_d ∫dζ_1/D(ρ_d(ζ_1))e^v_d ∫dζ_2/D(ρ_d(ζ_2)).

 
Upon substituting (<ref>) and (<ref>) into (<ref>) we obtain



    N_e = N {∫_-∞^∞ [ρ_d(ζ)]^2 exp[ v_d∫_μ^ζdξ/D(ρ_d(ξ))]dζ}^2 /∫_-∞^∞γ_f(ρ_d) [ρ_d(ζ)]^3 exp[ 2v_d∫_μ^ζdξ/D(ρ_d(ξ))]dζ.



In the expressions above, we explicitly account for the integration constant from the integral over [D(ρ_d)]^-1 by introducing the lower integration limit μ. Note that the contribution from the lower limit in the numerator is cancelled by the denominator and the expression is independent of μ. 




 §.§ Scaling of N_e in pulled, semi-pushed, and fully-pushed expansions


Clearly, if the integrals in Eq. (<ref>) converge then N_e ∼ N. As we now show, a different scaling appears if one of the integrals is divergent. This can be seen by considering the integral in the denominator. All factors in the integrand saturate at constant values or go to zero in the limit of ζ→ -∞. Therefore, if the integral fails to converge it can only do so in the limit ζ→ +∞. In this limit, ρ_d(ζ) decays exponentially (see Eq. (<ref>)). It is convenient to express the decay rate using the normalized velocity ν:



    k = v_d/2D(0)( 1 + √(1 - 1/ν^2)).


Introducing the exponential solution from above into Eq. (<ref>) and keeping only the dominant factors gives the asymptotic behavior of the effective population size:


    N_e ∝ N [ ∫^+∞ e^v_dζ/2D(0)( 1 - 3√(1 - 1/ν^2)) dζ]^-1.


Thus, the convergence of the integral above is determined by the sign of 1 - 3√(1 - 1/ν^2). The value of the normalized velocity at which the linear scaling breaks down then follows immediately:


    ν_c = 3/2√(2).




The divergence identified in Eq. (<ref>) is due to the region of the front where ρ_d(ζ) ≪ 1. In actual populations, this region is finite and no divergence occurs. 
We can account for this in our continuum theory by imposing a cutoff on the population density. 
An naive choice for the cutoff is ζ_c such that ρ_d(ζ_c) = 1/Na, where a is some length scale comparable to the size of a deme. Note that ζ_c is the position of the last unoccupied deme.
However, this naive cutoff fails to account for the large fluctuations in the density at the edge of the expansion.
The correct expressions were computed in Refs. <cit.> and are given by

 
ζ_c =
 1/k [ ln(Nw) + 3lnln(Nw) ],     ν= 1,

1/qln(Nw),     ν> 1,



where w∼ 1/k is the effective width of the front and q is given by Eq. (<ref>).


Substituting Eqs. (<ref>) into (<ref>) gives the formula for the effective population size from Eq. (12) in the main text, with the following values for the constants C_1, C_2, and C_3:


    C_1    = 1/2π^2 r_0,
    
    C_2    = γ_f(0) exp[ 2v_d ( ∫_μ^ζ_cdξ/D(ρ_d(ξ)) - ζ_c/D(0)) ]/[ ∫_-∞^∞ρ_d^2 exp[ v_d ∫_μ^ζ_cdξ/D(ρ_d(ξ))]dζ]^2,
    
    C_3   = ∫_-∞^∞γ_f(ρ_d)ρ_d^3 exp[ 2v_d ∫_μ^∞dζ/D(ρ_d(ζ))]dζ/{∫_-∞^∞ρ_d^2 exp[ v_d ∫_μ^∞dζ/D(ρ_d(ζ))]dζ}^2.








§ FRONT DIFFUSION
 

Demographic fluctuations, represented by the stochastic term in Eq. (<ref>), lead to a stochastic wandering of the front. In this section, we show that this wandering is diffusive and derive the general formula for the diffusion constant D_f. The derivation is an extension of the perturbation theory developed in Refs. <cit.> and <cit.>. The presentation closely follows the analogous calculation in Ref. <cit.> (SI, Sec. VI). We would also like to point the reader to a potentially relevant Ref. <cit.> that used an alternative method to obtain corrections to front velocity.




 §.§ Perturbation theory for demographic fluctuations


We begin with Eq. (<ref>) written in a more convenient form:


    ∂ρ/∂ t = x[ D(ρ)ρx] + r(ρ)ρ +1/√(N)Γ(ρ)η


where ρ=n/N is the normalized population density and Γ denotes the strength of the noise term:



    Γ(ρ) = √(γ_n(ρ)ρ).



We seek a solution to Eq. (<ref>) of the form


    ρ(t,x) = ρ_d(x-v_dt-ξ(t)) + Δρ (t, x-v_dt-ξ(t)),


where v_d and ρ_d are given by the solution of the deterministic equation (Γ=0). Substituting Eq. (<ref>) into Eq. (<ref>) results in the equation for Δρ


    ∂Δρ/∂ t - ℒ_fΔρ =  ρ'_d ξ' + 1/√(N)Γ(ρ_d)η,


where we introduced the operator ℒ_f:


    ℒ_f  =    ζ[ D(ρ_d) ζ] + [ v_d + D'(ρ_d) ρ_d' ]ζ
        + [ r_d(ρ_d) + r_d'(ρ_d)ρ_d + D'(ρ_d)ρ_d” + D”(ρ_d)(ρ_d')^2 ].



To solve (<ref>) we use the fact that ℒ_f has a pair of left and right eigenvectors with zero eigenvalue. 
We denote them by L_f and R_f, respectively.
We shown below that the existence of these eigenvectors is a consequence of the translational symmetry in Eq. (<ref>). 


Multiplying Eq. (<ref>) by L_f(ζ) and integrating over ζ gives


    t∫_-∞^+∞ L_fΔρ dζ  =  ∫_-∞^+∞ L_fρ'_dξ' dζ +  1/√(N)∫_-∞^+∞ L_fΓ(ρ_d)η dζ.


We now use the fact that ∫_-∞^+∞ L_f (ζ) Δρ(ζ) dζ=0. The projection of Δρ on L_f vanishes because we separated translations of ρ_d from shape fluctuations in Eq. (<ref>) are excluded from the fluctuations of the front shape and are instead included through ξ(t). Otherwise, according to Eq. (<ref>), ∫_-∞^∞ L_fΔρ dζ would perform an unconstrained random walk and grow arbitrarily large, which would violate the assumption that Δρ is small. After imposing ∫_-∞^+∞ L_fΔρ dζ=0, we obtain


    ξ'(t) = -1/√(N)∫_-∞^+∞L_f(ζ)Γ(ρ_d(ζ))η(t,ζ) dζ/∫_-∞^+∞L_f(ζ)ρ'_d(ζ)dζ.



From Eq. (<ref>) it is immediately clear that


    ⟨ξ'(t) ⟩ = 0.


Therefore, the translational motion of the front is stochastic and unbiased. To characterize this motion we determine the diffusion constant of the front D_f defined by


    D_f = Var{X_f^2}/2t = ⟨ξ^2⟩/2t,


where X_f is the position if the front. Substituting Eq. (<ref>) into Eq. (<ref>) gives the formal expression for the front diffusion constant:


    D_f = 1/2N∫_-∞^∞γ_n(ρ) [L_f(ζ)]^2 ρ(ζ) dζ/[ ∫_-∞^∞ L_f(ζ) ρ'(ζ) dζ]^2.





 §.§ The left eigenvector


To find L_f we first find R_f and then obtain L_f by transforming operator ℒ_f into Hermitian form.
It is easy to check that R_f(ζ) = ρ_d'(ζ).
This can be seen from the following heuristic argument. Since the unperturbed problem is translationally invariant, any shift in the deterministic front solution ρ_d(ζ) is also a solution. In particular, for an infinitesimal shift we have ρ_d(ζ + δζ) ≈ρ_d(ζ) + ρ_d'(ζ) δζ. Thus, we expect a perturbation Δρ∝ρ_d'(ζ) should leave the left hand side of Eq. (<ref>) unchanged. This can only happen if R_f(ζ) ∝ρ_d'(ζ). 


We now proceed to find the corresponding left eigenvector. This can be done by finding a function β(ζ) such that


    ℒ_f = e^β(ζ)ℋ_f e^-β(ζ),


and ℋ_f is a Hermitian operator.
Because ℋ_f is Hermitian, the left and right eigenvectors coincide and can be denoted by ψ. Then it can be easily shown that



    L_f(ζ) = e^-β(ζ)ψ(ζ),
    
    R_f(ζ) = e^β(ζ)ψ(ζ),


and


    L_f(ζ) = e^-2β(ζ)ρ_d'(ζ).


Substituting ψ(ζ) = e^-β(ζ)ρ_d'(ζ) into the eigenvalue problem ℋ_fψ = 0 gives the expression for β:



    β(ζ) = - v_d/2∫_μ^ζdξ/D[n_d(ξ)]   - 1/2ln[D[n_d(ζ)]].


Finally, upon substituting (<ref>) and (<ref>) into Eq. (<ref>) we arrive at the final expression for the diffusion constant



    D_f = 1/2N∫_-∞^∞γ_n(ρ_d) D^2(ρ_d(ζ)) ρ_d(ζ) [ρ_d'(ζ)]^2 exp[ 2v_d∫_μ^ζdξ/D(ρ_d(ξ))] d ζ/{∫_-∞^∞ D(ρ_d(ζ)) [ρ_d'(ζ)]^2  exp[ v_d∫_μ^ζdξ/D(ρ_d(ξ))] d ζ}^2.



As we have already shown, ρ_d(ζ) ∝ e^-kζ when ζ→ +∞. Therefore, ρ_d'(ζ) ∝ e^-kζ as well and the analysis of the integrals in Eq. (<ref>) is completely analogous to the one for N_e^-1, and the same scaling regimes apply.







§ EXACTLY SOLVABLE MODELS
 

Given importance of exact solutions, we present a few exactly solvable models of range expansions with density-dependent dispersal. Throughout this section we assume γ_n(ρ)=γ_n^0 (1 - ρ) and γ_f(ρ)=γ_f^0, which provide the correct continuum description of our simulations (see Sec. <ref>). Note that exact solutions can also be obtained for other forms of gammas (see Secs. VI and VIII of SI in Ref. <cit.>)


We begin with the model presented in the main main text and defined by the following equations:


    D(ρ)    = D_0 ( 1 +  Aρ),
    
    r(ρ)    = r_0 (1 - ρ).



Recently, exact solutions for the expansion velocity and steady-state profile for this model were found <cit.>. The exact solution shows that the wave is pulled for A<2 and pushed for A > 2. In the pulled regime, the solution is given by



    v    =  2 √(r_0D_0), 
    ρ(ζ)    ∼ e^√(r_0/D_0)ζ,


and in the pushed regime, the solution is given by


    v    = √(r_0 D_0)( √(A/2) + √(2/A)), 
    (1 - ρ)^1 + A/ρ   = e^√(r_0A/2D_0)ζ.



To the best of our knowledge, this is the only model in which a transition from a pulled to a pushed wave due to density-dependent dispersal can be demonstrated analytically. Using these results we can compute the effective population size and the front diffusion constant exactly in the fully-pushed regime. They are given by the following expressions:


    N_e = 6π N/γ_f^0√(2D_0/r_0A)A^2/(A + 2)(A + 4)2π/A,


and


    D_f = γ_n^0/10π N√(2D_0/r_0 A)A (A + 4) (3A + 4)/A + 1tan2π/A.



A separate class of exactly solvable models are given in Ref. <cit.> and are defined by


    D(ρ)    = D_0 ρ^q,
     
    r(ρ)    = r_0ρ^p ( 1 - ρ^k),


with


    k = p + q - 1.



In this case, since D(0) = 0, the waves are always fully-pushed. The exact solutions for the velocity and profile are given by the following formulas:


    v    = zD_0/λ_q
    ρ(ζ)    =
    ( 1 - e^ζ/λ_q)^1/q,    ζ < 0
    
    0,    ζ≥ 0,


where we used 


    λ_q    = √((q + 1)D_0/r_0),
    
    z    = √(q + 1/p + q).


In this case, the results for the effective population size and front diffusion constant are 


    N_e = λ_q N/γ_f^0[ B( q + 2/q -z, z ) ]^2 /B( q + 3/q - 2z, 2z ) .


and


    D_f = q^2 γ_n^0 λ_q/2NB( q + 3/q - 2z, 2 + 2z ) - B( q + 4/q - 2z, 2 + 2z )/[ B( 2 - qz/q, 2 + z ) ]^2.








§ COMPUTER SIMULATIONS
 

In this section, we explain the details of our computer simulations and the subsequent data analysis.

We simulated range expansions of two neutral genotypes in a one-dimensional habitat modeled by an array of patches separated by distance Δ x; the time was discretized in steps of duration Δ t. Thus, the abundance of each genotype was represented as n_i(t,x), where i∈{1,2} is the index of the genotype, and t and x are integer multiples of Δ t and Δ x. Each time step, we updated the abundance of both genotypes simultaneously by drawing from a multinomial distribution with N trials and probability p_i to sample genotype i. The values of p_i were chosen such that p_i N were equal to the expected abundances of the genotypes following dispersal and growth:


    p_i =  1 +  r( ñ )Δ t /N[ n_i(t) - m(n(t, x)) n_i(t, x) + m(n(t, x - Δ x))/2n_i(t, x - Δ x) + m(n(t, x + Δ x))/2n_i(t, x + Δ x) ],

 
where


    ñ =    [1-m(n(t,x))]n_1(t,x) + m(n(t,x-Δ x)/2n_1(t,x-Δ x) + m(n(t,x+Δ x))/2n_1(t,x+Δ x) 
        + [1-m(n(t,x))]n_2(t,x) + m(n(t,x-Δ x)/2n_2(t,x-Δ x) + m(n(t,x+Δ x))/2n_2(t,x+Δ x),


which is the expected total population density after dispersal, but before growth. For simplicity, we set both Δ x and Δ t to 1 in all simulations.




 §.§ Mapping simulation parameters to continuum model


To find the reaction-diffusion equation corresponding to our simulation algorithm we perform a Taylor expansion of the population density:


    n_i (t+Δ t, x)    = n_i (t, x) + Δ tn_it
    
    n_i(t, x ±Δ x)    = n_i(t, x) ±Δ xn_ix  + Δ x^2/2n_ix
    
    m [ n(t, x ±Δ x) ]    = m[n(t, x)] + m'[n(t, x)] ( ±Δ x nx + Δ x^2/2nx) + Δ x^2/2 m”[n(t, x)] ( nx)^2,

where we've kept terms up to first order in Δ t and second order in Δ x.


Substituting (<ref>) into Eq. (<ref>) and summing over the two genotypes gives the equation for the total population density:



    nt    = Δ x^2/2Δ t[ m(n)   + m'(n) n ] nx + Δ x^2/2Δ t[ m”(n) n + 2m'(n) ] (nx)^2  + r ( n ) n,


which can be rewritten as



    nt  = x[ Δ x^2 m(n)/2Δ t n ] + r (n) n.


In the continuum limit, when r(n)Δ t≪1 and kΔ x≪1, our model becomes equivalent to Eq. (<ref>) for the population density and to Eq. (<ref>) for the relative fraction of the two genotypes with



    D̃(n)    = Δ x^2 m(n)/2Δ t, 
    γ_n (n)    = 1-n/N/Δ t, 
    γ_f (n)    = 1/Δ xΔ t.



Note that D(n) = D̃(n) + D̃'(n)n as discussed in Sec. <ref>.




 §.§ Boundary and initial conditions

The most direct approach to simulating a range expansion is to use a stationary habitat, in which the range expansion proceeds from one end to the other. This approach is, however, inefficient because the computational time grows quadratically with the duration of the simulation. Instead, we took advantage of the fact that all population dynamics are localized near the expansion front and simulated only a region of 300 patches comoving with the expansion. Specifically, every simulation time step, we shifted the front backward if the total population inside the simulation array n_array exceeded 150N, i.e. half of the maximally possible population size. The magnitude of the shift was equal to ⌊ (n_array-150N)/N ⌋+1. The population density in the patches that were added ahead of the front was set to zero, and the number of the individuals moved outside the box was stored, so that we could compute the total number of individuals n_tot in the entire population including both inside and outside of the simulation array. 


Our choice of 300 patches in the simulation array was sufficient to ensure that at least one patch remained always unoccupied ahead of the expansion front and that the patches shifted outside the array were always at the carrying capacity. While neglecting the effect of patches deep inside the bulk is an approximation, its effect is relatively small if the size of the box is large compared to the front width. This can be seen from Eq.(<ref>), where the contribution of the bulk to the two integrals decays exponentially as e^-k_b | ζ |. In our case, 1/k_b is of the order of 10 patches and, thus, the error in determining N_e is negligible.


We initialized all simulations by leaving the right half of the array unoccupied and filling the left half to the carrying capacity. In each occupied patch, we determined the relative abundance of the two neutral genotypes by sampling from the binomial distribution with N trials and equal probabilities of choosing each of the genotypes.




 §.§ Computing front velocity

The velocity of the front was measured by fitting n_tot/N to vt+const. For this fit, we discarded the first 10% of the total simulation time (1000 generations for the shortest runs) 
to account for the transient dynamics. The length of the transient is the largest for pulled waves and is specified by the following result from Ref. <cit.>:


    v(t) = v_F[ 1 - 3/4r_0 t + 𝒪( t^-3/2) ].


Thus, discarding time points prior to t∼1/r_0 was sufficient to eliminate the transient dynamics in all of our simulations.




 §.§ Computing heterozygosity and the rate of its decay

For each time point, the heterozygosity h was computed as follows


    h(t) = 1/n_nonzero∑_x, n(t, x) ≠ 02n_1(t,x)n_2(t,x)/[n_1(t,x)+n_2(t,x)]^2,
 

where n_nonzero is the number of occupied patches. The average heterozygosity H was then obtained by averaging over independent simulation runs. To compute N_e, we fitted ln H to -t/N_e +const.


Because of the initial transient dynamics and large uncertainty in H for large t, due to the small number of replicas with H > 0, we restricted the analysis of H(t) to t∈(t_i,t_f).
The value of t_f was chosen such that at least 50 simulations had non-zero heterozygosity at t=t_f. The value of t_i was chosen to maximize the goodness of fit (R^2) between the fit to H∼ e^- t/N_e and the data subject to the constraint that t_f-t_i>1000. The latter constraint ensured that we had a sufficient number of uncorrelated data points for carrying out the fitting procedure to determine N_e.




 §.§ Comparison between simulations and theory


Table S1 shows a comparison between the theoretical prediction and the exponents obtained by fitting the data in Fig. 4 (main text). In the fully-pushed regime, the data are consistent with the linear scaling predicted by the theory. In the semi-pushed regime, the agreement is less good for the square root density-dependence and, to a lesser extent, for the linear density-dependence. This may reflect a larger contribution of the subleading corrections to N_e in these models. For pulled waves, the asymptotic log^3 N scaling sets in at very large values of N, as demonstrated in Ref. <cit.>. 
Nevertheless, as Figs. 4 and 5 show, our simulations are consistent with N_e ∼log^3 N, as predicted by theory.






