Radial Domany-Kinzel models with mutation and selection

We study the effect of spatial structure, genetic drift, mutation, and selective pressure on the evolutionary dynamics in a simplified model of asexual organisms colonizing a new territory. Under an appropriate coarse-graining, the evolutionary dynamics is related to the directed percolation processes that arise in voter models, the Domany-Kinzel (DK) model, contact process, and so on. We explore the differences between linear (flat front) expansions and the much less familiar radial (curved front) range expansions. For the radial expansion, we develop a generalized, off-lattice DK model that minimizes otherwise persistent lattice artifacts. With both simulations and analytical techniques, we study the survival probability of advantageous mutants, the spatial correlations between domains of neutral strains, and the dynamics of populations with deleterious mutations. “Inflation” at the frontier leads to striking differences between radial and linear expansions. For a colony with initial radius $R_0$ expanding at velocity $v$, significant genetic demixing, caused by local genetic drift, occurs only up to a finite time $t^{*}=R_0/v$, after which portions of the colony become causally disconnected due to the inflating perimeter of the expanding front. As a result, the effect of a selective advantage is amplified relative to genetic drift, increasing the survival probability of advantageous mutants. Inflation also modifies the underlying directed percolation transition, introducing novel scaling functions and modifications similar to a finite-size effect. Finally, we consider radial range expansions with deflating perimeters, as might arise from colonization initiated along the shores of an island.

Figure 1

A sketch of the phase diagram for mutation-selection balance in linear expansions with irreversible, deleterious mutations from green (light gray) to red (dark gray) cells occurring at rate ${\mu }_f$ at a fitness cost $s$ (see Appendix  for a biological motivation). The solid line between red (dark gray) and green (light gray) regions represents the directed percolation phase transition. We show as insets typical clusters generated from single green (light gray) cells in the active and inactive phases. The dashed line at the top represents the case of no mutations where the system falls into the compact directed percolation (CDP) universality class. The dynamics at criticality ($s={\mu }_f=0$ at the black dot) in this CDP regime is shown. The dotted line shows the position of the mean-field transition (i.e., for a well-mixed population). An analogous plot for inflating radial range expansions is shown in Fig. 11.

Figure 2

(a) A sketch of the Domany-Kinzel (DK) model implemented on a hexagonal lattice so each daughter cell has two potential parents. Periodic boundary conditions are employed, and each cell is evolved using input from the states of the two cells above it in the previous generation. (b) Illustration of the two update steps in our generalized DK model. The daughter shares the color of two parents of the same color. Otherwise, the daughter is green (light gray) with probability $p_G={(2-s)}^{-1}\approx 1/2+s/4$ if its parents are different colors, where $s\le 1$ is the selective advantage of the green strain. We set $p_R=1-p_G=\frac{1-s}{2-s}\approx 1/2-s/4$. A newly evolved cell is then allowed to mutate from green (light gray) to red (dark gray) with a forward mutation probability ${\mu }_f\ll 1$ or vice versa with a backward mutation probability ${\mu }_b\ll 1$.

Figure 3

An illustration of a radial expansion simulation on a hexagonal lattice with an initial state given by a “homeland” of seven cells. We add daughter cells one at a time, each time picking an unoccupied lattice site that is closest to the center of the homeland (the central purple cell). When two or more cells are at the same distance, we chose one of them at random. We implement the generalized Domany-Kinzel update rule that allows for more than two potential parent cells. Each daughter cell will be green (light gray) with a probability $p_G$ that will depend on the number of green (light gray) parents $n_G$, given by Eq. (1). Mutations can be implemented in a second update step as described in the text (also see Fig. 2).

Figure 4

(a) A radial DK model simulation on a hexagonal lattice. The structure factor of the lattice is shown in the right panel. The cell colors were assigned with neutral selective advantage $(s=0)$ and no mutations $({\mu }_f=0)$. The first 1500 cells generated by the algorithm are shown. The central shaded region is the homeland of radius $R_0\approx 10$ cell diameters, where green (light gray) and red (dark gray) cells are placed randomly with equal probability. (b) The same simulation performed on a Bennett model lattice with approximately equal numbers of small and large cells $(Q=0.5)$ and a radius ratio $\rho =0.769$. The structure factor in the right panel shows that the Bennett cell packing is isotropic.

Figure 5

Three examples of surfaces generated by the set of points $\mathbf{r}(t,\varphi )$ given by Eq. (12) with $0\le t\le 2$ and $-\pi \le \varphi <\pi$. We set the initial population circumference and the crossover time to one: $R_0=t^{*}=1$. The $\Theta =1/2$ surface is shaped like a bowl and has positive Gaussian curvature $K>0$. Conversely, the $\Theta =2$ surface is shaped like a trumpet and has $K<0$. Confining a population to grow along the $\Theta =1$ surface (with $K=0$) is another way of implementing the regular radial range expansion with linear inflation.

Figure 6

Nearly perfect collapse of the simulation data using the scaled variable $\xi =R_0\varphi /\sqrt{8D_r{t}_c\left(t\right)}$ for the heterozygosity with $H_0=1/2$, corresponding to random initial conditions, with equal probabilities for the two different colors. The initial colony radius is given by $R_0=50,100,200,300$. Linear range expansions correspond to the limit $R_0\to \infty$. All of the radial expansion simulations were fit with a single parameter $D_r$. The scaling variable is $\xi =x/\sqrt{8Dt}$ for linear expansions.

Figure 7

A calculation of the heterozygosity $H(\varphi ,t)$ at time $t=48\gg t^{*}$ on different lattices. $R_0=16$, $v=1$, ${\mu }_f=0.1$, and $s=0.57$ in the simulations. The length and time units are measured using the rescaling discussed in Sec. 2. Oscillations in the heterozygosity on regular lattices reflect the four- and sixfold rotational symmetries. These oscillations are absent for the amorphous Bennett model lattice with $\rho =0.769$ and $Q=0.5$.

Figure 8

Two examples of single sector dynamics with $R_0=50$ (in units of the average cell diameter) evolved over about 250 generations for two very different values of the key inflationary parameter $\kappa \approx s\sqrt{R_0/2v{\tau }_g}$ discussed in the text. In (a) a combination of selection bias and inflation dominates genetic drift after a short initial transient, with $\kappa \approx 2.4$. The initial sector width in radians at the edge of the homeland is ${\varphi }_0\approx 0.1$. In (b) $\kappa \approx 0.05\ll 1$, so genetic drift can influence the result expected from selection bias and inflation alone for an extended period. We set ${\varphi }_0\approx 0.2$ in this case. We see that the fate of the sector in (a) is decided much faster than the fate of the sector in (b), which still might eventually die out.

Figure 9

A plot of the crossover time ${\tau }_c$ that marks when genetic drift is no longer able to keep up with inflation. The time ${\tau }_c$ is given as a dimensionless time coordinate, related to the time $t$ via $\tau =vt/(R_0+vt)=vt/R\left(t\right)$. The dashed lines show the asymptotic approximations to ${\tau }_c$ given in Eq. (41).

Figure 10

Radial survival probabilities $S_r\left(\tau \right)$ for a green sector (like that in Fig. 8) of size ${\varphi }_0\approx 0.036$ initialized on the rim of a homeland with $R_0=300$, evolved over 1750 generations, for a variety of selective advantages $s$. The solid lines represent Eq. (45), the adiabatic approximation to $S_r\left(\tau \right)$. The integral in Eq. (45) is evaluated numerically. The survival probabilities were averaged over $8\times {10}^4$ independent runs. We fit the black line (the $s=0$ case) with the lattice parameter ${\tilde{a}}_r$, finding ${\tilde{a}}_r\approx 0.9$. The other lines use $\gamma \approx 1.2$ (see Appendix ).

Figure 11

The color map of the average fraction of green cells at long times ($t=1800$ generations) for a radial expansion with a homeland radius of $R_0=200$, averaged over 400 runs, for positive values of ${\mu }_f$ and $s$. The average fraction $\rho \left(t\right)$ approaches a steady state at long times so ${\rho }_{\infty }={\lim }_{t\to \infty }\langle f\left(t\right)\rangle ={\lim }_{t\to \infty }{\langle f(\mathbf{x},t)\rangle }_{\mathbf{x}}\approx \rho (t=1800)$. This approximation only breaks down extremely close to the phase boundary between the red (dark gray) and green (light gray) regions and the difference is not visible at the resolution of the plot. The dotted line shows the mean-field transition. The dashed line shows the approximate position of the directed percolation phase transition for a linear expansion, as shown in Fig. 1.

Figure 12

Linear and radial range expansions with ${\mu }_b=0$ and ${\mu }_f=0.1$ for various $s$ near the directed percolation phase transition. Green (light gray) organisms have a selective advantage $s$ but are subject to forward mutations to a less fit red (dark gray) strain. In the radial expansions (a) $R_0=50$ and the system evolved for $\sim$250 generations. The dark blue solid circles indicate the crossover time $t^{*}=50$. In (b) a linear range expansion of 300 cells evolves for 400 generations.

Figure 13

Data collapse for a radial range expansion with $R_0=2048$ for the average number of green cells over about 600 generations. We test many values for $\tau =s-s_c=\pm 0.01,0.02,...$ for ${\mu }_f=0.1$. The collapses are consistent with the directed percolation values $\alpha \approx 0.159464\left(6\right)$ and ${\nu }_{\parallel }\approx 1.733847\left(6\right)$. The critical selection parameter was approximately $s_c\approx 0.55\left(1\right)$.

Figure 14

We plot the average green cell density $\rho \left(t\right)\equiv \langle f\left(t\right)\rangle$ at long times $t$ versus a control parameter ${\mu }_f/s^2$ for an initially all-green population for radial and linear expansions. For radial expansions, $f\left(t\right)$ is averaged over ${10}^3$ runs at $t=2\times {10}^3$ generations. In the linear case, we averaged $f\left(t\right)$ over $2\times {10}^3$ runs, at time $t=2\times {10}^4$, for a system size of ${10}^3$. Away from the transition region (${\mu }_f/s^2\lesssim 0.2$), the time $t$ is long enough so the density $\rho \left(t\right)$ for both expansions is approximately equal to the steady-state density ${\rho }_{\infty }\equiv \langle f(t\to \infty )\rangle$, with an error smaller than the point size. The black line shows the predicted scaling of the steady-state density for small ${\mu }_f/s^2$. We varied $s$ between 0 and 1 and ${\mu }_f$ between 0 and $0.5$.

Figure 15

Data collapses for radial range expansion simulations with many values of $t^{*}=R_0/v$ (with ${\mu }_f=0.1$ and $s=0.547$). In (a) we plot the scaled green cell density $\rho \left(t\right)=\langle f\left(t\right)\rangle$ with an all-green initial condition for $R_0=8,16,32,...,2048$. The regimes $t\ll t^{*}$ and $t\gg t^{*}$ can be described by power laws. The $t^{-\alpha }$ decay for $t<t^{*}$ with $\alpha \approx 0.16$ is expected from directed percolation theory [22], but the exponent $\alpha \approx 0.26$ that arises for long times appears to be new. In (b) we plot the scaled sector survival probability $S_{\mathrm{ss}}\left(t\right)$ of a single green cell inoculated in an all-red homeland $R_0=16,32,...,2048$. The dashed line indicates the expected directed percolation result: a decay with slope $-\alpha \approx -0.16$. Note, however, that the survival probability levels off to a finite value for $t\gg t^{*}$ due to inflation.

Figure 16

Data collapse for radial expansion simulations with a single green cell in an otherwise red initial homeland perimeter (with ${\mu }_f=0.1$, ${\mu }_b=0$, and $s=0.547$) on log-log plots. (a) The rescaled mean-squared cluster spread $X_a^2\left(t\right)$. The dashed line for $t\ll t^{*}$ indicates the directed percolation power-law scaling with exponent: $2{\nu }_{\perp }/{\nu }_{\parallel }\approx 1.3$. The dashed line in the inflation-dominated $t\gg t^{*}$ regime agrees with a $R_a^2\left(t\right)\sim t^2$ behavior discussed in the text. (b) Plot of the rescaled average number of green cells $N_a\left(t\right)$. The expected directed percolation scaling for $t\ll t^{*}$ is shown with a dashed line, with $({\nu }_{\perp }-2\beta )/{\nu }_{\parallel }\approx 0.31$. For $t\gg t^{*}$, we find $N_a\left(t\right)\sim t^{1-\alpha }\approx t^{0.84}$, consistent with arguments incorporating inflation given in the text. The simulation consisted of $9.6\times {10}^4$ independent runs for each homeland with radii $R_0=10,30,50,100,...,2000$, evolved to a final population of size $R\left(t\right)=2800$.

Figure 17

An inoculation of two bacterial (Escherichia coli) strains on a Petri dish, labeled with blue (dark gray in printed image) and yellow (light gray) constitutively expressed fluorescent proteins. The strains are otherwise identical and, consequently, represent the $s={\mu }_f={\mu }_b=0$ case in our analysis. The initial inoculation has the shape of a partial ring to allow for nutrients to penetrate the ring center and allow for the bacteria to grow inward, creating a “deflating” radial range expansion. A regular “inflating” radial range expansion moves away from the ring. This single experiment illustrates both kinds of radial expansion discussed in this paper.

Figure 18

(a) A radial range expansion with a deflating population front, for $s={\mu }_f={\mu }_b=0$, starting from an initial population around a circular rim of radius $R_0=300$. We use the same Bennett model lattice described in Sec. 2, while running the update scheme in reverse. (b) Deflationary range expansions for a single green (light gray) sector with an initial size ${\varphi }_0=0.1$ radians for a large ($s=0.4$) and very small ($s=0.01$) selective advantage. In both cases, deflation forces the green (light gray) sector boundaries to collide at the center (when $t=t^{*}=R_0/v$).

Figure 19

A comparison of the long-time survival probabilities of green sectors for small selective advantages $w\propto s$ in range expansions with radial inflation [$S_r(t\to \infty )$, solid lines], radial deflation [$S_d(t\to t^{*})$, dashed lines], and for linear range expansions [$S_l(t\to \infty )$, dotted lines]. The probabilities are calculated using Eqs. (45), (61), and a similar result for a linear range expansion. In the forward and reverse radial expansions, the initial population radius is $R_0=100a_r$, where the effective lattice constant $a_r=1$ for convenience. The linear expansion assumes a total population size $L=2\pi R_0=200\pi a_r$. We used two different initial green sector sizes, ${\varphi }_0=0.01$ (upper set of curves) and ${\varphi }_0=0.002$ (lower set of curves).

Figure 20

A schematic of the quasispecies model in which each cell in a population has a sequence distributed in a fitness landscape with a single sharp peak in green (light gray) surrounded by a flat landscape in red (dark gray). As cells divide, each individual reproduces with a mutational error rate $\mu$ that moves it around in sequence space (black, thin arrows). Below a critical “error threshold” $(\mu <{\mu }_c)$, a finite fraction of the population has the master sequence [i.e., is a green (light gray) cell]. Cells deviating from this ideal sequence are red (dark gray). When $\mu >{\mu }_c$, the population fitness distribution becomes delocalized, and the fitness peak is depopulated.

Figure 21

The speed $w\left(s\right)$ of a green sector boundary at the frontier with various selection advantages $s$, invading an otherwise red population at the frontier. The solid red line is the best fit of the speed to a linear function for $s\le 0.1$ [$w\left(s\right)\approx 1.1s$]. The initial angular sector size is ${\varphi }_0\approx \pi /2$, and the initial homeland radius $R_0=300$. $w\left(s\right)$ is calculated by averaging $(\langle \varphi \left(t\right)\rangle -{\varphi }_0)/\ln [R\left(t\right)/R_0]$ over approximately ${10}^3$ generations, where $\langle \varphi \left(t\right)\rangle$ is the angular size of the green sector at time $t$, averaged over about $3\times {10}^4$ simulation runs.