Adaptive dynamics of unstable cancer populations: The canonical equation

Abstract In most instances of tumour development, genetic instability plays a role in allowing cancer cell populations to respond to selection barriers, such as physical constraints or immune responses, and rapidly adapt to an always changing environment. Modelling instability is a nontrivial task, since by definition evolving instability leads to changes in the underlying landscape. In this article, we explore mathematically a simple version of unstable tumour progression using the formalism of adaptive dynamics (AD) where selection and mutation are explicitly coupled. Using a set of basic fitness landscapes, the so‐called canonical equation for the evolution of genetic instability on a minimal scenario associated with a population of unstable cells is derived. We obtain explicit expressions for the evolution of mutation probabilities, and the implications of the model on further experimental studies and potential mutagenic therapies are discussed.


| INTRODUC TI ON
Cancer can be understood as the failure of those regulatory mechanisms that guarantee the maintenance of tissue and organ homoeostasis. Cooperative interactions along with extensive feedback signalling loops and replication checkpoints provide multiple paths to avoid the emergence of undesirable mutations or chromosomal abnormalities that can allow rogue cells to start proliferative growth.
In dynamical terms, what has to be avoided within multicellular organisms is any kind of individual cell Darwinian evolution (Gatenby & Brown, 2017;Greaves & Maley, 2012;Nowell, 1976).
It is generally acknowledged that genetic instability plays a key role in tumour progression and carcinogenesis (Hanahan & Weinberg, 2011). Unstable genomes result from the failure of molecular mechanisms responsible for the maintenance of genome integrity (Negrini, Gorgoulis, & Halazonetis, 2010). That cancer cells are unstable is fairly well illustrated by the observation of their karyotypes: in sharp contrast with healthy cells, cancer chromosomal arrangements reveal a wide degree of aneuploidy (Lengauer, Kinzler, & Vogelstein, 1998). Such high levels of mutational load deploy the potential to overcome selection barriers, as well as involve a rather uncommon process from multicellularity to reduced cellular complexity (Solé et al., 2014), giving place to a highly adaptive and heterogeneous population. Genetic instability acts as a driver as well as the search engine for disease progression. An important (and not always appreciated) consequence of instability is that, once unleashed, it can easily grow as the lack of proper repair can damage other components of the check-and-repair cellular network.
Despite increasing knowledge of the molecular basis of unstable tumorigenesis, there is still the need for understanding the role of instability on cancer evolution, namely discerning if it is a cause or a consequence of carcinogenesis, how does it evolve along tumour development, and what are the treatment strategies that arise from answering such questions. Many mathematical models have provided interesting points into this topic, with the introduction of relevant ideas such as the mutator phenotype (Loeb, 2001) and several multi-step models of mutation acquisition (see e.g., Komarova et al., 2002;Nowak et al., 2002) that have investigated the possible scenarios of correlation between instability and cancer progression.
The fact that genetic instability itself changes over cancer evolution makes it difficult to properly model its behaviour. Particular efforts, such as the computational models of Komarova, Sadovsky, &Wan, 2008 andDatta et al., 2013; have given interesting insight into understanding how a changing instability level affects, by means of modifying the probability of mutations, all kinds of replication and control mechanisms within the vast pathways towards cancer malignancy. Within this picture, instability cannot be taken as a parameter, but rather as an evolving phenotypic trait affected by the selective pressures of the tumour microenvironment. In this scope, the recent work by Asatryan and Komarova represents a further step for its proposal of an analytical approach where both instability and heterogeneity of cancer populations can be traced along time (Asatryan & Komarova, 2016). As a complementary point of view, we consider the need to include stochasticity in the process of acquiring either advantageous or deleterious mutations, together with considering instability as a trait evolving through changes within each single cell, compared to the idea of measuring it by following the competition dynamics between subpopulations with fixed mutation probabilities.
Here, we propose that the mean evolutionary paths of such stochastic process followed by unstable populations are describable by means of the framework of adaptive dynamics (AD) (Champagnat, Ferriere, & Ben Arous, 2001;Dieckmann & Law, 1996), which has been used in the study of cancer when focusing on niche construction (Gerlee & Anderson, 2015). AD models provide a powerful alternative to previous formal approaches by explicitly including replication, mutation and selection in a consistent way, allowing an exploration of the evolutionary dynamics of adaptive traits, while at the same time keeping a minimal, treatable model able to produce explicit expressions for trait evolution depending on a few parameters.
A central object in the AD framework is the so-called canonical equation. For a given quantitative phenotypic trait s, this equation describes the evolutionary trajectory for the mean trait value as where μ(⟨s⟩) is the probability under which mutant individuals are generated, σ 2 is the variance of the mutant distribution s′ derived from an individual with trait s, n the stationary population size and the last term in the right-hand side stands for the fitness gradient associated with the specific landscape at work. The standard formulation involves some strong assumptions on the mutation-selection process, and we will therefore review the mathematical process in order to understand up to which point the framework is suitable for our problem.
In the AD models, and in the work presented here, evolution takes place within a constant population context, where mutants appear and invade in a stepwise process, leading to a formalism for evaluating the trajectories of evolving trait values. This picture of cancer dynamics stems from classic work on ecological competition (Gatenby, 1995) where tumour cells act as invaders that cause the disruption of the local (tissue) ecology. These simplified models reveal how a proper formulation of competition can yield useful predictions (Gatenby, 1991). In this context, although the constant population falls short to describe the behaviour of some tumour growth processes, it is a much needed first approximation. Moreover, it can also be appropriate when dealing with some in vitro experiments involving long-term evolution of unstable cancer cell populations. We will go back to this at the end of the paper.
Understanding how instability becomes a driver of evolvability can give further insight about its role as a cancer hallmark, and might as well produce relevant steps towards contemplating genetic instability as potential target for treatment. Is it possible to formulate a canonical equation describing the time evolution of instability? The answer is affirmative and here we show how it can be obtained.

| P OPUL ATION DYNAMIC S
With the aim of obtaining a clear understanding of the questions proposed above, we look for a minimal model to implement the unstable evolutionary dynamics. Our goal is to consider the process of cancer progression, which involves a heterogeneous population of cells ( Figure 1a). In this population, cells are only characterized by their particular replication rate r i and mutation probability μ i . However, and in the eyes of the AD approach, this preliminary model uses a constant population approximation where mutation probabilities are small enough so that the dynamics remain in equilibrium in-between invasions. This approach, whose limitations will be later thoroughly discussed, is best described by means of a so-called Moran process (Moran, 1958).
A particularity of the Moran process-here coarse-grained into a continuous process, keeping in mind the long-term evolution of tumour progression-is that cells of type c i give birth by means of occupying other, randomly chosen cell sites at rate r i , so that the birth-death process is coupled into a single event (Figure 1c) that will eventually lead to selection towards cells with higher r i , thus producing a minimal environment where selection can take place.

Furthermore, mutation is introduced by considering that cells can
give birth to mutant offspring at probability μ i .
Mutations, however, do not occur as in quasispecies or replicatormutator models, where genomes mutate from one to another. In our model, a newly born mutant cell will have a modified mutation probability μ � = μ i + Δμ, where Δμ is taken from a continuous distribution

| S ELEC TI ON ON IN S TAB ILIT Y
As discussed above, a most common event during the process of tumorigenesis are mutations in oncogenes that usually result in increased levels of replication (Vogelstein & Kinzler, 2002), thus giving to instability a role in activating the paths towards higher replicative capacity. On the other hand, the same elevated levels of instability can trigger deleterious mutations in genes that are vital for correct cellular metabolism and functioning, eventually leading to reduced cell viability or death. This apparent trade-off supposes the existence of a clear coupling between replicative capacity, cell viability and mutation probability that sits at the basis of tumour replication, evolvability and adaptation. We hereby propose a minimal adaptive landscape that translates such coupling into replication rates being a function of instability, r(μ).

| Adaptive landscape
Within our scope of producing a minimal model we expect to describe evolutionary dynamics on an adaptive landscape containing a reduced, treatable set of components. Taking into account the previously mentioned trade-off, these follow from considering the effect of mutations enhancing malignant cell replication, provided that such mutations have not damaged any of the necessary machinery for cell viability. We start by considering that mutations on oncogenes can translate into a linear increase in replication rate, such that r(μ) = r 0 + N R δ R μ , with r 0 being the basal replication rate of normal cells, N R the number of oncogenes responsible for increased replication and δ R the mean effect on replication rate when mutating one of such genes. Following a linear approximation, we do not include a saturation term for the number of nonmutated oncogenes. This is actually consistent with early stages of tumour evolution, where only a small fraction of oncogenes has been affected and so N R can be kept as a constant. In this picture, we need as well to take into account the minimal genetic material needed for a cell to keep its basic functions. If we group such material into a number of house-keeping genes, N HK , the probability that none of them has been mutated is (1 −μ) N HK . Grouping both considerations together we obtain an analytical description of the coupling between replication and instability (Solé et al., 2014). This adaptive landscape is of course of qualitative nature, and realistic fitness landscapes for unstable tumour environments are still far from our knowledge. However, certain points can be made if we give values within realistic parameter ranges to our function. The number of both oncogenes and house-keeping genes have been widely assessed, and we take them to be about N R ≈ 140 (Vogelstein et al., 2013) and N HK ≈ 3804 (Eisenberg & Levanon, 2013), respectively. Interestingly enough, considering small replication effects for δ R , such experimental values produce an adaptive landscape (Figure 2) that has a positive gradient within the region of The Moran process rules associated with the model of a population of unstable cells competing for resources. We consider an idealized model of a heterogeneous cancer cell population (a) described by a well-mixed (mean field) model (b). Here cells occupy a given domain that is not explicit and each cell has a distinct phenotype described in particular by its intrinsic instability μ i . In the Moran process, when a cell replicates it occupies another cell's niche and produces an identical daughter (c) or a slightly different one due to a mutation event proportional to μ i , which can lead to an increase Δμ i of the instability levels (d) , so that our evolutionary trajectories will be bounded within a region of instability levels in accordance with those experimentally measured for tumour cells (Tomlinson, Novelli, & Bodmer, 1996).

| Distribution of new mutations
We have assessed so far what is the effect of instability in proliferation, thus coupling mutation and selection for mutation level.
Up next, we need to evaluate how does instability change during reproduction, so that we can finally compute the effects on replicative capacity of a mutated cell. As previously discussed, a broad range of mechanisms relates to variations in DNA replication fidelity.
Such variations, however, are hardly in the direction of increasing DNA stability, and in general account for an increase in the mutation probability of cancer cells due to accumulation of further tumoursuppressor or care-taker gene mutations (Vogelstein & Kinzler, 2002).
This trend of generating more unstable offspring is translated into a positively skewed distribution of mutants M(μ,Δμ). To keep the mathematical background of our model treatable, a Rayleigh distribution peaked at Δμ = 0 has been chosen 1 . Under this scheme, instability of a daughter cell is likely to be similar or slightly higher from its parent, controlled by a scale parameter σ 2 μ depicting the general size of mutational increases.

| ADAP TIVE DYNAMIC S
Adaptive dynamics is a set of techniques or a mathematical framework that models long-term phenotypic evolution of populations.
Several works by different authors cover a broad scope of possible applications, and we hereby focus on the work of Dieckmann and Law and others (Champagnat et al., 2001;Dieckmann & Law, 1996) and adapt it to our particular system. The main biological background behind the maths sits in considering the evolutionary step as a mutant appearing and invading in a population in ecological or dynamical equilibrium (Dieckmann & Law, 1996). Under this picture, the ecological and evolutionary time scales are considered to be uncoupled, so that the process of the mutant competing against the resident population, and eventually fixating in it, is considered instantaneous in the evolutionary process.
General AD literature (see e.g., Champagnat et al., 2001;Dieckmann & Law, 1996;Geritz et al., 1998)  The Rayleigh distribution is an asymmetric probability distribution defined for positive random variables (Forbes et al., 2010). We displaced so that its mode is zero, with shape accounting for asymmetric probabilities of possible mutation levels, with small forward mutations being the most common.

M(Δμ;σ
F I G U R E 2 Fitness landscape function associated with the evolutionary dynamics of unstable tumour cells. In (a) the full landscape, given by a replication rate r(μ) = (r 0 + N R δ R μ)(1 −μ) N HK , is plotted against the instability probability μ. Further discussion is focused on two limit cases representing initial linear progression of instability (b) and optimal mutation (c) domains The transition probabilities w(μ � |μ) describe the evolutionary step and contain the probability of the mutant with trait μ′ appearing (A) and fixating (ρ) in the population, so that w(μ � |μ) = A(μ,μ � )ρ(μ,μ � ). The probability that a mutant appears is A(μ,μ � ) = Nr(μ)μM(μ,μ � ) , the size of the population at equilibrium N, the probability of birth and mutation r(μ)μ and the probability that the mutant has mutation probability μ′ provided the parent cell had probability μ.
The probability ρ(μ,μ � ) that a mutant with fitness advantage r(μ � )∕r(μ) fixates in a population of N individuals has an analytical expression for the Moran model (Ewens, 2004) A common procedure of the AD framework is to expand ρ for small variations of the trait value under the assumption of large populations, assumptions that are not a restriction for our problem.
Under this view, the probability that the μ′ mutant fixates is zero for r(μ � ) ≤ r(μ) and for r(μ � ) > r(μ). Once a complete expression for the transition probabilities is build, we only need to recall how the evolution of the mean mutation probability can be written as so that, using the original master equation, we obtain where one recalls that a k (μ) = ∫ (μ � −μ) k w(μ � |μ)dμ � is the k-th jump moment. If the first jump moment were a linear function of μ, then ⟨a 1 (μ)⟩ = a 1 (⟨μ⟩) and the previous expression becomes directly treatable. 2 The evolutionary trajectory for the mean path (we cease denoting it by angle brackets) will therefore follow At this point, we recall that only fitter mutants can invade and so may eventually contribute to the exploration of the adaptive landscape. This translates onto the domain of the integration being restricted to μ � > μ , and so we integrate the positive part of our σ 2 μ skewed mutant distribution. These

| E VOLUTI ON OF IN S TAB ILIT Y
The canonical Equation 9 describes the evolution of instability in our model population depending on the population size N, the distribution of mutation jumps σ 2 μ and the product of the mutation probability and the gradient of the adaptive landscape μ ∂ μ r . For our model landscape (Equation 5), this turns out to be

Complicated analytical solutions for this equation might not
give best insight of the underlying dynamics. However, as a first test of our model we compare its numerical solution to averaged

Moran Process simulations (Figure 3). It is both relevant and useful
to understand the factors that cause deviations between computer experiments and our analytical approach, in order to further comprehend the approximations on which AD is build.
In terms of parameter range, these are mostly translated into the population being large enough, and mutation probabilities being proportionally small. The second is easily met for both healthy and cancerous human cells, but simulating full-size clinically detectable tumours (more than 10 8 cells (Bozic et al., 2013)) is of large computational cost, and keeping our model and It is interesting to understand the implications of such condition and how do they relate to the assumptions of the AD method and the limitations of our model. Between many approaches (see e.g., van Kampen, 1981;Kubo, Matsuo, & Kitahara, 1973), and without pretending to expose here a deep discussion on this aspect, one might consider extracting a Fokker-Planck equation from the Markov chain (6). By means of computing the evolution of the mean value ⟨μ⟩ for that equation, it can be seen that the dynamics of ⟨μ⟩ are only equivalent in both frameworks if either the selection gradient ∂r∕ ∂μ does not depend on μ (a linear adaptive landscape), or the population is strictly monomorphic on the trait μ. A detailed explanation of this considerations can be found in the Appendix section.
These conditions will impose a strong constraint when considering a realistic tumour environment, so that it will become necessary along our computational experiments and the discussion to assess the regions of validity of our model.
exercise minimal, we have used smaller populations, modelling smaller subclones or spatially segregated populations where drift comes into play. Such drift produces a nonmonomorphic population where evolution deviates from the gradient trajectory and so proceeds slightly slower than our estimate. As previously stated and discussed along with the Appendix, the high nonlinearity of our landscape ensures that ⟨a 1 (μ)⟩ = a 1 (⟨μ⟩) will be only valid up to a certain degree of approximation. It can be seen from Figure 3 that, still within this restricted range of validity, the canonical equation can capture the dynamics of instability up to a reasonable point.
A better understanding of the underlying dynamics can result from dividing the exploration of the landscape in well-behaved regions where simpler equations will arise.
On the one hand, in an initial phase of malignancy exploration for small values of μ, the shape of the adaptive landscape is dominated by the linear increase of mutated oncogenes, r(μ) = r 0 + N R δ R μ .

Within this region, dynamics of instability follow and the mean evolutionary trajectory is
It is remarkable to understand how, even in a linear adaptive landscape, the coupling between mutation and selection on unstable cells introduces a further nonlinearity that will account for exponential exploration of the space of instability and the consequent exponential increases in replication capacity.
Such results can be again compared to computer simulations of mutator-replicator cells (Figure 4). The smaller nonlinearity also ensures that AD remains a good approximation despite stochastic deviation.
Another interesting point is to understand the behaviour of the mean instability levels as the population approaches the landscape peak. This kind of behaviour is easily studied if one considers a simple landscape containing a peak, such as r(μ) = r 0 +δ R N R μ−δ HK N HK μ 2 , where the role of house-keeping genes is not considered totally deleterious but just reducing fitness quadratically with the mutation probability. This landscape has an optimal value at μ * = δ R N R ∕2δ HK N HK , and this peak is explored through By means of rewriting this trajectory as dμ∕dt = Aμ(B − Cμ) , with A = γNσ 2 μ , B = N R δ R and C = 2δ HK N HK , its solution simplifies to where c 1 ensures that μ(0)=μ 0 , the normal mutational probability of healthy cells. This trajectory saturates for long times at the expected result A/B=μ*, and can be again compared to computational experiments of replicating cells ( Figure 5).
The same deviation between simulations and the numerical fit is found in this case, with evolution proceeding slower than our estimate.
However, this minimal landscape approximation is able to capture the dynamical behaviour of our gene-related landscape model, mainly with an initial exponential growth followed by saturation around the peak, which can be proven to be an evolutionary stable strategy (Geritz et al., 1998).
Provided that the canonical equation has a nontrivial, singular point, as we found for μ * = δ R N R ∕2δ HK N HK , one can study the evolutionary stability of a quantitative trait. We can easily compute if this (14) dμ dt = γNσ 2 μ μ(N R δ R − 2δ HK N HK μ). F I G U R E 4 Exponential evolution of the mean mutation probability on a linear landscape: Moran process simulations (grey lines) of populations of 2,000 cells and the AD approximation (black curve), (σ = 0.01, μ 0 = 5 × 10 −6 ) singular mutation probability will be an evolutionary trap, that is, a strategy that no further mutants can invade, if which holds for our strategy: ∂ μμ r(μ) = −2δ HK N HK < 0 .

| D ISCUSS I ON
In the present work, we have discussed the implications of the Mounting evidence indicates that a successful approach to cancer therapy requires an explicit evolutionary perspective (Gatenby et al., 2009). One possible instance of this is provided by mutagenic therapies that have produced key results in the field of virology (Loeb et al., 1999). Would they be effective for cancer? Given some key analogies between RNA virus populations and unstable tumours (Solé & Deisboeck, 2004), this is an appealing possibility, although drug design or resistance mechanisms have yet to be assessed (Fox & Loeb, 2010 Regarding the later, our model allows to bring instability as the evolving trait, while providing potential insights, particularly before and beyond the optimal instability levels. The exponentially fast increase of small mutational loads indicates that reducing instability levels in hope for progression delay might result in rapid re-exploration of the mutator phenotype. On the other hand, pushing instability beyond optimal levels, even if a critical point is not trespassed (Solé & Deisboeck, 2004), might render tumour cells too unstable, and there exist relevant efforts towards using DNA repair inhibitors to produce critical instability levels (Helleday et al., 2008).
Our model differs from previous work in its simple analytical formulation, which do not depend on chosen parameter ranges, such as those of Datta et al., 2013 andAsatryan &Komarova, 2016; and they are thus qualitatively robust. On the one hand, this means that we are able to obtain analytical expressions for the exponential evolution and saturation of the mutation probability, which could eventually be used when studying in vitro long-term evolutionary experiments with cancer populations, using serial transfer methods similar to those performed on viruses (see e.g., Drake, 1993;Sanjuán et al., 2010;Solé et al., 1999) or bacterial populations (see e.g., Moxon et al., 1994;Sniegowski, Gerrish, & Lenski, 1997 F I G U R E 5 Mean mutation probability saturation at the fitness peak: Moran process simulations (grey lines) of populations of 2,000 cells and the AD approximation (black curve), (σ = 0.01, μ 0 = 10 −5 ) evolutionary time scales (Lambert et al., 2011), it would be worth exploring the evolution of instability of cancer cell cultures over many transfer generations (Langdon, 2004).
It remains an open question to analyse if any sort of saturating dynamics occur for both the fitness or the mutation probability when these experiments are performed on malignant cells. Furthermore, interesting theoretical approaches have been performed to infer the underlying adaptive landscape from the observable evolution of traits (Kryazhimskiy, Tkacik, & Plotkin, 2009). This seems a plausible point regarding how our model directly relates dynamics and landscape gradient and could therefore shed light onto understanding the evolutionary pressures underlying genetic instability at each stage of tumour progression.
On the other hand, while trying to produce a model that can be treated without the use of complex mathematical tools, we have been constrained to leaving aside many relevant considerations, the one we are most concerned with is the lack of growing population dynamics. This leaves our model interesting for the previously discussed specific confined experiments, while not yet complete when trying to study three-dimensional growing tumours. While studying modifications to our formalism, following the work of evolutionary game theory on growing populations (see e.g., Li et al., 2015;Melbinger, Cremer, & Frey, 2010), we have decided to present this basic model as it remains a first step into a comprehensible and qualitative insight for the dynamics of populations able to evolve their mutation probability.

ACK N OWLED G EM ENTS
The authors thank P. Ruiz and E. Beltran, as well as the CSL laboratory members for useful discussions. Special thanks to Oliver Law for his inspiring ideas. This work has been supported by the Botín

CO N FLI C T O F I NTE R E S T
None declared.

APPENDIX LI N E A R IT Y O F 〈a 1 ( μ )〉
As discussed when introducing the AD approach, we have considered that the first jump moment is linear in μ to obtain the canonical equation. Understanding the mathematical basis behind such condition can give further insight into the assumptions our model sits on.
Let us start again from a Master equation that describes the evolutionary biased random walk: the trait substitution sequence that we use to model evolution. This is We knew that w(μ, μ ′ ), the transition probability, was the probability of a mutant appearing times the probability of a mutant surviving. (17) dP(μ,t) dt = ∫ w(μ|μ � )P(μ � ,t) − w(μ � |μ)P(μ,t) dμ �