Survival Probabilities at Spherical Frontiers

Survival Probabilities at Spherical Frontiers

Maxim O. Lavrentovich Department of Physics and Astronomy, University of Pennsylvania, Philadelphia PA 19104 USA David R. Nelson Lyman Laboratory of Physics, Harvard University, Cambridge MA 02138 USA

Motivated by tumor growth and spatial population genetics, we study the interplay between evolutionary and spatial dynamics at the surfaces of three-dimensional, spherical range expansions. We consider range expansion radii that grow with an arbitrary power-law in time: , where is a growth exponent, is the initial radius, and is a characteristic time for the growth, to be affected by the inflating geometry. We vary the parameters and to capture a variety of possible growth regimes. Guided by recent results for two-dimensional inflating range expansions, we identify key dimensionless parameters that describe the survival probability of a mutant cell with a small selective advantage arising at the population frontier. Using analytical techniques, we calculate this probability for arbitrary . We compare our results to simulations of linearly inflating expansions ( spherical Fisher-Kolmogorov-Petrovsky-Piscunov waves) and treadmilling populations (, with cells in the interior removed by apoptosis or a similar process). We find that mutations at linearly inflating fronts have survival probabilities enhanced by factors of 100 or more relative to mutations at treadmilling population frontiers. We also discuss the special properties of “marginally inflating” expansions.

survival probability; genetic drift; range expansions; avascular tumor evolution; selection
journal: Journal of Theoretical Population Biology\biboptions


1 Introduction

Early tumor evolution is driven by rare driver mutations that sweep the prevascular tumor population at the frontier and push the growing cell mass further down the path toward metastasis. Hence, an understanding of the evolutionary dynamics governing the survival of such mutations is crucial in cancer prevention (Merlo et al., 2006; Vogelstein et al., 2013). One significant, largely unexplored aspect of this evolution is the effect of tumor geometry. An important in vitro model of cancer is the multicellular tumor with an approximately spherical shape, or “spheroid”. The spheroid captures many of the essential features of solid tumors in vivo and is a model for anti-cancer therapies (Kunz-Schughart, 1999; Santini and Rainaldi, 1999; Hirschhaeuser et al., 2010). Spheroids are especially useful for understanding small, avascular tumors. In the later stages of growth, in order for the tumor to survive, it requires a vascular system and undergoes angiogenesis (Weis and Cheresh, 2011). The growth then becomes more complicated, and more sophisticated modeling efforts are necessary (Shirinifard et al., 2009; Alarcón et al., 2005). We focus here on the earlier evolutionary dynamics of spheroidal range expansions in two and three dimensions. We assume that attractive cell-cell interactions keep such aggregates approximately spherical, i.e. that there is an effective surface tension, similar to that observed for yeast cell colonies (Nguyen et al., 2004). Although we are motivated by tumor evolution, our models are intended to be quite general. Two-dimensional and three-dimensional expansions may be realized in experiments, for example, using microbial or yeast populations in hard and soft agar, respectively (Korolev et al., 2012, 2011; Lavrentovich et al., 2013a).

We will be particularly interested in computing the survival probability of a mutation that occurs among the dividing cells at the surface of a spherical or circular population of initial radius , which may or may not increase in time. In general, the radius has a complicated time dependence, especially in tumor growth. At the early stages, cells divide everywhere inside the tumor, and the cluster radius grows exponentially in time. After the tumor reaches a size larger than a nutrient shielding length (Lavrentovich et al., 2013b), nutrients will no longer be able to diffuse into the tumor interior. This effect, combined with inward pressure from the surrounding non-cancerous tissue (Cheng et al., 2009; Montel et al., 2011, 2012), decreases the growth rate toward the center of the tumor. The radius then grows more slowly. We will model the growth generally as


where is the initial tumor radius, is a (possibly time-dependent) growth exponent, and is a characteristic time for the power-law growth in an inflating geometry. Both and may be tuned to model various growth regimes. For example, for a substantial portion of the growth in tumors, the radius grows linearly in time , and , where is the front speed (Brú et al., 2003). Linear growth and nutrient shielding are also present in microbial populations grown in Petri dishes (Korolev et al., 2012). Equation (1) can also model an exponential growth regime with rate if we let both , such that is held constant.

Figure 1: (Color online) A schematic of a treadmilling tumor. Due to nutrient shielding, cells divide in a thin green region at the frontier. In the red region, cells are in an arrested state and do not grow. In the necrotic core, cells undergo apoptosis and their contents are flushed out of the cluster, resulting in an overall volume loss. This volume loss can balance the gain of volume at the cluster periphery, resulting in a “treadmilling” effect and a cell mass with a constant radius (Cheng et al., 2009; Montel et al., 2011, 2012; Stott et al., 1999).

Eventually, apoptosis may be induced at the tumor center, creating a necrotic region (Cheng et al., 2009; Montel et al., 2011, 2012), illustrated in Fig. 1. The cells at the tumor periphery continue to divide relatively rapidly. Thus, a “treadmilling” effect is created, and the tumor experiences a rapid turnover of cells at its surface while remaining the same size, a situation we represent by a growth exponent in Eq. (1). We will show that the different growth regimes captured by varying have dramatically different consequences for the fate of mutations at the tumor frontier. We will focus on , , and , capturing, respectively, treadmilling, linear inflation, and an intriguing borderline growth regime.

The actively growing region in a tumor mass or a spherical microbial population can be quite thin, with a width of just a few cell diameters (Folkman and Hochberg, 1973; Lavrentovich et al., 2013b). In this case, genetic drift is strong and can locally fix the mutation at the population frontier. This local fixation creates a mutant “sector,” i.e., a region along the front that is entirely occupied by the mutant cells. Example sectors, marked in green, are shown in Fig. 2 for circular and spherical range expansions. Previous studies have focused on the deterministic movement of these mutant sectors: The sectors inflate or deflate due to a mutant selective advantage or disadvantage, respectively (Antal et al., 2013). However, genetic drift will introduce fluctuations in the sector motion at its boundaries that can drive the mutation to extinction, as illustrated on the left panels of Fig. 2(a) and Fig. 2(b). For example, in circular expansions [Fig. 2(a)], the sector has two boundaries which both perform random walks, as observed in microbial range expansions (Korolev et al., 2010, 2012). Selection introduces a bias to the sector boundary motion. Also, the increasing population radius will deterministically increase the distance between the boundaries. If the sector boundaries collide, the sector vanishes, and the mutation goes extinct [left panel of Fig. 2(a)].

Figure 2: Examples of simulated mutant clusters (green cells) in two- and three-dimensional range expansions (see Section 2), generated using two different values of the key dimensionless selection parameters defined in Eq. (2). (a) Circular range expansions with a uniform front with an initial radius average cell diameters and a single initial mutant green cell at the population frontier. (b) Spherical range expansions with uniform fronts and a single green cell at the initial population frontier with radius cell diameters.

Because the two boundaries of mutant sectors in two-dimensional expansions perform random walks, the distance between the boundaries performs a random walk as well. The mutant survival probability, then, is the probability that this distance never vanishes, i.e., that the random walk it performs never reaches the origin. Such a probability is a well-studied first-passage property of a random walk (Redner, 2001). Previous studies, such as Hallatschek and Nelson (2010); Korolev et al. (2010); Lavrentovich et al. (2013a), have exploited these known random walk results to calculate survival probabilities of mutations in two-dimensional populations. We will review some of these previous results for two-dimensional expansions and generalize them to arbitrary growth exponents . We will then use them to motivate our discussion of three-dimensional expansions.

In three-dimensional populations, the focus of the present work, the mutant sector can have a complex, branched shape, as shown in Fig. 2(b). Characterizing the boundary positions with just two random walks is impossible in this case. Although we can still treat mutant cell lineages using random walks [see, e.g., Cox and Griffeath (1986) or the chapter on voter models in Liggett (1985)], incorporating selection is more complicated than in the two-dimensional case (Bramson and Griffeath, 1981). Also, we know of no generalization to inflating frontiers. So, instead of mapping to random walks, we study the time evolution of the coarse-grained density of mutant cells along the population frontier of spherical range expansions. We then apply a field-theoretic analysis of the time evolution that allows us to treat genetic drift, selection, and inflation within a single theoretical framework.

By using random walk theory for two-dimensional range expansions and field-theoretic techniques for three-dimensional ones, we will show that the key dimensionless parameters for mutant survival are, respectively,


where is the characteristic time of the radius growth defined in Eq. (1), is the selective advantage, and is a generation time. The mutant survival will also depend on the initial number of mutant cells and the growth exponent . The mutant cluster shapes at different for linearly inflating frontiers [ in Eq. (1)] are illustrated in Figs. 2(a) and 2(b) for two- and three-dimensional expansions, respectively. For , the selective advantage is large enough to sustain a surviving mutant cluster after a short time, allowing it to survive indefinitely and (if no other mutant sectors are present) to sweep the population deterministically. When , however, the effects of genetic drift are important for a longer time, and can lead to an extinction of the cluster.

The remainder of the paper is organized as follows: We introduce the simulation model used to create two- and three-dimensional expansions in Section 2. We impose a global constraint that insures compact clusters, as a computationally efficient way of emulating the effect of surface tension. In Section 3 we review the deterministic dynamics of mutant sectors, exemplified by the average sector shapes corresponding to the right panels of Figs. 2(a) and 2(b). We also introduce the stochastic equation governing both the deterministic and stochastic dynamics. In Section 4 we review and extend results for two-dimensional expansions, as these provide insights into the three-dimensional case relevant to tumors. In Section 5 we calculate the survival probability of mutations in three dimensions, relegating details to A and B. Results for the case of with linear inflation are presented in C. The smooth circular and spherical fronts generated in our simulations facilitate comparison with analytical results. We make concluding remarks in Section 6.

2 Simulations

Our simulations will focus on populations with a single layer of growing cells. We also assume, to simplify the analysis, compact fronts where the population front closely approximates a uniform circle or sphere at all times. The latter approximation will be valid as long as front undulations and the selective advantage are small. (For discussions of how the dynamics change when these approximations are not valid, see Antal et al. (2013); Korolev et al. (2012); Kuhr et al. (2011)). As discussed in more detail by Lavrentovich et al. (2013a), these approximations have the advantage of allowing us to perform a “dimensional-reduction” by focussing on just the dynamics at the population frontier, which will be effectively one- and two-dimensional inflating geometries for two- and three-dimensional expansions, respectively.

All along these frontiers, cells will compete locally to divide into new territory and form the next generation of cells. We consider two types of cells: a wild-type red cell and a green cell with a selectively advantageous “driver” mutation. The mutant, “driver” green cells have a base growth rate that is greater than the red wild-type cell growth rate by a factor of , where is a selective advantage. The difference in the growth rates determines the probability of a green offspring cell occupying an available spot along the frontier. Specifically, if there are a total of “parent” cells, of which are green, competing to divide into an available spot, the probability of placing a green offspring in the spot is


New cells are placed one at a time in predefined locations to create expansions with uniform fronts, as described below. The probability of a particular cell color dividing into the new spot is determined by Eq. (3), with the probability of a red offspring given by .

Figure 3: Schematics of simulated radial range expansions in (a) two and (b) three dimensions. In (a), cells of two different sizes compete to divide into new territory. The possible cell locations are chosen in advance using an algorithm described in the main text and in Lavrentovich et al. (2013a). The algorithm is seeded with the two cell locations indicated by black crosses. The two cell sizes and their relative abundances are chosen to create an amorphous, isotropic cell packing. In (b), cells with identical diameters compete to divide. The cell locations are generated using the Bennett model (Bennett, 1972). In both (a) and (b), the virgin territory that is closest to the cell cluster center is settled first to create an effect similar to a surface tension. In both models, the green cells out-compete the red with probability given by Eq. (3).

We fill empty sites one at a time, starting with the site that is closest to a central reference point, as shown in Fig. 3. In one time step, the cells adjacent to the empty site compete to divide into the site. Repeating this process creates a growing cell cluster with a relatively smooth circular or spherical front, approximating the effect of a surface tension. Surface tension will suppress both front undulations and prominent bulges arising due to the faster growth of the mutant cells. Our methodology is reasonable provided surface tension forces only move daughter cells a short distance after they are born (i.e. about a cell diameter). The resulting compact cell cluster could also arise due to biological reasons. For example, cells may secrete a chemical that promotes cell divisions, which then are only feasible at the cluster surface where there is room to divide (Lavrentovich et al., 2013a). This smooth frontier approximation is not only more tractable analytically, but increases our simulations’ computational efficiency.

The set of empty sites into which cells may divide (i.e. the set of possible cell positions) is generated in advance. One possibility is to arrange the cell positions in a lattice. However, for outwardly growing circular range expansions, this approach produces strong artifacts of the lattice symmetry in the evolutionary dynamics (Lavrentovich et al., 2013a). We mitigate these effects by generating an amorphous packing of empty sites using a modified Bennett model (Rubinstein and Nelson, 1982). In this model, we first seed the packing with two adjacent sites, indicated by black crosses in Fig. 3(a). The center of this initial seed is our central reference point. Then, cell locations (sites) are generated one at a time and are placed as close as possible to the central reference point (without overlapping previously generated sites). Each site must also touch at least two previously placed sites. A smaller site is chosen over the larger with some specified probability. Lavrentovich et al. (2013a) found that isotropic, amorphous cell packings occurred if cells are chosen with a small-to-large radius ratio of , with the smaller one placed with 60% probability, choices we also make here. Finally, note that evolutionary dynamics do not depend on the cell size. A small cell may divide into an empty site that fits a large cell, or vice versa.

For spherical range expansions, we generate our set of possible cell locations using the original Bennett model (Bennett, 1972) with identical spherical cells with diameter . These cells are now placed closest to the central reference point such that each touches at least three previously placed cells. The initial seed is a tetrahedron of cells. The algorithm creates an amorphous packing without the need to vary the cell size as in the two-dimensional case. Circular and spherical expansions with relatively small population sizes are shown in Fig. 3. Note that a cell diameter is more complicated to define in the two-dimensional expansion because we use both small and large cells. So, in this case will denote an effective spacing between cells that is approximately equal to the average cell diameter.

To further mitigate lattice artifacts, we initialize the mutant green cells during each simulation run in one of at least 100 different, random locations, sampled randomly and uniformly along the population periphery. In two dimensions, all cells along the initial, circular population front have equal probability of being chosen as the mutant green cell. Thus, the mutant may be either a small cell or a large cell, depending on the location chosen. In three dimensions, we first perform a random rotation of the initial spherical population as described by Arvo (1992). We then seed the mutated, green cells at the north pole of the spherical population. This ensures a uniform, random sampling of the initial population front. We study initial fronts with just one green cell and fronts with a localized, compact clumps of many green cells.

This algorithm generates linearly inflating fronts with . It is possible to simulate treadmilling tumors with by adding a step that evolves the population backwards, toward the interior of the cluster: First, we evolve the range expansion forward just a few generations, generating a cluster of radius , where . Second, we go back to the cells at approximately a distance from the cluster center, and replace them with new offspring cells. This is done by evolving the population backwards, toward the central reference point. The cells at the population periphery (at distances from the center) divide toward the population interior and replace the cells generated during the forward sweep. The previously generated cells are replaced one at a time as the population evolves deeper into the interior. This replacement models both the volume loss due to apoptosis at the population center and pressure from the surrounding tissue, which prevents the tumor from growing outward. This backwards evolution is continued until a population front with radius is created. Finally, we go back to the cells at radius and repeat the first step, re-evolving the initial population radius outward up to a radius . We again let cells divide into sites previously occupied by cells generated during the backward sweep. The backwards and forward sweeps are repeated over and over to model a turn-over of cells at the tumor surface, taken to be the spherical shell of cells with fixed radius and thickness .

In both the treadmilling and expanding populations, our cells may divide as long as they are adjacent to either an empty site or, in the treadmilling case, a previously generated cell chosen to be replaced. So, the lifetime of a cell is entirely governed by the availability of adjacent empty sites. Once these have been filled, the cell may no longer divide.

3 Population genetics at compact population fronts

Before we analyze survival probabilities of mutations, it is instructive to first consider the deterministic dynamics of a mutation. In particular, we will now describe how a mutation with a selective advantage sweeps across the population frontier in a radial range expansion [corresponding to in Eq. (1)], ignoring sector extinction due to genetic drift (Antal et al., 2013; Korolev et al., 2012). A surviving mutant forms a sector as shown in Fig. 2 for two dimensions and schematically in Fig. 4 for three dimensions. We assume and/or a strong surface tension, so that we may ignore any bulges created by the green strain and to make contact with our simulations. The boundary of this sector is a genetic Fisher wave in which the mutant green strain “invades” the red wild-type region azimuthally along the frontier with some speed . This means that in a two-dimensional circular expansion, the two sector boundaries will move away from each other with a relative speed along the population frontier. The range expansion itself will grow outward with a velocity . Hence, the sector angle between the two boundaries (relative to the center of the range expansion) will increase as , where . Integrating this differential equation, we find that the angular sector size is given by


where is the initial sector size and is a rescaled sector boundary speed (Fisher wave speed) that depends on the selection coefficient and is influenced by the strength of the genetic drift. Equation (4) describes a sector with boundaries that form logarithmic spirals that spiral away from each other (Korolev et al., 2012).

In general, as in Eq. 4, independently of the strength of the genetic drift. For our simulations, , where is the distance the front moves in one generation and is the effective spacing between cells discussed in Sec. 2. The scaling is a consequence of strong genetic drift, and, as explored in more detail by Lavrentovich et al. (2013a), can be derived by mapping the sector boundaries to random walks (see Sec. 4). The mapping of sector boundaries to random walks allows us to generalize our analysis to sectors with bulges, as well. Indeed, we expect the green sector to bulge out because it will grow out radially faster than the red strain. However, provided the front has some effective line tension and does not roughen, the sector boundaries will still perform biased random walks with some effective rescaled bias speed in the presence of bulges. Hence, our analysis of the sector survival probability in Section 4 will apply for this case, as well. We will now turn to sector shapes in three dimensions. In three dimensions, mutations with a strong selective advantage form “logarithmic cone” sectors, i.e. solids of revolution formed by logarithmic spirals (Antal et al., 2013).

Figure 4: Schematic of the deterministic motion of a mutant sector arising at the surface of a spherical population of red cells. The green sector invades the red population with a selection-dependent lateral velocity , where is the radial front speed. The sector forms a logarithmic cone, formed by rotating a logarithmic spiral curve parameterized by the angle around the -axis. In our simulations, we ignore possible bulges in the population front and neglect the enhanced velocity of the green sector relative to the red domain expansion velocity . Hence, the green sector forms a spherical cap with some solid angle at the spherical population frontier.

We can compute the analogous angle for a spherical expansion by calculating the solid angle covered by the mutant sector at time , as shown in Fig. 4 (Antal et al., 2013). Suppose that the initial solid angle covered by the mutant sector is . The corresponding initial angle is . Then, the angle of the shape of the edge of the logarithmic cone is given by the relative angle :


where . We check this result with simulations in Fig. 5. We find that Eq. (5) accurately describes the average mutant sector shape, as shown by the linear fits through the data in the main plot. Note that when , the logarithmic cone sector boundary collapses, and the sector boundary approaches a constant solid angle. When many mutations with are present, these logarithmic cones of mutations can collide to form interesting limiting shapes discussed in detail by Antal et al. (2013) (see also Martens et al. (2011)). Upon using Eq. (5) to extract the rescaled lateral Fisher wave velocity , we find to an excellent approximation, as illustrated in the inset of Fig. 5. Note that this square-root scaling is markedly different from the circular expansion case. Three-dimensional inflating range expansions are thus consistent with the classic result that the Fisher wave-speed should approach when there is no genetic drift, where is a spatial diffusion coefficient (see Fisher (1937) and discussion below). We expect that the coefficient of the term depends on the noise (genetic drift) strength, which is fixed in our simulations. Note that for small is much bigger than the difference between the mutant and wild-type radial expansion velocities (of order ), thus justifying our neglect of the bulge in three dimensions. This scaling result is consistent with previous studies of noisy Fisher waves in higher dimensions (Moro, 2001; Hallatschek and Korolev, 2009).

Figure 5: Simulated mutant clusters in a spherical range expansion ( cell widths) with an initial solid angle steradians. We plot the angle delineating the edge of the conical shape formed by the mutant cluster edge, averaged over 38400 surviving sectors. After an initial transient, mutant sectors form logarithmic cones for any selective advantage , as illustrated by the linear fits through the data (black lines). In the inset, we plot the slope of these lines, . (Note that the case at the bottom is exceptional.) We find, to an excellent approximation, .

We now have to introduce genetic drift explicitly to take into account possible sector extinction. This will be done in Sec. 4 for two-dimensional expansions by mapping sector boundaries to random walks. However, to better understand our techniques for three-dimensional expansions in Sec. 5, we now develop some additional results for noisy Fisher waves in non-inflating range expansions. We define a coarse-grained green cell (mutant) fraction at positions along the front at time . We compute this fraction over box-shaped regions of volume , where is the box length or effective lattice spacing, and is the spatial dimension of the front, i.e., for a two-dimensional expansion and for a three-dimensional expansion. To compare the analytic results with our simulations where the fronts are monolayers of closely-packed cells, is chosen such that describes a region of the front occupied by one cell. For example, for a one-dimensional frontier in which the cells are collinear and touching each other, is a cell diameter. We assume that the characteristic size of a mutant cluster is small compared to the population front radius and that the front is uniform and flat. If the front is also not inflating (), the coarse-grained green cell fraction (upon exploiting dimensional reduction) obeys the Langevin equation


where is the selective advantage of the mutant, is a spatial diffusivity, is the generation time, is a genetic drift strength, and is a Gaussian, white noise [] interpreted in the Itô sense (see Section 4.3 in Gardiner (1985) or Korolev et al. (2010)). The strength of the genetic drift scales like , where the effective population size is approximately the number of organisms in each region that compete to divide at the population frontier. Equation (6) may also be derived using a stepping stone model (Korolev et al., 2010).

If we remove the genetic drift term in Eq. (6) by setting , we recover the Fisher-Kolmogorov-Petrovsky-Piscunov equation (Fisher, 1937; Korolev et al., 2010). This equation allows traveling-wave solutions that describe the deterministic sweep of a mutant sector (Korolev et al., 2012). These Fisher waves move at speed . The genetic drift induced by a non-zero can strongly influence these traveling-wave dynamics (Hallatschek and Korolev, 2009; Doering et al., 2003) and can extinguish an embryonic mutant sector. There are mathematical difficulties associated with interpreting Eq. (6) for when the noise term is present and when a continuous variable is used (Barton et al., 2002). However, we will always impose finite large and small distance cutoffs ( and , respectively) where necessary to get meaningful results, and we will test and confirm our analytic results with careful numerical simulations. So, although the sense in which Eq. (6) is mathematically sensible and may be used to compute averages and probabilities for arbitrary has been explored in the physics literature (Duty, 2000; Janssen, 2005), understanding this fascinating, subtle problem is not relevant for the calculations and regimes explored in our paper. Also, to understand the evolutionary dynamics of inflating circular and spherical range expansions in the presence of genetic drift () using Eq. (6), the equation has to be appropriately modified to take into account the inflating population fronts. We will make the appropriate modifications in Sec. 5, where we will also analyze the equation with field-theoretic techniques.

We will focus on an important, biologically relevant consequence of Eq. (6): What is the probability that a single green mutant cell establishes a surviving cluster at long times, given an initial green cell fraction at the population frontier? In the strong genetic drift limit, Doering et al. (2003) mapped Eq. (6) in one dimension to a reaction-diffusion model and found that the survival probability for the mutants is


where we integrate over all positions along the population frontier. Equation (7) is valid for non-inflating, large population frontiers. We expect Eq. (7) to hold for three-dimensional expansions, as well. This hypothesis was tested for a more sophisticated range expansion model developed by Pigolotti et al. (2013). We will show that Eq. (7) can be derived via a field theoretic technique in Section 5. Note that our choice of effective lattice spacing will allow us to relate the integral over the initial fraction to the initial number of mutant green cells :


where we may consider both one-dimensional and two-dimensional population fronts. Substituting Eq. (8) into Eq. (7) for , we see that the survival probability vanishes linearly with , , when . In contrast, we will find in Sections 4 and 5 that for inflating expansions in two and three dimensions, this probability instead approaches a non-zero value as .

The case, although not the focus of this paper, is interesting and subtle. In the deterministic regime, the boundaries of the sector in two dimensions form logarithmic spirals that spiral toward each other for , forming a tongue-like, or flower petal shape. Similarly, in three dimensions, the boundary of the logarithmic cone spirals inward instead of outward. Hence, we expect that the boundaries of any surviving mutant cluster will eventually be pinched off due to this deterministic motion for . Our theory predicts that the only way a mutant sector with can survive at long times is if it wraps all the way around the circumference of the population. Otherwise, the long time survival probability vanishes: . However, the eventual extinction of a sector might take a very long time, due to the logarithmically slow deterministic dynamics. The convergence of the survival probability to its long-time steady-state value is extremely slow for . Thus, we may not always be able to probe this long-time limit in our simulations, which can only evolve the range expansions up to some finite time. We briefly examine some of these issues in C.

4 Treadmilling and Expanding Populations in Two Dimensions

In two-dimensional, circular populations, mutant sector boundaries perform random walks along the population frontier (Korolev et al., 2010, 2012). Each mutant sector will have two boundaries, as shown in Fig. 2. Since each sector boundary performs a random walk, the sector angle between them, , performs a random walk, as well. We can write down a general stochastic differential equation for that describes its random walk on a circular frontier with a time-dependent radius given by Eq. (1):


where is a selection-dependent bias and is a diffusion constant. The function in Eq. (9) is a Gaussian white noise such that . We assume that a single sector boundary will only be able to move about a cell diameter per generation time due to cell rearrangement during division at the frontier. Hence, , where is the initial population radius. In addition, we assume that a green cell out-competing a red cell at the boundary can only shift the boundary by about a cell diameter. Then, since a green cell out-competes a red cell at a rate proportional to the difference between their respective growth rates, we have . A more detailed derivation of and is presented in Appendix C of Lavrentovich et al. (2013a). Note that this random walk model is consistent with the dynamics of mutant sector sweeps described in Sec. 3: if we set (so that ) and look at vanishing genetic drift , Eq. (9) yields the deterministic sector shape given by Eq. (4) in the previous section with .

4.1 Treadmilling and Marginally Inflating Fronts

When , Eq. (9) describes a drift-diffusion motion for the  variable with a constant drift and diffusion constant. At long times, there are only two possibilities for the fate of the sector: It either eventually vanishes (), or it wraps all the way around the population frontier (). The probability of the latter occurring is the survival probability. Using standard first-passage techniques [see, e.g., Eq. (2.3.8) in Redner (2001)], we find that this survival probability for our model in Sec. 2 is


where is the initial number of mutant cells at the frontier. We also set to describe expansions with compact frontiers that are a single cell wide. Let us now consider a single initial mutated cell, so . In this case, Eq. (10) resembles the classic Kimura formula for the fixation probability of a mutation with a selective advantage in a well-mixed population with initial frequency and total population size (see, e.g., Eq. in Crow and Kimura (1970)). Indeed, as shown by Maruyama (1970, 1974), the survival probability of a mutation with a selective advantage is the same in well-mixed populations and at non-inflating population fronts. When , we find a non-zero survival probability. Upon expanding Eq. (10) for small (with ), we have


Note that when in Eq. (10), the survival probability is exponentially suppressed,


for .

For , Eq. (9) has a time-dependent drift and diffusion term. This makes the survival probability analysis more challenging. However, neutral mutations with may be treated exactly, as discussed in Ali et al. (2013a, b); Ali and Grosskinsky (2010); Lavrentovich et al. (2013a). In the limit of a large population front in which the mutant sector stays small (), a neutral mutation will always eventually go extinct. The case is marginal, and a neutral mutation is nearly able to survive. We now move to a discussion of the case, for which inflation produces a striking enhancement of the survival probability. We will go back to the marginal case when we treat three-dimensional expansions in Section 5.

4.2 Inflating fronts

When , we may remove the time-dependence from the diffusion term in Eq. (9) by changing variables from the time to a unitless time given by


After this change of variables, Eq. (9) becomes


where we again have a Gaussian white noise such that . Equation (14) describes a simple diffusive motion for the sector boundary in the conformal coordinate , with a time-dependent bias. Note that when , then , also. However, if we send in Eq. (13), then . This means that the entire evolution of the sector, i.e., the interval , is compressed into the interval . The consequences of this compression are explored in more detail in Ali and Grosskinsky (2010); Lavrentovich et al. (2013a); Ali et al. (2013a, b).

The survival probability can be extracted from Eq. (14) by analyzing the first-passage probability that the sector angle vanishes (see Redner (2001) for a review of first-passage processes). Unlike the case, we will have to work in the limit of large population fronts so that we can ignore large . For , for example, we require that and . We will compare our results to simulations and to the exactly soluble case to check the approximation. We find a long-time survival probability that depends on two scaling variables and , given by


where we have again specialized to thin, compact population frontiers by setting and . Note that we can only make the replacement in Eq. (15) when .

The survival probability scaling function can be calculated using an adiabatic approximation. Although the survival probability depends on many parameters, i.e. , we find that the full time-dependent solution for in the adiabatic approximation depends on just two parameter combinations [see Eq. (15)] and on the growth exponent :


We can set (i.e. ) in Eq. (16) to find the ultimate survival probability
. Although the integral does not appear to be expressible in terms of elementary functions, we may numerically integrate it for arbitrary parameter values.

We now compare our analytic result, Eq. (16), to simulations for . Equation (16) is numerically integrated for to obtain the scaling function shown in Fig. 6. To determine and for our simulations (see Section 2), we first set our length units so that . Then, the effective cell spacing is extracted from an independent measurement in Lavrentovich et al. (2013a). In addition, simulations of the deterministic motion of the sector yield , as mentioned in Sec. 3. Finally, setting , we find and for our simulations using Eq. (15). This procedure is described in more detail in Lavrentovich et al. (2013a). The simulated survival probabilities are shown as data points in Fig. 6. We find that the simulation results match the predicted scaling form well. We also show the treadmilling result in Fig. 6 for a fixed to illustrate the large enhancement in the survival probability due to inflation.

Figure 6: The long-time (, ) survival probability scaling function of a green mutant sector [as in Fig. 2(a)] for circular range expansions with , as a function of the two dimensionless variables and displayed in Eq. (15). Clusters of contiguous mutants subtend an initial angle along the population front of an otherwise all-red population with initial radius . Different symbols lying on the same line correspond to data sets with different . In part (a), is tuned to yield a fixed . In part (b), is varied to find the probability at different at a fixed . The solid symbols on the bottom curve in part (b) correspond to neutral range expansions with . The lines show the analytic prediction using the adiabatic approximation in Eq. (16) for (accurate for and exact for ). In simulations (points), is estimated from the survival probability at a long (but finite) time generations. Hence, we expect some deviation of the data from the analytic prediction, as the survival probability in simulations converges slowly to for large and small . Solid orange lines show the survival probability for a treadmilling population front with and a single initial mutant cell. The orange line does go to a non-zero value as in part (a), but it is indistinguishable from zero in the plot. Note the greatly enhanced survival probability caused by inflation in part (a), as compared to treadmilling systems.

As we approach a neutral mutation in Fig. 6(a), we find a non-zero survival probability. This is markedly different from the infinite flat front result in Eq. (7), which vanishes as . Indeed, for a single initial mutant cell at a circular front with initial radius (), the limit is


It is instructive to compare this result to the much smaller survival probability of a mutation occurring at the surface of a two-dimensional treadmilling tumor. The survival probability as is much larger for an inflating rather than a treadmilling tumor, as shown in Fig. 6(a). Thus, for small selective advantages, inflation plays a much more important role in rescuing the mutation than one might expect from a naive analysis of a mutation sweeping a population with a finite size. We will show in the next section that there is an analogous scaling function , with appropriate three-dimensional analogues of the dimensionless scaling variables and . We will again find that the effects of a finite population front size are less important than inflation at small .

5 Treadmilling and Expanding Populations in Three Dimensions

Figure 7: Coordinate system for studying spherical expansions. At time when the expansion has radius , we can map any point on the spherical population front, , to a corresponding point on the initial population front with radius , along a radial line connecting to the center.

A convenient coordinate system for three-dimensional range expansions is the position on the initial spherical population front . Any positions on the spherical front of radius at any time can be traced backward to a unique on the initial front, as shown in Fig. 7. Thus, we can fully specify any position in the range expansion by a single position and a time . We then interpret the evolution of the range expansion as a process acting on the original population front with a superimposed inflation that locally dilates the front perpendicular to the local expansion direction. The corresponding Langevin equation for the green cell fraction with this inflationary effect reads


where is again a Gaussian white noise induced by genetic drift [as in Eq. (6), interpreted via the Itô calculus] and is a characteristic time. Equation (18) is a non-linear stochastic partial differential equation with a spatial dilation. Such equations are of interest not only in this biological context, but also in cosmology, interface growth, and other areas of condensed matter physics (Escudero, 2013). In the case of linear inflation , we have , where is the population front speed. In this case, Eq. (18) is a straightforward extension of the stochastic differential equation that describes circular range expansions (Lavrentovich et al., 2013a; Korolev et al., 2010). We now ignore the finite size of the sphere and extend our coordinate to the infinite two-dimensional space . This approximation will work as long as our sector sizes are small compared to the total surface area of the front. We move to the field-theoretic representation of Eq. (18), suitable for application of the methods of path integrals and statistical mechanics, and look for the ultimate survival probability . For a review of similar ideas from the perspective of population genetics, see de Vladar and Barton (2011). Readers unfamiliar with these formal developments may wish to pass on to our final conclusions in Sections 5.1 and 5.2 which have been subjected to extensive numerical checks.

Equation (18) can be converted to a field-theoretic action by employing standard methods, pioneered over 40 years ago by Martin et al. (1973), and extended to non-linear stochastic equations such as Eq. (6) and Eq. (18) by Janssen (1976) and de Dominicis (1976). We introduce a response function to account for the genetic drift term in Eq. (18). This procedure yields the response functional


where is the area taken up by a single cell at the initial spherical frontier with radius . Equation (19) allows us to treat inflationary dynamics on the sphere using the methods of statistical mechanics. The genetic drift term (proportional to ) is now on equal footing with the other contributions, making it easier to analyze than in the Langevin equation, Eq. (18). When such functionals are exponentiated, they determine the probability of particular path histories associated with the underlying stochastic differential equation (de Vladar and Barton, 2011; Martin et al., 1973; Janssen, 1976; de Dominicis, 1976).

The response functional in Eq. (19) can be used to calculate the survival probability for the initial patch of green cells. As shown in A, a variational principle applied to the exponentiated response functional, Eq. (19), yields a spatially uniform, mean field approximation to the response function . With this mean field approximation, our long-term survival probability for arbitrary growth exponents [see Eqs. (46)-(48) in A] can be obtained by a limiting procedure,


where satisfies


and where is the step function (, , , ).

Equation (21) may be solved for for an arbitrary . This is done in the Appendices. Substituting the solution for into Eq. (20) and taking the appropriate limits yields the survival probability


where is the incomplete gamma function (Gradshteyn and Ryzhik, 2007). We may now use Eq. (22) to identify two key dimensionless parameters and :


where is the initial number of mutant cells. Note the similarity between Eq. (23) and the analogous dimensionless variables for circular inflation in Eq. (15). The survival probability only depends on these two parameters and the growth exponent :


Equation (24) is the main result of our field theoretic approach. We will check that it approximates the survival probability well by comparing to simulations of treadmilling () and linearly inflating () population fronts. First, however, we probe the properties of Eq. (24) by looking at the survival probability for a neutral mutation (the limit). A careful analysis of this limit yields


Equation (25) tells us that is an interesting “marginal” case for which inflation is nearly able to rescue a neutral mutation. It is also possible to calculate the survival probability of a neutral mutation on the frontier of an exponentially growing population . If we substitute in Eq. (25) and take with a constant ratio , we find that . We will consider the three cases , , and in more detail in the next subsections.

5.1 Treadmilling fronts

When , the range expansion simply turns over cells at its surface, which we take for simplicity to be a one-cell-diameter thick spherical shell with radius . These dynamics are relevant for “treadmilling” tumors. If we ignore the finite size of the shell (a good approximation provided the mutant cluster is small, i.e., ), the survival probability can be calculated from the solution to Eq. (24) with . We find a survival probability


for a very large population front . This result is consistent with the results of Janssen (2005) and with the generalized formula of Doering et al. (2003) in Eq. (7). We expect corrections for both positive and negative when the finite size of the population front is taken into account, as in the treadmilling circle result in Eq. (10). We note that our result, Eq. (26), is also valid for range expansions with large frontiers of a fixed size that, instead of treadmilling, are growing in some direction by a single cell diameter per generation time . Since these fixed-size frontiers also do not experience inflation, we also expect them to obey our result, provided that the frontier is sufficiently large compared to the mutant cluster size.

Figure 8: A comparison of treadmilling and linearly inflating expansion survival probabilities for single cell, neutral mutations at different initial radii (in units of the cell diameter ). The survival probability for the treadmilling cases decays much faster with the time (measured in generations) and saturate at much smaller values. The plateaus at long times are consistent with , as in Eq. (27). The survival probability for an inflating front saturates earlier at a much higher value. The enhancement in the survival probability of inflating sectors for is at least 100-fold.

The neutral case is easy to analyze even for a finite size population front such as a treadmilling sphere. We know from studies of voter model dynamics (Cox and Griffeath, 1986; Liggett, 1985) that the descendants of one of the initial cells will eventually sweep the entire population. Hence, the survival probability of a single mutant cell out of a total of cells at the frontier is


We shall see in the next subsection that this treadmilling survival probability will be much smaller than the survival probability for an inflating tumor for large initial radii . This expectation is already evident upon comparing a treadmilling and an inflating front in our simulations. Note the dramatic qualitative difference in the survival probabilities with increasing time for in Fig. 8.

5.2 Linearly inflating fronts

When , our characteristic time is the time to inflate to twice the initial radius . It satisfies , where is the front speed. Then, similarly to the limit in the two-dimensional case [Eq. (15)] our dimensionless parameters in Eq. (23) reduce to


Using the limit of the general result in Eq. (24), we find


Note that when . Hence, since for , there is a jump discontinuity in the survival probability at the origin. This discontinuity occurs because we have assumed that our inflating population front is effectively infinite relative to the size of the inflating sector. Hence, when , a mutant sector will always be able to shrink back deterministically to a small enough size that it can be extinguished via genetic drift. This may happen very slowly, however, because the deterministic motion of the sector is logarithmic in time, as discussed in Section 3. Evaluation of the (small!) survival probability for negative for finite-sized fronts is a subject for future investigation (see also C).

It is instructive to study biologically relevant limits of Eq. (29). For example, if we let so that inflation becomes negligible, we recover the linear front result, consistent with Eq. (7) and Eq. (26):


Another important limit is of course (so that ). We find


Note that in our simulations, . Thus, for a single mutation and a small selection coefficient, we find, from Eq. (29), a linear increase in the survival probability as a function of :


where is Euler’s constant. In our simulations, the effective population size , so is on the order of unity as well. The inverse scaling of the zeroeth order term [] in Eq. (33) with the initial population radius is markedly different from the treadmilling result in Eq. (27). The inflation will rescue the mutation with much higher probability than the probability the mutation fixes at a treadmilling front due to the finite tumor size. Also, the logarithm in the first order term [] generates a diverging slope in as we send .

We check the approximation leading to our results with simulations in Fig. 9. There are no fitting parameters, except the genetic drift strength . We find , confirming our expectation that the genetic drift strength is on the order of unity in our simulations. The good agreement between the simulation results and the field theoretic solution justifies our approximation of an infinite population front. Indeed, from our scaling analysis in the previous section, we expect these particular finite size corrections to be small. We are inherently limited by our finite simulation time, which allows us to only simulate clusters which have grown to radii of at most 170 cell diameters.

Figure 9: The survival probability scaling function of a mutant sector for linearly inflating spherical range expansions. Clusters of mutants with initial cells start on the periphery of a population with initial radius and no other mutations. Different symbols lying on the same line correspond to data sets with different and , but a fixed in part (a) and a fixed in part (b). All range expansions were evolved up to a total cell cluster with a 170 cell diameter radius. The solid symbols in part (a) correspond to neutral range expansions with . The lines show the analytic prediction in Eq. (30). Note that at large values of and small values of in part (a), the survival probability calculated in the simulations has not yet converged to the steady-state value, so the data points lie above the analytic prediction (lines).

5.3 Marginally inflating fronts

We now examine the marginal case we identified previously in more detail. We do not have simulation results for this case, but we may examine the analytic results. Equation (24) yields a survival probability that reads


where is again the step function. This result also reduces to the linear front solution when . For small , we find the peculiar behavior


where we assume in the second line. This probability vanishes as , just as we found in Eq. (25). However, it does this logarithmically slowly. We plot this unusual marginal scaling function in Fig. 10.

Figure 10: The survival probability for a mutation living on the surface of a spherical expansion with a growth exponent [Eq. (1)]. The survival probability vanishes as , but logarithmically slowly.

6 Conclusions

We showed how geometry plays a crucial role in the evolutionary dynamics of spherical range expansions, just as it does for circular ones. We have calculated the long-time survival probability of a mutation at the spherical population frontier with either a fixed or growing radius which takes the form with arbitrary growth exponent for large compared to a characteristic time . Our primary focus has been on linearly inflating () and treadmilling () populations, both of which are relevant for tumor evolution. For linearly inflating fronts, we have , where is the front velocity, and is the time required for the front to be affected by the inflationary geometry. Inflating population fronts yield an enhanced survival probability of mutations at range expansion frontiers in two and three dimensions. We have shown how this survival probability scales with the population front radius when an isolated mutation occurs. The long time survival probability of a mutation decreases with increasing as in two dimensions and in three dimensions for linearly inflating range expansions, where is a typical cell diameter. The corresponding probabilities for treadmilling populations are and for two and three dimensions, respectively. Hence, inflation rescues mutations with a much higher probability than might be expected in populations with a fixed front size. For example, for spherical populations with thin population fronts and initial radii of ten cell diameters, , our analytic results for treadmilling populations and linearly inflating range expansions in Eq. (27) and Eq. (32), respectively, yield a 100-fold enhancement for inflating expansions: . We also verified this remarkably large enhancement in simulations, as shown in Fig. 8.

We have also gone beyond simple scaling and calculated the scaling function for the survival probability in the infinite time limit, , for spherical expansions. We applied field-theoretic techniques to calculate for arbitrary . On the other hand, using random walk theory, we have recaptured some previous results (Lavrentovich et al., 2013a) for the survival probability in circular range expansions, and generalized the results to an arbitrary [see Eq. (1)]. Using these two different techniques, we found that the scaling functions depend on key dimensionless parameters: a rescaled selection coefficient, , and a rescaled initial number of mutant cells, in two (2D) and three (3D) dimensions. We also looked at the fascinating case of marginal inflation, , where the population radius increase is almost able to rescue a neutral mutation from extinction via genetic drift. It would be interesting to apply the field-theoretic techniques described here to different evolutionary dynamics. For example, we may include more complicated, “mutualistic” interactions between cell types at tumor surfaces, which may occur when certain cell strains promote the growth of surrounding tumor cells (Marusyk et al., 2014).

The survival probability of a mutation in a three-dimensional range expansion captured by Eq. (24) has many biological applications. For example, we may evaluate a mutation survival probability for tumors, which have growth exponents that can decrease from to as the tumor approaches a “treadmilling” regime. Of course, in such tumors is time-dependent, and our results will hold only when varies slowly compared to the time it takes for the mutation’s fate to be decided. Then, if this variation in is slow, we could substitute the appropriate into Eq. (24) to find a mutation survival probability at various tumor growth stages. A precise calculation of how slow the variation of has to be and how the survival probability converges to its infinite-time limit for arbitrary would be an interesting subject for future work.

Our results may be tested experimentally in microbial populations, or, possibly, in various human cell lines (Debnath et al., 2003). For example, a spherical range expansion with might be created by growing yeast cell colonies in soft agar. Conversely, growing yeast cell pillars with a fixed cross-section, as described in Vulin et al. (2014), would generate an expansion with a fixed-size frontier and no inflation. Such expansions would obey our results in Eq. (26) for sufficiently large frontier sizes. Controlled cell cluster growth and the approach to a treadmilling regime may also be realized, for example, by confining tumor spheroids in hollow elastic microspheres (Alessandri et al., 2013). Again, we predict a strongly enhanced survival probability of mutations in the spherical expansions. This result could be checked using a fluorescing yeast strain, for example. If a circular or spherical range expansion has a negligible surface tension at the frontier, then the population front will roughen over time as the cells stochastically divide and push each other into virgin territory. This roughening can affect the evolutionary dynamics. Thus, another interesting avenue for future research, initiated by Anderson and Hallatschek (2014), is to calculate scaling functions for the survival probability of mutations on a rough population front with a non-zero selective advantage . It would also be interesting to study the regime in more detail. Since the convergence of the survival probability to its infinite time value is so slow in this regime (see C), a time-dependent theory would be necessary to find the survival probability even at very long times.

7 Acknowledgements

We thank James Glazier, Oscar Hallatschek, and Kirill Korolev for helpful discussions. One of us (DRN) would like to acknowledge helpful conversations about tumor growth with Suzanne Gaudet. This work was supported in part by the National Science Foundation (NSF) through grant DMR-1306367, and by the Harvard Materials Research Science and Engineering Center (grant DMR-1420570). Parts of this work were conducted during the workshop “Cooperation and the Evolution of Multicellularity” at the Kavli Institute for Theoretical Physics at Santa Barbara, supported by the NSF through grant PHY11-25915. Computational resources were provided by the Harvard Odyssey Research Computing group and the Center for Nanoscale Systems (CNS), a member of the National Nanotechnology Infrastructure Network (NSF grant ECS-0335765). CNS is part of Harvard University. MOL acknowledges support from NSF grant DMR-1262047.

Appendix A Response functional formalism

The response functional formalism allows us to calculate moments of the coarse-grained green cell fraction (with some position along the front at time ), such as , , etc., where the averages are over all possible realizations of the noise in the Langevin equation (18) obeyed by . Schematically, then, for any functional of , we may take an average , where we integrate over all possible noise instances , each of which has a probability weight


We perform the integration over by introducing a functional Dirac delta function. This allows us to change variables from to , ensuring that obeys Eq. (18), which we can write generally as


where is the force term and is the noise strength. Both of these are easily identifiable in Eq. (18). The integration over then becomes a Gaussian integral as follows: