Shape of population interfaces as an indicator of mutational instability in coexisting cell populations
SShape of population interfaces as an indicator ofmutational instability in coexisting cell populations
Daniel Castillo
Department of Physics & Astronomy, University of Tennessee, Knoxville, Tennessee37996, USAE-mail: [email protected]
Maxim O. Lavrentovich
Department of Physics & Astronomy, University of Tennessee, Knoxville, Tennessee37996, USAE-mail: [email protected]
Abstract.
Cellular populations such as avascular tumors and microbial biofilmsmay “invade” or grow into surrounding populations. The invading population isoften comprised of a heterogeneous mixture of cells with varying growth rates. Thepopulation may also exhibit mutational instabilities, such as a heavy deleteriousmutation load in a cancerous growth. We study the dynamics of a heterogeneous,mutating population competing with a surrounding homogeneous population, as onemight find in a cancerous invasion of healthy tissue. We find that the shape ofthe population interface serves as an indicator for the evolutionary dynamics within the heterogeneous population. In particular, invasion front undulations becomeenhanced when the invading population is near a mutational meltdown transition orwhen the surrounding “bystander” population is barely able to reinvade the mutatingpopulation. We characterize these interface undulations and the effective fitness of theheterogeneous population in one- and two-dimensional systems.
Keywords : invasion, interface roughening, mutational meltdown, spatial evolutionarydynamics
Submitted to:
Phys. Biol. a r X i v : . [ phy s i c s . b i o - ph ] S e p nterfaces as indicators of mutational instability
1. Introduction
Invasion and competitive exclusion is a common phenomenon in biology, with examplesspanning a wide range of length and time scales: An invasive land animal species maycompete with the species already present in the ecological habitat [1], microbial strainsmay compete and invade each other within a growing biofilm [2, 3], or virus strains maycompete for host resources [4]. Such competitions also exist within the tissues of variousorganisms, during development and in cancerous growth: a tumor which starts out asa small cluster of rapidly growing and mutating cells must compete with surroundinghealthy tissue [5]. In all of these examples, the spatial structure of the population mayhave a significant impact on the strain competition and evolution.Spatially-distributed populations are markedly different from their well-mixedcounterparts. Because local population sizes are small compared to the populationsize in a well-mixed test tube, genetic drift or small number fluctuations become moreimportant: Strains within spatially-distributed populations are more likely to locallyfix. Also, deleterious mutations more readily accumulate at leading edges of growingpopulations compared to well-mixed populations where natural selection would eliminatesuch variants [6, 7]. There may be mitigating factors that reduce this mutational load,however, such as the presence of an Allee effect due to strain cooperation, for example[8]. These considerations are particularly important for invading cancerous populationswhich exhibit genomic instability [9, 10] and are typically spatially heterogeneous,consisting of a wide distribution of strains [11, 12, 13, 14]. It is becoming increasinglyclear that spatial evolutionary models are necessary to understand the evolutionarydynamics of cancer cell populations [15, 16].The mutations that drive uncontrolled growth in cancerous populations are the so-called driver mutations . However, the majority of mutations are passenger mutations which have a neutral or slightly deleterious effect on the cancer cells. Such mutations areubiquitous in cancerous populations, although their importance for cancer progressionhas only recently been recognized [17]. Weakly deleterious passenger mutations canrapidly accumulate at the edges of spatially-distributed populations, and the combineddeleterious effect can lead to a cancer population collapse. Therefore, the elucidation ofthe impact of the passenger mutations may lead to new cancer therapies and a betterunderstanding of the efficacy of existing therapies [18, 19]. Indeed, an effective cancertreatment may involve increasing the mutation rate such that the passengers overwhelmthe drivers or increasing the deleterious effect of the passengers such that the driversare no longer able to sustain tumor growth. The accumulation of deleterious mutationsleading to population collapse is termed “mutational meltdown” [20, 21]. Already thereis evidence that cancer therapies may be developed that target passenger mutations toexpose vulnerabilities in cancer growth [22] and that cancer cell lines are particularlyvulnerable to mutational meltdown [23].In this paper, we develop a simple spatial model of the invasion of a cellularpopulation (i.e., a “bystander”) by another population (i.e., an “invader”) that is nterfaces as indicators of mutational instability nterfaces as indicators of mutational instability d = 1 + 1 and d = 2 + 1 evolutions, respectively. The +1 indicates the time dependence. Our focus here is on thecompetition between multiple strains within a population, so for simplicity we consider flat habitats which do not change shape as the population evolves. For d = 1 + 1 ,such a habitat may be a coast, a river bank, or a thin duct. For d = 2 + 1 , thestrains may be in a microbial population growing on a flat surface or in an epithelialtissue. Another possibility is that the population in which the strains compete is theleading edge of a range expansion. In this case, we assume that the population growthis confined to a thin region at the edge, which remains flat during the range expansion.This approximation will hold as long as there is a sufficient surface tension keeping thepopulation edge uniform which occurs in yeast cell colonies, for example [38]. However,if the population edge roughens over time, the roughening will generically change thegenetic sector motion [39], an analysis of which is beyond the scope of the current work.Note that for a cellular population at the edge of a range expansion, the +1 dimension represents the direction of the range expansion. So, in other words, for d = 1 + 1 , the strains we study may live along a thin, effectively one-dimensional edge ofa two-dimensional range expansion (e.g., a thin microbial colony grown on a Petri dish).For d = 2 + 1 evolution, the population may be the effectively two-dimensional, flatedge of a three-dimensional range expansion. A more realistic scenario is perhaps the d = 3 + 1 -dimensional case where a population embedded in three dimensions evolves intime with various strains within the population competing for the same space. Althoughwe do not study this case specifically here, the lower dimensional cases provide intuitionsfor considering this scenario. Also, if the geometry of the three-dimensional populationhas a large aspect ratio, then our one and two (spatial) dimensional models will describethe behavior of cross sections of the population. A similar kind of dimensional reductionwas recently employed for describing bacterial competition in three-dimensional channels[40]. Previous work has shown that range expansions develop frontiers with enhancedroughness when the population is near a phase transition in its evolutionary dynamics(e.g., at a mutational meltdown transition [39] or near the onset of mutualistic growth[41]). In this work, we consider invasion frontiers which are markedly different as theinvader grows into a surrounding population which may reinvade if the invader growthrate decreases. So, the invasion front speed v will depend on the relative growth ratesof the two populations and may vanish or change signs. In other words, the competitioninterfaces studied here have a variable speed, unlike a range expansion in which a nterfaces as indicators of mutational instability range shift , in which thepopulation growth is limited by the environment [42]. We will see that in the case of d = 2+1 -dimensional expansions, the average speed v of the invasion front will influencethe interface roughness.We will also show here that, like the range expansion, the invasion frontierdevelops an enhanced roughening at the mutational meltdown transition of the invaderpopulation. However, unlike a range expansion frontier, the roughening is more subtle,especially in the d = 2 +1 -dimensional case where the relative growth difference betweenthe invader and bystander populations (and consequent invasion front speed v ) alsoinfluences the roughening dynamics. The invasion frontier does not maintain a compactshape, and isolated pieces of the invading population may pinch off and migrate intothe surrounding “bystander” population, especially when the growth rates of the twopopulations are nearly equal. In this paper we will discuss these issues and connect theshape of the undulating frontier to the evolutionary dynamics of the invading population.The paper is organized as follows: In the next section, we present our lattice modelfor invasion by a mutating population for d = 1 + 1 and d = 2 + 1 -dimensional cases. InSection 3 we briefly review the nature of the mutational meltdown transition that mayoccur in the unstable invading population. In Section 4 we study the survival probabilityof the invading strain and construct phase diagrams characterizing whether or not theinvasion is successful as a function of the internal mutation rate µ and the selectiveadvantage of the various cellular strains. In Section 5, we analyze the rougheninginvasion front near the mutational meltdown transition for the invading population.Finally, we conclude with a discussion of the implications of our results in Section 6.
2. Model
We consider a simple lattice model, in the spirit of the Domany-Kinzel cellularautomaton [43, 44], of invasion of a stable population by a mutating invader consistingof two species, a fast-growing and a slow-growing strain into which the fast-growingone can mutate. We set the fast-growing strain growth rate to unity Γ f = 1 withoutloss of generality, so that time is measured in generation times τ g of the fast-growingstrain. The slow-growing strain within the invader population will have growth rate Γ s = 1 − s , where ≤ s < is a measure of the deleterious effect of the mutation.In a tumor or microbial colony, we know that the initial cluster of growing cells mayencounter other cells (e.g., surrounding healthy tissue or competing microbial strains).So, we have a third species representing the “bystander” population. The bystander willnot mutate, but will be able to displace the mutating population via cell division. Weset the bystander growth rate to Γ b = 1 + b − s , with b the selective advantage of thebystander over the slow-growing invader strain.The internal dynamics of the invading population (i.e., the mutation rate µ andselection parameter s ) will influence how the invader interacts with the bystander strain, nterfaces as indicators of mutational instability µ or s leading to an overall fitness decrease for the invading population,as might happen in a cancerous tissue during a course of therapy that increases thedeleterious mutation rate of the cancerous cells. We focus on the region between themutating population and the bystander, which we call the invasion front. As we willshow, when the invader is close to mutational meltdown, the invasion front develops anenhanced “roughness.” Figure 1. (a) Update rules for the bystander model for a population in d = 1 + 1 dimensions. Each generation is evolved by allowing for two cells from the previousgeneration to compete for an empty lattice site, as shown by the arrows. Theprobability of occupation by a cell of a type i = s, f, b is proportional to its growthrate Γ i , where s is the slow growing black strain, f is the fast growing red strain, and b is the yellow bystander. If a red (fast-growing) cell is placed in the empty spot, thenit in addition has a probability µ of mutating to the slower-growing black strain. (b)For a two-dimensional population ( d = 2 + 1) the generations are evolved on staggeredtriangular lattices, as shown. This time, three cells compete to divide into empty latticesites. Otherwise, the dynamics are the same as the d = 1 + 1 case. The specific lattice model rules are as follows: In both one- and two-dimensionalpopulation scenarios, we consider a three-strain model in which a single “bystander”strain (yellow cells in Fig. 1) grows in the presence of a fast-growing invading strain (redcells in Fig. 1) that can mutate to a more slowly growing strain (black cells in Fig. 1).These cells occupy a single lattice location, as shown in Fig. 1. During each generation(cell division) time τ g , the lattice of cells is regenerated by allowing for adjacent cellsto compete and divide into empty sites representing the next generation of cells. In arange expansion context, these empty sites would represent unoccupied territory at thefrontier. Alternatively, these updates can represent a turnover of cells due to birth anddeath within the population. After all empty sites in the next lattice have been filled(rows for d = 1 + 1 and sheets for d = 2 + 1 as shown in Fig. 1(a) and (b), respectively),the process can be repeated, generating a sequence of successive generations of the nterfaces as indicators of mutational instability µ .So, our model contains just three parameters: the deleterious effect of the mutation s , the selective advantage b of the bystander population over the slow-growing strain,and the mutation rate µ . We will be interested in the regime < Γ s ≤ Γ b ≤ Γ f ,for which ≤ s < and ≤ b ≤ s . In this case, the bystander either invades theslow-growing strain or is invaded by the fast-growing strain, as illustrated in Fig. 3 fora d = 1 + 1 simulation. Note that this reinvasion by the bystander population makesthe invasion frontier markedly different from, say, a range expansion frontier. In a rangeexpansion frontier, the range expansion always moves in one direction according to thegrowth rate of the total population. Here, the interface between the bystander andinvader can move in different directions or even remain, on average, stationary. We willsee that this aspect will be important when studying the roughness of the interface.Our parameterization allows us to tune the dynamics of the black/red mutatinginvader population separately. As we will analyze in the next section, the invader hasan internal “mutational meltdown” at which the fast-growing red strain is removed fromthe population due to mutation. This occurs for µ (cid:38) s in d = 1 + 1 dynamics and µ (cid:38) s ln s in d = 2 + 1 dynamics ( µ > s in well-mixed populations). Note that animportant limitation of our model is that we assume cells divide into adjacent spotson the lattice so that cell motility is essentially absent (apart from the short-range cellrearrangements occurring due to the cell division). This is a reasonable approximationfor certain microbial populations such as yeast cell colonies [29] or small, avascular solidtumors where cells primarily proliferate [45].
3. Mutational meltdown
Let us first focus on the invading population and perform a simple analysis of the internaldynamics. The invader consists of two strains: one fast-growing red strain and a secondslow-growing black strain into which the fast-growing strain mutates with rate µ per cellper generation. In the parameter space ( µ, s ) , we find two distinct phases [43]: In onephase, the fraction of the fast-growing strain remains positive after many generations, ρ f ( t → ∞ ) > ; we call this phase the active phase . In the other phase, called the absorbing or inactive phase , the fast-growing strain eventually completely dies out and ρ f ( t → ∞ ) = 0 . There is a line of continuous phase transitions ( µ ∗ , s ∗ ) which definesthe boundary between the two phases. Examples of these phases, and the critical region( µ ≈ µ ∗ , s ≈ s ∗ ), are shown in Fig. 2 where we have removed the bystander populationin order to see the internal dynamics of the invader.We can understand this transition in a well-mixed population (a mean-fieldapproach). Consider just the invading, mutating population. The fraction ρ f ≡ ρ f ( t ) of nterfaces as indicators of mutational instability Figure 2.
Simulated sectors of a black/red mutating population invading a bystanderpopulation (surrounding white area), initialized by a single red cell in a (a) one-dimensional and (b) two-dimensional population. The populations are evolved forabout 100 generations, with the time direction indicated. In (a), we indicate themotion of the invasion front (which in this case is a point) between the two populations.In (b), the invasion front would be the complicated boundary between the black/redpopulation and the white space at each time slice t along the vertical direction. Thephases of the internal dynamics of the invading population (inactive, critical, and activephases) are indicated. In the inactive phase, the red, fast-growing mutant is lost fromthe population over time. As the invader population transitions from the inactiveto the active phase in which the red strain is maintained, the invasion front exhibitsenhanced undulations. the fast-growing strain within the mutating population changes according to: dρ f dt = sρ f (1 − ρ f ) − µρ f , (1)which approaches ρ f ( t → ∞ ) = 1 − µ/s for µ < s , and ρ f ( t → ∞ ) = 0 for µ > s . The line µ = s is our set of critical points ( µ ∗ , s ∗ = µ ∗ ) . For a spatially distributed population,small number fluctuations or “genetic drift” dramatically alters the shape of the phaseboundary: The phase transition occurs for µ ∗ ∼ ( s ∗ ) in one-dimensional populations(such as at the edge of a growing biofilm [29]) and µ ∗ ∼ s ∗ ln( s ∗ ) for two-dimensionalpopulations [46]. This phase transition, called a “mutational meltdown,” is known tobelong to the directed percolation (DP) universality class [31]. For the particular latticemodel we consider here, this has been explicitly verified [43] by mapping the model tothe well-studied Domany-Kinzel cellular automaton [44]. The efficacy of this simplemodel has been demonstrated in a synthetic yeast strain, for which the parameters µ and s could be tuned over a broad range encompassing the DP phase transition [29].Note that when the population approaches the mutational meltdown transition, theslow-growing black strain within the population begins to take over. In the active phasein Fig. 2, the black strain makes finite-sized, small patches within the red population. nterfaces as indicators of mutational instability Figure 3. A d = 1 + 1 simulation of a red/black mutating population invading abystander yellow one. Here, the yellow strain grows faster than the black strain butslower than the fit red strain. The invasion front between the black/red population andthe bystander strain can be characterized by a random walk with alternating bias. Theyellow strain invades the black patches and is invaded by the red patches. The sizesof the red and black patches are controlled by the internal dynamics of the black/redstrain. We illustrate the characteristic sizes ξ ⊥ and lifetime ξ (cid:107) of the black patches. The two-species model may be generalized to include an arbitrary number ofpossible mutations, and such models have been shown to exhibit critical behaviorthat deviates from the DP universality class, but the loss of the fittest mutant in thepopulation is still well-described by DP [46]. The multi-species generalization has manyadditional interesting phenomena such as multi-critical behavior [47], which would allowfor interesting extensions of the work presented here. In this paper, for simplicity, weshall focus on the fittest mutant in an invading population with just two species. Thefittest strain could represent, for example, a driver mutation which has swept through acancerous tissue. The driver strain could then acquire deleterious mutations over timewith rate µ . We will focus here on just the initial loss of fitness, characterized by asingle mutation to the slower growth rate Γ s = 1 − s . nterfaces as indicators of mutational instability
4. Invasion probabilities
We now construct a phase diagram for successful invasion of the bystander strain bythe mutating invader (Fig. 4). We initialize a well-mixed population of equal parts ofthe mutating red and bystander yellow strains ( ρ b = ρ f = 1 / on the lattice, and thencalculate the density ρ m = ρ f + ρ s of the mutating red/black population at long times t . If ρ m → and the red/black population dies out at long times, then the “invasion”is unsuccessful. Otherwise, ρ m approaches a non-zero value and the invasion succeeds.The results are shown in Fig. 4 for d = 1 + 1 and d = 2 + 1 , where we see the two distinctphases. We also include the phase boundary µ = s − b for the well-mixed population(blacked dashed line in Fig. 4), which we derive in the next subsection. Note how faraway the well-mixed population transition line is from the actual transition in a spatialpopulation. The genetic drift associated with the spatial populations suppresses theinvasion by the red/black mutating population.We also know that as we approach this mutational meltdown transition at µ = µ ∗ fora fixed s [vertical dashed lines in Fig. 4], the characteristic sizes ξ ⊥ and the characteristiclifetimes ξ (cid:107) of black, slow-growing strain clusters diverge according to ξ ⊥ ∼ ∆ − ν ⊥ and ξ (cid:107) ∼ ∆ − ν (cid:107) , where < ∆ < is the distance from the phase transition in the ( µ, s ) planeand ν ⊥ and ν (cid:107) are critical exponents associated with directed percolation ( ν ⊥ ≈ . , ν (cid:107) ≈ . for d = 1 + 1 and ν ⊥ ≈ . , ν (cid:107) ≈ . for d = 2 + 1 [31]). We illustratethe sizes ξ ⊥ , (cid:107) in Fig. 3. The black patches will interact differently with the bystanderthan with the red patches of the fast-growing strain. As the patch sizes ξ (cid:107) , ⊥ diverge(when ∆ → ), they would have a more pronounced effect on the invasion dynamics.In particular, there will be larger regions over which either the yellow strain invades ablack patch, or a red patch invades the yellow bystander. This will increase the amountof “wiggliness" of the invasion front between the bystander and the mutating red/blackpopulation. We will see in the following that there is a significant enhancement of theroughness as we approach the mutational meltdown transition. To understand the behavior of this three-species model, we first briefly describe whathappens in a well-mixed (mean-field) context. Consider the time-evolution of thefractions ρ f , ρ s , ρ b of the fast-growing, slow-growing, and bystander strains, respectively.For a fixed total population size, we have ρ f + ρ s + ρ b = 1 . Given our growth rates Γ f = 1 , Γ s = 1 − s , and Γ b = 1 − s + b , we may define corresponding selection coefficientscharacterizing the competition between pairs of strains: s bs = Γ b − Γ s (Γ b + Γ s ) / b − s + bs fb = Γ f − Γ b (Γ f + Γ b ) / s − b )2 − s + b (2) s fs = Γ f − Γ s (Γ f + Γ s ) / s − s , nterfaces as indicators of mutational instability < s < and < b < s . In terms of these selectioncoefficients, the equations for the time-evolution of the bystander and fast-growing strainfractions ρ b,f ≡ ρ b,f ( t ) in a well-mixed population are (cid:40) ∂ t ρ b = s bs ρ b (1 − ρ b − ρ f ) − s fb ρ f ρ b ∂ t ρ f = s fs ρ f (1 − ρ b − ρ f ) + s fb ρ f ρ b − µρ f . (3)If ρ b = 0 , we recover the two-species dynamics of the invader population with a directedpercolation-like process between the fast-growing and slow-growing strains. We can alsoverify that there is no sensible stable fixed point where both the bystander populationand the invader coexist. Instead, if µ > s − b , then the bystander will sweep the totalpopulation and ρ b ( t ) → with increasing time t . Otherwise, if µ < s − b , the invasionby the mutating population is successful and we find ρ b ( t ) → over time. Moreover, if µ > s , we get a collapse of the fast-growing strain [ ρ f ( t ) → ], and then the bystanderstrain will win out since Γ s < Γ b . So, the mutational meltdown transition of the invaderpopulation occurs when µ = s in this mean field analysis.The mean-field analysis tells us that we should expect a critical surface in the ( µ, s, b ) parameter space given by µ = s − b separating a region of successful ( µ < s − b ) or failed ( µ > s − b ) invasion of the bystander strain by the mutating invader (whichitself may undergo a mutational meltdown when µ > s ). As we have already seen,the spatially-distributed populations also have this critical surface but the enhancedgenetic drift suppresses the phase space for successful invasion. To add the effects ofgenetic drift and the spatial distribution of the population to Eq. (3), we would have toincorporate a spatial diffusion term ∇ ρ b,f in each of the equations and stochastic noiseterms describing the birth/death dynamics (see [46] for a more detailed description).These additional terms significantly modify our mean field equations and introducedifferent phenomena, such as propagating waves (moving population interfaces) whichwe will analyze in Section 5.We can get a better approximation to the critical line for the d = 1 + 1 case thanthat given by the mean field theory by considering a single domain wall. We expect thetotal length along the domain wall to be split into sections of average length (cid:96) s wherethe slow-growing strain competes with the bystander and sections of average length (cid:96) f where the fast-growing strain competes with the bystander. During each generationtime τ g the domain wall can move by a single cell length a . So, in the fast-growingstrain segments, we expect the fast-growing clusters to out-compete the bystander andprotrude by an amount τ g (cid:96) f ( s − b ) /a , with τ g the generation time. Similarly, we expectthe slow-growing clusters to be out-competed by the bystander and recede by an amount τ g (cid:96) s b/a . At the phase transition, we expect these competitions to cancel each other outas the mutating invading population is barely able to invade the bystander in this case.Hence, we should have (cid:96) f ( s − b ) ≈ (cid:96) s b so that b ≈ (cid:96) f s(cid:96) s + (cid:96) f ≈ φ f s, (4)where φ f is the red-cell (fast-growing) fraction of the mutating invader population. We nterfaces as indicators of mutational instability Figure 4.
Phase diagrams for (a) the d = 1 + 1 case and (b) the d = 2 + 1 case,calculated by initializing a well-mixed population of the red and yellow strains andevolving the whole population for t ≈ generations. In (a) we use a one-dimensionalpopulation of L = 5000 cells and average over 256 runs of the evolution. In (b) wehave a two-dimensional population with L cells, where L = 500 . Here we average over40 evolution runs. In both cases we set s = 0 . . After evolving for generations,we calculate the red/black mutating fraction of the total population: ρ m = ρ f + ρ s .In each phase diagram, the black dashed line corresponds to the mean field prediction µ = s − b . The green and white dashed lines correspond to the improved predictions[see Eqs. (5 a ), (5 b ), (6 a ), (6 b )] for µ (cid:28) µ ∗ and µ ≈ µ ∗ , respectively, that take intoaccount the spatial structure of the population. We also indicate the line µ = µ ∗ along which we find a mutational meltdown transition within the invading red/blackpopulation. now can use Eq. (4) to predict the critical line in ( µ, b ) -space for a fixed s in two limitingcases: µ (cid:28) µ ∗ and µ ∼ µ ∗ , where µ ∗ is the critical value for µ for the specific fixed valueof s at which we get the mutational meltdown transition within the red/black invadingpopulation.To complete the derivation, we just need an estimate for the fraction of fast-growingstrain φ f . First, when µ (cid:28) µ ∗ , the invader population is in the “active” phase, and thepatches of black slow-growing strain are small and rarely collide, as shown in Fig. 2. In d = 1 + 1 -dimensions, the boundaries of these black patches are well-described by pairs nterfaces as indicators of mutational instability φ f ≈ − A µ/s [43, 48, 28]. The amplitude A is model-dependent, and we have A ≈ . for this Domany-Kinzel-type model,consistent with previous results [46]. As µ → µ ∗ , however, the fast-growing strainvanishes ( φ f → , and we have to make another approximation. From the randomwalk model, we expect that φ f vanishes when µ = µ ∗ ∼ s . Then, when µ ≈ µ ∗ , wewould be near a directed percolation (DP) phase transition where φ f serves as an orderparameter. The order parameter vanishes according to φ f ≈ A ( µ ∗ − µ ) β , where A isan amplitude that will depend on s and β ≈ . is a DP critical exponent [31]. Wemay now use Eq. (4) to make an estimate of the critical value of b for d = 1 + 1 : b = s (1 − A µ/s ) ( µ (cid:28) µ ∗ ) (5 a ) b = sA ( µ ∗ − µ ) β ( µ ≈ µ ∗ ) (5 b )where A and A can be calculated numerically from separate simulations of the two-species model. These improved estimates are plotted onto the phase diagrams in Fig.4(a) (green dashed line for the µ (cid:28) µ ∗ case and white dashed line for the µ ≈ µ ∗ case).For d = 2 + 1 -dimensional evolutions, the situation is more complicated becausethe patches of the invader strain no longer have a compact shape describable by asimple random walk [see Fig. 2(b)]. However, we expect that the bystander populationmay reinvade the invading mutating population when b > sφ f because, much like inthe d = 1 + 1 case, the average growth rate of the invader strain is approximately Γ ≈ φ f Γ f +(1 − φ f )Γ s = 1 − s + φ f s . The bystander strain has growth rate Γ b = 1 − s + b , sowe see that the growth rates are equal when b = φ f s . We now just need estimates for φ f for d = 2 + 1 . When µ (cid:28) µ ∗ , previous work [46] has shown that φ f ≈ − A µ ln( s/s ) /s ,with A ≈ . and s ≈ some model-dependent parameters. Conversely, when µ ≈ µ ∗ , we again find a DP transition with φ f ≈ A ( µ ∗ − µ ) β with critical exponent β ≈ . for d = 2 + 1 [31]. The corresponding estimates are for d = 2 + 1 : b = s [1 − A µ ln( s/s ) /s ] ( µ (cid:28) µ ∗ ) (6 a ) b = sA ( µ ∗ − µ ) β ( µ ≈ µ ∗ ) (6 b )These two approximations are plotted in Fig. 4(b), with µ (cid:28) µ ∗ in the green dashedline and µ ≈ µ ∗ in the white dashed line.The phase diagrams in Fig. 4 were constructed with simulations using mixed initialconditions; the first generation of cells on the lattice were populated by an even mixtureof fast-growing (red) cells and bystander (yellow) cells. These phase diagrams are heatmaps corresponding to the density ρ m of the mutating population (red/black strains)after many generations. Our simulations were performed for t ≈ generations, whichyields the steady state solution for the mutating population fraction ρ m for the vastmajority of points on the phase diagram in Fig. 4, except for points very near thephase transition line. Note that our improved mean field estimates based on directedpercolation and the random walk theory (white and green dashed lines, respectively)do a reasonable job of approximating the shape of the phase boundary, especially when µ ≈ and our system reduces to a simple competition between fast-growing red cells nterfaces as indicators of mutational instability µ ≈ µ ∗ , where we find the mutational meltdown transition of the invader populationwhich is in the directed percolation universality class. Figure 5.
Long-time survival probabilities P surv of clusters generated from a singlemutating red cell in a yellow bystander population for (a) d = 1 + 1 -dimensional rangeexpansions for s = 0 . (left panel) and s = 0 . (right panel), after t = 10 generationson a lattice with size L = 5000 cells averaged over 128 runs; and (b) d = 2 + 1 -dimensional range expansions for s = 0 . (left panel) and s = 0 . (right panel), after t = 2 × generations on a lattice with size L cells with L = 500 , averaged over 256runs. We show the expected transition shape near the µ ≈ µ ∗ DP transition in thewhite dashed line [Eqs. (5 b ),(6 b )]. The black dashed line is the transition position fora well-mixed population. The green dashed line is an improved mean-field estimate ofthe transition discussed in the main text [Eqs. (5 a ),(6 a )]. Another biologically interesting quantity to look at is the survival probability P surv of the progeny of a single red cell invader in a population of yellow bystander cells as t → ∞ . Such a survival probability would represent the probability of tumor invasion,for example, from a single mutated cell (i.e, a cell with a newly-acquired driver mutation)within an otherwise healthy population. If the bystander is replaced by the slow-growingstrain and we have a two-species evolution, then the evolution will be exactly the sameas a directed percolation with a “single-seed” initial condition. We would then have P surv ∝ ρ f due to the rapidity reversal symmetry of directed percolation [31]. In otherwords, the survival probability tracks the behavior of the fraction ρ f of the fast-growingstrain in a different simulation where the initial condition is a well-mixed population(or a population of just the mutating, fast-growing strain). In the three-species modelwe consider, there is no rapidity symmetry due to the presence of the bystander strain. nterfaces as indicators of mutational instability P surv vanishes on the same criticalsurface as the fraction ρ m (plotted in Fig. 4) because the invader strain will not beable to invade if the fast-growing strain is lost from the population. We show in Fig. 5the survival probability P surv at long times, which does indeed vanish at approximatelythe same place as ρ m in Fig. 4. So, the approximations we used to estimate where ρ m vanishes serve as good predictors of the transition of the survival probability, as well.We also show the phase boundary at a smaller values of s ( s = 0 . ) in the right panelsof Fig. 5. Note that our estimates work for the phase boundary in this case, also.
5. Roughening invasion fronts
We now study the shape of the interface between the mutating and bystanderpopulations. When either of the populations is invading the other, the invasionfront behaves as a noisy Fisher-Kolmogorov-Petrovsky-Piscounov wave [49, 50]. Mostprevious studies of such waves have focused on competition between two homogeneouspopulations or the range expansion of a population into virgin territory. The noise playsa crucial role here [51], strongly modifying, for example, the wave speed. Also, in the(exactly soluble [52]) d = 1 + 1 case, there is a diffusive wandering of the front aroundits average position.For d = 2 + 1 , the situation is more complicated, but generally the noisy wave frontwill have a characteristic roughening. This roughening falls in the Kardar-Parisi-Zhang(KPZ) universality class [53], although observing the predicted scaling behavior of thisclass is challenging for noisy Fisher waves [54, 55]. For example, for the KPZ class ofinterfaces, the characteristic size σ w of the interface should grow as σ w ∝ t / . However,a basic analysis of the noisy Fisher waves [56] is more consistent with σ w ∝ t . , whichis also what we observe in our model. Although this apparent discrepancy has beenexplained, the proper recovery of the KPZ exponents requires a deeper analysis outsidethe scope of the current paper [54]. So, for our simulations, we will find consistencywith previous analyses of noisy Fisher waves and leave the more extensive analysis of theinterface shape scaling for future work. Also, as the speed of the invasion goes to zero,we expect a transition to a different, “voter-model” [57] interface coarsening behavioras both the mutating and the bystander populations become stable and do not invadeeach other (on average). The interface roughens in a different way in this “critical” case,with the characteristic size σ w of the interface increasing diffusively as σ w ∝ t / . Wewill observe such crossovers in our simulation results.We discuss these issues in more detail below and show that our model exhibits arange of behaviors depending on the invasion speed and the proximity of the mutatingpopulation to the meltdown transition. These invasion waves are examples of “pulled”wavefronts [58], which are driven by the growth (invasion) at the leading edge of thewave. Various aspects of such wave fronts are reviewed in, e.g., Ref. [59]. We shall seein the following that adding mutations to one of the populations significantly modifiesthe expected pulled front wave behavior and, in the d = 1 + 1 case, introduces a super- nterfaces as indicators of mutational instability d = 2 + 1 case presents an even richer setof behaviors depending on the mutation rate and relative fitness of the mutating andbystander populations. Our purpose here will not be the particular value of scalingexponents, but rather general features of the roughening dynamics such as a change inroughening due to the internal evolutionary dynamics of the invading strain. d = 1 + 1 -dimensional invasions In d = 1 + 1 , domain walls between the bystander and the invading populations canbe characterized by a random walk with alternating bias (when Γ s < Γ b < Γ f ) asthe bystander will invade the slow-growing species and be invaded by the fast-growingspecies within the mutating population. As we approach the mutational meltdowntransition, the average size of clusters of the slow-growing strain will diverge as expectedfrom the directed percolation transition. In Fig. 6 we see a comparison of two domainwalls for d = 1 + 1 . At the bottom of the figure, Fig. 6(b), we see a domain wallwhere the mutating red/black population is far from the two-species phase transition.In this case, the black patches in the population are small and do not influence themotion of the invasion front much. Conversely, in Fig. 6(a), we see a domain wallwith the mutating population near a mutational meltdown. In this case, there is anenhancement of the “roughness” of the domain wall as the large black patches createmore regions of alternating bias in the domain wall between the yellow bystander andthe red/black invading population.To obtain a more quantitative estimate of this roughening effect, we set up asimulation with initial conditions that include a sharp boundary between the bystanderand the mutating population: the bystander occupies lattice sites i ≤ L/ , and all otherlattice sites i > L/ are occupied by the mutating population (taken to be all red, fast-growing cells initially). We then track the position x ( t ) of the invasion front over time.We measure the roughness of the front by calculating the variance of the position: (cid:104) [ w ( t )] (cid:105) = (cid:104) [ x ( t ) − (cid:104) x ( t ) (cid:105) ] (cid:105) = (cid:104) [ x ( t )] (cid:105) − (cid:104) x ( t ) (cid:105) , (7)where we average over sufficient runs to ensure convergence. In the case of a domainwall between just two strains, perhaps with a difference in growth rates, the domainwall performs a biased random walk [43]. Therefore, we may expect that our position x ( t ) also performs a diffusive motion in time. The number fluctuations at the boundaryintroduce a stochasticity to the motion, while the difference in growth rates providesa deterministic bias. So, for a boundary between two strains, we expect the variance σ w ( t ) ≡ (cid:112) (cid:104) [ w ( t )] (cid:105) to satisfy σ w ( t ) ≈ (cid:112) D w t, (8)with D w a diffusion coefficient for the domain wall. Indeed, x ( t ) itself should perform abiased random walk and we may extract the diffusion constant D w from a time series ofthe position x ( t ) performing a time average of the squared displacements of the interface[60]. We did this for the simulations shown in Fig. 6. We see that when the population is nterfaces as indicators of mutational instability Figure 6.
A picture of the domain wall in a d = 1 + 1 -dimensional evolution betweenthe invading black/red population and the bystander yellow population (a) near thetwo-species (DP) phase transition ( µ = 0 . , b = 0 . , s = 0 . ) and (b) far fromthe mutational meltdown DP phase transition ( µ = 0 . , b = 0 . , s = 0 . ). We seethat the roughness of the domain wall becomes enhanced for (a) with a domain walldiffusion coefficient D w ≈ . as compared to D w ≈ . for (b) (calculated using thetime-averaged squared displacements of the wandering interface in the correspondingfigures [60]). The width of the population is L = 300 cells, and we show the evolutionover generations. This long evolution time allows for an observation of the domainwall wandering. However, as a result, the cells are compressed along the time directionin this figure. near a mutational meltdown at µ ≈ µ ∗ [Fig. 6(a)], the observed diffusivity is much largerthan for a population far away from this transition [Fig. 6(b), with µ (cid:28) µ ∗ ]. However,a proper measurement of D w requires an ensemble averaging over many simulation runsand also a longer time series.We shall see in the following that a more detailed analysis of the boundary motionwill show that x ( t ) actually performs a super-diffusive motion near the mutationalmeltdown µ ≈ µ ∗ , with displacements satisfying σ w ( t ) ∝ t ν , with ν > / . Super-diffusivity is not uncommon in spatial population dynamics: In a range expansion, forexample, the roughness of the expansion front may contribute to the motion of sectors ofstrains, introducing super-diffusivity to the sector boundary motion [27]. However, thissuper-diffusivity depends on the conditions of the growth, and a diffusive motion oftenserves as a reasonable approximation [61, 62]. Moreover, if we are just thinking aboutpopulations living in a fixed one-dimensional geometry, then we expect sector boundarymotion to be diffusive. We will find diffusive motion of our sector boundaries everywherein the ( s, b, µ ) parameter space, except near the mutational meltdown µ ≈ µ ∗ where thesector motion becomes super-diffusive. nterfaces as indicators of mutational instability Figure 7.
For varying values of s = 0 . (top row), . (middle), and . (bottom), weshow how, as we move along the critical line starting from ( µ, b ) = (0 , s ) towards thetwo-species critical point at ( µ, b ) = ( µ ∗ , , the domain walls separating a bystanderpopulation from a mutating one acquire super-diffusive behavior. The phase diagramson the left illustrate where we calculate the roughness exponent ν ( t ) on the right. Twolimiting values of the exponent are indicated with dashed lines in each plot: a diffusive ν = 0 . (lower lines) and a super-diffusive, directed percolation value ν ≈ . (upper lines). The phase diagrams were created from simulations with L = 5000 , t = 10 generations, and averaged over runs. The exponent curves on the rightwere created from simulations with L = 5 × , t = 10 generations, and runs. Let us now analyze the dynamics in more detail. For a domain wall or invasionfront between our mutating, heterogeneous invader population and the homogeneousbystander, the slow- and fast-growing strain patches of the invader will interactdifferently with the bystander. We can analyze how this impacts the domain wall motionby studying the standard deviation σ w ( t ) = (cid:112) (cid:104) [ w ( t )] (cid:105) , averaged over an ensemble ofsimulation runs. We sample our evolved population at times t = t i and calculate theeffective exponent associated with the interface width: ν ( t = t i ) ≡ ln[ σ w ( t i ) /σ w ( t i − )]ln[ t i /t i − ] , (9) nterfaces as indicators of mutational instability Figure 8.
The diffusion coefficient, D w , of the boundary between healthy andcancerous populations, measured for various values of µ approaching µ c with s = 0 . and b adjusted so that the mutating and bystander populations remain relativelyneutral. We see that as µ → µ ∗ ≈ . , D w diverges as the domain wall becomessuper-diffusive. The coefficients were calculated from simulations with L = 3 × , t ≈ × generations, and averaged over N = 1280 runs. where we choose t i /t i − ≈ . The effective exponent ν ( t ) approaches a limiting value atlong times. Moreover, any super-diffusive enhancement to the roughness would be seenas a limiting value ν > / . The exponent is plotted for various values of ( s, b, µ ) inFig. 7. We find that there is an enhanced, super-diffusive roughness ( ν > / wheneverthe mutating population is close to the mutational meltdown (DP) transition at µ = µ ∗ [along the vertical line in the phase diagram in Fig. 4(a)]. The enhanced value of ν nearthe DP transition may be understood by considering the limiting case b = 0 . In thiscase, the bystander population and slow-growing strain within the mutating populationwill grow at the same rate, so then an initial condition with a single red fast-growingcell in a population of yellow bystander cells will expand as it would in a standardDP process with a single seed initial condition. Hence, the standard deviation σ w ( t ) scales with the DP dynamical critical exponent: σ w ( t ) ∼ t ν DP , with ν DP = 1 /z ≈ . for d = 1 + 1 [31]. This exponent is indicated with the upper dashed line in Fig. 7.Introducing a non-zero b > should not change the situation much; we would onlyexpect a difference in the bias of the domain wall motion.Away from the DP transition, the invasion front has a diffusive behavior, with σ w ( t ) = √ D w t . The diffusion constant D w may be measured and serves as a goodindicator of the mutational meltdown transition because D w should diverge as µ → µ ∗ for fixed b and s . This is illustrated in Fig. 8 for s = 0 . and values of b along thephase transition boundary. In this d = 1 + 1 -dimensional case, the value of b , accordingto our simplified analysis, does not change the wandering behavior of the domain wallsas it only serves to change the domain wall bias . This hypothesis is consistent withthe data shown in Fig. 7, where the red squares and purple crosses, despite havingvery different b values, both exhibit super-diffusive exponents ν ( t ) > / at long timesbecause both points are near the mutational meltdown transition at µ = µ ∗ . We do seesmall b -dependent differences, but these may be due to the finite time of our simulations. nterfaces as indicators of mutational instability s in Fig. 7, that the exponent ν ( t ) has not yet saturated to its long-time, limiting value within the simulation time.In any case, it is clear that we find super-diffusive motion whenever µ approaches µ ∗ in the various panels of Fig. 7. As we shall see in the next subsection, the situationchanges dramatically for the d = 2 + 1 -dimensional case where the interface betweenthe bystander and mutating populations will behave differently depending on the valueof b . d = 2 + 1 -dimensional invasions Figure 9.
Snapshots at t = 8192 generations of an initially circular cluster of thefast-growing strain (initial diameter of 400 cells). In (a), the internal parameters ( µ, s ) of the cancerous population are set such that there is an average bias of the interface sothe growth speed | v | > . In (b), we have | v | = 0 , and thus we are on the critical line ofthe 3-species phase diagram. In this case, the interface between populations dissolves,thus our characterization of the interface roughness becomes more complicated in the d = 2 + 1 -dimensional case than it was for the d = 1 + 1 -dimensional case, where theboundary was a single point performing a biased (super)diffusive random walk. For a two-dimensional population, the invasion front is no longer a simple point, butrather an undulating line between the bystander and the mutating red/black population.Moreover, this line can thicken as pieces of the invader population pinch off and migrateinto the bystander population due to rearrangements induced by the cell division. Thisdissolution of the front is more prominent when the bystander and the invader haveapproximately the same growth rates. An example of the complicated frontier shapesare shown in Figs. 9 and 10. In Fig. 9 an initially circular mutating population getsreinvaded by the yellow bystander population in (a) and is approximately neutral withrespect to yellow population in (b). We see that in (b) the initially circular interfacedissolves. In (a), the dissolution is less prominent, but still has an effect. This differencein roughness properties between Fig. 9(a) and (b) indicates that the overall growth speed v of the interface between the bystander and the invader will influence the interfaceroughness. This is in marked contrast to the d = 1 + 1 -dimensional case where theaverage difference in growth rates between the invader and bystander only changes the nterfaces as indicators of mutational instability d = 2 + 1 , we will haveto consider the dynamics at particular values of the average interface speed v . Figure 10.
A comparison of the invasion front (in a population with a total of L = 100 cells) along the critical line (equal growths for the bystander and invader)at t ≈ generations, starting from an initially uniform, flat interface. In (a), anon-mutating invader is in neutral competition with the bystander ( Γ f = Γ b = 1 and µ = 0 ) and the interface remains overall stationary and dissolves over time. In (b), weare just below the mutational meltdown transition with µ ≈ . (for s = 0 . ) forthe invader. We set b = 0 so that the interface is stationary. We see dissolving of theinterface, but there is an enhanced roughness due to the mutational meltdown. Theroughness may be quantified directly for these simulation snapshots via the width σ w as defined in Eq. (11). Note that width (given in cell diameters) in (b) is twice that in(a). We begin with some qualitative observations of the dynamics. For relativelyneutral populations with an average interface speed v = 0 , we again expect to see anenhanced roughening of the interface over time when the invader population approachesa mutational meltdown, much like in the d = 1 + 1 -dimensional case. We can see theenhanced roughening qualitatively in Fig. 10 for this special case where the invader andbystander have the same growth rate [i.e., we are on the phase transition boundary inFig. 4(b) and v = 0 ]. In Fig. 10(a) we have a non-mutating invader and in Fig. 10(b) wehave an invader near a mutational meltdown transition. The frontier is more undulatedin Fig. 10(b) near the meltdown transition. The increased undulation may be quantifiedby studying the average interface width σ w , which we now describe.In order to partially mitigate the effects of the “fuzzing out” of the interface, wequantify the roughening by looking at the average location of the interface at each time t during the evolution. To do this, we set up a coordinate system where we orient alinear population interface along the x -direction of our lattice and we let x i representthe zigzagged columns of our hexagonal lattice along this direction, as shown in Fig. 11.Then, for each column x i , we define the interface location by averaging over all N u locations of red/black cells within a certain range: y ( x i , t ) = 1 N u y max (cid:88) y = y min y ( x i , t ) , (10) nterfaces as indicators of mutational instability y min = y min ( x i , t ) is the location of a black/red cell on the lattice at the point x i (and time t ) such that all cells with y < y min are also black or red. Similarly, y max isthe location of the red/black cell for which all cells with y > y max are all yellow. Thescheme is illustrated in Fig. 11.Using the average location y ( x i , t ) allows us to define an interface width σ w ( t ) byaveraging over all x i along the interface: σ w ( t ) = (cid:118)(cid:117)(cid:117)(cid:116)(cid:42) L L (cid:88) x i =0 (cid:2) y ( x i , t ) − y ( t ) (cid:3) (cid:43) , (11)where y ( t ) = 1 L L (cid:88) x i =0 y ( x i , t ) . (12)The angular brackets in Eq. (11) indicate an ensemble average over many populationevolutions. However, we may also use σ w ( t ) as an indicator of the front roughness for asingle snapshot of a population at a particular time, as shown in Fig. 10. Examples ofthe calculated σ w ( t ) (averaged over many simulation runs) for various values of selectionparameter b and mutation rate µ are shown in Fig. 12. For example, in the case where theinvader and bystander populations are relatively neutral and there are no mutations, theroughening of the interface illustrated in Fig. 10(a) is shown with blue circles (connectedby a dashed line) in Fig. 12. The interface in Fig. 10(b) approximately corresponds tothe red squares in Fig. 12. Note that, as expected, σ w ( t ) increases significantly faster intime for the latter case compared to the former. Figure 11.
Schematic for finding an average location of the interface between theyellow bystander and red/black mutating populations. The interface runs along the x direction. We identify columns x i in the hexagonal lattice as shown with the bluezigzagged line. At each column x i , the average position y is calculated by averagingover all red/black cell locations between the red/black cell which is the furthest intothe mutating region [at y min ( x i ) ] and the black/red cell which is the furthest into thebystander population [at y max ( x i ) ]. nterfaces as indicators of mutational instability σ w in other ways, includingestimating the interface position using the location y min or y max (see Fig. 11).Alternatively, one might use the difference y max − y min as a measure of the “fuzziness” ofthe interface, which we might also expect to roughen near a mutational meltdown. Wehave verified that using other definitions of the interface roughness does not change thelong-time scaling properties of the interface roughness or the relative enhancement ofthe roughness near a mutational meltdown. It would be interesting, however, to moresystematically study the consequences of using alternative definitions of the roughness.We will now focus our quantitative analysis on the v = 0 case of a stationary(on average) interface, since it is along the critical line where we find a predictableroughening effect. We will then take a closer look at the cases | v | > where eitherthe mutating population or the bystander has an overall selective advantage. Thisintroduces complications as the roughening behavior of a moving front is different froma stationary one. Indeed, whenever | v | > , the invasion becomes a noisy Fisher wavewhich has its own particular roughening properties. We shall see that a non-zero velocity v will suppress the interface roughness at long times, but signatures of the rougheningdue to mutational meltdown persist at shorter times. v = 0 Along the 3-species critical surface, where theinvader and bystander are relatively neutral, we expect to see an enhancement of theinterface roughening as we approach the mutational meltdown transition µ → µ ∗ forthe invader population [the bottom terminal end of the phase boundary in Fig. 4(b)].To quantify the roughening, we can calculate the effective exponent ν ( t ) [see Eq. (9)]from the interface width σ w ( t ) defined in Eq. (11). Without a bias, we expect that theinterface coarsening should be described by voter-model-like dynamics [63] because theinvader and bystander populations divide into each other without a surface tension. Wegenerally expect a diffusive behavior σ w ∝ √ t in this case.In Fig. 13 we see an enhanced roughening as µ → µ ∗ as we move along the phasetransition boundary ( v = 0 ): The limiting value ν of the exponent increases as wemove along the phase transition line towards the mutational meltdown at µ = µ ∗ .Interestingly, near mutational meltdown, the width σ w seems to grow approximatelydiffusively with σ w ∝ t . (red squares in Fig. 13), whereas the non-mutating case µ = 0 coarsens according to the power law σ w ∝ t . (blue circles in Fig. 13). We mighthave expected larger values for these exponents as the non-mutating case should beclosest to the voter model dynamics where interfaces dissolve diffusively, similarly to thedynamics of σ w in the d = 1 + 1 case away from the meltdown transition [63]. However,generalizations of the voter model can yield different results for interface coarsening anddetermining the value of the exponent ν can be subtle [64]. Another possibility is that ν is suppressed due to our particular choice of lattice update rules. It would be interestingto study the behavior with simulations with overlapping generations (independentlydividing cells).Although the behavior for d = 2 + 1 is different from the d = 1 + 1 case where nterfaces as indicators of mutational instability Figure 12.
The interface width σ w [see Eq. (11)] in units of cell diameters of a d = 2 + 1 -dimensional invasion, starting from an initially flat interface between themutating and bystander population ( σ w = 0 ) cells long and with s = 0 . forvarious values of selection parameter b and mutation rate µ . [Note that it is helpful toconsult the phase diagram in the top panel of Fig. 13 for identifying the locations ofthese points in the ( µ, b ) plane.] The interface is evolved for t = 10 generations, andwe average over runs. Lines connect the points to guide the eye. Note that thevalue of b strongly influences the behavior of σ w , as seen by comparing the red squaresand the purple crosses, both of which have the mutating population near meltdown ( µ ≈ µ ∗ ) . In general, we find a suppressed roughness when the mutating and bystanderpopulations are not relatively neutral (compare blue dashed line to orange diamondsand purple crosses). Otherwise, for (on average) stationary interfaces, we see theenhanced roughness due to population meltdown (green triangles and red squares).The smaller plot shows the roughness at short times. the domain wall roughening was clearly super-diffusive near mutational meltdown anddiffusive away from it (see Fig. 7), we also find here that the mutational meltdownenhances the interface undulations by modifying the exponent ν associated with theinterface width σ w ∝ t ν , increasing ν from a value of approximately 0.4 to 0.5. A moredramatic difference is found when we move away from the v = 0 critical line and haveeither the mutating invader population or the bystander grow with an overall selectiveadvantage. We then find a moving Fisher wave with a suppressed exponent ν , as wewill see in the next section. | v | > The comparison between | v | > and v ≈ dynamics can be seen prominently if we consider an initially disc-like population of theinvader strain. Then, any non-zero velocity will either shrink or grow the initial disc.An example of an v < evolution is shown in Fig. 9(a) where the bystander strainreinvades the invader, which eventually dies out. Conversely, when v ≈ , we can seein Fig. 9(b) that the boundary between the invader and bystander gradually dissolves. nterfaces as indicators of mutational instability Figure 13.
Interface roughening exponents ν ( t ) are calculated on the right plots fordifferent combinations of ( b, µ ) indicated on the phase diagrams on the left, for varyingvalues of s = 0 . (top row), . (middle), and . (bottom). As we move along thecritical line (blue circles, green triangles, and red squares), we show the enhancementof the boundary roughness [from σ w ∝ t . to σ w ∝ t . ]. Away from the critical line(purple crosses and orange diamonds), we see the effects of Fisher wave dynamics. Hereeither the mutating population (orange diamonds) or the bystander (purple crosses)has a selective advantage, and the moving interface has a suppressed roughness atlong times, approaching σ w ( t ) ∝ t . (bottom dashed lines in plots on the right),consistent with previous Fisher wave simulation results [56]. The phase diagrams havethe same simulation parameters as in Fig. 5. The exponent curves on the right useinterfaces that are initially cells long, and we average over 160 runs. nterfaces as indicators of mutational instability | v | > different from the critical line: oneof the populations (either the mutating population or the bystander) becomes unstable and will deterministically shrink in the presence of the other population.Let us consider first the simplest case when µ = 0 and we have an interfacebetween a (non-mutating) fast-growing red strain and the yellow bystander. The orangediamond point data in Fig. 13 show what happens in this case. The interface behavesas a noisy Fisher-Kolmogorov-Petrovsky-Piskunov wave [49, 50] describing the invasionof the bystander. Without fluctuations (in the mean field limit), these waves admitstationary shapes and we have no roughening over time. However, fluctuations preventthe formation of stationary wave fronts for the d = 1+1 and d = 2+1 -dimensional cases.For d = 2+1 , previous simulations [56] show that the interface width is expected to growas t ν with ν ≈ . . This coarsening is consistent with our results for σ w , as the orangediamond data points in the right panels of Fig. 13 approach the ν ≈ . limiting valueat long times, indicated by the lower dashed line. The time until convergence, however,is quite long as the effective exponent ν ( t ) continues to decrease over the course of theentire simulation run time.The case of a non-mutating invader is interesting for d = 2 + 1 because we wouldnaively expect our system to fall into the KPZ universality class. The average interfaceposition y ( x t , t ) could be interpreted as a kind of “height function” and the interfacewidth σ w should scale like σ w ∝ t / at early times, consistent with d = 1+1 -dimensionalKPZ dynamics. A broad class of systems fall into this universality class (see [36] fora review) as the KPZ equation includes the most relevant nonlinearity associated withlateral growth. However, we see here that the behavior is more subtle as the fuzzingout of the interface will contribute to the measured roughness. This complication inmeasuring the interface roughness was discussed and analyzed in previous work [54].Our focus here, however, is not the particular exponent associated with the rougheningbut rather the effects of adding mutations. We will see that adding mutations doesenhance the roughness, but only at short/intermediate times while the fast-growing,mutating strain maintains a significant fraction within the population.Consider the portion of the phase diagram where the bystander can reinvade themutating population due to fitness loss at a non-zero mutation rate µ (purple crosses inFig. 13). Here, the evolution of our system begins at first as biased competition betweentwo species (between fast-growing and bystander species) but as the fast-growing cellsmutate and die off, the bystander population begins to reinvade the slow-growing species,and eventually we should find a Fisher wave of the bystander invading the less fit,mutating population. On the right side of Fig. 13 we see that the exponent ν ( t ) forthe purple crosses at first is enhanced as we would expect near mutational meltdown( µ ≈ µ ∗ ). At later times, however, once a Fisher wave is established, the exponenteventually dips down and is consistent with a Fisher-wave-like coarsening [56]. One cantrack this especially easily in the s = 0 . case (top row of Fig. 13) where we see thatthe purple cross data points follow the critical roughening points (red squares) and thentransition to a slower roughening more consistent with a regular Fisher wave (orange nterfaces as indicators of mutational instability σ w ( t ) for this case is also shown in Fig. 12. One sees herethat at times t < (smaller plot), the purple cross data points have a larger width σ w ( t ) due to the mutational meltdown dynamics, but σ w ( t ) then crosses over to smallervalues for longer times when the Fisher wave behavior dominates.
6. Conclusion
We have now analyzed a simple model of invasion of a stable, homogeneous populationby a population acquiring deleterious mutations at a rate µ . We examined this invasionin both one- and two-dimensions as a function of the mutation rate µ , the selectiveadvantage s of the fast-growing strain within the mutating population, and the selectiveadvantage b of the bystander population. We have shown that the effectively small localpopulation sizes (compared to a well-mixed population) suppress the probability thatthe invasion succeeds. This suppression can be understood by analyzing the motionof the boundary between the mutating population and the bystander population it isinvading. We find a reasonable estimate of the phase transition position in the ( µ, b, s ) phase space, as shown in Figs. 4,7, and 13. Our model assumed that cell motility withinour population is suppressed, with the only cell rearrangements occurring due to celldivision and local competition for space. It would be interesting to consider the effectsof a spatial motility as it has been demonstrated that some of the expected features ofspatial dynamics, such as spatial heterogeneity and local fixation of strains is partiallymitigated by increased cell motility [65].Next, we considered the properties of the invasion front and showed that this frontundulates more when the mutating population is near the meltdown transition at whichit loses the fittest strain. For d = 1 + 1 dimensions, this transition is well-characterizedby the directed percolation universality class, and we used properties of this class tounderstand the enhancement of the roughening. In the future, it would be interestingto compare our results to experiments. One possibility is to use microbial populationssuch as bacteria or yeast where one may design strains with varying ( µ, b, s ) . Anotherpossibility would be to examine such invasions in cancers. For instance, it would beinteresting to monitor the edges of a tumor over time as it either grows or shrinks. Wepredict that if the tumor begins losing fitness due to accumulated deleterious mutations(during treatment, for example), then we should be able to observe this transition to“mutational meltdown” as a roughening of the tumor edges.For d = 2 + 1 dimensions, we find a range of behaviors for the roughening interface.When the speed of the invasion front approaches zero, the interface roughens moresignificantly due voter-model-like coarsening. We also find an enhancement of theroughening exponent as the mutating population approaches meltdown. On the otherhand, when either the mutating or the bystander population has a selective advantageand the population interface develops an overall velocity, the roughening is suppressed,and we find roughening exponents consistent with those observed for noisy Fisher wavesat long times. Therefore, the long-time behavior of the population interface roughness nterfaces as indicators of mutational instability σ w ∝ t ν . For long times t , the interface undulation size will eventuallysaturate due to the finite system size L , and we might expect a general scaling form σ w = t ν f ( t/L β ) , with f ( x ) a scaling function and β a new critical exponent. The scalingproperties of this saturation should also depend on the proximity to the mutationalmeltdown transition. It would also be interesting to consider a d = 3 + 1 -dimensionalevolution such as the invasion of surrounding tissue by a compact cluster of cancerouscells. In this case, the invasion front would be an entire surface which could also pinchoff and coarsen. Previous simulations of the noisy Fisher wave dynamics suggest thatthe situation in this case is similar to the d = 2 + 1 case considered here [56]. We wouldagain expect to find some enhancement of the interface roughening when a mutatinginvader is near a mutational meltdown transition.Interestingly, increased roughening is typically an indicator of more malignantcancerous growths, and the roughness of tumor edges has been a useful prognosticindicator in a wide variety of cancers [66]. Also, in general, increased heterogeneityresults in a worse clinical prognosis [67]. While our results point to the possibility of anopposite correlation, our model does not take into account tumor vasculature or cancercell motility. Conversely, most of the clinical studies focus on more mature tumors whichhave developed a vasculature. Hence, we expect our model to be relevant for early, smallavascular tumors or regions of larger tumors lacking vasculature. These small tumorsare not easily detected as they are typically just a few millimeters in size. Nevertheless,small spheroidal avascular tumors are good in vitro models for early cancer growth [68].It would thus be interesting to study the edges of such cultured tumors under a largemutational load. We may also verify some of our results in microbial range expansions(e.g., in yeast cell colonies grown on Petri dishes) where there is little cell motility.A promising experimental realization of a d = 2 + 1 -dimensional expansion may be agrowing cylindrical “pillar” of yeast cells, as realized in Ref. [69]. Acknowledgments
M.O.L. thanks B. Weinstein for helpful discussions. Computational support wasprovided by the University of Tennessee and Oak Ridge National Laboratory’s NationalInstitute for Computational Sciences. M.O.L. acknowledges partial funding from theNeutron Sciences Directorate (Oak Ridge National Laboratory), sponsored by the USDepartment of Energy, Office of Basic Energy Sciences. nterfaces as indicators of mutational instability References [1] Petren K and Case T J 1996 An experimental demonstration of exploitation competition in anongoing invasion
Ecology Proc. Natl. Acad.Sci. U.S.A
Nat. Rev. Microbiol. Virology
Trends Cell Biol. (10) 776–788[6] Willi Y, Fracassetti M, Zoller S and Buskirk J V 2018 Accumulation of mutational load at theedges of a species range Mol. Biol. Evol. Genetics bioRxiv [9] Cahill D P, Kinzler K W, Vogelstein B and Lengauer C 1999 Genetic instability and darwinianselection in tumours
Trends Cell Biol. M57–M60[10] Negrini S, Gorgoulis V G and Halazonetis T D 2010 Genomic instability - an evolving hallmark ofcancer
Nat. Rev. Mol. Cell Biol. Biochim. Biophys.Acta Rev. Cancer
Breast Cancer Res. R48[13] Yoo J, Chong S, Lim C, Heo M and Hwang I G 2019 Assessment of spatial tumor heterogeneityusing CT growth patterns estimated by tumor tracking on 3D CT volumetry of multiplepulmonary metastatic nodules
PLoS ONE e0220550[14] Sottoriva A, Kang H, Ma Z, Graham T A, Salomon M, Zhao J, Marjoram P, Siegmund K, PressM, Shibata D and Curtis C 2015 A Big Bang model of human colorectal tumor growth Nat.Genet. Bull. Math. Biol. Nat. Rev.Cancer Proc. Natl. Acad. Sci. U.S.A.
Cancer Res. Proc. Natl. Acad. Sci. U.S.A.
Mutat. Res. Evolution Nature nterfaces as indicators of mutational instability [23] Zhang Y, Li Y, Li T, Shen X, Zhu T, Tao Y, Li X, Wang D, Ma Q, Hu Z, Liu J, Ruan J, CaiJ, Wang H Y and Lu X 2019 Genetic load and potential mutational meltdown in cancer cellpopulations Mol. Biol. Evol. Radiology
Breast Cancer Res. Treat. New J. Phys. Proc. Natl. Acad. Sci. U.S.A.
Evolution Biophys. J.
Nat.Commun. Adv. Phys. Biophys. J. Phys. Rev. Lett. Phys. Rev. E Phys. Rev. E Phys. Rep.
PLOS Comput. Biol. e1007243[38] Nguyen B, Upadhyaya A, van Oudenaarden A and Brenner M P 2004 Elastic instability in growingyeast colonies Biphys. J. New J. Phys. Phys. Rev. E Phys. Rev. Lett.
PLoS Genet. (9) e1007450[43] Lavrentovich M O, Korolev K S and Nelson D R 2013 Radial Domany-Kinzel models with mutationand selection Phys. Rev. E Phys. Rev. Lett. (4) 311–314[45] Giese A, Bjerkvig R, Berens M E and Westphal M 2003 Cost of migration: invasion of malignantgliomas and implications for treatment J. Clin. Oncol. J. Stat. Mech. Theory Exp.
P05027 nterfaces as indicators of mutational instability [47] Täuber U C, Howard M J and Hinrichsen H 1998 Multicritical behavior in coupled directedpercolation processes Phys. Rev. Lett. Phys.Biol. Ann. Eugen. (4) 355–369[50] Kolmogorov A, Petrovsky N and Piscounov N 1937 Étude de l’équation de la diffusion aveccroissance de la quantité de matière et son application a un problème biologique Moscow Univ.Math. Bull. Phys. Rev. Lett.
J. Stat. Phys. Phys. Rev. Lett. Phys. Rev. Lett. Phys. Rev. Lett. Phys.Rev. Lett. Ann. Prob. Math. Biosci. Phys. Rep.
Biophys. J. PLOS Comput. Biol. e1005866[62] Korolev K S, Xavier J B, Nelson D R and Foster K R 2011 A quantitative test of populationgenetics using spatiogenetic patterns in bacterial colonies Am. Nat.
Phys. Rev. Lett. Phys. Rev. E Nature
Insights Imaging Oncotarget Neoplasia Biophys. J.106