Nucleation of antagonistic organisms and cellular competitions on curved, inflating substrates

TL;DR

This paper investigates how curvature, inflation, and cellular interactions influence the survival and nucleation of antagonistic organisms on curved, inflating substrates, revealing critical sizes and the impact of motility.

Contribution

It introduces a model for the nucleation and survival of antagonistic organisms on curved, inflating surfaces, highlighting the effects of spatial diffusion, line tension, and curvature.

Findings

Critical cluster size depends on selective advantage and interaction strength.

Higher motility reduces survival probability when antagonistic interactions are present.

Curvature and inflation increase the survival chances of mutant cells at the frontier.

Abstract

We consider the dynamics of spatially-distributed, diffusing populations of organisms with antagonistic interactions. These interactions are found on many length scales, ranging from kilometer-scale animal range dynamics with selection against hybrids to micron-scale interactions between poison-secreting microbial populations. We find that the dynamical line tension at the interface between antagonistic organisms suppresses survival probabilities of small clonal clusters: the line tension introduces a critical cluster size that an organism with a selective advantage must achieve before deterministically spreading through the population. We calculate the survival probability as a function of selective advantage $\delta$ and antagonistic interaction strength $\sigma$. Unlike a simple Darwinian selective advantage, the survival probability depends strongly on the spatial diffusion constant $D_s$ of the strains when $\sigma>0$, with suppressed survival when both species are more motile. Finally, we study the survival probability of a single mutant cell at the frontier of a growing spherical cluster of cells, such as the surface of an avascular spherical tumor. Both the inflation and curvature of the frontier significantly enhance the survival probability by changing the critical size of the nucleating cell cluster.

Equations

$$\{\begin{array}{l}{\Gamma }_y={\Gamma }_y^0+\alpha (1-f)\\ {\Gamma }_b={\Gamma }_b^0+\beta f\end{array},$$

$${\partial }_{t}f=D_s{\nabla }^{2}f+\frac{f(1-f)}{⟨{\Gamma }^0⟩}\left[\frac{\delta }{2}+\sigma \left(2f-1\right)\right]+\bar{\eta },$$

$$⟨\bar{\eta }(\mathbf{x},t)\bar{\eta }({\mathbf{x}}^{\prime },t^{\prime })⟩=D_{g}f(1-f)\delta (\mathbf{x}-{\mathbf{x}}^{\prime })\delta (t-t^{\prime }),$$

$$P_{\infty }^s$$

$$\approx 1-\exp \left[-\frac{\delta }{D_g}\int {\mathrm{d}}^2\mathbf{x}f(\mathbf{x},t=0)\right]\approx \frac{A_0\delta }{D_g},$$

$${\partial }_t\varphi$$

$$⟨\eta (\mathbf{x},t)\eta ({\mathbf{x}}^{\prime },t)⟩$$

$$E[\varphi ]\equiv \int {\mathrm{d}}^2\mathbf{x}\left[\frac{1}{2}(\nabla \varphi {)}^2+\frac{V[\varphi ]}{D_s}\right].$$

$$\beta \equiv \frac{2D_s}{D_g},$$

$$S(\mathbf{k},t)=\sum _{\mathbf{x}}⟨\varphi (\mathbf{y},t)\varphi (\mathbf{y}+\mathbf{x},t){⟩}_{\mathbf{y}}{e}^{i\mathbf{k}\cdot \mathbf{x}},$$

$${\rho }_{int.}$$

$$\propto -{\nabla }^{2}S(\mathbf{x}=0,t)\propto \sum _{\mathbf{k}}|\mathbf{k}{|}^{2}S(\mathbf{k},t),$$

$$N_S(t)\equiv \frac{1}{N_{\mathbf{k}}}\sum _{\mathbf{k}}|\mathbf{k}{|}^{2}S(\mathbf{k},t)\propto {\rho }_{\mathrm{int}.}(t).$$

$${\nabla }^2\varphi \approx \frac{{\partial }^2{\varphi }^{\ast }}{\partial z^2}=\frac{1}{D_s}\frac{dV}{d\varphi }[\varphi ={\varphi }^{\ast }],$$

$${\varphi }^{\ast }(z)=-\tanh \left[\sqrt{\frac{\sigma }{4D_s}}z\right],z=r-r_d(t).$$

$$E[r_d]=2\pi \gamma r_d-\pi c{r}_d^2,$$

$$\frac{d{r}_d}{dt}=-\frac{D_s}{r_d}+\frac{\delta }{2}\sqrt{\frac{D_s}{\sigma }}.$$

$$P_{\infty }^s=\Omega e^{-\beta E[r_d^{\ast }]}=\Omega e^{-\frac{8\pi D_s\sigma }{3D_g\delta }},$$

$$P_{\infty }^s\propto \Gamma =\frac{\beta c}{2\pi }{e}^{-\beta \frac{\pi {\gamma }^2}{c}}\approx \frac{\delta }{3\pi D_g}{e}^{-\frac{8\pi D_s\sigma }{3D_g\delta }},$$

$$P_{\infty }^s=\frac{A_0\delta }{D_g}\exp \left[-\frac{8\pi D_s\sigma (1+\delta /{\delta }_0)}{3D_g\delta }\right],$$

$$R(t)=R_0(1+t/t^{\ast }),$$

$$P_{\infty }^s\approx \frac{2A_0}{t^{\ast }{D}_g}\exp \left[-\frac{4\pi t^{\ast }\left(2\sqrt{D_s\sigma }-t^{\ast }\delta \right)}{3D_g}\right],$$

$$P_{\infty }^s(\sigma =0,\delta ,t^{\ast })$$

$$\approx \left[\frac{e^{-t^{\ast }\delta /2}}{(t^{\ast }\delta /2{)}^{2}G(t^{\ast }\delta /2)}\right]\frac{\delta A_0}{D_g},$$

$$r_d^{\ast }(t^{\ast })=\left[\frac{2\pi t^{\ast }\delta }{9}{e}^{\frac{2\pi t^{\ast }\delta }{9}}{E}_{1/3}\left(\frac{2\pi t^{\ast }\delta }{9}\right)\right]\sqrt{\frac{4D_s\sigma }{{\delta }^2}},$$

$$\{\begin{array}{l}{\bar{r}}_d^{\ast }=R_0\arctan \left[\frac{r_d^{\ast }}{R_0}\right]<r_d^{\ast }\\ \bar{E}=2\pi R_0\left[\gamma \sin \left[\frac{{\bar{r}}_d}{R_0}\right]-2c{R}_0{\sin }^2\left[\frac{{\bar{r}}_d}{2R_0}\right]\right]\end{array},$$

$$P_{\infty }^s(\sigma ,\delta ,t^{\ast })\approx \frac{A_0}{D_g}[p_1(t^{\ast }\delta /2)]\delta \times$$

$$\exp \left[-\frac{4\pi R_0^2\delta }{3D_g}\left[\frac{\frac{4D_s\sigma }{(R_0\delta {)}^2}{p}_2(t^{\ast }\delta )+1}{\sqrt{1+\frac{4D_s\sigma }{(R_0\delta {)}^2}{\left[p_2(t^{\ast }\delta )\right]}^2}}-1\right]\right],$$

Full text

I Introduction

Antagonistic interactions between species, genotypes, or different organism types within a population are ubiquitous across the tree of life and span many length and time scales. For example, animal species such as the fire-bellied and yellow-bellied toad can interbreed and create selectively disadvantageous hybrids. The selection against hybrids can stabilize the ranges of these toads, creating narrow interfaces or “clines” between the species Szymura and Barton (1986). On much smaller length scales, many Gram-negative bacteria have a type-VI secretion system which allows them to inject nearby cells with “effector” molecules that can poison incompatible strains Russell et al. (2014). Genetically encoded toxin-antitoxin cassettes confer protection on the poisoners. Like the selection against hybrids, this antagonistic interaction stabilizes linear interfaces between spatial domains or circular boundaries surrounding genetically homogeneous patches of organisms McNally et al. (2017). For a review of social interactions between microorganisms in two dimensions, see Ref. Yanni et al. (2019).

Distinct clones within cell tissues may also have antagonistic interactions. For example, cancers have a particularly rich ecology Merlo et al. (2006) due to their heterogeneous genetic composition Sottoriva et al. (2015). Competitive interactions such as amensalism, in which a strain can inhibit the growth of nearby strains while remaining unaffected, is likely present amongst the heterogeneous cells of a tumor. Such an interaction may lead to the stabilization of distinct clones in solid tumors Marusyk and Polyak (2010). Understanding such ecological effects within the tumor is crucial for developing better cancer treatment options Korolev et al. (2014). Similar antagonistic effects are relevant for species with artificial gene constructs, such as the CRISPR-Cas9 technology introduced to contain certain pest populations McFarlane et al. (2018). In this case, the effective antagonistic interaction between the engineered and wild-type pest strains provides a mechanism for control of the wild-type population Tanaka et al. (2017). We model two-dimensional populations of two antagonistic strains, relevant for all of these scenarios. In this paper, we study not only flat, two-dimensional populations (as one might find in the animal range scenario), but also thin, effectively two-dimensional frontiers of growing spherical clusters of cells. The latter evolutionary dynamics has implications for antagonism near the frontiers of avascular tumors.

We focus on the dynamics of two species, labelled blue and yellow, with growth rates $\Gamma_{b}$ and $\Gamma_{y}$, respectively. Consider a monolayer or thin layer of individuals at carrying capacity, such that the fraction of yellow cells is $f$ and the fraction of blue is $1-f$. To model the antagonism, we suppose that the growth rates depend on the local fraction of the other species. Hence,

[TABLE]

where $\Gamma_{y,b}^{0}$ are base growth rates and the parameters $\alpha$ and $\beta$ are the amplitudes of the fraction-dependent contributions for the yellow and blue strain, respectively Korolev and Nelson (2011). We will take $\alpha,\beta\leq 0$, so that the two strains have mutually antagonistic interactions and suffer a growth penalty when either type grows next to the other. As discussed below, these interactions create a dynamical “line tension” $\gamma$ between cell types, as shown in Fig. 1, where yellow and blue strains demix more readily as we make $\alpha,\beta$ more negative. Here we use new coordinates $\delta=\alpha-\beta$ and $\sigma=-(\alpha+\beta)/2$. The line tension $\gamma$ is a monotonically increasing function of $\sigma$, $\gamma\equiv\gamma(\sigma)$, in the third quadrant of the phase diagram in Fig. 1 Lavrentovich and Nelson (2014).

The line tension $\gamma$ is generated by the spatial evolutionary dynamics between strains and describes a tendency of interfaces between domains of yellow and blue strains to minimize their lengths. As we shall see, we can interpret this tendency as an “energetic penalty” for creating a larger interface, analogous to a surface tension term in a free energy describing phase ordering. The presence of the dynamical line tension will allow us to make a connection between the evolutionary dynamics and phase-ordering dynamics of nucleation and growth. These dynamics become especially interesting on the curved, inflating frontiers of growing tissues such as microspheroidal tumors because the underlying space on which the nucleation and growth dynamics plays out is now curved and time-dependent! Also, note that unlike typical surface tensions, the dynamical line tension $\gamma$ does not depend on mechanical forces such as those driving microbial colony shapes Nguyen et al. (2004); Kayser et al. (2019).

We also discuss briefly the $\alpha,\beta>0$ case, i.e., the first quadrant of Fig. 1. Here, we have mutualism where each cell type benefits from the presence of the other type. In this case, there is a possibility of a “mutualistic phase” in which the two species remain mixed as the population evolves (shown schematically as a wedge centered on the line $\alpha=\beta$ in the first quadrant of the phase diagram in Fig. 1). In this case, the mutualistic phase has to combat number fluctuations embodied in genetic drift, which can drive the population locally to fixation of one or the other strain. Although mutualism prevails in the entire first quadrant within mean field theory, the region of stable mutualism is reduced to a narrow wedge by these fluctuations. This mutualistic regime has been observed in yeast range expansions Müller et al. (2014) and has been the subject of previous theoretical work Korolev and Nelson (2011); Lavrentovich and Nelson (2014).

We will focus on thin populations and monolayers of individuals in particular. This scenario could describe a thin layer of microorganisms on a Petri dish, a geographical range of animals, or the frontier of a growing cancerous cell mass. In the latter case, we assume the growth to be avascular, with the tumor receiving nutrients only via diffusion from the surrounding medium. Then, only the cells close to the curved frontier will divide and this interface becomes in effect a locally two-dimensional population. Such a growth modality has been observed in microspheroidal tumor cultures Zacharaki et al. (2004).

Consider the local fraction $f\equiv f(\mathbf{x},t)$ of yellow cells at some position $\mathbf{x}=(x,y)$ in the population at time $t$. To derive a spatial dynamics for our strains, we employ a “stepping stone model” Kimura and Weiss (1964) where we assume that at each position $\mathbf{x}$ the population can be modelled as well-mixed with some characteristic genetic drift strength $D_{g}\approx a^{2}/(\tau_{g}N_{\mathrm{eff}})$, where $N_{\mathrm{eff}}$ is a local population size and $a$ is a characteristic length over which we can treat the population as well-mixed. In our case, we will be concerned primarily with monolayers of cells for which $a$ is a cell diameter and $N_{\mathrm{eff}}\approx 1$. All of the dynamics associated with antagonism and selection embodied by the growth rates in Eq. (1) will occur within these well-mixed populations. Then, we assume that cells exhibit some small diffusive motion which leads to a cell exchange between adjacent positions. This means that the fraction $f$ will have a spatial diffusion with a characteristic spatial diffusion constant $D_{s}$.

Combining the well-mixed dynamics with the spatial diffusion leads to a simple time evolution for the spatial fraction when $\delta$ and $\sigma$ are small compared to the average base growth rate $\langle\Gamma^{0}\rangle\equiv(\Gamma_{y}^{0}+\Gamma_{b}^{0})/2$ :

[TABLE]

where $\sigma=-(\alpha+\beta)/2$ is an antagonistic interaction strength ($\alpha,\beta<0$), $\delta=\Gamma_{y}^{0}-\Gamma_{b}^{0}+\alpha-\beta$ is a selective advantage of the yellow cells, and $\bar{\eta}\equiv\bar{\eta}(\mathbf{x},t)$ is a stochastic noise due to individual cell births and deaths in each local well-mixed population. Since only $\delta$ is influenced by the difference in base growth rates, we will set, without loss of generality, $\Gamma_{b}^{0}=\Gamma_{y}^{0}=\tau_{g}^{-1}$ with $\tau_{g}$ a generation time. The nonlinear noise is interpreted in the Itô sense, has a zero average $\langle\bar{\eta}\rangle=0$, and has correlations bilinear in the organism fractions $f(\mathbf{x},t)$ and $1-f(\mathbf{x},t)$,

[TABLE]

where $D_{g}$ is the genetic drift strength mentioned previously. Details of such stochastic “stepping stone models” Kimura and Weiss (1964), with a focus on one spatial dimension, are reviewed in Ref. Korolev et al. (2010). We shall also see that the dynamics embodied in Eqs. (2) and (3) can be interpreted as a noisy evolution in an effective potential $V[f]$, which is sketched in Fig. 1 for various $\alpha$ and $\beta$ parameter regimes of interest. This potential has two local minima at $f=0$ and $f=1$. The dashed lines bounding the third quadrant in Fig. 1 are limits of metastability for $V[f]$.

As is typical of non-equilibrium systems, the initial condition for the populations plays a crucial role. Unlike equilibrium systems which evolve to a unique state fully specified by a free energy, the nature of the initial distribution of strains within the population can dramatically influence the subsequent evolutionary dynamics. In this paper we will focus on two kinds of initial conditions: (1) a random mixture of yellow and blue cells with equal concentrations [$\langle f(\mathbf{x},t=0)\rangle_{\mathbf{x}}=1/2$, where $\langle\cdot\rangle_{\mathbf{x}}$ represents a spatial average] which coarsens over time as illustrated in Fig. 1, creating genetic “sectors” whose typical size grows in time, and (2) a dilute suspension of yellow cells with a selective advantage $\delta>0$ in an otherwise all-blue population, as shown in Fig. 2(a). The first initial condition is analogous to spinodal decomposition Bray (2002), and would arise from an initial condition of genetic mixing in the mutualistic wedge in the first quadrant of Fig. 1 with $\alpha\approx\beta>0$, if the population were suddenly moved to an environment with $\alpha\approx\beta<0$.

In the first scenario, strains can compete with each other locally and the mixture coarsens over time. This behavior is typical in microbial studies where such an initial condition may be generated by inoculating, say, a well-mixed planar array of bacteria without flagella on a Petri dish McNally et al. (2017). We will first study the coarsening properties of the domains in this scenario. In this case, we assume cells that are killed off by antagonism give up their volume by bursting, with contents that are either absorbed into the agar gel, or provide nutrients to the antagonists. The phase diagram describing the long-time dynamics as a function of $\sigma$ and $\delta$ for a random mixture of blue and yellow cells in equal proportions is sketched in Fig. 1. When $\delta>0$, the yellow cells will eventually sweep the population. When $\delta<0$, the blue cells are victorious. We will study the coarsening dynamics near the line $\delta=0$, and study the effect of increasing $\sigma$.

In the second scenario, yellow cells (representing, e.g., rare mutants within the population), if they survive genetic drift, can out-compete neighboring blue cells when their selective advantage $\delta>0$. If the resultant yellow cluster reaches a certain critical cluster size, it can spread through the blue population, as shown in Fig. 2(a). We will connect such a dynamics to a nucleation and growth process where the blue population represents a metastable state and the yellow cells correspond to a stable, nucleating minority phase. This connection is possible because Eq. (2), apart from the noise correlations in Eq. (3), has an identical form to the stochastic partial differential equation describing non-conserved, model A dynamics at a first-order phase transition Bray (2002); Hohenberg and Halperin (1977). We shall now briefly comment on the difference in the noise.

In the typical (model A) nucleation and growth scenario, such as in a kinetic Ising model in two dimensions, or near a two-dimensional liquid-vapor phase transition, the noise comes from thermal fluctuations throughout the system and has a simple correlation: The nonlinearity $f(1-f)$ in Eq. (3) is replaced by a constant proportional to the temperature. In contrast, for our antagonistic evolutionary dynamics, the noise has nonlinear correlations and acts only close to interfaces between “phases” of $f=1$ (yellow cells) and $f=0$ (blue cells). However, when interfaces between the two phases are sharp (the so-called “thin-wall limit”), thermal noise in model A primarily acts to distort the position of the interface: the low-energy modes excited by thermal noise are translations or undulations of the nucleating droplet boundary Voloshin (1985); Günther et al. (1980), as shown schematically in Fig. 3(b). In one spatial dimension, it can be shown rigorously that the interface dynamics behaves like a random walk biased by the asymmetry in potential well depth Brassesco and Buttà (1998). Hence, it seems reasonable to conjecture that the model A scenario and the antagonistic dynamics are analogous in the thin-wall limit for single “droplets” of the stable phase (yellow cells) surrounded by the metastable phase (blue cells). Specifically, the dynamics described by Eq. (2) should reduce, in this thin-wall case, to the dynamics of a fluctuating, sharp interface. Thus, for the problems of interest to us here, it seems plausible that we can evaluate the nonlinear noise correlations in Eq. (3) at the middle of the interface where $f=1/2$, a conjecture we check with our numerical simulations below.

It is also worth noting that in the case $\sigma=0$ and variable selective advantage $\delta$, the survival probability $P^{s}_{\infty}$ at long times $(t\rightarrow\infty)$ of a single yellow cell in a sea of blue is given by the celebrated Kimura formula, which, although originally derived for well-mixed populations, is remarkably independent of the spatial structure of the population under rather general conditions Maruyama (1974). For a single cell in a population of $N\gg 1$ blue cells, this formula reads:

[TABLE]

where $\tau_{g}$ is the generation time. The first line is the original formula of Kimura Crow and Kimura (1970), and the second line is its generalization to two dimensions with $\sigma=0$ (no antagonism). The first approximate equality only holds for $N$ so large that $\tau_{g}N\delta\gg 1$. Note the remarkable independence of this result on the precise spatial arrangement of yellow cells at $t=0$, $f(\mathbf{x},0)$. The initial area occupied in the case of an isolated yellow cluster is $\int\mathrm{d}^{2}\mathbf{x}\,f(\mathbf{x},t=0)=A_{0}$ ($A_{0}\approx a^{2}$ for a single cell with diameter $a$). The last approximate equality in Eq. (4) requires a small selective advantage, $\delta\ll D_{g}/A_{0}$, which will be our regime of interest. Because $D_{g}\sim a^{2}/(\tau_{g}N_{\mathrm{eff}})$, we require $\delta\ll\tau_{g}^{-1}$ for the special case of a single mutant cell and $N_{\mathrm{eff}}=1$. Note also that this formula is independent of the spatial diffusion coefficient $D_{s}$. This is a general phenomenon and has been checked for a wide range of systems Maruyama (1974); Doering et al. (2003); Pigolotti et al. (2013).

We plot Eq. (4) on one of the panels of Fig. 2(b), where it clearly fits our simulation data (described below) for $\sigma=0$ (with $\tau_{g}D_{g}/A_{0}=1.00\pm 0.01$ as a fitting parameter). We shall see that when we include an antagonistic interaction, $\sigma>0$, the survival probability of a single mutant cell (or small patch of cells) as a function of selective advantage $\delta$ is strongly suppressed, as illustrated in Fig. 2(b). Moreover, the formula for $P^{s}_{\infty}$ is qualitatively different for $\sigma>0$ and depends explicitly on the spatial diffusion constant $D_{s}$ and the spatial structure of the population.

The paper is organized as follows: In Section II, we describe our simulations and analytic approaches. In Section III, we consider how domains of yellow and blue cells coarsen over time for $\delta=0$ and $\sigma\geq 0$, i.e., right along the line of first-order transitions in the third quadrant of Fig. 1. We show that increasing $\sigma$ increases the line tension between yellow and blue cell clusters, as illustrated in the population snapshots in Fig. 1. In the Section IV, we analyze the probability that a single yellow cell surrounded by blue cells is able to create a cluster that survives and eventually expands at long times. In Section V, we extend our analysis of the survival probability to range expansions with inflating spherical frontiers, such as the surface of an avascular solid tumor. We conclude with some discussion of our results and prospects for future work in Section VI.

II Simulations and Analysis

To begin, it is convenient to make some connections to models describing the relaxation of an (Ising-like) magnetic spin system to equilibrium. To do this, we change variables from $f$ to $\phi=2f-1$, so that the yellow and blue strains correspond to $\phi=\pm 1$, respectively. With this variable change, the equation of motion in Eq. (2) is analogous to the dynamics of an Ising-like system well below the critical temperature $T_{c}$, with a coarse-grained “magnetism” $\phi$ evolving in a potential $V[\phi]$ which exhibits two minima at $\phi=\pm 1$ (e.g., the two magnetizations of an Ising spin) in such a way that the total magnetization of the system is non-conserved (model A) Hohenberg and Halperin (1977). Note that, although competing cell populations are clearly a form of “active matter” Marchetti et al. (2013), the dynamics of motile bacterial and Janus particles that motivate many of these studies do not involve cell divisions, and hence are better described by variants of model B dynamics, with a conserved order parameter Hohenberg and Halperin (1977). In our evolutionary dynamics, model A time evolution means that the total fraction of yellow and blue cells is not conserved, and instead varies in time. The equation of motion in Eq. (2) is now conveniently rewritten as

[TABLE]

where $\eta=2\bar{\eta}$ is again a nonlinear noise and $V[\phi]=(\phi+1)^{2}\left[2\delta(\phi-2)+3\sigma(\phi-1)^{2}\right]/24$ has a double-well structure for $\sigma>\delta/2>0$, as shown in Fig. 3(a) and in Fig. 1. In this regime, the blue state becomes metastable. The double-well structure of the potential allows us to naturally connect the antagonistic dynamics to nucleation and growth, a connection that was made previously in the ecological context by Rouhani and Barton Rouhani and Barton (1987) and in the context of gene drives Tanaka et al. (2017). Although it is tempting to interpret the limit $\sigma\rightarrow 0$ in Fig. 3 in terms of an Ising model in a magnetic field just below its critical point, this is precisely where the unusual noise correlations of Eq. (3) become important, and produce the distinctly non-Ising model features in the first quadrant of the “phase diagram” shown in Fig. 1.

How shall we deal with the peculiar noise correlations $\langle\eta(\mathbf{x},t)\eta(\mathbf{x}^{\prime},t^{\prime})\rangle=D_{g}(1-\phi^{2})\,\delta(\mathbf{x}-\mathbf{x}^{\prime})\delta(t-t^{\prime})$ in Eq. (6), which differs from conventional model A dynamics Hohenberg and Halperin (1977), even deep in the third quadrant? Although the noise $\eta$ here is manifestly non-thermal and provided by the intrinsic birth-death dynamics of the individual cells, it will nevertheless be helpful for us to think about Eq. (5) as the dynamics of a field $\phi$ as it relaxes to a minimum of the following “free energy”:

[TABLE]

The field $\eta(\mathbf{x},t)$, then, can be thought of as a thermal noise that acts only at interfaces between domains of $\phi=\pm 1$. At these interfaces, we have $\phi\approx 0$ and the noise correlations in Eq. (6) are approximately independent of $\phi$. Here, we can connect the noise to an inverse temperature via a formula analogous to the fluctuation-dissipation theorem:

[TABLE]

so that near domain interfaces, we have $\langle\eta(\mathbf{x},t)\eta(\mathbf{x}^{\prime},t^{\prime})\rangle\approx 2\beta^{-1}D_{s}\delta(\mathbf{x}-\mathbf{x}^{\prime})\delta(t-t^{\prime})$ Hohenberg and Halperin (1977); Täuber (2014). We shall test this approximation further in Section IV. We may now draw analogies between relaxation dynamics of a magnetic spin system and the evolutionary dynamics.

We also implement simulations of monolayers of cells, growing either in a two-dimensional triangular lattice or at the edge of an inflating spherical clump of cells. The simulation approaches will be similar to the ones employed in Refs. Lavrentovich and Nelson (2014, 2015), which we review briefly here. We evolve cells one cell division at a time, either in a triangular lattice or on the surface of a spherical cluster of cells generated using the Bennett hard sphere packing algorithm Bennett (1972). During each division, three adjacent cells compete to divide into the new spot, as shown in Fig. 4(a). It is convenient to parameterize as follows in order to have the simulations best approximate our choice of parameters in Eqs. (1,2): If two cells are blue and one is yellow, then the probability the cell is yellow is $p_{y}=1/3+\bar{\alpha}$, with $\bar{\alpha}\equiv\delta/9-2\sigma/27$. Otherwise, if we have two yellow and one blue cell, $p_{y}=2/3-\bar{\beta}$, where $\bar{\beta}\equiv-\delta/9-2\sigma/27$, where $\delta$ and $\sigma$ are in units of $\tau_{g}^{-1}$. All other possibilities follow from probability conservation. Because we assume no spontaneous mutations on the time scale of our simulations, the daughter cell that appears in a pocket formed by three adjacent cells of the same color will have the same color. Some examples are illustrated in Fig. 4(a). For the triangular lattice, we keep dividing cells to fill in a staggered copy of the lattice, representing the next generation of cells. For the spherical cluster growth, a generation time $\tau_{g}$ passes when the total spherical cluster grows by a full cell diameter $a$ .

An important quantity we use to characterize the coarsening dynamics in the next section is the structure factor, defined as

[TABLE]

with $\mathbf{k}=(k_{x},k_{y})$ a 2-dimensional wavevector and $\langle\ldots\rangle_{\mathbf{y}}$ is a spatial average and an average over many simulation runs. For the triangular lattice with each cell a diameter $a$ away from its six nearest neighbors, $\mathbf{k}$ ranges over a hexagonal Brillouin zone as shown in Fig. 4(b). If we have an initially isotropic distribution of organism types, we then expect that $S(\mathbf{k},t)$ will also be isotropic during the evolution of the population in time and hence only dependent on $k\equiv|\mathbf{k}|$ at least close to the origin, where we expect details of the lattice structure to drop out. We may use $S(k,t)$ to estimate the characteristic size $\xi(t)$ of strain domains at time $t$. As $S(k,t)$ is expected to evolve as in a system with non-conserved dynamics, we anticipate that it develops a peak at $k=0$ [see Fig. 4(b)] which will grow and sharpen over time Bray (2002). Hence, we can estimate a domain size $\xi(t)$ from the half-width at half maximum of the structure factor peak at $k=0$, $k_{\mathrm{HWHM}}$, where $k_{\mathrm{HWHM}}(t)\approx 2\pi/\xi(t)$.

We can also estimate the average interface density $\rho_{\mathrm{int.}}$ for a population of $N$ cells using the structure factor. If the cells are close-packed on some lattice with cell-diameter spacing $a$, then the density of interfaces $\rho_{\mathrm{int}.}$ at time $t$ is given by

[TABLE]

where $\langle\ldots\rangle$ is an average over many simulation runs and we sum over all cells at positions $\mathbf{x}$ and their $z$ nearest neighbors at relative positions $\bm{\delta}$. Note that the summand in the first line of Eq. (10) is only non-zero when $\bm{\delta}$ spans an interface between organism types. In the second line, we take the continuum limit $|\bm{\delta}|\rightarrow 0$ and use a finite-difference approximation to the Laplacian, which can be written in Fourier space as a sum over wave-vectors $\mathbf{k}$ in the Brillouin zone (BZ) in Fig. 4(b). A convenient quantity, then, measurable in either experiments or simulations in which we uniformly sample $N_{\mathbf{k}}$ different wave-vectors (over the Brillouin zone, say) reads:

[TABLE]

The constant of proportionality will depend on the particular lattice or microscopic implementation of the simulated evolution. We now proceed to the results of our simulation and analytic approaches.

III Domain Coarsening

We first consider a well-mixed initial condition with equal proportions of blue and yellow organisms and ask how the domains of blue and yellow coarsen over time with no selective advantages $(\delta=0)$ but variable antagonism $\sigma$ (as in Fig. 1). The structure factor is indeed isotropic at long wavelengths, as shown in Fig. 4(b). When $\sigma=0$, the dynamics of the system reduce to the voter model Lavrentovich and Nelson (2014), at which we get a peculiar “line-tension-less” coarsening. In fact, the full structure factor for the voter model can be calculated exactly and matches our results for $\sigma=0$ Scheucher and Spohn (1988). In this case, although the largest domain size in the system grows as $\sqrt{t}$, the density of interfaces decays with a much different (and slower) $1/\ln t$ behavior. The latter slow decay is due to the presence of very diffuse domains in the system. The behavior of $S(k,t)$ at larger wavevectors $k$ probes the interface topology. Specifically, for sharp interfaces, we expect according to Porod’s law that $S(k,t)\sim 1/k^{3}$ for large $k$ (for a two-dimensional system) Bray (2002). This is indeed what happens for $\sigma>0$ when we have antagonistic interactions, as shown in Fig. 5(b). In the voter model case $\sigma=0$ in Fig. 5(a), we see that $S(k,t)\sim 1/k^{2}$ for large $k$, instead, consistent with the absence of sharp interfaces between domains and suggesting that the contribution to $S(k,t)$ for large $k$ is dominated by small, isolated cell clusters, as can be seen in the $\sigma=0$ panel of Fig. 1.

To detect this change in the character of the domain interfaces between $\sigma=0$ and $\sigma>0$, we look at both the interface density $\rho_{\mathrm{int.}}$ and the typical domain size $\xi$ [as judged from $N_{S}(t)$ in Eq. (11) and the inverse of the half-width at half-maximum $k_{\mathrm{HWHM}}$ of the structure factor peak at $k=0$, respectively] as functions of time $t$. We observe in Fig. 6 the marked difference between these two quantities. The characteristic domain size $\xi(t)$ always increases as $\xi\propto\sqrt{t}$, with the constant of proportionality only weakly dependent on the antagonistic interaction strength $\sigma$, as indicated by the substantial overlap between the various lines in Fig. 6(a). In sharp contrast, the interface density $\rho_{\mathrm{int.}}$ decay depends very strongly on $\sigma$. As $\sigma$ is increased, the line tension at the boundary increases and the interfaces between the domains sharpen. This leads to a transition from a $1/\ln t$ decay for $\sigma=0$ to a much faster $1/\sqrt{t}$ decay for $\sigma>0$, as shown in Fig. 6(b). This increasing line tension with $\sigma$ also has a pronounced effect on the survival probability of organism genotypes, as we will see in the next section.

IV Single Cluster Dynamics

Let us now consider the fate of a single yellow cell in a (very large) blue cell population, representing some rare mutation, for example. When $\sigma=0$ and the yellow cell has a small selective advantage $\delta>0$, then the probability that it sweeps the population is given by the formula in Eq. (4). However, as soon as $\sigma>0$ and there are antagonistic interactions, the effective line tension discussed in the previous section strongly suppresses the survival probability, as is evident in Fig. 2(b). This suppression can be understood by noting that for $\sigma>\delta/2$, the potential shown in Fig. 3 develops a barrier for transitioning between the blue state ($\phi=-1$) and the more favorable yellow type $(\phi=1)$. Similar to conventional nucleation theory, the yellow cell must eventually generate a sufficiently large cluster of offspring to overcome this barrier.

We first consider the deterministic dynamics of a yellow cell cluster for $\sigma>\delta/2$. That is, consider a “droplet” solution to the dynamical equation Eq. (5) in the absence of noise $(\eta=0)$. If the droplet interface width $\sim\sqrt{D_{s}/\sigma}$ is much smaller than the droplet radius $r_{d}(t)$ (see Fig. 3), then we can look for solutions $\phi(r,t)\equiv\phi^{*}[r-r_{d}(t)]$, in which the cluster maintains a stationary interface profile $\phi^{*}(z)$ and has a time-dependent radius $r_{d}(t)$. Under the assumption that the interface profile does not change as the yellow cell cluster expands or shrinks, the profile $\phi^{*}(z)$ can be estimated from the stationary solution to Eq. (5):

[TABLE]

which, under the appropriate boundary conditions [$\phi(z\rightarrow-\infty)=1$ and $\phi(z\rightarrow+\infty)=-1$] is solved by

[TABLE]

This solution holds provided the variation of $\phi^{*}(z)$ is localized around $z=r-r_{d}=0$ and the cluster radius $r_{d}$ is large compared to the interfacial width: $r_{d}\gg\sqrt{D_{s}/\sigma}$. In this thin-wall limit, we can use Eq. (13) to recast our free energy in Eq. (7) as an energy for a disc-shaped cluster

[TABLE]

with a line tension $\gamma=(2/3)\sqrt{\sigma/D_{s}}$ and a condensation energy $c=\delta/(3D_{s})$. Moreover, the dynamics in Eq. (5) becomes Bray (2002)

[TABLE]

There are two important aspects of the Eq. (15) dynamics. First, when $\delta=0$ (so that the competing organisms are selectively neutral), the droplet nevertheless shrinks and we get the solution $r_{d}(t)=\sqrt{(r_{d}^{0})^{2}-2D_{s}t}$, with $r_{d}^{0}$ the initial cluster radius. This solution suggests that curved clusters will have radii that shrink as $\sqrt{D_{s}t}$ as we would expect from the non-conserved, Ising-like Allen-Cahn (model A) dynamics Allen and Cahn (1979). Note that the dynamics is independent of $\sigma$, consistent with the results for the domain size $\xi(t)$ we saw in Fig. 6(a). Second, there is a stationary radius $r_{d}^{*}$ at which $dr_{d}/dt=0$. This is the critical droplet size: If the initial $r_{d}(t=0)$ is larger than this size, then the two-dimensional droplet grows systematically larger with time, with eventually $r_{d}(t)\approx(\delta/2)\sqrt{D_{s}/\sigma}\,t$, where the interface can be regarded as a pushed genetic wave Szymura and Barton (1986); Tanaka et al. (2017). Otherwise, the droplet shrinks. Upon setting $dr_{d}/dt=0$, we find $r_{d}^{*}=\gamma/c=\sqrt{4D_{s}\sigma/\delta^{2}}$. Therefore, in order for a yellow cell cluster to survive, genetic drift [i.e., the fluctuations embodied in the interfacial noise of Eq. (6)] must overcome the line tension to inflate the droplet to the critical size $r_{d}^{*}$, after which it may grow deterministically. This is the rate-limiting step in the nucleation and growth process for an organism with selective advantage $\delta$, and we conjecture that the probability of this step for this nonequilibrium dynamical system, i.e. the survival probability, is given by an Arrhenius factor:

[TABLE]

where $E[r_{d}^{*}]=4\pi\sigma/(3\delta)$ is the energy of the critical droplet, $\beta$ is the inverse “temperature” in Eq. (8), and $\Omega$ is a more subtle prefactor which we now describe.

The prefactor $\Omega$ will have two contributions. First, note that there is nothing in Eq. (16) that specifies an initial condition. In the usual (homogeneous) nucleation and growth dynamics, thermal fluctuations are responsible for initially nucleating the stable state, so the initial condition could be just the uniform metastable (blue) phase. On the other hand, for our evolutionary dynamics with the nonlinear noise that vanishes away from interfaces, all nucleating yellow cell clusters must come from yellow cells already present in the initial condition. Hence, the specific initial condition (a single yellow cell, say), will influence the prefactor $\Omega$. Secondly, it is known that shape fluctuations around the critical droplet shape (reminiscent of entropy effects in equilibrium systems) will contribute to $\Omega$ Langer (1969); Voloshin (1985). In the thin-wall limit, these fluctuations are undulations of the interface position as shown schematically in Fig. 3(b).

In the thin-wall limit $r_{d}\gg\sqrt{D_{s}/\sigma}$, the critical droplet can be thought of as a closed loop with a spectrum of fluctuations which can be derived from the “free energy” in Eq. (7). The fluctuations of this loop are driven by thermal noise in the case of conventional nucleation and growth. For our evolutionary dynamics, the interface shape fluctuates due to the nonlinear noise in Eq. (3). However, because the dynamics at the interface dominates in this thin-wall limit, we expect that the usual nucleation and growth analysis will work well, provided we evaluate the noise at the middle of the interface where $\phi=0$ (or $f=1/2$). Another possibility, similar to the approach for deriving Eq. (15) from Eq. (5), is to average the nonlinear noise correlations over the stationary interface profile $\phi^{*}$ given by Eq. (13). This alternative procedure does not change our results in any substantial way (compared to evaluating the noise at $\phi=0$) and amounts to a simple rescaling of the genetic diffusion constant: $D_{g}\rightarrow 4D_{g}/5$. As we use $D_{g}$ as a fitting parameter when we check our analysis with simulations, the two approaches give identical results.

Once we approximate the two-dimensional droplet of genetic material with a selective advantage $\delta$ as a fluctuating loop, we can exploit the results of Voloshin, Voloshin (1985) who showed that in the thermal system (i.e., an Ising model below $T_{c}$), the nucleation rate (including the prefactor and the Arrhenius term) in two dimensions is completely determined by the temperature, and the fluctuation-renormalized condensation energy $c$ and line tension $\gamma$. The result for the nucleation rate $\Gamma$, which is necessarily proportional to the survival probability $P_{\infty}^{s}$ in our problem, is

[TABLE]

where we have substituted in for the inverse temperature $\beta$ using Eq. (8) and used the results for the line tension and condensation energy $\gamma=(2/3)\sqrt{\sigma/D_{s}}$ and $c=\delta/(3D_{s})$, respectively. As $\Gamma$ is a probability per unit area, the constant of proportionality must be an area related to the initial condition for our population genetics problem, which is in this case a single cell. We thus expect that this factor is some constant close to $a^{2}$. As we shall see, this assumption allows standard formulas from nucleation theory to be matched smoothly onto the Kimura results for fixation probabilities when $\sigma=\gamma=0$, and leads to an understanding of fixation probabilities in a much larger domain.

Although for a selective advantage such that $\delta>2\sigma$ we no longer have a potential barrier and a critical droplet size [recall that $|\delta|=2\sigma$ defines the limit of metastability for the potential energy displayed in Fig. 3(a)], we can still use the $\sigma=0$ limit to set the constant of proportionality between $\Gamma$ and $P_{\infty}^{s}$. In addition, we expect corrections to the nucleation formula in Eq. (17) in the exponential term, correcting the effective free energy associated with the critical droplet (to account for deviations from the thin wall limit, possible pinching off of the interface, deviations from the $\delta\ll\sigma$ limit, etc.) Voloshin (1985). With these considerations in mind, we conjecture the following formula for the survival probability of a single yellow cell in a two-dimensional population of blue cells:

[TABLE]

where $\delta/\delta_{0}$ is a correction term that takes into account behavior away from the $\delta\ll\sigma$ limit. Note the important dependence on the spatial diffusion constant $D_{s}$ when $\sigma>0$, in contrast to the Kimura result for a single cell in Eq. (4), which is independent of $D_{s}$.

We now test this result via computer simulations. The $D_{g}/A_{0}$ part of the prefactor can be set by fitting to $\sigma=0$ data; we find $D_{g}/A_{0}=1.00\pm 0.01$ in units of inverse generation time $\tau_{g}^{-1}$. This result is plausible since both $\tau_{g}D_{g}$ and $A_{0}$ are expected to be of order $a^{2}$. We then determine just two dimensionless fitting parameters: $D_{s}/D_{g}=0.0280\pm 0.0002$ and $\tau_{g}\delta_{0}=0.054\pm 0.002$, which can be found by fitting $P_{\infty}^{s}$ as a function of both $\sigma$ and $\delta$ for the data with $\sigma/\delta>1/2$. The fit indicates that the spatial diffusion coefficient $D_{s}$ is small compared to the constant $D_{g}$ controlling genetic drift, consistent with the small effective population size embodied in our simulation approach. The cells we model on a triangular simulation lattice only diffuse via cell division displacements and do not have an independent motility. With just two fitting parameters, we find good agreement between theory and simulation for $P_{\infty}^{s}$ as a function of both $\sigma$ and $\delta$, as shown in Fig. 7. Note that apart from the $\delta/\delta_{0}$ correction, we get excellent data collapse of $P_{\infty}^{s}e^{\sigma/\sigma_{0}}/\delta$ [where $\sigma_{0}\equiv 3D_{g}\delta_{0}/(8\pi D_{s})$] versus $\sigma/\delta$, as shown in Fig. 7(c), where we also indicate the best fit.

Our fit to nucleation theory, shown in Fig. 7(b,c), reveals that the survival probability, apart from some corrections to scaling at large $\delta$, decays exponentially with the ratio $\sigma/\delta$. This prediction is verified in the simulations. To get this good agreement, it was necessary to run the simulations until either the yellow cell cluster died out completely or spread through a large enough fraction of the population so that ultimate fixation is inevitable; this fraction was typically taken to be a quarter or a half of our simulation area. These simulation times can be quite long because the survival probability $P^{s}(t)$ typically decays as a power-law in time $t$. This slow decay can be understood for small $\sigma$ and $\delta$ from the behavior at the voter model point $\sigma=\delta=0$. Here, standard results Bramson and Griffeath (1980) show that the survival probability in two spatial dimensions for a yellow cell cluster spawned from a single initial yellow cell decays according to $P^{s}(t)\sim\ln t/(\pi t)$, consistent with our simulation results from early to intermediate times before $P^{s}(t)$ saturates to its limiting value of $\lim_{t\rightarrow\infty}P^{s}(t)=P_{\infty}^{s}$ at long times. Some examples of the behavior of $P^{s}(t)$ as a function of time $t$ are shown in Fig. 8(c).

An interesting aspect of the evolutionary dynamics studied here is the fixation time, i.e., the time necessary for the yellow cell cluster to definitively sweep the population. A rough estimate of the fixation time $t_{f}$ for a surviving cluster is $t_{f}=t_{\mathrm{nucl}}+(L-r_{d}^{*})/v_{d}$, where $L$ is the linear extent of the system, $v_{d}=(\delta/2)\sqrt{D_{s}/\sigma}$ is the yellow droplet growth velocity (after it becomes larger than the critical size), and $t_{\mathrm{nucl}}$ is the “nucleation time” it takes for the yellow cell cluster to grow from a single cell to the critical size. The quantity $t_{\mathrm{nucl}}$ is challenging to estimate because, unlike the late-time droplet growth governed by the deterministic velocity $v_{d}$, the early-time yellow cluster dynamics are strongly influenced by the nonlinear noise $\eta$ in Eq. (5) that is ultimately responsible for overcoming the line tension induced by the antagonism $\sigma>0$. A proper analysis would examine the evolution of the initial condition to the critical droplet shape under the stochastic evolution of Eq. (5). Here we will only make the approximate estimate given below.

The motion of the yellow-blue interface in the presence of noise requires a delicate analysis, and various regimes are possible depending on the magnitude of the genetic drift $D_{g}$ Panja (2004). We can make a crude estimate by assuming that a surviving subcritical yellow cluster propagates into the blue territory as a very noisy Fisher wave, driven by the genetic drift, spatial diffusion, and selective advantage $\delta>0$. The velocity of such a wave scales according to $v_{\mathrm{F}}\sim D_{s}\sqrt{\delta/D_{g}}$ Hallatschek and Korolev (2009), which yields an estimated $t_{\mathrm{nucl}}\sim r_{d}^{*}/v_{\mathrm{F}}\sim\sqrt{D_{g}\sigma/(D_{s}\delta)}\,/\delta$. This approximation yields, for example, a reasonable $t_{\mathrm{nucl}}\approx 500$ generation times for the evolution in Fig. 2(a). However, the antagonism $\sigma$ will surely modify the growth speed of the droplet and a proper analysis of the early-time dynamics would take into account the metastability for $\sigma>\delta/2$. The interplay of the antagonism and the genetic drift noise is subtle; understanding the detailed interface motion is an active area of current research, even in one spatial dimension Birzu et al. (2018). A more detailed analysis of the fixation time is outside of the scope of this paper. We will now analyze new features for survival probabilities on inflating spherical frontiers.

V Survival at Inflating Frontiers

Spatial evolutionary dynamics can also play out on the surface of growing tissue or at the edge of a three-dimensional microbial colony. One such example is the surface of an avascular tumor Lavrentovich and Nelson (2015). Tumors may contain a heterogeneous distribution of strains which prefer different microenvironments and have a complicated ecology Merlo et al. (2006); Marusyk and Polyak (2010); Korolev et al. (2014). It is thus important to consider evolutionary interactions between cancer strains that go beyond a simple selective advantage. One possibility is to include a frequency-dependent selection, as we do for our growth rates in Eq. (1). Such interactions are increasingly thought to be relevant for cancer strains Korolev et al. (2014). For example, evolutionary game theory approaches suggest that a frequency-dependent selection may describe the effect of the exchange of diffusible goods, relevant for cancer ecology Altrock et al. (2015). Thus, it is interesting to consider our antagonistic interactions in this context, as they represent a regime of frequency-dependent selection (the third quadrant in Fig. 1).

With this motivation, consider a spherical clump of cells (e.g., a solid tumor) with a radius $R(t)$ that grows linearly in time:

[TABLE]

where $R_{0}$ is the initial radius and $t^{*}=R_{0}/v$ is a time scale associated with the speed $v$ of the growing tumor frontier. This simple microspheroid model captures many of the essential features of real tumors, especially in the early avascular stages of tumor development Kunz-Schughart (1999). We will assume for simplicity that the frontier of the spherical tumor remains smooth, and neglect the interesting enhancement of genetic boundary wandering associated with rough interfaces Hallatschek et al. (2007); Lavrentovich and Nelson (2014). In simulations, this constraint can be implemented by requiring the cell closest to the center of the spherical cluster to divide first. This condition mimics an effective “surface tension” for the mass of cells, a plausible assumption since cell masses of microbes and cell tissues are known to behave as (viscoelastic) fluids under many circumstances Nguyen et al. (2004); Foty et al. (1994). We assume that nutrients are available primarily to cells at the frontier of the growing three-dimensional tumor, so that only these cells can undergo cell divisions.

We again consider a single mutant yellow cell appearing at the frontier of a blue spherical cluster with an initial radius $R_{0}$. The yellow variant grows with some selective advantage $\delta$ over the surrounding (blue) tumor cells, as would be the case for a “driver” mutation at a two-dimensional tumor surface Pon and Marra (2015). This mutated strain might prefer its own special microenvironment and try to poison or attack its competitors, who may, in turn, produce secretions harmful to the mutant. We thus assume the curved two-dimensional blue/yellow cell interface is associated with some antagnostic interaction with a characteristic $\sigma$. An important question from the point of view of the tumor ecology is: What is the probability that the yellow cell is able to sweep the population? As illustrated in Fig. 8(a,b), this quantity will depend strongly on $\sigma$, as discussed above for the flat two-dimensional environments. However, in this curved setting with growth, inflation at the tumor frontier enhances the survival probability $P^{s}$, as investigated previously with theory and simulations for the $\sigma=0$ case as a function of $\delta$ Lavrentovich and Nelson (2015).

When $\sigma>0$, we again expect a nucleation and growth scenario, with a line tension on the curved two-dimensional tumor surface playing a key role. However, there are two important modifications to the theory in a flat environment: First, we have the inflationary effect due to the growing spherical frontier, with $R(t)$ increasing linearly with time according to Eq. (19). This factor will be especially important at late times $t\gg t^{*}$. Second, the curvature of the frontier will also play a role, especially if the two-dimensional critical radius $r_{d}^{*}$ of a curved yellow domain at the frontier is comparable to $R_{0}$ Gómez et al. (2015); Horsley et al. (2018). A more detailed discussion of these factors is given in Appendix A.

The simulation results for the survival probability $P_{s}$ for a single yellow cell on an initially all blue spherical cluster of cells (after about 200 generations) is shown in Fig. 9. Note how smaller initial radii $R_{0}$ yield larger survival probabilities, as expected due to the increased inflationary effect. We also see a decay of the probability with increasing $\sigma$ and decreasing $\delta$ in Fig. 9. It is possible to derive a theoretical prediction that has similar trends. This is done in detail in Appendix A. A key feature of the calculation is that the spatial diffusion coefficient plays a significant role, unlike the analysis for a simple selective advantage Lavrentovich and Nelson (2015). We find that for $t^{*}\delta\ll 1$, the long-time survival probability $P_{\infty}^{s}\equiv P_{\infty}^{s}(\sigma,\delta,t^{*})$ is given by

[TABLE]

where $A_{0}\approx a^{2}$ is the initial area occupied by the yellow cell. A detailed comparison of our theoretical prediction with simulation results, along with a more detailed formula valid for larger values of $t^{*}\delta$ is given in Appendix A.

As discussed in more detail in the Appendix, our theoretical prediction captures the essential features of the simulated survival probability, even at finite times. We see the strong enhancement due to both the inflationary effect and the curved spherical surface. To find better approximations, it would be necessary to study in more detail the interplay between the evolutionary dynamics and the inflation. It would also be helpful to consider the time-dependence of the survival probability. This would be an interesting topic for future research.

VI Conclusion

We have explored how antagonism plays a crucial role in the spatial evolutionary dynamics of a population, leading to effects unanticipated in a simple theory of natural selection. As originally appreciated by Barton and collaborators Szymura and Barton (1986); Rouhani and Barton (1987), the primary effect of antagonistic interactions is to introduce a dynamical line tension between interfaces of antagonistic species. This line tension modifies the genetic sectoring phenomenon, creating more compact sectors with more well-defined boundaries. In contrast, genetic sectors are much more diffuse when the strains are neutral, as at the voter model point $\alpha=\beta=0$ in Fig. 1. The transition between these two different coarsening modalities may be easily observed by tracking either the density of interfaces between species [Fig. 6(b)], or by looking at the large-$k$ behavior of the structure factor associated with the identities of the species.

Antagonism modifies the tail of the structure factor from $S(k)\sim 1/k^{2}$ to $S(k)\sim 1/k^{3}$ [Fig. 5(a) compared to (b)]. The structure factor is experimentally accessible, as demonstrated recently in microbial colony experiments McNally et al. (2017). However, the $1/k^{3}$ behavior may be challenging to observe, as any sources of noise will obscure this tail and the structure factor will exhibit an equilibrium $1/k^{2}$ behavior, instead Brown et al. (1998). Thus far, only a $1/k^{2}$ tail has been observed in microbial colony experiments with antagonism McNally et al. (2017). In this previous work of McNally et al., the coarsening of genetic domains is compared to the Ising model below the critical temperature. We have shown that such a connection is incomplete due to the nonlinear noise that acts only at domain interfaces in the population genetics problem. The noise becomes especially important as the antagonistic strength vanishes ($\sigma\rightarrow 0$) near the voter model point $\delta=\sigma=0$. Moreover, varying $\sigma$ significantly alters the time-dependence of the genetic domain interface density. Hence, it would be interesting to study, for example, the structure factor tail and the interface density in an experiment which systematically varies the antagonistic interaction strength $\sigma$. We note that, independent of the degree of antagonism, the typical sector size of the strains should grow as $\sqrt{t}$, as is indeed observed in the microbial experiments.

We have also examined the survival probability of single mutant (yellow) cells in the presence of antagonism. Antagonism strongly suppresses the survival probability with an exponential factor which can be interpreted as an Arrhenius factor $e^{-\beta E^{*}}$ associated with the formation of a “critical cluster” of the yellow strain, where $E^{*}$ is a kind of free energy of the critical cluster and $\beta$ is an effective genetic “temperature” given by Eq. (8). An important aspect of this suppression is that the factor $\beta E^{*}$ is proportional to the spatial diffusion coefficient $D_{s}$, characterizing the motility of the strain. This means that, unlike the survival probability of a mutant with a simple selective advantage and no antagonism [Eq. (4)] which is independent of $D_{s}$, antagonistic interactions introduce a very strong dependence on motility. Thus, we expect that introducing a higher motility to the mutant strain will significantly suppress its survival probability. Consequently, it would be interesting to explore whether strains with antagonistic interactions evolve ways such that genetically similar offspring “stick together” and thus be better able to survive.

Finally, we considered the survival probability of the mutants on the surface of inflating spherical clusters of cells, as one might find in avascular tumor growth. We showed that both the curvature of the spherical cluster and the inflationary effect due to the growth enhance the survival probability of the mutant, as seen in Fig. 9. We also developed a theory of this enhancement, which focussed on effects due to the Gaussian curvature of the sphere for small times and inflationary effects for large times. In the future, it would be interesting to build a dynamical theory which would treat the full time evolution of the mutant organisms. One possibility would be to use the field-theoretic methods employed in Ref. Lavrentovich and Nelson (2015) for strains with no antagonism. We should emphasize that our statistical analysis of number fluctuations for these survival probabilities, in both flat and curved environments, required replacing nonlinear noise correlations such as Eq. (3) by their values at the interface by setting $f=1/2$ [or $\phi=0$ in Eq. (6)]. Although this approximation produced reasonable agreement with our numerical simulations and is mathematically well-motivated in the thin-wall limit, future work should investigate its validity. In addition, some spatial populations have rough, undulating frontiers which can modify these evolutionary dynamics Hallatschek et al. (2007); Lavrentovich and Nelson (2014). We expect frontier roughness to further enhance the survival probability of mutants, because a mutant at the frontier may be able to “grow around” any nearby antagonistic strains by creating a protrusion into the surrounding empty space.

Acknowledgements.

M.O.L. thanks B. T. Weinstein and E. M. Horsley for helpful discussions. D.R.N. benefited from discussions with R. Benzi. Computational support was provided by the University of Tennessee and Oak Ridge National Laboratory’s Joint Institute for Computational Sciences. M.O.L. gratefully acknowledges partial funding from the Neutron Sciences Directorate (Oak Ridge National Laboratory), sponsored by the U.S. Department of Energy, Office of Basic Energy Sciences. Work by D. R. N. was supported by the National Science Foundation, primarily through Grant No. DMR-1608501, and through the Harvard Materials Science and Engineering Center, through Grant No. DMR-1420570.

Appendix A Calculating the Survival Probability with Inflation

In this Appendix we give additional details related to calculating the long-time survival probability $P_{s}^{\infty}(\sigma,\delta)$ of a single mutant yellow cell with selective advantage $\delta$ evolving in a sea of blue cells at the surface of a growing spherical cluster of cells. A convenient coordinate system for describing these dynamics is one in which each point $\mathbf{y}$ on the population frontier at time $t$ is traced back to a reference position $\mathbf{x}=R_{0}\mathbf{y}/R(t)$ on the initial frontier. In this coordinate system, the dynamics evolves according to Eq. (2) but with the altered spatial and genetic diffusion constants $D_{s,g}\rightarrow D_{s,g}/(1+t/t^{*})^{2}$. A detailed analysis for $\sigma=0$, the case of zero antagonism, gives us an inflationary analog of the Kimura formula (see Lavrentovich and Nelson (2015) for details) for the survival probability of a single yellow cell on an initially all blue cell spherical frontier:

[TABLE]

with $G(x)=\int_{x}^{\infty}z^{-2}e^{-z}\mathrm{d}z$ an incomplete gamma function. Note that Eq. (21) reduces to Eq. (4) in the $t^{*}\delta\rightarrow\infty$ limit of a planar frontier. It can be checked that the inflating survival probability at long times is always larger than the non-inflating one: $P^{s}_{\infty}(0,\delta,t^{*})>P^{s}_{\infty}(0,\delta,\infty)$ for any $t^{*}<\infty$. As expected, the inflating environment on the surface of the sphere always enhances survival. Moreover, note that the survival probability no longer vanishes as the selective advantage $\delta\rightarrow 0$. Indeed, inflation can rescue neutral mutations from extinction, as discussed in more detail in Lavrentovich and Nelson (2015).

Note that the formula in Eq. (21) does not include the spatial diffusion coefficient $D_{s}$, much like the original Kimura formula in Eq. (4). The effects of inflation (within the approximations discussed in Ref. Lavrentovich and Nelson (2015)) are entirely captured by the increase in local population size and subsequent decrease in the genetic drift coefficient $D_{g}$. When we include antagonistic interactions, however, we know that the analysis will have to change substantially as the dynamical line tension $\gamma$ introduces an explicit dependence on $D_{s}$. So, to proceed in the case $\sigma>0$, we now consider the nucleation and growth scenario on the surface of an inflating sphere.

We first return to the deterministic dynamics of a cluster of yellow cells. We expect that the overall spherical growth will deterministically inflate the initial droplet area by a factor of $(1+t/t^{*})^{2}$, while inflating its perimeter by a factor of $1+t/t^{*}$. The deterministic dynamics in Eq. (15) then changes and we find that $r_{d}(t)$ has a modified critical radius $r_{d}^{*}(t^{*})$, where $r_{d}(t)\rightarrow\infty$ if $r_{d}(t=0)>r_{d}^{*}(t^{*})$ (and $r_{d}(t)\rightarrow 0$, otherwise), which reads

[TABLE]

where $E_{1/3}(x)=\int_{1}^{\infty}\mathrm{d}z\,e^{-zx}/z^{1/3}$ is an exponential integral function. The expression correctly reduces to the non-inflating result when the curvature effects vanish in the limit $t^{*}\rightarrow\infty$: $r_{d}^{*}(t^{*}\rightarrow\infty)=\sqrt{4D_{s}\sigma/\delta^{2}}$. Moreover, we have $r_{d}^{*}(t^{*})<r_{d}^{*}(t^{*}\rightarrow\infty)$ for any $t^{*}$, so that inflation always makes it easier for the yellow cluster to achieve a critical size (inflation protects against genetic drift).

So far we have only been concerned with the effects of inflation at the frontier. However, the spherical frontier also has a nonzero Gaussian curvature. Consequently, the flat space formulas for the perimeter and area of a circle must be modified. For a fixed perimeter of a nucleating yellow cluster, the area enclosed will be larger on a spherical surface than on a flat surface due to the positive curvature. On a spherical surface of radius $R_{0}$, this geometric effect modifies the energy $E\equiv E[r_{d}]$ [see Eq. (14)] and the critical radius $r_{d}^{*}$ Gómez et al. (2015) as follows:

[TABLE]

where the bars indicate the quantities on the sphere. Here, the droplet radius $\bar{r}_{d}$ is measured along a great circle of the sphere surface from the droplet center. Note that the critical radius $\bar{r}^{*}_{d}$ on the sphere is always smaller than the flat surface value $r_{d}^{*}$. Hence, both the inflationary expansion and the curved droplet at the frontier will enhance the yellow sector survival probability. We see from the simulation results in Fig. 10(a) that the survival probability is in fact significantly enhanced for smaller $R_{0}$ for all values of $\sigma$ and $\delta$.

If the fate of the yellow sector is decided (by a successful nucleation event or extinction) early compared to the inflationary time scale $t^{*}$, then we may use the non-inflating droplet energy $\bar{E}$ on a spherical surface given in Eq. (23) to calculate the survival probability, since we expect $R(t)\approx R_{0}$ at early times. This approximation is reasonable provided the genetic drift is sufficiently strong, which we expect to be the case for a thin layer of growing cells at the population frontier, as in our simulations. The idea here is that surviving yellow cell clusters will reach size $\bar{r}_{d}^{*}$ early when $R(t)\approx R_{0}$ and then deterministically spread through the population as both inflation and the selective advantage $\delta$ take over. Thus, evaluating $\bar{E}$ at $\bar{r}_{d}=\bar{r}_{d}^{*}$ and matching to the (inflating) Kimura formula in Eq. (21) at $\sigma=0$, we find that the inflationary survival probability with antagonism is approximately

[TABLE]

where we found it convenient to define two new functions $p_{1}(x)\equiv e^{-x}[x^{2}G(x)]^{-1}$ and $p_{2}(x)\equiv 2\pi xe^{2\pi x/9}E_{1/3}(2\pi x/9)/9$ in terms of the incomplete gamma function $G(x)$ and the exponential integral function $E_{1/3}(x)$. Note that the effects of the curved front are negligible when $\sigma D_{s}\ll(R_{0}\delta)^{2}$. Our simulations are in the regime where $\sigma D_{s}$ is comparable to $(R_{0}\delta)^{2}$. An approximation to Eq. (24) for $\delta t^{*}\ll 1$ is given by Eq. (20) in the main text.

Unlike the flat, two-dimensional population simulations, our analysis here is limited because we only had the resources to simulate our spherical clusters out to a maximum radius of 200 cell diameters. Although this is a rather large cluster (approximately 40 million cells in the Bennett model of closely packed hard spheres Bennett (1972)), the survival probability may still be decaying and we expect in general simulated probabilities larger than the prediction in Eq. (24). Nevertheless, we can get a sense of the behavior of the survival probability by proceeding as in the flat case. We first use the $\sigma=0$ data to fit the parameter $\tau_{g}D_{g}/A_{0}\approx 0.983\pm 0.008$. Since our Bennett packing model simulation is implemented such that $v=a/\tau_{g}=1$, we then just need the fitting parameters $D_{g}=33\pm 4$ and $D_{s}=0.43\pm 0.06$ (in units where $a=\tau_{g}=1$) to fully specify $P_{\infty}^{s}$ in Eq. (24). The results are shown in Fig. 10(b). Our fitted genetic drift magnitude $D_{g}$ is large compared to the flat population, probably due to the incompatibility of comparing $P_{\infty}^{s}$ to a simulated survival probability at a finite time. In addition, we do generally expect different values for our parameters compared to the flat populations because the Bennett cluster algorithm generates a disordered sphere packing. The inset to Fig. 10 shows the relative error between the theoretical $P_{\infty}^{s}$ and the simulated $P^{s}$ for the $R_{0}=30a$ case. The theoretical $P_{\infty}^{s}$ is generally less than the simulated $P^{s}$, as we might expect from the finite time of the simulation. Note that the theory does not match the simulations as well for small $\delta$ and for small $\sigma$, as can also be seen in the inset in Fig. 10. This anomaly arises because we are outside the region of metastability for small $\sigma$. Also, at small values of $\delta$, the critical yellow cluster size becomes large and the fate of the yellow cell cluster will no longer by decided at early times $t\ll t^{*}$, so the approximations described above break down.