Genetic demixing and evolution in linear stepping stone models

Results for mutation, selection, genetic drift, and migration in a one-dimensional continuous population are reviewed and extended. 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 first 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 micro-organisms undergoing a range expansion on a Petri dish. 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 effects of mutation, however, are not significantly affected by the segregation. Although fixation in the neutral well-mixed (or “zero-dimensional”) model occurs exponentially in time, it occurs only algebraically fast in the one-dimensional model. An unusual sublinear increase is also found 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 different. The relatively unexplored problem of evolutionary dynamics at the edge of an expanding circular colony is studied as well. Also reviewed are how the observed patterns of genetic diversity can be used for statistical inference and the differences are highlighted between the well-mixed and one-dimensional models. Although the focus is on two alleles or variants, $q$-allele Potts-like models of gene segregation are considered as well. Most of the analytical results are checked with simulations and could be tested against recent spatial experiments on range expansions of inoculations of Escherichia coli and Saccharomyces cerevisiae.

Figure 1

(Color online) Spatial segregation in an expanding microbial population. Different colors label different 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-defined monoallelic domains. The colony is of the order $1\mathrm{cm}$ in height. Details of the experiment are presented in 22.

Figure 2

(Color online) An illustration of the two models of a growing front. (a) and (c) Illustration of the model with a rough undulating front, which is a natural result of an unconstrained two-dimensional growth. (b) and (d) Illustration of the model with a flat front, which is constrained to have no lateral undulations to simplify the analytical analysis. The blank hexagons represent empty sites, and different colors of the occupied hexagons represent different alleles. (a) The model of an undulating population front. The highlighted hexagon is a randomly chosen cell that can reproduce and deposit an identical offspring 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 filled only by an offspring of one of its nearest neighbors in the previous generation (shown with arrows). Both (a) and (b) show the effects of genetic drift (sampling error) when, e.g., the second from the left cell in the bottom row leaves no offspring. Such events lead to coarsening seen in (c) and (d). (c) and (d) 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 flat, whereas in (c) it is rough. This roughness affects 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 by 22, domain walls are expected to wander more vigorously in (c) than in (d). Despite the apparent differences, both models exhibit the same qualitative behavior.

Figure 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.

Figure 4

(Color online) The stationary distribution $P(\infty ,f)$ in the presence of selection, mutation, and genetic drift [see Eqs. (19, 20)]. The dotted line shows $P(\infty ,f)$ for $s={\mathfrak{D}}_g$ and ${\mu }_{12}={\mu }_{21}=10{\mathfrak{D}}_g$, which corresponds to weak genetic drift, ${\mu }_{12}∕D_g⪢1$. The solid line shows $P(\infty ,f)$ for $s={\mathfrak{D}}_g$ and ${\mu }_{12}={\mu }_{21}=0.1{\mathfrak{D}}_g$, which corresponds to strong genetic drift, ${\mu }_{12}∕D_g⪡1$. Note the difference 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 to ${\mathfrak{D}}_g∕2$ and $P(\infty ,f)$ diverges at $f=0$ and 1. Also note that the effect 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.

Figure 5

(Color online) Solutions of Eq. (33) at various times given $H(0,x)=\frac{1}{2}$. Time and distance are measured in units such that $D_s=1$ and $D_g=1$. Time increases from the top curves to the bottom curves. Note the statistical reflection symmetry, $H(t,x)=H(t,-x)$.

Figure 6

(Color online) The number of monoallelic domain boundaries as a function of time in the undulating-front model. The simulation of 100 demes averaged over 100 runs is plotted (dots), and the theoretically expected decay of the number of boundaries as $t^{-2∕3}$ [see 52] is plotted (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 dots are the simulation data and the solid line is the expected slope according to 52.

Figure 7

(Color online) The number of monoallelic domain boundaries as a function of time in the flat-front model. The simulation of 100 demes averaged over 100 runs is plotted (dots), and the theoretically expected decay of the number of boundaries as $t^{-1∕2}$ (see 24) is plotted (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 dots are the simulation data and the solid line is the expected slope according to 7 and 24 and Eq. (38).

Figure 8

(Color online) A single run of the flat-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 [defined in the sense of Eq. (C1) but without taking the limit $L\to \infty$] for three times measured in generations: $t=0$, 1024, and 4096. 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 fluctuates 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{\tau }_g$ lies completely below $1∕2$ because, at this time, the relative fraction of the alleles deviates significantly from the initial 50:50 ratio due to genetic drift.

Figure 9

(Color online) The effects of coarse graining on the time evolution of the spatial heterozygosity $H(t,l)$ averaged over 100 realizations of the flat-front model with a 1000 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)$. (b) Qualitative agreement with the solution of stepping stone model displayed in Fig. 5, and (a) outside the region around $l=0$. Note that, unlike the calculation presented in Fig. 5, this simulation was run sufficiently long to show the effects of the boundary conditions.

Figure 10

(Color online) Comparison between the analytical prediction for the limiting shape of $H(\infty ,\overline{x})$ and the simulations of the flat-front model. The continuous curve (black) is formed by the data points representing $H(t,x)$ for several times between $t=2\times {10}^4$ and $4\times {10}^5$ plotted in the rescaled coordinates $\overline{x}$ and the circles 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\times {10}^5$ generations with averaging over 500 realizations. At $t=0$, each site is assigned either allele one or allele two with equal probability.

Figure 11

(Color online) Genetic drift in a finite population. At $t=0$, each site is assigned either allele one or allele two with equal probability. (a) The total fraction of allele one $\mathfrak{f}\left(t\right)$ vs time in four single runs of the neutral model with a flat front. Here $L=1000$, and there are no mutations. (b) The average standard deviation of the frequency of allele one $\Delta \left(t\right)$ is obtained from 200 realizations of the simulations described in (a). The solid line shows the best power-law fit, 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.

Figure 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(\infty ,x)$ and the simulations with ${\tilde{\mu }}_{12}={\tilde{\mu }}_{21}={10}^{-4}$. The black dots represent the results of the simulation and the circles represent the best fit of theoretical result given by Eq. (48) to the data. Here only $D_g$ is a fitting parameter; the values of $D_s$, ${\mu }_{12}$, and ${\mu }_{21}$ follow from the correspondence between the discrete and continuum models. The data are obtained in a simulation of 3200 individuals for $2\times {10}^5$ generations with averaging over 100 realizations. At $t=0$, each site is assigned either allele one or allele two with equal probability.

Figure 13

(Color online) The effective extinction rate $\alpha$ vs $\tilde{s}$ in the limit of weak selection. The 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 support $\alpha \propto {\tilde{s}}^2$. The values of $\alpha$ are obtained from graphs like the one shown in the inset. Inset: $\ln (1-F)$ vs $t$ for $\tilde{s}=0.12$. The circles are the actual data points and the line is the best least-squares linear fit. The simulation confirms exponentially fast fixation. 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.

Figure 14

(Color online) Spatial segregation in an expanding circular bacterial colony of E. coli. Different colors label different 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-defined monoallelic domains, which coalesce at early times but seem to stop coalescing in the final stages of the experiment. Note that the boundaries between the domains are biased to move away from each other due to inflation, in addition to their diffusive random-walk-like motion. Details of the experiment are presented by 22.

Figure 15

(Color online) Solutions of Eq. (63) with $\sigma =1$ at various rescaled times $\mathsf{t}$ given random initial conditions on the circle bounding the homeland, $H(0,\varphi )=\frac{1}{2}$. Time increases from the top curves to the bottom curves. Note that there is no observable difference between $H(500,\varphi )$ and $H(1000,\varphi )$ because $H(\mathsf{t},\phi )$ reaches a nontrivial limit shape as $\mathsf{t}\to \infty$.

Figure 16

A plot of $K\left(\sigma \right)$ for Eq. (71) from a numerical solution of Eq. (63).

Figure 17

(Color online) An illustration of the backward-in-time 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 differ in at least one nucleotide, but in experiments one often finds organisms that have identical sequences. We trace the spatial diffusion 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 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. The ancestral process shown here satisfies the infinite site 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 shown in the text.

Figure 18

(Color online) Illustration of spinodal-decomposition-like genetic demixing in a one-dimensional population. (a) Initially well-mixed population with colors labeling different 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.

Figure 19

(Color online) Genetic drift during a linear range expansion in the flat-front model with three alleles. (a) The genetic composition of the population $[{\mathfrak{f}}_1\left(t\right),{\mathfrak{f}}_2\left(t\right),{\mathfrak{f}}_3\left(t\right)]$ projected on the plane ${\sum }_{i=1}^3{\mathfrak{f}}_i\left(t\right)=1$ in a single run of the neutral three-allele model with a flat front. The population is finite, $L=1000$, and there are no mutations. (b) The average standard deviation of the frequency of allele one $\Delta \left(t\right)$ is obtained from 200 realizations of the simulations described in (a). The solid line shows the best power-law fit, 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).

Figure 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 $[{\mathfrak{f}}_1\left(t\right),{\mathfrak{f}}_2\left(t\right),{\mathfrak{f}}_3\left(t\right)]$ projected on the plane ${\sum }_{i=1}^3{\mathfrak{f}}_i\left(t\right)=1$ in a single run of the neutral three-allele model with an undulating front. The population is finite, $L=1000$, and there are no mutations. (b) The average standard deviation of the frequency of allele one $\Delta \left(t\right)$ is obtained from 200 realizations of the simulations described in (a). The solid line shows the best power-law fit 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).