Genetic Demixing and Evolutionary Forces in the One-Dimensional Stepping Stone Model
K.S. Korolev∗
Department of Physics and FAS Center for Systems Biology, Harvard University, Cambridge, Massachusetts 02138, USA

Mikkel Avlund
Department of Physics and FAS Center for Systems Biology, Harvard University, Cambridge, Massachusetts 02138, USA and Niels Bohr Institute, University of Copenhagen, Denmark

arXiv:0904.4625v2 [q-bio.PE] 13 Apr 2011

Oskar Hallatschek
Max Planck Research Group for Biological Physics and Evolutionary Dynamics, Max Planck Institute for Dynamics & Self-Organization (MPIDS), G¨ttingen, Germany o

David R. Nelson†
Department of Physics and FAS Center for Systems Biology, Harvard University, Cambridge, Massachusetts 02138, USA (Dated: April 14, 2011)
We review and extend results for mutation, selection, genetic drift, and migration in a onedimensional continuous population. The population is described by a continuous limit of the stepping stone model, which leads to the stochastic Fisher-Kolmogorov-Petrovsky-Piscounov equation with additional terms describing mutations. Although the stepping stone model was ﬁrst proposed for population genetics, it is closely related to “voter models” of interest in nonequilibrium statistical mechanics. The stepping stone model can also be regarded as an approximation to the dynamics of a thin layer of actively growing pioneers at the frontier of a colony of microorganisms undergoing a range expansion on a Petri dish. We ﬁnd that the population tends to segregate into monoallelic domains. This segregation slows down genetic drift and selection because these two evolutionary forces can only act at the boundaries between the domains; the eﬀects of mutation, however, are not signiﬁcantly aﬀected by the segregation. Although ﬁxation in the neutral wellmixed (or “zero dimensional”) model occurs exponentially in time, it occurs only algebraically fast in the one-dimensional model. We also ﬁnd an unusual sublinear increase in the variance of the spatially averaged allele frequency with time. If selection is weak, selective sweeps occur exponentially fast in both well-mixed and one-dimensional populations, but the time constants are diﬀerent. The relatively unexplored problem of evolutionary dynamics at the edge of an expanding circular colony is studied as well. We also brieﬂy review how the observed patterns of genetic diversity can be used for statistical inference, and highlight the diﬀerences between the well-mixed and one-dimensional models. Although we focus on two alleles or variants, q-allele Potts-like models of gene segregation are considered as well. Most of our analytical results are checked with simulations, and could be tested against recent spatial experiments on range expansions oﬀ inoculations of Escherichia coli and Saccharomyces cerevisiae.

Keywords: stepping stone model, stochastic Fisher-Kolmogorov-Petrovsky-Piscounov equation, selective sweep, voter model, Eden model.

Contents
I. Introduction II. Population Genetics In Well-Mixed Populations III. One-Dimensional Stepping Stone Model IV. Neutral Model Without Mutations 2 4 8 9

V. Neutral Model With Mutations VI. Selection VII. Simulations VIII. Inﬂation IX. Genetic Inference X. Conclusions Acknowledgments

11 12 14 17 19 22 23 23

∗ Electronic

address: papers.korolev@gmail.com † Electronic address: nelson@physics.harvard.edu

A. The Itˆ calculus o

2
B. Solution of the Neutral Model Without Mutations C. Average domain density from the spatial heterozygosity H(t, x) D. Inﬁnite Alleles Model E. A model with several neutral alleles F. Connection with the voter model and one-dimensional reaction kinetics References

24 25 25 26 28 29

I. INTRODUCTION

The quantitative theory of evolution is an important open problem. The theory is necessary to determine the history of species migrations, and it could shed light on the origin and development of life. Moreover, a better understanding of the evolutionary dynamics could help control epidemics (Murray, 2003), ﬁght diseases with an evolutionary character such as cancer and acquired immune deﬁciency syndrome (Nowak, 2006), and guide the engineering of artiﬁcial evolution for practical applications (Bar-Yam, 2005; Poli et al., 2008). Most of the current understanding of evolutionary dynamics comes from population genetics, a scientiﬁc discipline that studies how evolutionary forces shape the genetic diversity of populations. The majority of theoretical models and experiments in population genetics study only one or a few well-mixed populations, i.e. populations without spatial structure, where every individual is equally likely to interact with any other individual inside the same population. Microorganisms growing and evolving in a well-mixed liquid culture provide an important example. While nonspatial models are often easier to analyze than spatial ones, they do miss what can be essential features of natural populations. In nature, organisms often occupy areas that are much larger than the square of the dispersal distance, i.e. the distance typically traveled by an individual in one generation. This causes two main problems for well-mixedpopulation models. First, well-mixed-population models underestimate the role of genetic drift (ﬂuctuations due to the discreteness of the number of individuals). The diﬀerence arises because the organisms can only interact with their neighbors, and the number of neighbors within the dispersal distance is much smaller than the total number of organisms in the entire population. Second, wellmixed-population models neglect the spatial structure of the population that can be created by external factors or by internal dynamics. Such spatial structures often exist, and, as we show in this paper, they can signiﬁcantly aﬀect evolutionary processes in the population. Well-mixed-population models are particularly inadequate when applied to expanding populations. Expansions are very common in biology. Species spread to new territories from the locations where they ﬁrst evolved.

FIG. 1 (Color online) Spatial segregation in an expanding microbial population. Diﬀerent colors label diﬀerent alleles. The Petri dish was inoculated with a well-mixed population occupying a narrow horizontal linear region between the arrows, which show the direction of the growth. As this population expands, it segregates into well deﬁned monoallelic domains. The colony is of the order 1 cm in height. Details of the experiment are presented in Hallatschek et al. (2007).

Expansions also occur because of environmental changes such as the global warming and the glacial cycles or due to sudden long distance migrations to new habitats. Even though well-mixed-population models can account for the growing number of individuals (population size), these models do not capture the fact that the newly settled areas are colonized by the oﬀspring of only a small number of individuals at the expanding front. Since the ancestral population is small, the genetic drift is strong. As a result, neutral genetic diversity decreases with the distance from the origin of the expansion. This reduction in genetic diversity, which is often called “the founder eﬀect” (Mayr, 1942), has been observed in humans (Ramachandran et al., 2005; Templeton, 2002) and many other species. For example, the founder eﬀect in the population waves following the receding glaciers is believed to be responsible for the reduced genetic diversity in high latitude regions compared to equatorial ones (Hewitt, 1996). The spreading of Escherichia coli (E. coli) and Saccharomyces cerevisiae (S. cerevisiae) on Petri dishes has been investigated in recent experiments by Hallatschek et al. (2007). In these experiments, microbes grown in the dark carried one of two selectively neutral alleles, diﬀering only in a gene encoding for proteins with two distinct ﬂuorescence spectra. Figure 1 shows the expansion of an initially well-mixed 50 : 50 population of E. coli into two unoccupied half planes initiated by a razor blade inoculation with cells grown up in liquid culture. The distinctive feature illustrated by the typical experiment in Fig. 1 is that the population does not remain well-mixed; instead, it segregates into well-deﬁned domains. The segregation occurs because the strong genetic drift associated with reduced population size facilitates ﬁxation of one of the two alleles at the front. Analogous phenomena should also occur in a nonexpanding one-dimensional population because its dynamics is similar to the dynamics of the front of a grow-

3 ing population. The front of a population wave and a literally one-dimensional habitat are not exactly equivalent because the contour of the front undergoes undulations while a one-dimensional habitat has a ﬁxed linear shape. Nevertheless, both are eﬀectively one-dimensional and should deviate from the predictions of well-mixedpopulation models in similar ways. The advantage of a ﬂat one-dimensional habitat is that it is easier to analyze. In addition, although most species live in eﬀectively two dimensional habitats, a quasi one-dimensional habitat could describe a bank of a river, a sea coast, and a slope of a linear mountain range. To study the dynamics of a population analytically, we adopt the stepping stone model proposed by Kimura and Weiss (Kimura and Weiss, 1964). This model considers many well-mixed populations, demes, located on a spatial lattice. Each deme is subject to mutation, selection, genetic drift, and short range migration between neighboring demes. In the limit of weak evolutionary forces and large number of demes, the stepping stone model is equivalent to the continuous models proposed by Wright (1943) and Mal´cot (1955) and is dee scribed by the stochastic Fisher-Kolmogorov-PetrovskyPiscounov equation (Fisher, 1937; Kolmogorov et al., 1937) with additional terms representing mutation. On the other hand, when each deme contains only one organism, the model is analogous to the Eden model (Saito and M¨ller-Krumbhaar, 1995) used to u describe the growth of interfaces and the voter model (Cox and Griﬀeath, 1986) discussed in Appendix F. We also performed numerical simulations to better understand the relationship between the experiments in Hallatschek et al. (2007) and our analytical results. An illustrative simulation (with periodic boundary conditions) is shown in Fig. 2, which also shows the diﬀerence between a growing population front with undulations and a literally one-dimensional habitat advancing uniformly in time. Figure 3 shows qualitative agreement between the experiments and the simpliﬁed row-by-row growth model that we studied analytically. In this paper, we ﬁrst focus on the spatial segregation due to genetic drift and its eﬀect on the dynamics of a linear one-dimensional population. We ﬁnd that segregation of two neutral alleles has two stages. During the ﬁrst stage, distinguishable domains emerge from the well-mixed population. During the second stage, domain boundaries diﬀuse and annihilate upon collision. As a result, some of the domains vanish whereas others grow. We show how our calculations might be used to extract the diﬀusion constant and the eﬀective population size from experiments like those in Hallatschek et al. (2007), and discuss how well the model describes the behavior of microbes. A detailed comparison (beyond the qualitative agreement we ﬁnd with the main features) would require more extensive and precise experiments; we hope such experiments will be carried out in the future. The spatial segregation dramatically changes the eﬀects of ge-

FIG. 2 (Color online) An illustration of the two models of a growing front. (a) and (c) illustrate the model with a rough, undulating front, which is a natural result of an unconstrained two-dimensional growth. (b) and (d) illustrate the model with a ﬂat front, which is constrained to have no lateral undulations to simplify the analytical analysis. The blank hexagons represent empty sites, and diﬀerent colors of the occupied hexagons represent diﬀerent alleles. (a) The model of an undulating population front. The highlighted hexagon is a randomly chosen cell that can reproduce and deposit an identical oﬀspring in any of its four empty nearest neighbor sites (shown with arrows) with equal probability. (b) The model of a one-dimensional habitat, where each row represents a generation. Thus, each row is completed before moving on to the next one, so an empty site can be ﬁlled only by an oﬀspring of one of its nearest neighbors in the previous generation (shown with arrows). Both (a) and (b) show the eﬀects of genetic drift (sampling error) when, e.g., the second from the left cell in the bottom row leaves no oﬀspring. Such events lead to coarsening seen in (c) and (d). (c) and (d) are single simulation runs for models in (a) and (b) respectively. A population of 100 cells was wrapped around a cylinder to illustrate periodic boundary conditions used in this paper. Note that in (d) the front is ﬂat whereas in (c) it is rough. This roughness aﬀects some aspects of the shapes of the monoallelic domains shown in (c): A domain boundary followed from its lowest point to its highest point always goes up in (d), but, in (c), it sometimes goes down. As discussed in Hallatschek et al. (2007), domain walls are expected to wander more vigorously in (c) than in (d). Despite the apparent diﬀerences, both models exhibit the same qualitative behavior.

netic drift and selection on the population compared to the predictions of well-mixed-population models. For the neutral model without mutation, we ﬁnd that local diversity or “heterozygosity” decays as t−1/2 , and the standard deviation of the global fraction of an allele grows subdiﬀusively as t1/4 . The evolutionary dynamics during a radial expansion (see Fig. 14) is studied as well. In this case, migration and genetic drift slowly weaken as the circumference grows. As the result, the domains boundaries eventually stop coalescing leading to a ﬁnite number of domains in the long time limit. We ﬁnd that

4 correlation functions. In Sec. IV and Sec. V we solve these equations for zero and nonzero mutation rates respectively. While the neutral stepping stone model has been treated before, we derive some new results and use a diﬀerent technique that can be easily extended to radially expanding populations. The eﬀects of selection are considered in Sec. VI, and in Sec. VII we test our analytical results with simulations. In Sec. VIII, evolutionary dynamics during a radial range expansion is analyzed, and Sec. IX deals with genetic inference. Various details are relegated to Appendices A–F. In Appendix E, we indicate how some of the 2-state (i.e., “2-allele”) results can be generalized for the Potts-model-like nonequilibrium dynamics of q-alleles with q ≥ 3.
FIG. 3 (Color online) Qualitative comparison of a gene segregation experiment from a linear inoculation (inset) and the simulation of a one-dimensional habitat. The experiment is analogous to the one depicted in Fig. 1. II. POPULATION GENETICS IN WELL-MIXED POPULATIONS

this ﬁnal number of domains grows as a square root of the initial radius of the colony. We also study the dynamics in the presence of weak selection and ﬁnd that it diﬀers markedly from that of a well-mixed population. Because of the spatial segregation into domains, selection acts only near domain boundaries, which constitute only a small fraction of the population. Hence, extinction of a deleterious allele proceeds much more slowly in one-dimensional populations than in well-mixed populations. Unlike genetic drift and selection, the eﬀects of mutation in the spatial model are essentially the same as in the well-mixed-population model, but the spatial model gives a more accurate description of the population and accounts for the spatial correlations. Finally, we discuss how one can estimate important model parameters by sampling and sequencing DNA from organisms in a natural population. The diﬀerences between spatial and nonspatial models used for genetic inference are highlighted. A substantial fraction of our results for the neutral dynamics in a one-dimensional habitat has been derived previously in population genetics (Barton et al., 2002; Kimura and Weiss, 1964; Mal´cot, 1975), ecole ogy (Houchmandzadeh and Vallade, 2003), and nonequilibrium statistical mechanics (Bramson and Lebowitz, 1991; Cox and Griﬀeath, 1986). Here, we present a single self-contained derivation of these earlier results in a novel context of expanding populations in two dimensions and in a language familiar to physicists, with future microbial tests of the theory in mind. Our new results are primarily conﬁned to the analysis of radial expansions and natural selection. This paper is organized as follows. First, we review classical results for well-mixed populations in Sec. II. We then introduce the one-dimensional stepping stone model in Sec. III and derive the equations of motion for spatial

Well-mixed-population models are relevant to microorganisms vigorously shaken in a test tube, but they do not describe spatial phenomena. Indeed, if cells visit all parts of the test tube during a cell division time, they live in an eﬀectively zero-dimensional habitat. Nevertheless, well-mixed-population models can serve as a useful reference point to which spatial models can be compared. Nonspatial populations also provide a simple context to introduce genetic drift, mutation, and selection; and the stepping stone model presented in Sec. III uses a well-mixed-population model to describe the dynamics of allele frequencies within the demes. This section summarizes the classical results of nonspatial population genetics, which are primarily due to Wright, Fisher, Haldane, and Kimura; the books by Hartl and Clark (1989), and Crow and Kimura (1970) provide a good introduction to the subject and refer to the original literature, which is too extensive to be discussed here; see also Blythe and McKane (2007) for a recent review written for physicists. To simplify the discussion and to make a direct connection with the experiments in Hallatschek et al. (2007), we consider two alleles in a population of N haploid organisms, i.e. organisms with a single set of chromosomes. 1 The two-allele approximation may seem very restrictive, but many of our results can be generalized to an arbitrary integer number of q ≥ 3 alleles. In addition, a two-allele model can be used to describe the dynamics of an allele of interest (with or without a selective advantage) when all other alleles have the same ﬁtness. We assume that each of the individuals in the population can die, give

1

The theory of haploid organisms also describes the dynamics of genes in cellular organelles like mitochondria and chloroplast and on certain sex chromosomes like Y-chromosome in Homo sapiens. For N diploid organisms, the theory is essentially the same under certain assumptions, provided one focuses on the dynamics of 2N gene copies in each generation; see Hartl and Clark (1989).

5 birth (divide), and mutate. The details of this birth and death process are species dependent, but the dynamics on time scales larger than the generation time τg is believed to be universal provided N is large. This universal dynamics is often referred to as the diﬀusion or continuous approximation. Two simple models are commonly used to illustrate the continuous approximation: the WrightFisher model and the Moran model. Here, we use the latter because it more closely resembles microbes with overlapping generations. First, we consider the Moran model without selection and mutation. During a time step, two individuals are randomly selected with replacement from the population. The ﬁrst individual is chosen to reproduce, and the second one to die; thus, the total number of the organisms is conserved. If the “frequency” of allele one (i.e., the fractional number of individuals with genotype one) ˜ ˜ at time step t is f (t), then, at the next time step, it is f + 1/N with probability f (1 − f ), f − 1/N with probability f (1 − f ), and f with probability f 2 + (1 − f )2 . ˜ The expectation value and variance of f (t + 1) are then given by, ˜ ˜ f (t + 1) = f (t), (1) 1989) lead to an equation identical to Eq. (3), but with a diﬀerent numerical coeﬃcient in Eq. (4). Equation (3) is subject to absorbing boundary conditions 2 at f = 0 and f = 1 because, if one of the alleles is lost, it cannot appear again in the absence of mutation. Therefore the population eventually becomes ﬁxed at one of the absorbing states. We calculate the rate of the ﬁxation by considering the average heterozygosity of the population H(t) ≡ h(t) = 2f (t)[1 − f (t)] , (5)

which is the (averaged over realizations) probability that two randomly selected individuals have diﬀerent alleles. When the population is close to the ﬁxation (f ≈ 0 or f ≈ 1), the heterozygosity is close to zero. The equations of motion for F (t) ≡ f (t) and H(t) follow from Eq. (3) by multiplying both sides with f or h, integrating over f , and eliminating the derivatives with respect to f via integration by parts. The results are dF (t) = 0, dt dH(t) = −Dg H(t). dt

(6)

˜ ˜ [f (t + 1) − f (t + 1) ]2

˜ ˜ 2f (t)[1 − f (t)] = , N2

(2)

(7)

where angular brackets represent average with respect to the random choice of individuals for reproduction and death. Because only one of N organisms gives birth in a ˜ Moran time step, t measures time in fractional generation time, τg /N . ˜ Equations (1) and (2) imply that f (t) performs an unbiased random walk in the space of allele frequencies. In the continuum limit, this random walk can be described by the following Fokker-Planck equation with a frequency dependent diﬀusion coeﬃcient (Crow and Kimura, 1970; Hartl and Clark, 1989) ∂P (t, f ) Dg ∂ 2 [f (1 − f )P (t, f )] , = ∂t 2 ∂f 2

Equations (6) and (7) imply that, while the average frequencies of these neutral alleles do not change F = f = f (t = 0) ≡ F0 , the population reaches ﬁxation exponentially fast, H(t) = H(0)e−Dg t = F0 (1 − F0 )e−Dg t . The average heterozygosity is closely related to the variance of f (t), the fraction of the ﬁrst allele, 1 V (t) = (f (t) − f (t) )2 = F (t)[1 − F (t)] − H(t). (8) 2 Thus, even if a population starts with zero variance, the ﬂuctuations grow until the variance reaches its maximum value of F0 (1 − F0 ), which corresponds to a population ﬁxed to allele one with probability F0 and to allele two with probability 1−F0 . Note that, for small t, V (t) grows linearly with time, but, at large times, the variance approaches its limiting value exponentially fast. The linear growth of variance at small times also follows from the Fokker-Planck equation because, at small times, Eq. (3) can be approximated by a diﬀusion equation with a constant diﬀusivity. Next, we generalize Eq. (3) to account for mutations. In the Moran model, mutation is included at the end of a

(3)

where P (t, f ) is the probability density function for f at time t measured in generations, and Dg is the genetic diﬀusion constant. Here, t is the time measured in generations; as discussed above, N Moran time steps constitute a generation time τg . Thus, in the Moran model, we have 2 . N τg

Dg =

(4)
2

Alternative reproduction schemes, such as Wright-Fisher sampling, (Crow and Kimura, 1970; Hartl and Clark,

Since Eq. 3 is singular at the boundaries, we require limf →0,1 f (1 − f )P (t, f ) = 0. See (Risken, 1989) and (Kimura, 1955) for a more detailed discussion.

6 time step by allowing the oﬀspring to mutate with probability µ12 from allele one to allele two and with proba˜ bility µ21 from allele two to allele one. If the frequency ˜ ˜ ˜ of allele one at time step t is f (t), then, at the next time ˜ step, the expectation value of f (t + 1) is given by ˜ ˜ f (t + 1) = f (t) + ˜ ˜ µ21 [1 − f (t)] − µ12 f (t) ˜ ˜ , N From Eqs. (14) and (8), we see that, when the population size is large enough, i.e. Dg ≪ (µ12 + µ21 ), H(∞) ≈ 2F (∞)[1 − F (∞)], the stationary value of the heterozygosity is consistent with f (t) ≈ F (∞). Thus V (∞) ≈ 0, and the ﬂuctuations of f (t) are negligible. In the oppo+µ site limit, H(∞) ∝ µ12Dg 21 is signiﬁcantly smaller, which suggests that most of the time the population is ﬁxed to one of the alleles, and mutations lead to rare transitions between states with f = 0 and f = 1. Consequently, the stationary distribution is dominated by the regions around f = 0 and f = 1, as one can see in Fig. 4. Our interpretation of Eq. (14) is consistent with a more rigorous analytical and numerical analysis by Duty (Duty, 2000). Finally, we introduce Darwinian natural selection, which is usually related to the diﬀerence in the reproduction or survival probability of the organisms. In the continuous time limit considered here, both mechanisms of selection lead to the same dynamics; therefore, we only consider selection due to diﬀerent growth rates. In the Moran model, a growth rate diﬀerence is embodied in modiﬁed probabilities of reproduction: the individual with allele one is chosen to reproduce not with probabilw1 ity f but with probability w1 f +w2f(1−f ) , where w1 and w2 are the ﬁtnesses (i.e., growth rates) of alleles one and two respectively. In the absence of mutations, this modiﬁcation results in ˜ ˜ f (t)[1 − f (t)](w1 − w2 ) ˜ ˜ . f (t + 1) = f (t) + ˜) + w2 [1 − f (t)]} ˜ N {w1 f (t

(9)

˜ and the variance of f (t + 1) is given by Eq. (2) to the leading order in the mutation rates and the inverse population size. ˜ Since the expectation value of f (t) changes with time, mutation leads to an f -dependent drift term in the Fokker-Planck equation. Upon recalling that N Moran time steps equal one generation time, we have ∂ ∂P (t, f ) =− {[µ21 − (µ12 + µ21 )f ]P (t, f )} ∂t ∂f Dg ∂ 2 [f (1 − f )P (t, f )] , + 2 ∂f 2

(10)

where µ12 ≡ µ12 τg and µ12 ≡ µ21 τg are the mutation ˜ −1 ˜ −1 rates per generation. Because the alleles can mutate into each other, the probability ﬂux through the boundaries must be zero, so Eq. (10) has reﬂecting boundary conditions, and a nontrivial stationary solution for P (t, f ) exists. While the stationary distribution can be obtained easily, see Eq. (20), Fig. 4, and Crow and Kimura (1970), it is suﬃcient to analyze the moments F (t) and H(t) introduced above. This will also allow us to make a direct comparison with the corresponding solutions of the one-dimensional stepping stone model in Sec. V. The equations of motion for F (t) and H(t) are obtained from Eq. (10) in the same way as for the absence of mutation. The results are dF (t) = µ21 − (µ12 + µ21 )F (t), dt dH(t) = − (Dg + 2µ12 + 2µ21 )H(t) dt + 2[µ21 + (µ12 − µ21 )F (t)].

(15)

When selection is weak, that is |w1 − w2 | ≪ w1 + w2 , Eq. (15) reduces to s ˜ ˜ ˜ ˜ ˜ f (t + 1) = f (t) + f (t)[1 − f (t)], N

(16)

(11)

(12)

Since these equations are linear diﬀerential equations with constant coeﬃcients, the equilibrium is approached exponentially fast. The stationary solutions, which are obtained in the limit t → ∞, are given below F (∞) = µ21 , µ12 + µ21 (13)

where s = 2(w1 −w2 )/(w1 +w2 ) is the selective advantage ˜ of allele one, which has to be much smaller than one for the approximation to hold. When s > 0, allele one is ad˜ vantageous; for s < 0, it is deleterious. In the following, ˜ we assume that allele one is advantageous because one can always relabel the alleles to satisfy this condition. Similar to the case of mutations without selection, the ˜ variance of f (t + 1) is given by Eq. (2) to the leading order in s and N −1 , and the corresponding Fokker-Planck ˜ equation acquires an f -dependent drift term due to selection: ∂ ∂P (t, f ) =−s [f (1 − f )P (t, f )] ∂t ∂f Dg ∂ 2 [f (1 − f )P (t, f )] , + 2 ∂f 2

(17)

H(∞) =

2F (∞)[1 − F (∞)] 1+
Dg 2(µ12 +µ21 )

.

(14)

where s = sτg . The equations for F and H are not ˜ −1 as useful as before because the hierarchy of the moment

7 equations does not close. Nevertheless, Eq. (17) can be easily analyzed in two limits. When the population size is large (Dg ≪ s), ﬂuctuations are not important, and dF (t) ≈ sF (t)[1 − F (t)]. Upon setting F0 ≡ F (0), dt we have 1 1+
1−F0 −st F0 e

6 5 4

P (∞, f )

F (t) ≈

,

(18)

3 2 1 0 0

so the selective sweep is exponentially fast. When the ﬂuctuations dominate the dynamics, the selection slightly increases the odds of ﬁxation of the advantageous allele, but does not signiﬁcantly aﬀect the rate of ﬁxation. For a detailed analysis of Eq. (17) see Crow and Kimura (1970). In the continuous limit, the population genetics of a well-mixed eﬀectively zero-dimensional population with genetic drift, selection, and mutation is summarized by the following Fokker-Planck (or forward Kolmogorov) equation: ∂P (t, f ) ∂ =−s [f (1 − f )P (t, f )] ∂t ∂f ∂ − {[µ21 − (µ12 + µ21 )f ]P (t, f )} ∂f Dg ∂ 2 [f (1 − f )P (t, f )] . + 2 ∂f 2

0.2

0.4

f

0.6

0.8

1

(19)

The stationary distribution for Eq. (19) is reached exponentially fast and takes the following form (Crow and Kimura, 1970; Duty, 2000) P (∞, f ) = Ce2sf /Dg f 2µ21 /Dg −1 (1 − f )2µ12 /Dg −1 , (20) where C is the normalization constant chosen to 1 set 0 P (∞, f )df = 1. This stationary distribution is plotted in Fig. 4 for both strong and weak genetic drift. Although the formulation in terms of a Fokker-Plank equation is appropriate for nonspatial models, an alternative formulation via a stochastic diﬀerential equation can be generalized to spatial models more easily. Equation (19) is equivalent to df (t) =sf (t)[1 − f (t)] + µ21 − (µ12 + µ21 )f (t) dt + Dg f (t)[1 − f (t)]Γ(t) (Itˆ), o

FIG. 4 (Color online) The stationary distribution P (∞, f ) in the presence of selection, mutation, and genetic drift, see Eqs. (19) and (20). The blue dotted line shows P (∞, f ) for s = Dg , and µ12 = µ21 = 10Dg , which corresponds to weak genetic drift, µ12 /Dg ≫ 1. The red solid line shows P (∞, f ) for s = Dg , and µ12 = µ21 = 0.1Dg , which corresponds to strong genetic drift, µ12 /Dg ≪ 1. Note the diﬀerence in curvature between the two cases and the fact that the distribution is dominated by the central region in the weak genetic drift limit, but by the tails in the opposite limit. The transition between these two regimes occurs when the mutation rates equal Dg /2, and P (∞, f ) diverges at f = 0 and f = 1. Also note that the eﬀect of natural selection is to bias the distribution toward f = 1. As s increases, the maximum of the distribution shifts to the right for weak genetic drift, and the right tail of the distribution becomes more prominent for strong genetic drift.

(21)

Γ(t1 )Γ(t2 ) = δ(t1 − t2 ),

(22)

where Γ(t) is a white, zero mean Gaussian noise, and δ(t) is Dirac’s delta-function; to get the correct Fokker-Planck Equation (19), one must use Itˆ’s prescription to deﬁne o how Eq. (21) steps the dynamics forward in time. This

interpretation of the noise term ensures that f (t) depends only on Γ(t′ ) with t′ < t as it is appropriate for population genetics. Itˆ’s prescription is adopted throughout o the paper, and a brief introduction to the Itˆ calculus o is given in Appendix A [see also (Duty, 2000; Gardiner, 1985; Risken, 1989)]. In Sec. III, we use Eq. (21) to formulate the stepping stone model in one dimension. Well-mixed-population-models do not describe migration and subdivision of natural populations (Hartl and Clark, 1989). To remedy this deﬁciency, two common approaches exist: to assume a uniformly populated spatial habitat with free diﬀusion or to assume a patchy habitat with a prescribed pattern of limited migration between the patches. The former is the subject of this paper, and can be regarded as the continuum limit of the stepping stone model (Kimura and Weiss, 1964), see Sec. III. The simplest variant of the latter approach is known as the island model (Wright, 1931). The island model assumes that all patches or islands have the same number of organisms and populations in every patch obey well-mixed-population dynamics. The migration occurs between any two patches with equal probability, so, in some sense, this is a mean ﬁeld or inﬁnite-dimensional model. The island model

8 successfully predicts that the organisms are more likely to be related locally than globally, but most of its predictions are similar to those of well-mixed-population models because the migration does not account for spatial structure. In the limit of an inﬁnitely large number of islands, the eﬀect of migration in and out of any patch is equivalent to an eﬀective mutation rate; however, this is not the case in a one-dimensional model considered below.
III. ONE-DIMENSIONAL STEPPING STONE MODEL

∂f ∂2f =Ds 2 + sf (1 − f ) + µ21 − (µ12 + µ21 )f ∂t ∂x + Dg f (1 − f )Γ (Itˆ), o

(25)

Γ(t1 , x1 )Γ(t2 , x2 ) = δ(t1 − t2 )δ(x1 − x2 ),

(26)

In Sec. II, we formulated a model to describe genetic drift, mutation, and selection in an eﬀectively zero-dimensional habitat. For a one-dimensional population considered in this section, we extend the model to account for short range migrations during every generation. Migration is usually modeled either as exchange of individuals between neighboring island populations (demes) (Kimura and Weiss, 1964; Wright, 1931) or as dispersal of oﬀspring or adults within a continuous population (Mal´cot, 1975; Nagylaki, 1974; Wright, e 1943). Although the ﬁrst approach was developed to model patchy populations, it can be used to describe continuous populations if the deme sizes are much smaller than the whole population, and spatial variations are gradual. In this limit, both migration models should give essentially the same results. Here, we adopt the ﬁrst approach because it is conceptually simpler. To specify the one-dimensional stepping stone model, we consider an inﬁnite set of demes arranged on a line. Neighboring demes are separated by distance a and indexed by an integer l = −∞, ..., −1, 0, 1, ..., ∞. Each deme has N organisms (but the total population size is inﬁnite), and the frequency of allele one in deme l is fl (t). Migration occurs only between nearest neighbors, and, every generation, a deme exchanges mN/2 individuals ˜ with its right neighbor and mN/2 individuals with its ˜ left neighbor. We assume that the exchange fraction m ˜ is much smaller than one, and that the individuals of both allelic types are equally likely to be exchanged. Thus, in one generation, fl changes by m(fl−1 + fl+1 − 2fl )/2 ˜ due to migration. The variance of fl grows due to randomness in the exchange process, but this increase is negligible compared to the genetic drift within an island. In the continuous time limit, fl (t) obeys the following generalization of Eq. (21): dfl m = (fl−1 + fl+1 − 2fl ) + sfl (1 − fl ) + µ21 dt 2 − (µ12 + µ21 )fl + Dg fl (1 − fl )Γl (Itˆ), o (24)

where the spatial and genetic diﬀusion constants are Ds = ma2 /2 and Dg = aDg = (2a)/(τg N ) respectively. Thus the continuous time and space limit the stepping stone model is described by the stochastic Fisher-Kolmogorov-Petrovsky-Piscounov equation (Fisher, 1937; Kolmogorov et al., 1937) with additional terms describing mutation. Similar to the analysis of the well-mixed-population model discussed in Sec. II, we use equal-time correlation functions of f (t, x) to characterize the dynamics of the stepping stone model. The spatial versions of the average frequency and heterozygosity are deﬁned as follows: F (t, x) = f (t, x) , (27)

H(t, x1 , x2 ) = f (t, x1 )[1−f (t, x2 )] + f (t, x2 )[1−f (t, x1 )] . (28) The equation of motion for F (t, x) depends on H(t, x, x), and is readily derived by averaging Eq. (25), which gives ∂F ∂2F s = Ds 2 + µ21 − (µ12 + µ21 )F + H(t, x, x). (29) ∂t ∂x 2 The dynamics of H(t, x1 , x2 ) is obtained by diﬀerentiating Eq. (28) with respect to t and then eliminating ∂f ∂t with the help of Eq. (25). Note, Itˆ’s formula (see Apo pendix A) must be used to diﬀerentiate Eq. (28) correctly. The result is ∂ ∂2 ∂2 + 2 H(t, x1 , x2 ) H(t, x1 , x2 ) =Ds ∂t ∂x2 ∂x2 1 − Dg H(t, x1 , x2 )δ(x1 − x2 ) − 2(µ12 + µ21 )H(t, x1 , x2 ) + 2µ21 + (µ12 − µ21 )[F (t, x1 ) + F (t, x2 )] s + [H(t, x1 , x1 ) + H(t, x2 , x2 )] 2 − s 2f (t, x1 )[1 − f (t, x1 )]f (t, x2 ) − s 2f (t, x2 )[1 − f (t, x2 )]f (t, x1 ) . (30) Equations (29) and (30) agree with the ones derived in Nagylaki (1978) in the limit of no mutations considered there.

(23)

where m = and δl1 l2 is Kronecker’s delta. We can also write Eq. (23) in the continuous space limit by introducing a spatial coordinate x = la,

Γl1 (t1 )Γl2 (t2 ) = δl1 l2 δ(t1 − t2 ), mτg ˜ −1

9 From Eq. (30), one can see that the hierarchy of the moment equations does not close unless selection is absent. Similar to the well-mixed case, the correlation functions for neutral models with and without mutations can be found analytically, see Secs. V and IV, but diﬀerent methods are required to analyze the dynamics in the presence of selection, see Sec. VI. To simplify the analysis, we consider well-mixed, spatially homogeneous initial conditions. Then F is only a function of t, and H is a function of t and x = x1 − x2 . With these simplifying assumptions, the equations of motion for F (t) and H(t, x) take the following form: dF (t) s = µ21 − (µ12 + µ21 )F (t) + H(t, 0), dt 2 alleles, Eq. (33) is valid for an arbitrary number of spatially diﬀusing neutral alleles. See Appendix E for a more detailed discussion of the q-allele problem, with q ≥ 3. To better understand the microbiology experiments on neutral alleles by Hallatschek et al. (2007), we consider uncorrelated initial conditions F (0) = F0 and H(0, x) = H0 , where F0 is the fraction of allele one and H0 = 2F0 (1 − F0 ), which is the heterozygosity of a well-mixed population with the frequency of allele one equal to F0 . For these initial conditions, Eq. (33) is solved in Appendix B. The results are
t − 8D
x2 s (t−t′ )

(31)

H(t, x) =H0 − Dg

dt′ e

0

H(t′ , 0) , 8πDs (t − t′ ) (34)

∂2 ∂ H(t, x) =2Ds 2 H(t, x) − Dg H(t, 0)δ(x) ∂t ∂x − 2(µ12 + µ21 )H(t, x) + 2µ21 + 2(µ12 − µ21 )F (t, x) + sH(t, 0) − 2s 2f (t, 0)[1 − f (t, 0)]f (t, x) . (32)

IV. NEUTRAL MODEL WITHOUT MUTATIONS

where erfc(y) is the complementary error function. The spatial heterozygosity at vanishing separation, H(t, 0), is particularly interesting because it indicates the degree of spatial segregation: if H(t, 0) ≪ 1, then, locally, the demes are ﬁxed to one of the two alle2 les. From Eq. (35), one can see that, for t ≫ 8Ds /Dg ,
2 πDg t 8Ds −1/2

H(t, 0) = H0 erfc 



 2 2 Dg t Dg t  e 8Ds , 8Ds

(35)

We start the analysis of the one-dimensional stepping stone model by considering neutral alleles that do not mutate. In practice, this means N 2 µ12 , N 2 µ21 ≪ 1 ˜ ˜ and N 2 s ≪ 1 (as we show below). Although these as˜ sumptions are not always realistic, they help to clarify the role of genetic drift in a spatial context. In addition, neglecting mutations is a good approximation on time scales shorter than the waiting times for the mutations µ−1 and µ−1 . Under these assumptions, F does 12 21 not change, F (t) = F0 , and Eq. (32) reads ∂2 ∂ H(t, x) = 2Ds 2 H(t, x) − Dg H(t, 0)δ(x). ∂t ∂x

H(t, 0) = H0

+ O(t−3/2 ),

(36)

(33)

Equation (33) can also be derived by tracing the ancestral lineages of organisms backward in time. The average spatial heterozygosity H(t, x) is the average probability of sampling two diﬀerent individuals chosen at time t from demes separated by distance x. As we trace the lineages of the two sampled organisms backward in time, the lineages diﬀuse in space due to migration and, when they are at the same point, they have a chance to coalesce, in which case the sampled organisms must be identical because they have a common ancestor. Such a coalescence event changes the probability of being diﬀerent from H(t, 0) to 0, and acts like a sink at x = 0. The ﬁrst term on the right hand side of Eq. (33) describes the diﬀusion, and the second term describes the coalescence. Since this argument is valid for an arbitrary number of

which means that at long times one of the alleles reaches ﬁxation locally. Therefore we see that the spatial model we are considering is consistent with the experiments by Hallatschek et al. (2007) (see Fig. 1) because it predicts the formation of domains (regions of local ﬁxation). Thus, similar to the well-mixed model considered in Sec. II, one of the alleles reaches ﬁxation lo2 cally with ﬁxation time τf = 8Ds /(πDg ) ∼ N 2 . But not only is this ﬁxation time proportional to N 2 , instead of N , the functional form of heterozygosity decay is different: instead of a rapid exponential decay, the spatial model shows a slow algebraic decay of local heterozygosity. These results agree with the previous works on population genetics by Mal´cot (Mal´cot, 1975) and Nagye e laki (Nagylaki, 1974). Local ﬁxation and t−1/2 decay of local heterozygosity have also been found in the voter model (Cox and Griﬀeath, 1986), which corresponds to the stepping stone model with N = 1 (see Appendix F). The characteristic demixing time can also be estimated by the following scaling argument: the characteristic population size at time t in the coarsening process √ is Nch (t) ∼ n0 tDs , where the population density n0 ∼ N/a. Upon recalling that the ﬁxation time in zero dimensions is τf ∼ Nch (τf )τg , and solving self-consistently 2 2 for τf , we have τf ∼ Ds N 2 τg /a2 ∼ Ds /Dg ∼ N 2 τg .

10 walks, and the average domain size can be easily calculated (Hallatschek and Nelson, 2010). Equations (35) and (38) suggest that the processes driven by the genetic drift slow down with time because the logarithmic time derivatives of H(t, 0) and ℓ tend to zero as time goes to inﬁnity. In the annihilating random walk picture of Hallatschek and Nelson (2010), annihilations become rarer and rarer as the coarsening progresses. A more direct measure of genetic drift, which is also interesting from the biological point of view, is the ﬂuctuations of the total fraction of, say, the ﬁrst allele f(t) in a ﬁnite population of length L. We deﬁne f(t) as
FIG. 5 (Color online) Solutions of Eq. (33) at various times, given H(0, x) = 1 . Time and distance are measured in units 2 such that Ds = 1 and Dg = 1. Time increases from the top curves to the bottom curves. Note the statistical reﬂection symmetry, H(t, x) = H(t, −x).

f(t) =

1 L

L

f (t, x)dx,
0

(39)

Another important characteristic of H(t, x) is the length scale over which H(t, x) changes from its minimum to its maximum values. Figure 5 plots Eq. (34) and shows the spatial variation of H(t, x) at diﬀerent times. One can see that the spatial heterozygosity is reduced near the origin due to the local ﬁxation, but H(t, x) rises to its initial value H0 at large x, where the alleles remain un2 correlated. After the domains form, i.e. for t ≫ 8Ds /Dg , this change from H(t, 0) to H(0, x) happens on a length scale that is set by the average size of the domains ℓ, √ which is proportional to the diﬀusion length 2Ds t, as follows from Eq. (34). Since this characteristic length scale changes with time, it is convenient to rescale dis√ tances: x = x/ 2Ds t. Upon using Eq. (36) to simplify ¯ Eq. (34), we see that H(t, x) approaches a nontrivial limit in terms of x as time goes to inﬁnity: ¯ 1 π
1 0

and compute its variance ν(t) to characterize its ﬂuctuations. Upon integrating Eq. (25)over x with s = µ12 = µ21 = 0, we obtain the equation of motion for f: df 1 = dt L
L 0

Dg f (t, x)[1 − f (t, x)]Γ(t, x)dx,

(40)

where the spatial diﬀusion term vanishes after integration by parts provided periodic or Newman boundary conditions are imposed. Upon noting that f = F = const (the Itˆ interpretation of the noise Γ(t, x) is crucial o here) and deﬁning ν = f2 − f 2 , (41)

H(t, x) −→ H0 ¯

t→∞

1−

dζ ζ(1 − ζ)

e−

|¯|2 x 4ζ

,

(37)

d one ﬁnds immediately that dν = dt f2 . To evaluate dt the time derivative, we use the rules of the Itˆ calculus, o sketched in Appendix A, and ﬁnd

which agrees with the known results for the voter model (Cox and Griﬀeath, 1986). A more precise evaluation of the domain density and hence an average domain size ℓ(t) can be obtained from H(t, x), as shown in Appendix C. From Eq. (C4), we know that ℓ(t) = Dg4Ds , so using Eq. (36) we see H(t,0) that √ 2πDs t , 2f0 (1 − f0 )

dν(t) 1 = 2 dt L ×

L 0 0

L

Dg f (t, x1 )[1 − f (t, x1 )]

Dg f (t, x2 )[1 − f (t, x2 )]δ(x1 − x2 )dx1 dx2 , (42)

where the delta function comes from averaging over the noise and using Eq. (26). From Eq. 42, it follows that Dg 2L
t 0

ℓ=

(38)

ν(t) =

H(t′ , 0)dt′ ,

(43)

which is consistent with the analysis of Hallatschek and Nelson (2010). Note that the genetic diﬀusion constant Dg ∼ 1/N drops out because, at large times, the only dynamics left is the diﬀusive motion of the domain walls. With neutral alleles, these boundaries behave as annihilating random

where we assume ν(0) = 0. Hence, we know ν(t) exactly because H(t, 0) is given by Eq. (35). 2 For small times, t ≪ 8Ds /Dg , the variance grows linearly with time. For large times, we can use the asymptotic expansion of H(t, 0) given by Eq. (36) to calculate ν(t). The result is

11 the two alleles. We assume, as before, statistically homogeneous initial conditions and note that the dynamics of the one and two-point correlation functions is then given by dF (t) = µ21 − (µ12 + µ21 )F (t), dt ∂ ∂2 H(t, x) =2Ds 2 H(t, x) − Dg H(t, x)δ(x) ∂t ∂x − 2(µ12 + µ21 )H(t, x) + 2µ21 + 2(µ12 − µ21 )F (t), (46)

√ H0 8Ds t ν(t) = √ +O πL

Ds Dg L

.

(44)

Equation (44) is consistent with Bramson and Lebowitz (1991), and we immediately conclude that the standard deviation ∆(t) = ν(t) grows as t1/4 for large times! This important result is generalized for the ﬂat-front and undulating-front models with q-alleles in Appendix E by approximating the dynamics of the domain boundaries by annihilating random walks. Thus, f performs a subdiﬀusive random walk, and genetic drift of the global frequency f(t) becomes weaker with time. Equation (44) is L2 valid only for t ≪ Ds because it relies on Eq. (36), which is valid for an inﬁnite population, and should break down at times that are long enough for a domain boundary to diﬀuse from one end of the population to the other. Using Equations (8) and (43), one can also calculate the behavL ior of the global heterozygosity H(t) = L−1 0 H(t, x)dx, i.e. the probability to sample two diﬀerent alleles from the population regardless of their spatial locations: Dg t H(t′ , 0)dt′ L 0 √ Ds H0 32Ds t +O = √ πL Dg L
2 8Ds /Dg

(47)

2F0 (1 − F0 ) − H(t) =

(45) ,

where F (t) ≡ f (t, x) is independent of x. The equation of motion for F in the spatial model is exactly the same as Eq. (11), which describes the well-mixedpopulation model. Therefore F relaxes to its equilibrium µ21 value, F (∞) = µ12 +µ21 , [see Eq. (13)] exponentially fast with time constant (µ12 + µ21 )−1 ≫ τg . The similarity to the nonspatial model is not surprising because neutral mutations occur equally likely at any point within the population, regardless of its spatial structure. The dynamics of H(t, x) is, however, more complicated because both mutation and genetic drift determine the behavior of the spatial heterozygosity. The stationary solution of Eq. (47) reads 2µ12 µ21  e 1− (µ12 + µ21 )2 1+4 
−
µ12 +µ21 Ds

where the second equality requires ≫ t ≪ L2 Ds for the reasons mentioned above. In the opposite
L limit t ≫ Ds , the global heterozygosity H(t, x) obeys zero-dimensional dynamics of a well-mixed population with an eﬀective Dg = Dg /L as shown in Nagylaki (1974). The local heterozygosity and average domain size can be obtained from experiments on microbial spreading like the one shown in Fig. 1. If the data are suﬃciently precise, Eqs. (36) and (38) could be used to extract Ds and Dg from the experiments. Since Dg ∼ 1/N , extracting Dg from experimental data determines the eﬀective deme size for the equivalent stepping stone model. Ds can be obtained from the diﬀusion of individual domain boundaries or ν(t). These two parameters completely determine the neutral dynamics without mutation and play an important role when selection or mutation is present.
2

H(∞, x) =

|x|

Ds (µ12 +µ21 ) 2 Dg



proaches 2F (∞)[1 − F (∞)]. Thus mutations cause the frequencies of allele one to eventually become uncorrelated at large separations. At shorter distances, however, there are correlations, and H(∞, x) < H(∞, ∞) for all x < ∞. Note, in particular, that H(∞, 0) = 2F (∞)[1 − F (∞)] < H(∞, ∞). Dg 1+ 1√ 4
Ds (µ12 +µ21 )

(48) Equation (48) agrees with the solution by Kimura and Weiss (Kimura and Weiss, 1964), which was obtained in the discrete space and time limit. One can see Ds that, for x ≫ µ12 +µ21 , the spatial heterozygosity ap-

.

(49)

V. NEUTRAL MODEL WITH MUTATIONS

While on short time scales mutation can be neglected, it is the long time scales and the patterns of genetic diversity created by mutations that are of particular interest in population genetics. Noticeable mutations also arise in microbiology experiments like those in Fig. 1, especially if mutation rates are enhanced by DNA damaging chemicals or radiation. In this section, we extend the results of Sec. IV by allowing for nonzero mutation rates between

Note also that, for small mutation rates, the heterozygosity is proportional to µ12 + µ21 [see Eq. (14)] in a wellmixed population, but the local heterozygosity in a one√ dimensional population is proportional to µ12 + µ21 whenever τg (µ12 + µ21 ) ≪ a2 /(N 2 Ds τg ), which is a consequence of weaker genetic drift in one dimension. 3

3

In population genetics, population structure and spatial correlations are often reported via Fst = H(∞, 0)/H(∞, ∞), which can readily be obtained from Eq. (49).

12 When H(∞, 0) ≪ H(∞, ∞), the population is segregated into domains of diﬀerent allelic types. Upon invoking Eq. (C4), we obtain the following average domain size: 2Ds (µ12 + µ21 )2 Dg µ12 µ21
1 3

2f (t, x1 )[1 − f (t, x1 )]f (t, x2 ) ≈ 2F (t)[1 − F (t)][1 − 2F (t)] − [1 − 2F (t)]H(t, x) + H(t, 0)F (t), and

(54)

ℓ=

1+

1 4

Dg Ds (µ12 + µ21 )

(50) f (t, x1 )f (t, x2 )f (t, x3 ) ≈

2 Ds (µ12 + µ21 ) 2 . ≈ 2µ12 µ21

This result, together with Eq. (13), can be used to extract the mutation rates from experimental data. We can also determine how fast H(t, x) reaches its stationary value. Since the heterozygosity cannot be in equilibrium unless the frequency of the alleles has equilibrated, we assume, for simplicity, that F (0) equals its stationary value. Then, the deviation of the spatial heterozygosity from its long time equilibrium ˜ value H(t, x) = H(t, x) − H(∞, x) obeys the following equation: ∂2 ˜ ∂ ˜ ˜ ˜ H = 2Ds 2 H − Dg Hδ(x) − 2(µ12 + µ21 )H. (51) ∂t ∂x Equation (51) can be further simpliﬁed by the change ˜ ˆ of variables H = e−2(µ12 +µ21 )t H, which leads to ∂2 ˆ ∂ ˆ ˆ H = 2Ds 2 H − Dg Hδ(x). ∂t ∂x

f (t, x1 )f (t, x2 ) f (t, x2 )f (t, x3 ) , F assuming x1 ≤ x2 ≤ x3 . (55)

(52)

Since Eq. (52) is identical to Eq. (33), we conclude that, at long times, the diﬀerence between H(t, x) and the staˆ D ˆ tionary solution decays as C D2st e−2(µ12 +µ21 )t , where C
g

The ﬁrst scheme is a simple factorization approximation; the second scheme, which assumes small ﬂuctuations, is due to Nagylaki (1978); and the third scheme, which provides a good approximation for some diﬀusion limited reactions, was proposed in Lin (1991). Unfortunately, none of the schemes describe the behavior of the system correctly. The progress can be made, however, for some initial conditions in two limiting cases of strong selection sD2s ≫ 1 and weak selection sD2s ≪ 1 Note that we Dg Dg now use the term weak selection in a diﬀerent sense than in Sec. II. For the rest of this section, we include spatial diﬀusion and genetic drift, but neglect mutations, which is justiﬁed on short time scales. First, let us consider the initial condition f (0, x) = 1 − θ(x), where θ(x) is the Heaviside step function. This initial condition speciﬁes just one domain boundary, which, for any positive s, undergoes Brownian motion with a drift to the right. This is a good description of an expansion of a new advantageous mutant spreading through the population. In the strong selection 2 limit (s ≫ Dg /Ds ), Fisher (1937) found that the sharp boundary above broadens to a width of order Ds /s, and the velocity of the genetic wave is given by vs = 2 sDs . (56)

is a constant. Thus, apart from an algebraic prefactor (and a nontrivial spatial dependence), the dynamics of H(t, x) is essentially the same as in the well-mixed case. In this section, we considered a model with only two alleles; however, in many circumstances, an inﬁnite alleles model is more appropriate. The inﬁnite alleles model is brieﬂy discussed in Appendix D. Some results forqalleles, 2 < q < ∞, are discussed in Appendix E.
VI. SELECTION

When, in contrast, selection is weak compared to genetic drift, it was recently found that the velocity is given by (Doering et al., 2003; Hallatschek and K. S. Korolev, 2009) 2sDs . Dg

vw =

(57)

Unlike the neutral models with spatial diﬀusion and mutation discussed above, the one-dimensional stepping stone model with selection are diﬃcult to treat analytically because the hierarchy of moment equations does not close. We brieﬂy examined three closure schemes: 2f (t, x1 )[1 − f (t, x1 )]f (t, x2 ) ≈ H(t, 0)F (t), (53)

When the population contains multiple domains, the domain walls bordering a favorable genetic variant (“allele one”) expand to engulf the regions occupied by the more deleterious allele. Another interesting initial condition is f (0, x) = F0 = const., i.e. the population is initially well-mixed. This scenario, for example, describes the quasi-onedimensional strip of pioneers advancing at the front of a two-dimensional population wave that originated from

13 a well-mixed ancestral population and is propagating in the region where one of the alleles has higher ﬁtness (see Hallatschek and Nelson (2010)). If selection is strong enough, then allele one (“the preferred variant”) takes over the population before spatial correlations have time to appear. To see this, note from Eq. (18) that allele one wins locally on the time scales of s−1 , but, from Eq. (35), the time for spatial correla2 tions to appear is on the order of Ds /Dg , which is much larger than s−1 in the strong selection limit. Thus the behaviors of one-dimensional and well-mixed populations 2 2 are similar when s ≫ Dg /Ds = a2 /(τg N 2 Ds ). In the limit of weak selection, however, spatial correlations appear before allele two is eliminated. Qualitatively, we can divide the selective sweep into two stages. During the ﬁrst stage, the eﬀects of selection are negligible and spatial segregation occurs as described in Sec. IV. During the second stage, the domains of allele two shrink at each end with wall velocity vw given by Eq. (57), and the stochastic motion of the domain boundaries can be neglected. The crossover time between the stages occurs when the diﬀusive displacement of the walls is of the same order as their deterministic displacement, i.e. √ when vw t = 2Ds t. Thus the crossover time τ ∗ is on the 2 2 order of Ds /vw , which can be expressed as Dg /(s2 Ds ) with the help of Eq. (57). Then, from Eq. (38), the average domain size at the crossover ℓ∗ is on the order of Dg /(sH0 ). The dynamics during the second stage depends on the probability distribution of domains of size η, Pd (η) at time τ ∗ . For annihilating random walkers, which are a good approximation to domain boundaries during the ﬁrst stage, Bramson and Griﬀeath (1980) proved that Pd (η) has exponential tail for large η of the ′ ∗ form e−γ η/ℓ , where γ ′ is a number of the order 1. Since each domain shrinks with velocity 2vw , the fraction of allele two can be expressed as
∞ 0 − λsDDs t 2
g 2

Finally, we address a slightly diﬀerent, but equally important, question: What is the probability psurv that a few copies of the advantageous allele survive and establish a growing domain? Doering et al. (2003) solved this problem exactly: 2s Dg
+∞

psurv = 1 − exp −

f (0, x)dx .
−∞

(59)

Surprisingly, the survival probability does not depend on the diﬀusion constant. For a small initial number of advantageous alleles, we can qualitatively explain this result in the limit of weak selection by the following argument. Initially, the dynamics is almost neutral, and the probability of survival within a small interval of length ∆x is proportional to the relative fraction of the advantageous allele in this interval, f (0, x)dx/∆x, because every organism has approximately the same probability to reach ﬁxation. Once a domain of size ∆x is formed, its survival probability equals the probability that the two biased random walks performed by the domain boundaries never meet, which is 1 − exp(−v∆x/Ds ) ≈ v∆x/Ds , for small ∆x and s, see (Hallatschek and Nelson, 2010; Redner, 2001). Then, using Eq. (57), the survival probability is v∆x/Ds f (0, x)dx/∆x = 2s f (0, x)dx/Dg ≈ 1 − exp −2s f (0, x)dx Dg ]. Note that, even though psurv does not depend on Ds , the expression for the survival in one dimension does not reduce to its analog in well-mixed populations (Crow and Kimura, 1970; Doering et al., 2003) 1 − e−2sf (0)/Dg , 1 − e−2s/Dg

psurv =

(60)

1 − F (t) ∝

ηe−

γ(η+2vw t) ℓ∗

dη ∝ e

,

(58)

where λ is a number of order 1. From Eq. (58) it follows that, as in the well-mixed case, the selective sweep is exponentially fast, but the time constant of this process is proportional to s−2 rather than s−1 . The analysis leading to Eq. (58) can be generalized to an arbitrary initial probability distribution, Pd (η), provided the dynamics is dominated by selection and genetic drift. For example, if a population initially in equilibrium with respect to mutations and genetic drift (see Sec. V) is aﬀected by an abrupt environmental change that makes allele one advantageous, then the shift to the new equilibrium occurs exponentially fast with a time constant 1/2 sDs µ12 µ21 proportional to Dg (µ12 +µ21 )3/2 , assuming Pd (η) has exponential tail of the form e−γ η/ℓ , where γ ′′ is a constant of order unity, and ℓ is given by Eq. (50).
′′

unless one assumes s/Dg ≫ 1. We make two important observations based on the results of this section. First, the temporal dynamics of the one-dimensional stepping stone model with selection can depend strongly on the initial conditions. Second, the results in the weak selection limit are sometimes related to the results in the strong selection limit or in well-mixed-population models by a parameter substitu2 tion, e.g. s → s2 Ds /Dg , at least up to a numerical factor. The second observation suggests that, while data can be naively ﬁtted to a well-mixed-population model, the ﬁt in fact gives the “renormalized” value of s instead of the “bare” one. Although the one-dimensional stepping stone model provides a reasonable approximation to neutral genetic demixing at a linear front of an expanding twodimensional microbial colony (see Fig.3), there are additional subtleties associated with the dimensional reduction from two to one dimensions when one of the alleles is more ﬁt. Apart form the undulations of the front mentioned in the caption of Fig. 2, the front develops additional structure because favorable sectors bulge out

14

Number of boundaries

ahead of their less ﬁt neighbors. In the limit of very small genetic drift, there are kink singularities where favorable and unfavorable domains meet. Nevertheless, the basic picture of domain boundaries engulﬁng unfavorable sectors at a constant velocity is still valid. See Hallatschek and Nelson (2010) for further details.
VII. SIMULATIONS

100
(d(t)-d(0)) 2

100 10 1 0.1
pe =4 /3

slo

t

-2

/3

10

1

10

100

t (generations)

In Secs. III, IV, V, and VI, we reviewed and extended the theoretical analysis of the one-dimensional stepping stone model. This model, while of great theoretical interest, relies on a restrictive set of assumptions including large deme sizes and slow diﬀusive migration. The recent experiments by Hallatschek et al. (2007), on the other hand, were carried out with bacterial fronts that were only a monolayer thick; therefore, demes consisted of only a few microbes. Moreover, depending on the microorganism, nearby demes can exchange a signiﬁcant fraction of cells in each generation. In this section, we discuss numerical simulations not subject to these restrictions and compare them with the theoretical predictions. We simulate L organisms arranged on a line and labeled by an integer l, l = 1, 2, ..., L. Each organism can be either of allelic type 1 or allelic type 2. During even generations, the oﬀspring at site l comes from an organism at either site l − 1 or site l, whereas, during odd generations it comes from either site l or l + 1. The simulations embody the process illustrated in Fig. 2, laid out on a triangular lattice in space and time. Periodic boundary conditions are imposed at the left and right ends of the population. Let 12 → 2 refer to the event that the oﬀspring has allelic type 2, while one of its possible parents has allelic type 1 and the other has allelic type 2, etc. The transition probabilities, which depend on the states of the possible ancestors, are then given by 11 → 1 : 11 → 2 : 22 → 1 : 22 → 2 : 12 → 1 : 12 → 2 : 21 → 1 : 21 → 2 : 1 − µ12 , ˜ µ12 , ˜ µ21 , ˜ 1 − µ21 , ˜ 1 2 1 2 1 2 1 2 s ˜ 4 s ˜ + 4 s ˜ + 4 s ˜ + 4 +

1

0

50 t (generations)

100

FIG. 6 (Color online) The number of monoallelic domain boundaries as a function of time in the undulatingfront model. The simulation of 100 demes averaged over 100 runs is plotted in blue (dots), and the theoretically expected decay of the number of boundaries as t−2/3 [see Saito and M¨ller-Krumbhaar (1995)] is plotted u in red (solid line). Note that the agreement between the theory and the simulations is not expected during the transitory regime at early times. At t = 0, each site is assigned either allele one or allele two with equal probability. Inset: Log-log plot of the mean square displacement of the domain boundaries as a function of time in the same set of simulations as the main plot. The blue dots are the simulation data, and the solid red line is the expected slope according to Saito and M¨ller-Krumbhaar (1995). u

1 2 1 s ˜ µ12 + ˜ − 2 4 1 (1 − µ12 ) + ˜ 2 1 s ˜ µ12 + ˜ − 2 4 (1 − µ12 ) + ˜

−

s ˜ µ21 , ˜ 4

(61)

(1 − µ21 ), ˜ − s ˜ µ21 , ˜ 4

(1 − µ21 ). ˜

The event 12 → 2 can happen either if allele one was selected for reproduction (probability 1/2 + s/4) and mu˜ tated (probability µ12 ), or if allele two was selected for ˜ reproduction (probability 1/2 − s/4) and did not mu˜ tate (probability 1 − µ21 ). Other transition probabilities ˜

are obtained analogously. Thus, the system we simulate is very similar to the voter model (Cox and Griﬀeath, 1986), which equivalent to population genetics models with N = 1 (see Appendix F). However, to make the calculation faster, we use discrete generations rather than exponentially distributed waiting times until reproduction. We found no signiﬁcant diﬀerences in dynamics between the voter model and the model used here. First, we simulate the neutral model without mutations. To illustrate the similarities and diﬀerences between the stepping stone model and the undulating-front model, we also simulate a linear population wave in a two-dimensional habitat; both models are displayed in Fig. 2. Our model with an undulating front is the same as in Saito and M¨ller-Krumbhaar (1995), but we use a u triangular grid instead of a square one. Figures 6 and 7 show how the average number of domain boundaries decreases with time; the insets show the mean square displacements of the boundaries as a function of time. The simulations conﬁrm that the domain move diﬀusively for the one-dimensional model with a ﬂat front and superdiffusively for the undulating front (Eden) model, in agreement with Saito and M¨ ller-Krumbhaar (1995). Thereu fore, the number of domain boundaries decays faster for the undulating-front model compared to the ﬂat-front model, which, as we show below, most closely tracks the prediction of the one-dimensional stepping stone model.

15

Number of boundaries

100
(d(t)-d(0)) 2

50

10

slo

pe

=1

1

10

t

-1/ 2

0.1

1

10

100

t (generations)

0

50 t (generations)

100

FIG. 7 (Color online) The number of monoallelic domain boundaries as a function of time in the ﬂat-front model. The simulation of 100 demes averaged over 100 runs is plotted in blue (dots), and the theoretically expected decay of the number of boundaries as t−1/2 (see Hallatschek and Nelson (2010)) is plotted in red (solid line). Note that the agreement between the theory and the simulations is not expected during the transitory regime at early times. At t = 0, each site is assigned either allele one or allele two with equal probability. Inset: Log-log plot of the mean square displacement of the domain boundaries as a function of time in the same set of simulations as the main plot. The blue dots are the simulation data, and the solid red line is the expected slope according to Hallatschek and Nelson (2010), Bramson and Lebowitz (1991), and Eq. (38).

FIG. 8 (Color online) A single run of the ﬂat-front model with 1000 organisms. At t = 0, each site is assigned either allele one or allele two with equal probability. The spatially averaged heterozygosity [deﬁned in the sense of Eq. (C1) but without taking the limit L → ∞] for three times measured in generations: t = 0 shown in cyan, t = 1024 shown in red, and t = 4096 shown in dark blue. The separation l is the shortest distance between two points around the cylinder, and we take the clockwise direction to be positive. At inoculation, the heterozygosity ﬂuctuates around 1/2 since the two alleles have equal probabilities of occupying any site. The only exception is the site l = 0, where the heterozygosity is zero automatically because we only allow a single microorganism per site. After 1024 generations, short range correlations are clearly visible, and, after 4096 generations, one can relate the abrupt changes in the slope of H(t, l) to the sizes of the sectors in the population (not shown). The wiggles are eliminated when averaged over many realizations, as shown in Fig. 9. Note that the curve for t = 4096τg lies completely below 1/2 because, at this time, the relative fraction of the alleles deviates signiﬁcantly from the initial 50 : 50 ratio due to genetic drift.

A single run of a simulation is shown in Fig. 8, and the spatial heterozygosity averaged over many realizations is shown in Fig. 9. Figure 10 shows that H(t, x) for large t is described well by the limiting shape of spatial heterozygosity given by Eq. (37). To properly represent H(t, 0), we artiﬁcially deﬁne demes of a larger size by grouping M neighboring individuals into one deme. From a theoretical point of view, this procedure is similar to the formation of Kadanoﬀ block spins as in renormalization group methods (Goldenfeld, 1992; Wilson and Kogut, 1974) whereas, from the point of view of population genetics, this procedure is similar to the methods of collecting data from a dispersed natural population. In ﬁeld studies, scientists do not typically sample every single individual; instead, they often divide the habitat into patches and sample a representative number of individuals from those patches. To summarize, we keep the dynamics of the simulation exactly the same, but deﬁne the spatial heterozygosity on the demes of size M rather than on single individuals, see Fig. 9. We found that the local heterozygosity has the form H(t, 0) = β(M )t−1/2 , as predicted by our analysis of the stepping stone model, for all M studied (1 ≤ M ≤ 64). From Eq. (36), we expect that β ∝ M −2 ; this expectation is also conﬁrmed by our simulations. As discussed in Sec. IV, the total fraction of allele one f(t) ﬂuctuates in an unusual way with time. Fig-

FIG. 9 (Color online) The eﬀects of coarse-graining on the time evolution of the spatial heterozygosity H(t, l) averaged over 100 realizations of the ﬂat-front model with a thousand organisms. At t = 0, each site is assigned either allele one or allele two with equal probability. (a) Each deme hosts only one organism. Consequently, the heterozygosity at l = 0 is zero at all times. (b) The same set of simulations, but the organisms have now been grouped into demes of size 5 for the purpose of calculating H(t, l). Part b shows qualitative agreement with the solution of stepping stone model displayed in Fig. 5, as does part a outside the region around l = 0. Note that, unlike the calculation presented in Fig. 5, this simulation was run suﬃciently long to show the eﬀects of the boundary conditions.

ure 11a shows examples of these remarkable variablestep-length random walks. The ﬂuctuations of f(t) obey Eq. (44) and grow subdiﬀusively, as shown in Fig. 11b. We also ﬁnd good agreement between the theory and the simulations in the presence of mutation for all values of M studied. The stationary average heterozygos-

16

FIG. 10 (Color online) Comparison between the analytical prediction for the limiting shape of H(∞, x) and the simula¯ tions of the ﬂat-front model. The continuous curve (black) is formed by the data points representing H(t, x) for several times between t = 2·104 and t = 4·105 plotted in the rescaled coordinates x, and the circles (red) represent the theoretical ¯ prediction of the limiting shape of the average spatial heterozygosity, see Eq. (37). The data are obtained in a simulation of 3200 individuals for 4 · 105 generations with averaging over 500 realizations. At t = 0, each site is assigned either allele one or allele two with equal probability.

FIG. 12 (Color online)Equilibrium between mutation and genetic drift in the absence of selection. Comparison between the analytical prediction for the steady state heterozygosity H(∞, x) and the simulations with µ12 = µ21 = 10−4 . ˜ ˜ The black dots represent the results of the simulation, and the red circles represent the best ﬁt of theoretical result given by Eq. (48) to the data. Here, only Dg is a ﬁtting parameter; the values of Ds , µ12 , and µ21 follow from the correspondence between the discrete and continuum models. The data are obtained in a simulation of 3200 individuals for 2 · 105 generations with averaging over 100 realizations. At t = 0, each site is assigned either allele one or allele two with equal probability.

FIG. 11 (Color online) Genetic drift in a ﬁnite population. At t = 0, each site is assigned either allele one or allele two with equal probability. (a) The total fraction of allele one f(t) versus time in four single runs of the neutral model with a ﬂat front. Here, L = 1000, and there are no mutations. (b) The average standard deviation of the frequency of allele one ∆(t), shown in blue, is obtained from 200 realizations of the simulations described in a. The red solid line shows the best power-law ﬁt, and the slope is close to the exponent expected from Eq. (44). The gray area encloses the points within one standard deviation from the mean. FIG. 13 The eﬀective extinction rate α versus s in the limit ˜ of weak selection. The red dots are the data from the simulation, and the black line has the slope equal to 2, which is the expected slope from Eq. (58). The data supports α ∝ s2 .The ˜ values of α are obtained from graphs like the one shown in the inset. Inset: ln(1 − F ) versus t for s = 0.12. The green ˜ circles are the actual data points, and the blue line is the best least squares linear ﬁt. The simulation conﬁrms exponentially fast ﬁxation. The data are obtained in a simulation of 1600 individuals for 6000 generations with averaging over 100 realizations. At t = 0, each site is assigned either allele one or allele two with equal probability.

ity for M = 1 is shown in Fig. 12. Finally, we studied selective sweeps in an initially well-mixed population. Figure 13a conﬁrms the prediction from Eq. (58) that F ∝ (1 − e−αt ), and Fig. 13b conﬁrms the result of Sec. VI that, for strong genetic drift, the eﬀective extinction rate α is proportional to s2 . Numerical results for three neutral alleles are presented in Appendix E.

17 oculated with microorganisms. A radial range expansion of E. coli is illustrated in Fig. 14, which highlights the eﬀects of genetic drift at the front. The growth of a circular colony lengthens the front, thereby increasing characteristic local length scales. This inﬂation is speciﬁed by the dependence of the radius of the colony R on time t. Here, we assume R(t) = R0 + vt, which corresponds to a colony expanding with a constant velocity v from an initial radius R0 . The velocity of the expansion has been found constant in the experiments by Hallatschek et al. (2007) and in the theoretical studies of the two dimensional Fisher equation (Murray, 2003), provided the width of the front is much smaller than its length, a condition necessary for a one-dimensional model to hold. To highlight the eﬀects of inﬂation, we consider the simplest version of the one-dimensional stepping stone model without mutations and natural selection. In a circular geometry, it is convenient to use the angle ϕ = x/R(t) instead of x to reference positions along the front. Then, the equation of motion for H(t, x) takes a form analogous to Eq. (33): ∂ ∂2 Dg 2Ds H(t, ϕ)− H(t, ϕ) = H(t, 0)δ(ϕ), 2 ∂ϕ2 ∂t (R0 + vt) R0 + vt (62) where the factors of R0 + vt have been introduced to account for the inﬂation. Like Eq. (33), Eq. (62) is valid for an arbitrary number of neutral alleles, and can be understood by tracing two lineages backward in time. The time dependence of the coeﬃcients in front of the diffusion and coalescence terms accounts for the fact that, as the colony grows, the same sizes in the ϕ-space correspond to diﬀerent sizes in the x-space, where the diﬀusion and coalescence terms have their familiar, time independent form as in Eq. (33). When reexpressed in terms of t and x = (R0 + vt)ϕ, Eq. (62) contains an advection term describing the deterministic decrease of the separation between the lineages as they go back to the initial radius R0 . Equation (62) is deﬁned on a bounded domain ϕ ∈ [−π, π] with periodic boundary conditions. Nevertheless, we can approximate the problem well by considering an unbounded domain ϕ ∈ (−∞, ∞), provided two diametrically opposite lineages are suﬃciently unlikely to coalesce. From Eq. (62), we see that diﬀusion eﬀectively stops after a characteristic time R0 /v, so our approximation of an unbounded domain should be valid if the distance traveled by the lineages during this time is small compared to the radius of the colony: R0 Ds /v ≪ R0 or Ds /v ≪ R0 ; this corresponds to a regime with many sectors as we show later. One can also test the goodness of the approximation by evaluating H(0, π)−H(t, π), which is expected to be small if the approximation is valid. To simplify the analysis, we make Eq. (62) dimensionless in terms of the new variables t and φ such 2 that t = tDs /Dg and ϕ = φDs /(Dg R0 ). The equation of motion for H(t, φ) then reads

FIG. 14 (Color online) Spatial segregation in an expanding circular bacterial colony of E. coli. Diﬀerent colors label diﬀerent alleles. The Petri dish was inoculated with a well-mixed population occupying the circled region of the colony, leading to many small domains in the central “homeland.” As this population expands (shown with arrows), it segregates into well deﬁned monoallelic domains, which coalesce at early times, but seem to stop coalescing in the ﬁnal stages of the experiment. Note that the boundaries between the domains are biased to move away from each other due to inﬂation, in addition to their diﬀusive random-walk-like motion. Details of the experiment are presented in Hallatschek et al. (2007).

VIII. INFLATION

Throughout this review we have focused on the evolutionary forces acting at a linear (ﬂat or undulating) front, whose total length (averaged over the undulations) does not change in time. In this section, we explore the changes in the evolutionary dynamics caused by a constant increase of the total front length, for example, at the edge of an expanding circular colony; see Fig. 14. We now show that this increase, which we term “inﬂation,” in an analogy with cosmology (Guth, 1981), slows down genetic drift and natural selection at the front. Models of both linear and circular fronts are relevant to biology. Linear fronts describe the essential features of the dynamics when the eﬀects of curvature and changes in the front length are negligible, or when the spreading is limited by some geographical barriers, say, a receding glacier between two parallel rivers. If one focuses on genetic markers in Homo sapiens (e.g. in mitochondria DNA), dynamics of a linear front also resembles the abrupt settlement by pioneers via a “land run” in 1889 across the border of Oklahoma (Gibson, 1965). Circular fronts are more appropriate for modeling an initial colonization by a small number of pioneers arriving in the interior of a large, spatially homogeneous habitat. Semicircular fronts are relevant to colonization after landing on a coast line. The circular scenario is often realized in microbiological experiments when a Petri dish is in-

18
0.6

∂ ∂2 1 2 H(t, φ) − H(t, φ) = H(t, 0)δ(φ), ∂t (1 + σt)2 ∂φ2 1 + σt (63) 2 where the dimensionless parameter σ = vDs /(R0 Dg ) is proportional to the ratio of two characteristic time scales 2 in the problem: the local ﬁxation time τf ∼ Ds /Dg in the model of a linear front and the time in which the colony doubles its initial radius. Upon assuming φ ∈ (−∞, ∞), we obtain the exact solution of Eq. (63) for the initial condition H(0, ϕ) = H0 by a generalization of the method presented in Appendix B:
t

0.5 0.4 0.3 0.2 0.1 0 t=0 t=1 t=5 t=10 t=100 t=500 t=1000 −4 −2 0

H(t, φ)

φ

2

4

H(t, φ) =H0 −

dt′
0

H(t′ , 0) 1 + σt′

φ2 (1 + σt)(1 + σt′ ) × exp − , 8(t − t′ ) where √ H0 1 + σt π H(t, 0) =H0 − √ − arcsin 2 2πσ + H0 × erf − H0 + H0 π(1 + σt) exp 8σ 1 + σt 8σ 1 + σt 8σ 1 √ 8σ

(1 + σt)(1 + σt′ ) 8π(t − t′ )

(64)

FIG. 15 (Color online) Solutions of Eq. (63) with σ = 1 at various rescaled times t, given random initial conditions 1 on the circle bounding the homeland, H(0, φ) = 2 . Time increases from the top curves to the bottom curves. Note that there is no observable diﬀerence between H(500, φ) and H(1000, φ) because H(t, ϕ) reaches a nontrivial limitshape as t → ∞.

1 √ 1 + σt

− erf

π(1 + σt) t/8 e −1 8σ 1 + σt t/8 t ′ −t′ /8 1 e dt e arcsin √ 128πσ 1 + σt′ 0 (65)

domains, should decay as t−1 , and the shape of H(t, ϕ) should approach a nontrivial limit. From the exact solution (64), (65), we compute the average angular size of the domains ℓϕ (t) and the average number of the domains N (t) = 2π/ℓϕ (t). Similar quantities were calculated by Hallatschek and Nelson (2010) in the approximation of random walking domain boundaries, which appropriate when t > τf . Although we cannot use Eq. (C4) because of the inﬂation, Eq. (C2) remains valid and takes the following form . ∂H(t, +0) ∂ϕ
−1

ℓϕ (t) =

.

(66)

The behavior of H(t, φ) is shown in Fig. 15. Similar to a linear front, the local heterozygosity H(t, 0) vanishes for large times, and the characteristic angular length scale over which H(t, ϕ) changes from 0 to H0 increases. Thus, Eq. (62) predicts the formation and growth of the domains shown in Fig. 14. However, there are two important diﬀerences that distinguish radial expansions from linear ones. First, H(t, 0) tends to zero as t−1 rather than as t−1/2 . Second, the curve H(t, ϕ) approaches a nontrivial limit-shape, unlike the case with a linear front, where the analogous curve widens indeﬁnitely. Following Hallatschek et al. (2007), we can qualitatively understand this behavior by noticing that the diffusion and lineage coalescence eﬀectively stop after the characteristic time R0 /v. After this point, the number of domain boundaries and the angular width of the domains remain approximately constant. Therefore, H(t, 0), which is proportional to the fraction of the circumference occupied by the boundaries between the

By integrating Eq. (62) over ϕ in the neighborhood of 0, we express ∂H(t,+0) in terms of H(t, 0) and obtain ∂ϕ ℓϕ (t) = 4Ds , Dg (R0 + vt)H(t, 0) (67)

which approaches a constant at long times. This limit can be computed analytically, with the results H0 v + H0 Dg 2πH0 v + H0 Dg R0 v 2πDs
−1

ℓφ (∞) =

,

(68)

N (∞) =

2πR0 v . Ds

(69)

Equation (69) implies two things. First, by measuring N (∞) as a function of the initial homeland radius R0 ,

19 one can estimate both Ds and Dg for a microbial population, which could potentially be easier than the experiments with linear fronts that we proposed in Sec. IV. Second, if all individuals in the founding population are distinguishable (H0 = 1), then each of the ﬁnal sectors must originate from a single ancestor. Hence, N (∞) for H0 = 1 gives the average number of ancestors of the genetically segregated population at the periphery, which contains most of the organisms. This number is remarkably small. Figure 14, where H0 = 1/2, has about 20 domains, so, since N (∞) ∝ H0 , the segregated part of the population descended from only about 40 ancestors, a tiny fraction of about 20, 000 founding cells. Although some of these cells are trapped in the interior of the homeland, a large number of them are piled in a ring at the edge of the homeland within minutes of inoculation, as the carrier ﬂuid dries out (Hallatschek et al., 2007). We can further quantify the amount of genetic drift in the population by the variance ν(t) of the total fraction of allele one f(t). For simplicity, we assume a population with only two alleles. For several alleles, the global heterozygosity H(t) is more appropriate and can be easily obtained from our expressions for ν(t) because H(t) = H0 −2ν(t). We compute H(t) and thus ν(t) by integrating Eq. (62) over ϕ; the result is Dg 4π
t 0
2

10

10

1

K(σ)

10

0

10

−1

10 −4 10

−2

10

−2

σ

10

0

10

2

FIG. 16 A plot of K(σ) for Eq. (71) from a numerical solution of Eq. (63).

IX. GENETIC INFERENCE

ν(t) =

H(t′ , 0) ′ dt . R0 + vt

(70)

We are mostly interested in the long time limit ν(∞), which is approached asymptotically as t−1 . This limit can be expressed as D s H0 K(σ), 4πR0 Dg

ν(∞) =
∞

(71)

where K(σ) = 0 H ∗ (t′ , 0)/(1 + σt′ )dt′ , and H ∗ (t, φ) is the solution of Eq. (63) for H0 = 1. The dependence of K on σ is shown in Fig. 16. In the limit of large Dg (approximated by the voter model, see Appendix F), an analytical expression for ν(∞) is given by Eq. (F7). Even though inﬂation slows down lineage diﬀusion and coalescence, genetic drift can still cause large ﬂuctuations in the relative frequency of the alleles. These ﬂuctuations are particularly important for any organism that undergoes spatial colonizations followed by almost complete extinctions (such life cycles are common both in nature and in the laboratory). For such organisms, ν(∞) or H(∞) characterizes the eﬀective genetic drift, which is much larger than the one predicted by well-mixedpopulation models. We also note that the eﬀects of genetic drift could be more pronounced at a circular front because natural selection is less eﬃcient in the presence of inﬂation: Domains of deleterious alleles persist longer because the contraction of the domain due to Darwinian selection must also be able to overcome its natural expansion due to inﬂation.

So far, we have focused on forward-in-time dynamics, while trying to calculate the patterns of genetic diversity from simple models of evolutionary dynamics. However, it is often necessary to reverse the question: Given the observed genetic diversity, how do we infer the recent history of the population and estimate important parameters like the mutation rates and the eﬀective population size? This question is particularly important because the current state of genetic diversity is often the only clue to the past. Fortunately, genetic inference can be very powerful because the diﬀerences among the genomes of individuals contain valuable information about the evolution of the population, and these diﬀerences can now be easily measured via DNA sequencing. For example, genetic inference has been used to determine the time and origin of the recent expansion of Homo sapiens (Ramachandran et al., 2005) and to test whether Homo sapiens and Homo neanderthalensis used to interbreed (Nordborg, 1998). Genetic inference is a well-developed subject, which becomes rather technical when one wants to incorporate biological details and use advanced statistical tools. In this section, we address some of the basic questions in genetic inference and highlight the diﬀerences between the spatial and nonspatial models. The results for well-mixedpopulation models presented here are usually attributed to Kingman (1982); we refer the interested reader to the book by Wakeley (2008) for a lucid introduction. In a typical study, n organisms are sampled from the population, and parts of their genomes are sequenced (see Fig. 17). A genetic sequence, say . . . ACTGAA. . . , is an ordered string of letters taken from a four-letter alphabet: A, T, C, and G, where the letter stand for the nucleotides:

20 adenine, thymine, cytosine, and guanine respectively. For haploid organisms considered here, an oﬀspring inherits its sequences from the parent with possibly a few mutations (but no recombination). While a wide range of mutations is possible, we consider only point mutations, i.e. substitutions of one letter for another. Moreover, we assume that every new point mutation occurs at a new site (position) in the genome (Kimura, 1969). Because the mutation rate per site µ is very small and the total number of sites Lg (i.e. the length of the sequenced section of the genome) is large, most of the mutations occur at diﬀerent positions along the sequence, so this inﬁnite sites approximation is reasonable on time scales shorter than µ−1 . For simplicity, we neglect the dependence of the mutation rates on the position within the genome as well as on the type of substitution, i.e. all 12 possible substitutions are assumed to occur at the same rate µ. We further assume that all genetic variation is neutral (Kimura, 1983). In principle, the complete data set of n sequences of length Lg can be and often is used to estimate parameters in the model. However, we can understand the basic principles of genetic inference by considering two simple summary statistics: the average number of pairwise differences Π and the number of segregating sites S, i.e the number of sites that do not have identical nucleotides in at least two sequences in the sample. The former is intimately related to the average heterozygosity H and the latter illustrates the use of genealogical trees in genetic inference, see Fig. 17. We ﬁrst consider Π, which is deﬁned as the expected number of diﬀerent sites in two randomly selected sequences. In a ﬁnite population, any two sequences have a common ancestor, so, as we trace them backward in time, their lineages must coalesce. Let us denote the average time it takes two lineages to coalesce by T2 . Then the expected number of pairwise diﬀerences is given by Π = 2µLg T2 , (72) pled. Therefore, we introduce U2 (τ, x1 , x2 ) as the probability that two lineages have not coalesced and are at positions x1 and x2 respectively at reverse time τ . Then, U2 (τ ) is given by
L L

FIG. 17 (Color online) An illustration of the backward-intime dynamics of the ancestral lineages in a one-dimensional habitat with periodic boundary conditions. Five organisms (i.e. n = 5) are sampled from the much larger population at t = 0 and their DNA is sequenced. We do not display the sites that are identical for all organisms, which are usually the majority of the sequenced sites, i.e. only the segregating sites are shown. For illustration purposes, we also assumed that all samples diﬀer in at least one nucleotide, but in experiments one often ﬁnds organisms that have identical sequences. We trace the spatial diﬀusion and coalescence of the lineages backward in time until they merge into a single lineage of the common ancestor of the whole sample. The coalescence events are denoted by red circles, and the mutations are denoted by arrows and the resulting mutated sequences. Note that lineages may cross without coalescing as shown in the top left corner of the ﬁgure. The ancestral process shown here satisﬁes the inﬁnite sites model and illustrates the fact that the more genetically similar the lineages are the more likely they are to have a common ancestor in the recent past. Most genetic inference methods rely on this relationship as we show in this section.

where the factor of 2 accounts for the fact that mutations occur in both lineages. Since genetic inference deals with backward-in-time dynamics, it is convenient to use reverse time τ = −t. To calculate T2 , we introduce the persistence probability U2 (τ ), the probability that two lineages sampled at t = 0 have not coalesced between t = 0 and t = −τ . Because U2 (τ ) is the cumulative probability distribution function for the coalescence times, the desired probability density function is −dU2 (τ )/dτ , and T2 can be calculated from the following equation
∞ 0

U2 (τ ) =
0 0

U2 (τ, x1 , x2 )dx1 dx2 ,

(74)

T2 =

τ −

dU2 (τ ) dτ = dτ

∞ 0

U2 (τ )dτ.

(73)

For a one-dimensional population, we have to take into account the positions where the organisms are sam-

where L is the length of the habitat. The time evolution of the persistence probability and the average heterozygosity are intimately related both in well-mixed and spatial populations because both quantities describe the fate of two lineages traced backward in time. In fact, the equation of motion for H(t) is identical to that of U2 (τ ), and the same is true for H(t, x1 , x2 ) and U2 (τ, x1 , x2 ). For example, in the well-mixedpopulation model considered in Sec. II, H(t) and U2 (τ )

21 change only due to coalescence events, and each coalescent event changes both quantities from their current values to zero. Thus, analogously to Eq. (7) the equation of motion for U2 (τ ) reads d U2 (τ ) = −Dg U2 (τ ), dτ for two lineages to meet for the ﬁrst time.4 Note that T2 (x1 , x1 ) is identical to T2 in a well-mixed population, provided we take the eﬀective population size to be the total size of the spatially extended population: L/Dg = D−1 L/a = (N τg /2)L/a, where L/a is the g number of demes (a is the distance between neighboring demes). Note that the distribution of the coalescence times is highly skewed, and the average coalescence time does not characterize the distribution well: most of the time coalescence occurs very fast compared to T2 (x1 , x1 ), but in rare cases lineages persist for times much longer than T2 (x1 , x1 ) (Charlesworth et al., 2003). The average number of pairwise diﬀerences for the whole data set is obtained by averaging over the spatial positions of the samples xj , j = 1, 2, . . . , n:
n j1 −1 j2 =1

(75)

with the initial conditions U2 (0) = 1. For the onedimensional stepping stone model, we obtain that ∂ ∂2 ∂2 + 2 U2 (τ, x1 , x2 ) U2 (τ, x1 , x2 ) =Ds ∂τ ∂x2 ∂x2 (76) 1 − Dg U2 (τ, x1 , x2 )δ(x1 − x2 ) in analogy with Eq. (30) for µ12 = µ21 = s = 0. The initial condition is U2 (0, x1 , x2 ) = δ(x1 − x0 )δ(x2 − x0 ), 2 1 where x0 and x0 are the positions of the ﬁrst and second 1 2 lineages at the time of sampling. For the well-mixed case, we integrate both sides of Eq. (75) with respect to τ from zero (U2 (0) = 1) to inﬁnity (U2 (∞) = 0) and use Eq. (73) to ﬁnd that T2 = D−1 . g Then, from Eq. (72), we obtain the average number of pairwise diﬀerences 2µLg . Dg

Π1d =

4µLg n(n − 1) j

T2 (xj1 , xj2 ),

(79)

1 =2

Πwell−mixed =

(77)

The mutation rate µ can often be measured experimentally (Araten et al., 2005; Drake, 1991), so Eq. (77) and the knowledge of Πwell−mixed can be used to estimate the eﬀective population size encoded in Dg [see Eq. (4)]. For the one-dimensional stepping stone model, one has to specify the spatial boundary conditions for Eq. (76). Since a lineage can neither go outside the habitat nor disappear at its edge, reﬂecting boundary conditions should be used. With these boundary conditions, Eq. (76) has been analyzed by Wilkins and Wakeley (2002). Here, we assume periodic boundary conditions, which are appropriate for a population living on a coast line of an island. These boundary conditions are simpler because they ensure translational invariance: the average coalescence time for two lineages sampled at x0 and x0 can only 1 2 depend on |x0 − x0 |, but not on x0 and x0 separately. 1 2 1 2 Following (Wilkinson-Herbots, 1998), we solve Eq. (76) with periodic boundary conditions by the Fourier transform in the positions and the Laplace transform in reverse time. The result is |x0 − x0 |(L − |x0 − x0 |) L 2 1 2 + 1 , Dg 4Ds

where the factor n(n − 1)/2 accounts for the total number of diﬀerent ways to pair up the sequences. Given the mutation rate µ and a suﬃciently large sample size n, one can use Eqs. (79) and (78) to estimate Dg and Ds . Both parameters can be estimated because Π1d , unlike Πwell−mixed, depends on the spatial positions of the samples as well as on the properties of the population. Thus, one can generate independent equations to estimate Dg and Ds by using diﬀerent subsets of the samples; for example, Dg can be estimated from the samples taken from the same point, and Ds can be estimated from the remaining samples. Note, however, that Π1d depends only on µ/Dg and µ/Ds , so at most two parameters can be estimated from the data; similar considerations hold for the well-mixed case as well. The average number of pairwise diﬀerences is relatively easy to compute in both spatial and nonspatial models because it depends on the history of only two lineages. For the same reason, Π does not illuminate the underlying tree-like genealogy of the sample (see Fig.17), and a diﬀerent statistic is needed for that purpose. Under the inﬁnite site assumption, a given site is either monomorphic, i.e. all samples have the same nucleotide at this site, or polymorphic, i.e. two diﬀerent nucleotides are found: one is ancestral and the other is due to a mutation. Only the polymorphic sites contain information about the underlying genealogy, and the frequencies of mutations at each site are often used for genetic inference. Here, we consider a simpler summary statistic, the

T2 (x0 , x0 ) = 1 2

(78)

4

where the ﬁrst term on the right hand side is the average coalescence time for two lineage sampled at the same point, and the second term is the average time

The average time to the ﬁrst encounter of two random walks on an interval with periodic boundary conditions is equal to the average survival time of a single random walk with twice the diﬀusion constant on the same interval, but with absorbing boundary conditions. This equivalent problem can be solved by a standard method; see, e.g., (Redner, 2001).

22 number of segregating sites S, i.e. the expected number of polymorphic sites in the sample. As we go backward in time, the number of lineages decreases due to coalescence events from n to n−1, to n− 2, etc. until it eventually reaches 1; we consider only pairwise coalescence events assuming the population size is suﬃciently large so that the coalescence of more than two lineages at one time is unlikely. Let Tj be the average time when j lineages are present. Then the expected number of polymorphic sites is given by
n

Since Tj is the time during which the number of lineages changes from j to j − 1, it follows from Eq. (83) that Tj ≈ τ (j − 1) − τ (j) ≈ − dτ (j) 1 , = dj πDs j 3 (84)

S = µLg
j1 =2

jTj ,

(80)

where the factors of j account for the fact that mutations can occur in any of the j lineages during the time interval Tj . For a well-mixed population, we compute Tj by noticing that any of j(j − 1)/2 distinct pairs of lineages can coalesce next, and, from Eq. (75), each pair has a constant coalescence rate of Dg . Hence, 2 , j(j − 1)Dg

Tj = and from Eq. (80) 2µLg Dg
n−1 j1

(81)

Swell−mixed =

1 1 2µLg , ln(n) + γ − ≈ j Dg 2n =1

(82) where the approximation is valid for large sample sizes n (Gradshteyn and Ryzhik, 1980), and γ is the Euler constant. For a one-dimensional population, one could try to generalize the approach used to calculate T2 for multiple lineages, but this method seems prohibitively difﬁcult for large n. However, we can qualitatively understand the eﬀects of spatial structure by considering n lineages sampled uniformly in x from the habitat. While analyzing a related problem of annihilating random walks (see Appendix F), Doering and ben Avraham (1988), and Zhong and D. ben-Avraham (1995) showed that for a generic uniform spatial distribution of the samples the number of surviving lineages j at reverse time τ decays as 1 j(τ ) ∼ √ 2πDs τ

where τ (j) is the inverse function of j(τ ) used in Eq. (83). Equation (84) is only valid for intermediate j values: 1 ≪ j ≪ n because of the similar restrictions on Eq. (83). Upon comparing this one-dimensional result for Tj to Eq. (81), we see that the well-mixed model overestimates the contribution to the number of segregating sites from the recent part of genealogy with a large number of lineages. This should also be true for the initial stage j ≈ n, where Eq. (83) is not valid, because faster coalescence results from the fact that a lineages has to travel only about L/n to meet its neighbor. Other statistics that rely on the relative duration of periods with j lineages should be aﬀected in a similar way. This is particularly important when the deviations of the observed genealogical data from the predictions based on Eq. (81) are used to infer past evolutionary events, such as a selective pressure or geographic isolation (Wakeley, 2008), because some of these deviations could be due to the spatial structure of the population rather than external or internal perturbations. In summary, the classic theory of genetic inference can be extended to spatial populations. These extensions are not only more accurate and realistic than assuming the well-mixed-population dynamics, but also can be used to obtain information about the migration within the habitat (Wilkins and Wakeley, 2002). As spatially resolved genetic data sets become more readily available, better statistical tools based on spatial population genetics models will be needed.
X. CONCLUSIONS

(83)

for intermediate times, when, on one hand, the time is suﬃciently small for any lineage to diﬀuse across the whole habitat (τ ≪ L2 /Ds ), but, on the other hand, the time is suﬃciently large for neighboring lineages to −1 −1 coalesce (τ ≪ Ln−1 Dg + L2 n−2 Ds ).

Fluctuations due to sampling error during reproduction signiﬁcantly aﬀect the evolutionary dynamics of quasi one-dimensional populations, e.g. two-dimensional populations undergoing range expansions. These ﬂuctuations lead to the genetic demixing illustrated in Fig. 18, where an initially well-mixed population of alleles “phase separates” into monoallelic domains. The transition is somewhat analogous to spinodal decomposition in physics and material science (Scheucher and Spohn, 1988), but is also markedly diﬀerent. In particular, unlike conventional demixing phase transitions in ﬁnitetemperature statistical mechanics, genetic demixing occurs only in a low number of spatial dimensions d (d ≤ 2) (Duty, 2000; Scheucher and Spohn, 1988). The dependence of genetic demixing on the number of spatial dimensions d is illustrated by the decay of local heterozygosity in the absence of selection and mutation. For long times, the functional form of the decay is given by (Duty, 2000),

23 in s in the well-mixed-population model. The eﬀects of mutation do not diﬀer as dramatically in spatial and nonspatial models, but the stepping stone model reveals nontrivial spatial correlations and predicts a diﬀerent value for √ local steady state heterozygosity, proporthe tional to µ12 + µ21 for small mutation rates, compared to µ12 + µ21 in the well-mixed-population model. The evolutionary dynamics of spatial models also depends on the geometry of the expansion. For radial expansions, we found that the number of domains approaches a ﬁnite limit, which is, up to an additive constant, proportional to the square root of the initial radius of the colony R0 . Our main conclusion is that the data from natural populations may not always conform to the predictions of the well-mixed-population model and, even when it does, the estimated parameters from the model may not be biologically meaningful. The spatial model contains an important additional parameter, the spatial diﬀusion constant parameter Ds , which enters into many of the predictions. For example, the timescale for local ﬁxa2 tion is given by Ds /Dg rather than N τg [see Eq. (36)] and, for small selective advantage, s is sometimes re2 placed by s2 Ds /Dg , see Sec. VI. Moreover, as we saw in Sec. VII, the timescale of ﬁxation depends on the partitioning of the population by the experimenter into measurement sites. Thus care must be taken when interpreting the data from the natural populations. Finally, well-mixed-population models and experiments without spatial resolution do not account for spatial correlations, which contain important information about the population (see Secs. V and IX).

FIG. 18 (Color online) Illustration of spinodaldecomposition-like genetic demixing in a one-dimensional population. (a) Initially well-mixed population with red and green colors labeling diﬀerent genotypes. (b) The same population several generations later. The frequency of one of the alleles is now oscillating between 0 and 1 because the population segregates into monoallelic domains.

Note that d = 2 is the critical dimension. Here, we have shown that the one-dimensional stepping stone model has very diﬀerent dynamics compared to the standard well-mixed-population models used in population genetics. Most of the diﬀerences arise because, in the spatial model, populations segregate into monoallelic domains. As a result, genetic drift and selection can only act at the boundaries of the domains, which slows down the dynamics of the model. In particular, we found that, in the neutral model without mutation, ﬁxation occurs exponentially fast in a well-mixed population, but the decay of heterozygosity is algebraic in the spatial model. Genetic drift in the population as a whole becomes weaker with time, as spatial diﬀusion causes the eﬀective population size to increase. For a linear onedimensional model, we also found that the standard deviation of the total fraction of one of the alleles (in the absence of selection and mutation) increases subdiﬀusively as t1/4 . Selective sweeps also occur more slowly in the 2 spatial model: for weak selection, s ≪ Dg /Ds , we found that the time constant of the selective sweep is quadratic in s in the linear spatial model, but it is only linear

 −2t/(N τg ) e    (τf /t)1/2 H(t, 0) ∼  1/ ln(t)    const

d = 0, d = 1, d = 2, d > 2. (85)

Acknowledgments

M. A. is grateful for ﬁnancial support from the Danish National Research Foundation through Center for Models of Life. Overall support for this project was provided by the National Science Foundation, through Grant DMR-0654191, National Institute of General Medical Sciences Grant GM068763 of the National Centers for Systems Biology, and by the Harvard Materials Research Science and Engineering Center through DMR-0820484.

Appendix A: The Itˆ calculus o

In this appendix, we brieﬂy discuss the Itˆ calculus. o Our presentation relies on Risken (1989) and Gardiner (1985), which can be consulted for a more extensive presentation. For simplicity, we only consider nonspatial stochastic diﬀerential equations, but the results can be extended to spatial problems straightforwardly. Let us analyze the following stochastic diﬀerential equation, which includes Eq. (21) as a special case, dψ = ω(ψ) + g(ψ)Γ(t), dt

(A1)

24 where Γ(t) satisﬁes Eq. (22), and ω(ψ) and g(ψ) are arbitrary continuously diﬀerentiable functions. From the point of view of ordinary calculus, Eq. (A1) is not welldeﬁned because Γ(t) is discontinuous at every point. One way to circumvent this problem is to use discrete time steps of inﬁnitesimal length δt rather than continuous time. Then, Eq. (A1) takes the following form: ψ(t + δt) − ψ(t) = ω[ψ(t)] + g[ψ(t)]Γ(t). δt

d u[ψ(t)] =u′ [ψ(t)]ω[ψ(t)] + u′ [ψ(t)]g[ψ(t)]Γ(t) dt 1 + u′′ [ψ(t)]g 2 [ψ(t)], 2

(A4)

(A2)

However, this is not the only way to interpret Eq. (A1). For example, an alternative way to go from the continuous to a discrete description is to write Eq. (A1) as, ψ(t + δt) − ψ(t) ψ(t) + ψ(t + δt) =ω δt 2 (A3) ψ(t) + ψ(t + δt) Γ(t). +g 2 In fact, there is an inﬁnite number of ways to interpret Eq. (A1), depending on the relative weight of ψ(t) and ψ(t + δt) inside the arguments of the functions on the right hand side of the equation. The two most commonly used interpretations are Itˆ’s and Stratonovich’s o prescriptions. The former corresponds to Eq. (A2), and the latter to Eq. (A3). In physics, Stratonovich’s prescription is commonly used because Γ(t) is usually an approximation to a thermal noise with small but ﬁnite correlation time; therefore, the argument of g(·) should be an average value of ψ over the time that the correlations persist. In population genetics, on the other hand, Ito’s prescription is appropriate because a random change of the allele frequencies depends only on the genetic composition of the population prior to the change. Without the stochastic term, Eqs. (A2) and (A3) would yield the same results provided δt is suﬃciently small, but the stochastic terms remain diﬀerent even in the limit δt → 0. An easy way to see this difference is to average Eqs. (A2) and (A3) with respect o to the nondiﬀerentiable noise function Γ(t). Itˆ’s prescription gives ψ(t + δt) − ψ(t) = ω[ψ(t)] δt because g[ψ(t)]Γ(t) = g[ψ(t)] Γ(t) = 0 due to the independence of ψ(t) and Γ(t). A similar simpliﬁcation, however, cannot be applied to Stratonovich’s prescription because, generically, the stochastic term depends on ψ(t + δt), which is not independent of Γ(t). Because of the aforementioned ambiguity in interpreting stochastic diﬀerential equations with multiplicative noise, care must be taken while diﬀerentiating stochastic variables. While the rules of ordinary calculus apply to Stratonovich’s prescription, special rules of the Itˆ calcuo lus are required for Itˆ’s prescription when tracking the o evolution of a composite function u[ψ(t)] of the stochastic variable obeying Eq. (A1). In this paper, we use Itˆ’s o formula, namely (Gardiner, 1985; Risken, 1989)

where u(ψ) is a twice continuously diﬀerentiable function, and the primes now indicate diﬀerentiation with respect to ψ. The last term is the crucial addition due to the Itˆ o calculus. We conclude this discussion with an illustration of how Eq. (A4) can be used by deriving Eqs. (6) and (7) from Eq. (21) assuming s = 0 and µ12 = µ21 = 0. Thus, we start from the following equation of motion for f (t) df (t) = dt Dg f (t)[1 − f (t)]Γ(t) (Itˆ). o (A5)

Thus, ψ(t) = f (t), ω[ψ(t)] = 0, and g[ψ(t)] = Dg ψ(t)[1 − ψ(t)]. Since F (t) = f (t) , we obtain Eq. (6) by averaging Eq. (A5). For H(t) = h(t) = 2f (t)[1−f (t)] , we use Eq. (A4) with u[ψ(t)] = 2ψ(t)[1− ψ(t)] to obtain the equation of motion for h(t) dh(t) =0 + 2[1 − 2f (t)] Dg f (t)[1 − f (t)]Γ(t) dt 1 + (−4)Dg f (t)[1 − f (t)] (Itˆ). o 2

(A6)

Upon averaging Eq. (A6) with the rules described above, we obtain Eq. (7).
Appendix B: Solution of the Neutral Model Without Mutations

In this appendix, we solve Eq. (33) subject to the initial condition H(0, x) = H0 . It is advantageous to ﬁrst solve a simpler equation: ∂ ∂2 H = 2Ds 2 H − b(t)δ(x), ∂t ∂x (B1)

where b(t) is an arbitrary function of time. Equation (B1) is a standard diﬀusion equation with a sink term, and it can be readily solved in the Fourier domain. The result is
t ′ b(t ′

H(t, x) = H0 −

dt
0

)e

− 8D

x2 s (t−t′ )

8πDs (t − t′ )

.

(B2)

Note the convolution of b(t′ ) with the diﬀusion propagator. Now, we impose a self-consistency condition b(t) = Dg H(t, 0), which leads to
t

H(t, 0) = H0 − Dg

dt′

H(t′ , 0) 8πDs (t − t′ )

.

(B3)

0

25 This is Abel’s integral equation of the second kind, canonically written as
x

This relation is analogous to the one derived in ben-Avraham (1998). We can further simplify Eq. (C2) by observing that Dg H(t, 0) ∂H(t, +0) , = ∂x 4Ds (C3)

y(x) + λ
a

y(t)dt √ = g(x), x−t

(B4)

where g(x) is a known function. The general solution of Eq. (B4) given in Polyanin and Manzhirov (1998) reads y(x) = G(x) + πλ2
a x x

eπλ

2

(x−t)

G(t)dt,

where (B5)

which follows from integrating Eq. (32) or Eq. (D1) with respect to x from −ǫ to ǫ, 0 < ǫ ≪ 1, and noticing that H(t, x) is an even function of x. The ﬁnal result then reads ℓ(t) = 4Ds . Dg H(t, 0) (C4)

G(x) = g(x) − λ

a

g(t)dt √ . x−t

Equations (34) and (35) follow from Eqs. (B2), (B3), (B4), and (B5). For radial expansions considered in Sec. VIII, one can solve the equation of motion for H(t, φ) by following the same set of steps.
Appendix C: Average domain density from the spatial heterozygosity H(t, x)

It should be emphasized that this result is only valid in the limit of very large domain sizes compared to the boundary regions, which means H(t, 0) ≪ 1. Therefore the leading term in H(t, 0) is suﬃcient at this level of approximation. Note that Eqs. C2 and C4 are valid in the presence of genetic drift, migration, selection, and mutation. For radially expanding populations subject to inﬂation, Eq. C2 remains valid, but Eq. ( C4) is replaced by Eq. (66).
Appendix D: Inﬁnite Alleles Model

In this appendix, we derive the relationship between the spatial heterozygosity, H(t, x), and the average domain density nd (t). From nd , we can get a domain size by deﬁning ℓ ≡ n−1 . The result for the domain dend sity is valid for an arbitrary number of alleles, so in this appendix we use a broader deﬁnition of H(x, t) as the average probability of sampling at time t two diﬀerent alleles from two demes distance x apart. We assume that the domains have formed, and they are on average much larger than the boundary regions. Let h(t, x1 , x2 ) equal to 1 if both x1 and x2 are occupied by organisms in diﬀerent allelic state and 0 otherwise. To compute ℓ, we use an alternative deﬁnition of H(t, x) with ensemble average replaced by space average: 1 L→∞ L
L

H(t, x) = lim

h(t, ξ, ξ + x)dξ,
0

(C1)

where we assume periodic boundary conditions. Let us 1 L compute H(t, x+δx)−H(t, x) = limL→∞ L 0 [h(t, ξ, ξ + δx)−h(t, ξ, ξ)]dξ for δx small compared to typical domain size, but large compared to the deme spacing a. To do so, we expand both sides in δx. At the lowest order in δx, each domain boundary contributes δx to the right ∂ hand side; therefore, ∂x H(t, +0) equals the density of the domain boundaries. Upon deﬁning the average domain size ℓ(t) as the inverse of the domain boundary density, we obtain the following relationship: ∂H(t, +0) ∂x
−1

In this appendix, we extend the analysis of the stepping stone model with mutations presented in Sec. V to the inﬁnite alleles model. The inﬁnite alleles model assumes that every new mutation creates a new allele, which is a good approximation for genes encoded by a large number of nucleotides because the number of all possible mutations is much larger than the number of all possible back mutations (Hartl and Clark, 1989). The equation of motion for H(t, x), which we interpret as the average probability of sampling two diﬀerent alleles from demes x apart, can be derived by following two lineages backward in time, as done in Sec. IV. In the presence of mutation, the right hand side of Eq. (33) should contain an additional term describing the rate of increase of H(t, x) due to mutations in both of the lineages. Because, in the inﬁnite alleles model, a mutation changes the probability that the organisms have diﬀerent alleles from H to 1, that is by 1 − H, the new term is 2µ(1 − H), where µ is the mutation rate that is assumed to be the same for all types of mutations. Thus, Eq. (33) becomes ∂2 ∂ H = 2Ds 2 H + 2µ(1 − H) − Dg Hδ(x) ∂t ∂x

(D1)

for the inﬁnite alleles model [compare Eq. 47]. The stationary solution of Eq. (D1) is given by e− √
µ Ds

ℓ(t) =

.

(C2)

H(∞, x) = 1 −

|x|

1+

1 √Dg 4 Ds µ

.

(D2)

26 At large separations, H(∞, x) approaches one, which is consistent with the inﬁnite number of alleles. LoD cally, H(∞, 0) = (1 + 1 √Dg µ )−1 , and if H(∞, 0) ≪ 1 4 s the population is segregated into domains containing only one allelic type. The average size of such domains follows from Eq. (C4): 4Ds ℓ= Dg 1 Dg 1+ √ 4 Ds µ ≈ Ds , µ (D3)

∂Fij (t, x1 , x2 ) ∂2 ∂2 Fij (t, x1 , x2 ) =Ds 2 + ∂x2 ∂t ∂x1 2 + Dg δ(x1 − x2 )[δij Fi (t, x1 ) − Fij (t, x1 , x2 )] (E2)
q q

+
i′ =1

[µi′ i Fi′ j (t, x1 , x2 )
j ′ =1

where the last equality follows from the assumption that H(∞, 0) ≪ 1. The approach to the stationary state can be either obtained by methods of Appendix B or by ˆ the change of variables H(t, x) = H(∞, t)+ e−2µt H(t, x), which reduces Eq. (D1) to Eq. (33). The result is that the ˜ ˜ slowest decaying mode vanishes as Ct−1/2 e−2µt , where C is a constant. The inﬁnite allele model and Eq. (D1) has been analyzed before by Mal´cot (1975) and Nagylaki (1974), e who calculated the stationary solution and the long time approach to the equilibrium. Our results are consistent with their ﬁndings.
Appendix E: A model with several neutral alleles

+ µj ′ j Fij ′ (t, x1 , x2 )], where δij is Kronecker’s delta, which is zero if i = j and one otherwise. Thus, for a generic mutation matrix µij one has to solve a system of coupled linear partial diﬀerential equations. For simplicity and the ease of comparison with the other results in this paper, let us assume spatial homogeneity and identical mutation rates between any two alleles, µi=j = µ/q. Under these assumptions, Eq. (E2) can be simpliﬁed by introducing averaged spatial heterozygosity
q q q

H(t, x) =
i=1 j=1 j=i

Fij (t, 0, x) =
i=1

A model with q neutral alleles is an intermediate case between the two-alleles model that we focus on in this paper and the inﬁnite alleles model discussed in Appendix D. The q-alleles model is also analogous to nonequilibrium q-state Potts models. In this appendix, we brieﬂy outline how the q-alleles model can be formulated and solved in the language of one and two-point correlation functions, compare our analytical predictions to simulations, and extend Eq. (44) to the undulatingfront model. To specify the q-alleles model, we let fi (t, x) be the frequency of allele i at time t and position x; these quanq tities satisfy i=1 fi (t, x) = 1. The spatial diﬀusion and coalescence probability of two lineages are still characterized by Ds and Dg respectively. Intra-allelic mutations are described by the mutation matrix µij , which is the probability of allele i mutating into allele j if i = j. q When i = j, we let µii = − j=1, j=i µij to describe the outﬂow of alleles from allelic state i due to mutations. The dynamics of the q-alleles model can be analyzed by considering one-point correlation functions Fi (t, x) = fi (t, x) and two-point correlation functions Fij (t, x1 , x2 ) = fi (t, x1 )fj (t, x2 ) . Fi (t, x) is the probability to ﬁnd allele i at position x at time t, and Fij (t, x1 , x2 ) is the probability to simultaneously ﬁnd at time t allele i at position x1 and allele j at position x2 . The evolution equations for these correlation functions are obtained by tracing one and two lineages backward in time; the results are ∂Fi (t, x) ∂ 2 Fi (t, x) µji Fj (t, x), + = Ds ∂t ∂x2 j=1
q

fi (t, 0)[1 − fi (t, x)] ,

(E3) which is the probability to sample two diﬀerent alleles at time t distance x apart. The equation of motion for H(t, x) can be derived both from Eq. (E2) and, more simply, by tracing two lineages backward in time: ∂2 ∂ H = 2Ds 2 H + 2µ ∂t ∂x q−1 − H − Dg Hδ(x). (E4) q

(E1)

Note that Eq. (E4) agrees with Eq. (D1) in the limit q → ∞ and with Eq. (47) for µ12 = µ21 = µ/2. Since Eq. (E4) has the same functional form as Eq. (D1), the methods of Appendix D can be used to solve for H(t, x). In the absence of mutations, Eq. (E4) is identical to Eq. (33), as we brieﬂy mentioned in Sec. IV. However, q-alleles models with diﬀerent q may have slightly different dynamics due to q-dependent initial conditions: for example, an initially well-mixed population is represented by H(0, x) = H0 = 1 − 1/q. Thus the results of Sec. IV apply to the q-alleles model, provided appropriate initial conditions are used. In particular, we expect the standard deviation of fi (t), the total frequency of allele i, in a ﬁnite population to grow as t1/4 . This is indeed conﬁrmed by our simulations shown in Fig. 19. Spatial correlations in the nonequilibrium q-state Potts model have recently been analyzed by Masser and benAvraham (Masser and ben-Avraham, 2000), who also found that two-point correlation functions obey the same q-independent equation of motion. Finally, one can obtain the behavior of the standard deviation of the total frequency of allele one, ∆(t), in

27

FIG. 19 (Color online) Genetic drift during a linear range expansion in the ﬂat-front model with three alleles. (a) The genetic composition of the population [f1 (t), f2 (t), f3 (t)] pro3 jected on the plane i=1 fi (t) = 1 in a single run of the neutral 3-alleles model with a ﬂat front. The population is ﬁnite, L = 1000, and there are no mutations. (b) The average standard deviation of the frequency of allele one ∆(t), shown in blue, is obtained from 200 realizations of the simulations described in a. The red solid line shows the best power-law ﬁt, and the slope is close to the exponent expected from Eq. (44). The gray area encloses the points within one standard deviation from the mean. At t = 0, each site is assigned either allele one or allele two with equal probability, which corresponds to the center of the triangle in (a).

the undulating-front model by the following scaling argument. We consider ∆(t) at large times after monoallelic domains have formed. Let Nd (t) be the number of domains consisting of allele one and dk (t), k = 1, 2, ..., Nd(t) be lengths of these domains. Then, ∆(t) is given by   2

∆(t) =

1 L2

Nd (t)

Nd (t)

k=1

dk (t) −

k=1

We simplify Eq. (E5) by making an approximation that Nd (t) and dk (t) for k = 1, 2, ..., Nd(t) are independent random variables, which gives ∆2 (t) ≈ 1 Nd (t) [ d2 (t) − d1 (t) 2 ] 1 L2 2 2 + d1 (t) [ Nd (t) − Nd (t) 2 ] ,

dk (t) 

.

(E5)

(E6)

FIG. 20 (Color online) Genetic drift during a linear range expansion in the undulating-front model with three alleles. (a) The genetic composition of the population [f1 (t), f2 (t), f3 (t)] projected on the plane 3 fi (t) = 1 i=1 in a single run of the neutral 3-alleles model with an undulating front. The population is ﬁnite, L = 1000, and there are no mutations. (b) The average standard deviation of the frequency of allele one ∆(t), shown in blue, is obtained from 200 realizations of the simulations described in a. The red solid line shows the best power-law ﬁt, and the slope is close to the exponent expected from Eq. (E8). The gray area encloses the points within one standard deviation from the mean. At t = 0, each site is assigned either allele one or allele two with equal probability, which corresponds to the center of the triangle in (a).

where we used the fact that di (t) are identically distributed. By using ﬁrst passage time analysis discussed 2 in Redner (2001), one can show that Nd (t) − 2 2 2 Nd (t) ∝ Nd (t) [ d1 (t) − d1 (t) ]/ d1 (t) 2 ∝ 2 L[ d1 (t) − d1 (t) 2 ]/ d1 (t) 3 . Thus ∆2 (t) ∝ 1 [ d2 (t) − d1 (t) 2 ]. 1 L d1 (t) (E7)

tζ/2 t1/3 ∆(t) ∝ √ ∝ √ , L L

(E8)

Upon recalling, that, in the undulating-front model, d1 (t) ∝ tζ , and d1 (t) 2 ∝ t2ζ , we conclude that

where, in the last proportionality, we used ζ = 2/3 from Saito and M¨ ller-Krumbhaar (1995). Equau tion (E8) is in good agreement with the simulations of the undulating-front model shown in Fig. 20.

28
Appendix F: Connection with the voter model and one-dimensional reaction kinetics

The stepping stone model with only one organism per island or “deme,” N = 1, has been extensively studied in probability theory (Durrett, 1988; Liggett, 2004) ´ and nonequilibrium statistical mechanics (Odor, 2004), where it is known as the voter model. The model typically considers a set of voters on a hypercubic lattice in d-dimensions. Each voter holds one of the q possible opinions about an issue (corresponding to q alleles in population genetics), and, at a certain rate, each voter reconsiders the issue, and adopts the opinion of a randomly chosen nearest neighbor. The voter model can be mapped onto the dynamics of the q-state Potts model at zero temperature. In one and two dimensions, opinions in the voter model coarsen spatially with time, and the model approaches one of the q absorbing states, in which all the voters have the same opinion (Cox and Griﬀeath, 1986; Duty, 2000). In higher dimensions, the voters still form cluster of opinions, but these clusters stop growing after reaching a certain limiting size. Selection and mutation are typically not considered in voter models. The voter model can be solved exactly by tracing the history of opinion adoptions backward in time (Cox and Griﬀeath, 1986; Duty, 2000; Scheucher and Spohn, 1988). The opinion of a given individual performs a random walk as we follow the opinion from its current holder to its ultimate ancestor. With this observation, we can easily understand how the behavior of the voter model depends on the number of spatial dimensions. In one and two dimensions, a pair of random walks always meet (Redner, 2001), so the histories of opinion adoptions starting from two diﬀerent voters will eventually converge to a single voter as we trace them backward in time. Therefore, any two voters should have the same opinion after a suﬃciently long time has elapsed. In higher dimensions, however, there is a ﬁnite probability that two random walkers never meet (Redner, 2001); therefore, the voters never agree, and an absorbing state is never reached. Another important property of the voter model is that the dynamics occurs only at the boundaries between the opinion clusters; inside a cluster the opinions cannot change because every voter has the same opinion as its nearest neighbors. This property is particularly useful in one spatial dimension, where it allows us to map the dynamics of the voter model to the one-dimensional diﬀusion-limited chemical kinetics of point particles. We identify each domain wall with a particle performing a random walk due to opinion changes at the boundary. When two particles meet, they react with two possible outcomes. They annihilate (A + A → 0) if the ﬂanking domains have the same opinion or coalesce (A + A → A) otherwise. If the initial state is uncorrelated, the annihilation occurs with probability 1/(q − 1), and the coalescence with probability (q − 2)/(q − 1). In onedimension, this reaction diﬀusion system has been an-

alyzed by Masser and ben-Avraham (2000), who found that the density of the domain walls decays as t−1/2 in agreement with Eq. (38). A related model of annihilating random walks for radial and linear range expansions was solved by Hallatschek and Nelson (2010). It is not surprising that the voter model and the stepping stone model have the same long time behavior in one dimension. At long times, most of the voters belong to large domains; therefore, we do not aﬀect the system by combining neighboring sites into larger coarse-grained demes, as in Sec. VII. For these large demes, the equations of motion of the stepping stone model are valid, so the two models are equivalent in the long time limit. The voter and the stepping stone models are also equivalent in the small-Ds limit of very slow migration. In this case, each deme reaches ﬁxation much faster than it sends out or accepts new migrants; hence, the stepping stone model reduces to the voter model with one voter representing an entire deme. We can further illustrate the connection between the stepping stone model and the voter model by calculating the probability that two voters l sites apart have diﬀerent opinions. This probability is analogous to the average spatial heterozygosity, so we call it H(t, l). The equation of motion for H(t, l) is obtained by following the histories of opinion adoptions backward in time. Since H(t, l) changes only due to the diﬀusion of the history traces, the equation of motion reads d H(t, l) = [H(t, l − 1) + H(t, l + 1) − 2H(t, l)], (F1) dt where we measure time in such units that the rate of opinion adoption is set to unity. While Eq. (F1) can be solved exactly (Houchmandzadeh and Vallade, 2003), it is more instructive to go to the continuum limit, in which the equation of motion for H(t, x) takes the following form ∂ ∂2 H(t, x) = 2Ds 2 H(t, x), ∂t ∂x

(F2)

where Ds denotes the spatial diﬀusion constant as in Eq. (33). Upon comparing Eqs. (33) and (F2), one might naively conclude that the voter model corresponds to Dg = 0 limit of the stepping stone model; in fact, the opposite is true: the voter model corresponds to the limit Dg = ∞. Qualitatively, one can see this from the fact that Dg ∝ N −1 , so, as the deme size N approaches its lowest value of 1, we expect Dg to increase. On more rigorous grounds, we should note that the role of the delta function in Eq. (33) is to enforce a boundary condition at x = 0, provided one considers H(t, x) only for x > 0. This boundary condition is derived in Appendix C and is given by Eq. (C3). The corresponding boundary condition for Eq. (F2) is H(t, 0) = 0 because the probability

29 of one voter having two diﬀerent opinions is zero. We indeed recover H(t, 0) = 0 by letting Dg → ∞ in Eq. (C3). One can solve Eq. (F2) for the initial condition H(0, x) = H0 by the Laplace transform in time or a self-similar ansatz; the solution reads |x| √ 8Ds t
References Araten, D., D. Golde, R. Zhang, H. Thaler, L. Gargiulo, R. Notaro, and L. Luzzatto, 2005, Cancer Res. 65, 8111. Bar-Yam, Y., 2005, Making Things Work: Solving Complex Problems in a Complex World (Knowledge Press). Barton, N., F. Depaulis, and A. Etheridge, 2002, Theor. Popul. Biol. 61, 31. ben-Avraham, D., 1998, Phys. Rev. Lett. 81, 4756. Blythe, R., and A. McKane, 2007, J. Stat. Mech 7018, 1. Bramson, M., and D. Griﬀeath, 1980, Ann. Probab. 8, 183. Bramson, M., and J. Lebowitz, 1991, J. Stat. Phys. 62, 297. Charlesworth, B., D. Charlesworth, and N. Barton, 2003, Annu. Rev. Ecol. Evol. Syst. 34, 99. Cox, J., and D. Griﬀeath, 1986, Ann. Probab. 14, 347. Crow, J., and M. Kimura, 1970, An Introduction to Population Genetics Theory (Harper & Row, New York). Doering, C., and D. ben Avraham, 1988, Phys. Rev. A 38, 3035. Doering, C., C. Mueller, and P. Smereka, 2003, Physica A 325, 243. Drake, J., 1991, Proc. Natl. Acad. Sci. USA 88, 7160. Durrett, R., 1988, Lecture Notes on Particle systems and Percolation (Wadsworth Publishing Company). Duty, T. L., 2000, PhD thesis, (The University of British Columbia). Fisher, R., 1937, Ann. Eugenics 7, 353. Gardiner, C., 1985, Handbook of Stochastic Methods (Springer, New York). Gibson, A., 1965, Oklahoma: A History of Five Centuries (Harlow, Norman). Goldenfeld, N., 1992, Lectures on Phase Transitions and the Renormalization Group (Westview Press). Gradshteyn, I., and I. Ryzhik, 1980, Table of Integrals, Series, and Products (Academic Press, New York). Guth, A., 1981, Phys. Rev. D 23, 347. Hallatschek, O., P. Hersen, S. Ramanathan, and D. R. Nelson, 2007, Proc. Natl. Acad. Sci. USA 104, 19926. Hallatschek, O., and K. S. Korolev, 2009, Phys. Rev. Lett. 103, 108103. Hallatschek, O., and D. R. Nelson, 2010, Evolution 64, 193. Hartl, D., and A. Clark, 1989, Principles of population genetics (Sinauer Associates, Sunderland). Hewitt, G., 1996, Biol. J. Linn. Soc. 58, 247. Hinrichsen, H., and M. Howard, 1999, Eur. Phys. J. B 7(4), 635. Houchmandzadeh, B., and M. Vallade, 2003, Phys. Rev. E 68, 61912. Kimura, M., 1955, Proc. Natl. Acad. Sci. USA 41, 144. Kimura, M., 1969, Genetics 61, 893. Kimura, M., 1983, The Neutral Theory of Molecular Evolution (Cambridge University Press). Kimura, M., and G. Weiss, 1964, Genetics 49, 561. Kingman, J., 1982, J. Appl. Prob. , 27. Kolmogorov, A., N. Petrovsky, and N. Piscounov, 1937, Moscow University Bulletin of Mathematics 1, 1. Liggett, T., 2004, Interacting Particle Systems (Springer Verlag). Lin, J., 1991, Phys. Rev. A 44, 6706. Mal´cot, G., 1955, Cold Springs Harbor Symp. Quant. Biol e 20, 52. Mal´cot, G., 1975, Theor. Popul. Biol. 8, 212. e Masser, T., and D. ben-Avraham, 2000, Phys. Lett. A 275,

H(t, x) = H0 erf

.

(F3)

We can now compute the average size of the domains with the help of Eq. (C2). As we expect, the result is given by Eq. (38) because the long time limits of the stepping stone model and the voter model agree. The Dg = ∞ approximation is particularly valuable for circular fronts undergoing inﬂation because the exact solution of the stepping stone model in this case [Eqs. (64) and (65)] is rather unwieldy. The equation of motion for H(t, ϕ) in the voter model with inﬂation is given by ∂ ∂2 2Ds H(t, ϕ). H(t, ϕ) = 2 ∂ϕ2 ∂t (R0 + vt)

(F4)

One can solve Eq. (F4) in the Fourier domain, and compute the nontrivial limit-shape as t → ∞. The result reads R0 v 8Ds

H(∞, ϕ) = H0 erf

|ϕ|

.

(F5)

With the help of the angular version of Eq. (C2) [see Eq. (66)], we calculate the ﬁnal number of sectors: 2πR0 v , Ds

N (∞) = H0

(F6)

which agrees with Eq. (69) in the limit Dg = ∞. In the same limit, we can also obtain an analytical expression for the long time variance ν(∞): 1 4π
π −π

ν(∞) =

[H0 − H(∞, ϕ)]dϕ = H0

2Ds , (F7) π 3 R0 v

where we used the relationship between the variance ν(t) and the global heterozygosity H(t, ϕ) given by the spatial generalization of Eq. (8). Finally, we note that the mapping to a one-dimensional reaction-diﬀusion system of particles could be generalized to account for super-diﬀusive boundaries in the undulating-front model, for example, by considering continuous time L´vy ﬂights instead of random walks; e see (Hinrichsen and Howard, 1999). In the chemical kinetics picture, one can also account for mutations by introducing a birth process 0 → 2A and for natural selection by imposing an attraction between the particles ﬂanking domains of the deleterious allele.

30
382. Mayr, E., 1942, Systematics and the Origin of Species from the Viewpoint of a Zoologist (Columbia University, New York). Murray, J., 2003, Mathematical Biology (Springer). Nagylaki, T., 1974, Proc. Natl. Acad. Sci. USA 71, 2932. Nagylaki, T., 1978, Proc. Natl. Acad. Sci. USA 75, 423. Nordborg, M., 1998, The Am. J. of Hum. Genet. 63, 1237. Nowak, M., 2006, Evolutionary Dynamics: Exploring the Equations of Life (Harvard University, Cambridge MA). ´ Odor, G., 2004, Rev. Mod. Phys. 76, 633. Poli, R., W. Langdon, and N. McPhee, 2008, A Field Guide to Genetic Programming (Lulu Enterprises, UK Ltd). Polyanin, A., and A. Manzhirov, 1998, Handbook of Integral Equations (CRC Press, Boca Raton). Ramachandran, S., O. Deshpande, C. Roseman, N. Rosenberg, M. Feldman, and L. Cavalli-Sforza, 2005, Proc. Natl. Acad. Sci. USA 102, 15942. Redner, S., 2001, A Guide to First Passage Processes (Cambridge University Press). Risken, H., 1989, The Fokker-Planck equation: Methods of Solution and Applications (Springer, Berlin and Heidelberg). Saito, Y., and H. M¨ller-Krumbhaar, 1995, Phys. Rev. Lett. u 74, 4325. Scheucher, M., and H. Spohn, 1988, J. Stat. Phys. 53, 279. Templeton, A., 2002, Nature 416, 45. Wakeley, J., 2008, Coalescent Theory: An Introduction (Roberts & Company Publishers). Wilkins, J., and J. Wakeley, 2002, Genetics 161, 873. Wilkinson-Herbots, H., 1998, J. Math. Biol. 37, 535. Wilson, K., and J. Kogut, 1974, Phys. Rep. 12, 75. Wright, S., 1931, Genetics 16, 97. Wright, S., 1943, Genetics 28, 114. Zhong, D., and D. ben-Avraham, 1995, Phys. Lett. A 209, 333.