Abstract

Dynamics and bifurcations in a simple quasispecies model of tumorigenesis

Vanessa Castillo, J. Tomás Lázaro Josep Sardanyés,

1. Universitat Politècnica de Catalunya, Barcelona, Spain

2. ICREA-Complex Systems Lab, Department of Experimental and Health Sciences, Universitat Pompeu Fabra, Dr. Aiguader 88, 08003 Barcelona, Spain

3. Institut de Biologia Evolutiva (CSIC-Universitat Pompeu Fabra), Passeig Maritim de la Barceloneta 37, 08003 Barcelona, Spain

Corresponding authors. E-mail: jose.tomas.lazaro@upc.edu; josep.sardanes@upf.edu.

## Abstract

Cancer is a complex disease and thus is complicated to model. However, simple models that describe the main processes involved in tumoral dynamics, e.g., competition and mutation, can give us clues about cancer behaviour, at least qualitatively, also allowing us to make predictions. Here we analyze a simplified quasispecies mathematical model given by differential equations describing the time behaviour of tumor cells populations with different levels of genomic instability. We find the equilibrium points, also characterizing their stability and bifurcations focusing on replication and mutation rates. We identify a transcritical bifurcation at increasing mutation rates of the tumor cells population. Such a bifurcation involves an scenario with dominance of healthy cells and impairment of tumor populations. Finally, we characterize the transient times for this scenario, showing that a slight increase beyond the critical mutation rate may be enough to have a fast response towards the desired state (i.e., low tumor populations) during directed mutagenic therapies.

## 1 Introduction

Cancer progression is commonly viewed as a cellular microevolutionary process [1, 2]. Genomic instability, which seems to be a common trait in most types of cancer [3], is a key factor responsible for tumor progression since allows a Darwinian exploratory process required to overcome selection barriers. By displaying either high levels of mutation or chromosomal aberrations, cancer cells can generate a progeny of highly diverse phenotypes able to evade such barriers [4]. Genomic instability refers to an increased tendency of alterations in the genome during the life cycle of cells. Normal cells display a very low rate of mutation ( changes per nucleotide and replication cycle). Hence, it has been proposed that the spontaneous mutation rate in normal cells is not sufficient to account for the large number of mutations found in human cancers. Indeed, studies of mutation frequencies in microbial populations, in both experimentally induced stress and clinical cases, reveal that mutations that inactivate mismatch repair genes result in times the background mutation rate [5, 6, 7]. Also, unstable tumors exhibiting the so-called mutator phenotype [4] have rates that are at least two orders of magnitude higher than in normal cells [8, 9]. This difference leads to cumulative mutations and increased levels of genetic change associated to further failures in genome maintenance mechanisms [10]. The amount of instability is, however, limited by too high levels of instability, which have been suggested to exist in tumor progression [3], thus indicating that thresholds for instability must exist. In fact, many anti-cancer therapies take indirectly advantage of increased genomic instability, as is the case of mitotic spindle alteration by taxol or DNA damage by radiation or by alkilating agents.

The mutator phenotype is the result of mutations in genes which are responsible of preserving genomic stability, e.g., BRCA1 (breast cancer 1), BLM (bloom syndrome protein), ATM (ataxia telangiectasia mutated) or gene protein P53 (which is involved in multitude of cellular pathways involved in cell cycle control, DNA repair, among others). This mutator phenotype undergoes increases in mutation rates and can accelerate genetic evolution in cancer cells that can ultimately drive to tumor progression [4]. As mentioned, genomic instability is a major driving force in tumorigenesis. Tumorigenesis can be viewed as a process of cellular microevolution in which individual preneoplastic or tumor cells acquire mutations that can increase proliferative capacity and thus confer a selective advantage in terms of growth speed. The rate of replication of tumor cells can increase due to mutations in both tumor suppressor genes, e.g., APC (adenomatous polyposis coli) or P53; and oncogenes, e.g., RAS (rat sarcoma) or SRC. Tumor suppressor genes protect cells from one step on the path to cancer and oncogenes are genes that, when mutated, have the potential to cause cancer. In terms of population dynamics, alterations in both types of genes drive to neoplastic process through increases in cancer cells numbers. Mutations in replication-related genes that confer an increase of fitness and thus a selective advantage are named driver mutations [11]. This evolutionary process allows tumor cells to escape the restrictions that limit growth of normal cells, such as the constraints imposed by the immune system, adverse metabolic conditions or cell cycle checkpoints.

The iterative process of mutation and selection underlying tumor growth and evolution promotes the generation of a diverse pool of tumor cells carrying different mutations and chromosomal abnormalities. In this sense, it has been suggested that the high mutational capacity of tumor cells, together with an increase of proliferation rates, may generate a highly diverse population of tumor cells similar to a quasispecies [12, 13]. A quasispecies is a “cloud” of genetically related genomes around the so-called master sequence, which is at the mutation-selection equilibrium [14, 15]. Due to the heterogeneous population structure, selection does not act on a single mutant but on the quasispecies as a whole. The most prominent examples of a quasispecies are given by RNA viruses (e.g. Hepatitis C virus [16], vesicular stomatitis virus [17], the human immunodeficiency virus type 1 [18], among others).

An important concept in quasispecies theory is the so-called error threshold [14, 15]. The error threshold is a phenomenon that involves the loss of information at high mutation rates. According to Eigen’s original formulation, a quasispecies can remain at equilibrium despite high mutation rates, but the surpass of the critical mutation rate will upset this balance since the master sequence itself disappears and its genetic information is lost due to the accumulation of errors. It has been suggested that many RNA viruses replicate near the error threshold [19]. Another important concept in quasispecies theory is lethal mutagenesis. As a difference from the error threshold (which is a shift in sequence space), lethal mutagenesis is a process of demographic extinction due to an unbearable number of mutations [20, 21]. Most basically, it requires that deleterious mutations are happening often enough that the population cannot maintain itself, but it is otherwise no different from any other extinction process in which fitness is not great enough for one generation of individuals to fully replace themselves in the next generation. In simple words, increased mutagenesis could impair the maintainance of a quasispecies due to the crossing of the error threshold or due to lethal mutagenesis.

Quasispecies theory has provided a population-based framework for understanding RNA viral evolution [22]. These viruses replicate at extremely high mutation rates and exhibit significant genetic diversity. This diversity allows a viral population to rapidly adapt to dynamic environments and evolve resistance to vaccines and antiviral drugs. As we previously mentioned, several features have been suggested to be shared between RNA viruses and tumors [23], at least qualitatively. One is the presence of high levels of heterogeneity, both at genotype and phenotype levels. Typically, cancer cells suffer mutations affecting cell communication, growth and apoptosis (i.e., programmed cell death). Accordingly, escape from the immune system (and other selection barriers) operates in both RNA viruses and tumors. Viruses use antigenic diversity whereas tumors evade the immune system by loosing their antigens through mutation, or making use of antigenic modulation and/or tumor-induced immune suppression. Even more, similarly to RNA viruses, genetic instability in cancer cells will have detrimental effects on cells fitness, since most random mutations are likely to be harmful. As indicated by Cahill et al. [3], the best chance of cure advanced cancers might be a result of tumor genetic instability: cancer cells are more sensitive to stress-inducing agents. In this sense, possible therapies increasing mutation of tumor cells could push this populations towards the error threshold or induce lethal mutagenesis. This is the topic that we will address in this work by using a mathematical model describing the dynamics of competition between different cell populations with different levels of genomic instability. Specifically, we will find equilibrium points, characterizing their stability, focusing on possible changes in the stability as a function of the mutation rate and the fitness of cells populations. We will also investigate the transient times of the system to reach a given equilibrium state.

Each of the mathematical results will be given in terms of the mutation rates or the replication fidelity of the cells. This is important to understand how tumor populations behave. As we previously mentioned, unstable cells may reach an error threshold and start loosing genetic information until its population decreases or, even more, disappears. We will see how mutation rates affect this error threshold. Furthermore, we will see how they affect to the velocity of the tumor population to reach an equilibrium point, which is interesting since mutation rates can be modified through drugs (similarly to RNA viruses [24]) or radiotherapy.

## 2 The model

In this section we introduce the model by Solé and Deisboeck [25] , which describes the competitive behaviour between cell populations with different levels of genomic instability. The model is given by the following set of differential equations:

 ⎧⎪ ⎪⎨⎪ ⎪⎩˙x0=f0Qx0−Φ(x0,x1,x2)x0,˙x1=f0(1−Q)x0+f1Q′x1−Φ(x0,x1,x2)x1,˙xi2=f1(1−Q′)q′ix1+∑nj=1fj2μijxj2−Φ(x0,x1,x2)xi2.

Variable is the fraction of cells with anomalous growth but no genetic instability; is the fraction of cells derived from by mutation that allows genetic instability; and for , is the fraction of tumor mutant cells that can be generated from due to mutation (see Fig. 1). Thus, , and denote the fraction of each population and, therefore, the sum of all these variables must be exactly one, i.e., . Moreover, the probability of mutation from to can be denoted by and is given by . Notice that . Coefficients denote the mutation rate (or probability) from to . The probability of error-free replication of is represented by . Additionally, cross-mutations connect the different subclones of through the term , being the growth rate of each population and the growth rate of each subpopulation of . These and are known as the fitness of each population. The greater is the fitness of a population, the greater is its replicative rate. Finally, the term is the average fitness of the population vector , i.e., . is also known as the “constant population constraint” and it ensures that the population remains constant, also introducing competition between the three populations of cells. Although the model does not explicitly consider environmental constraints, such as blood supply, hypoxia or acidosis, they can be considered as implicitly introduced through the set [12]. Population mutates to with a probability of . In the same way, mutates in different -sequences with a probability . In both cases, , being and the copying fidelity during replication for and respectively. So, in this model, mutations from to and from to have not been considered. It is also worth noticing that the model with is the two-variable quasispecies model [26].
The set of equations (2) can also be written in a matricial way:

 ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝˙x0˙x1˙x12⋮˙xn2⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝f0Q0f0(1−Q)f1Q′0f1(1−Q′)q′1⋮⋮0f1(1−Q′)q′n[]ccc0⋯00⋯0MμDf2⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝x0x1x12⋮xn2⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠−Φ(x)⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝x0x1x12⋮xn2⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (1)

where and

So, if we set and as the following matrices

 Unknown environment '%

then is a Markov matrix by columns. This means that . This kind of matrix appears in mathematical models in biology [29], economics (e.g., the Markov Switching Multifractal asset pricing model [30]), telephone networks (the Viterbi algorithm for error corrections [31]) or ”rankings” as the PageRank algorithm from Google [32, 33]. Therefore, the system can be rewritten as:

 ˙x=MDfx−Φ(x)x. (2)

Our approach to this problem consists on assuming behaves as an average variable [25]. As a consequence, only two different mutation rates are involved in such simplified system: and . Hence, the set of equations is given by:

 ⎧⎪⎨⎪⎩˙x0=f0Qx0−Φ(x0,x1,x2)x0,\par˙x1=f0(1−Q)x0+f1Q′x1−Φ(x0,x1,x2)x1,\par˙x2=f1(1−Q′)x1+f2x2−Φ(x0,x1,x2)x2. (3)

In this simplified system, the effect of mutations can be represented by means of a directed graph as shown in Fig. 1. From now on, let us assume that we always start with a population entirely composed by stable cells, i.e., and . Notice that for this particular case every point of the trajectory always verifies .

The model by Solé and Deisboeck [25] considered , where was a decreasing function such that , indicating that the speed of fitness decays as mutation increases, and where was a competition term. The dependence of on and was introduced to take into account the deleterious effects of high genetic instability on the fitness of the unstable population . In our work we will consider as a constant, being independent of and . This assumption allows the analysis of a more general scenario where population can produce by mutation, and where can have different fitness properties depending on the selected values of both and . For instance, we can analyze the case of deleterious mutations from to with , the case of neutral mutations from to when , or the case of increased replication of with . This latter case would correspond to mutations in driver genes that might confer a selective advantage to the tumor populations .

## 3 Results and discussion

### 1 Equilibrium points

The system (3) has three different fixed points. Namely, a first one showing the total dominance of the unstable population, i.e., and go extinct and then . A second one where goes extinct and coexists with , and a third possible fixed point, where the three populations coexist. We can not consider a fixed point with , and because, as shown in Fig. 1, cells only mutate in one direction and it is not possible to have such scenario. Furthermore, the trivial equilibrium is never achieved because .

Let us seek for these three possible fixed points. Then we have to find the solutions of the following system:

 ⎛⎜⎝000⎞⎟⎠=⎛⎜⎝f0Q00f0(1−Q)f1Q′00f1(1−Q′)f2⎞⎟⎠⎛⎜⎝x0x1x2⎞⎟⎠−Φ(x)⎛⎜⎝x0x1x2⎞⎟⎠, (4)
• If we consider and , from the system (4) we obtain , which has two possible solutions or . As we require the sum of the three variables to be 1, then the solution is .

• Let us consider now and then solve the obtained system of equations:

 Unknown environment '% (5)

Its solution is and . Notice that, as and are ratios of population, and must take values from to . Then, from their expressions and remembering that , we get the conditions: and .

• Finally, for the last fixed point we consider . Then, from the first equation of the system we have . Hence the new system has this form:

 {f0(1−Q)x0+(f1Q′−f0Q)x1=0,f1(1−Q′)x1+(f2−f0Q)x2=0, (6)

under the constraint . Then its solution is as follows:

 x∗0=(f0Q−f1Q′)(f0Q−f2)φ(f0,f1,f2,Q,Q′),
 x∗1=f0(1−Q)(f0Q−f2)φ(f0,f1,f2,Q,Q′),

and

 x∗2=f1(1−Q′)f0(1−Q)φ(f0,f1,f2,Q,Q′).

where

 φ(f0,f1,f2,Q,Q′)=(f0Q−f1Q′)(f0Q−f2)+f0(1−Q)(f0Q−f2)+(f1(1−Q′)f0(1−Q).
• missingmissingmissing3.1.

We can consider as particular case. If these two fitness parameters are equal, the system (3) only has two possible fixed points. As we previously mentioned, this case would correspond to the production of neutral mutants from population to .

From now on, we provide numerical results1 of the system (3). To compute the solutions of this system we have used the Taylor software. Taylor is an ODE solver generator that reads a system of ODEs and it outputs an ANSI C routine that performs a single step of the numerical integration of the ODEs using the Taylor method. Each step of integration chooses the step and the order in an adaptive way trying to keep the local error below a given threshold, and to minimize the global computational effort [27].

Numerically, we integrate the solution of system (3) with initial condition . The next Lemma ensures that any point of this orbit satisfies . It is used as an accuracy control while integrating the ODE system.

###### Lemma 3.1

is a first integral of system (3), that is .

• From the equations (3), we know that

 Extra close brace or missing open brace

Since the initial condition is and , and the assertion follows.

As a particular case we consider , and . Then we make to take values from to and to take values from to both with a step of . We can not start with or because this way or would become extinct. So, we compute the analytical result for all the possible fixed points and, with the iterative method, we integrate the ODE until the distance between the result of the iterate and one of the fixed points is less than an error tolerance previously set, in our case . We also fix an upper limit for the time taken by the system to reach a pre-established distance to the fixed point, therefore we consider it does not reach any fixed point if this upper limit is surpassed.

These computations are displayed in Figure 2, which shows the fixed points that are reached by the system depending on the parameters and . Here each fixed point is indicated with a different colour. Red indicates that the fixed point is ; blue indicates that the fixed point is ; and green indicates the fixed point where all of subpopulations coexist. In addition, black lines indicate that none of the fixed points were reached by the system in the given time. Notice that in our analysis we are tuning and . The decrease of these pair of parameters is qualitatively equivalent to the increase of mutation rates and , respectively, since and . For instance, going from to can be achieved increasing .

Figure 3 shows the density of each population at the equilibrium point in the same parameter space of Fig. 2. This analysis allows us to characterize the regions of this parameter space where the density of stable cells is high, and the malignant population remains low. As Solé and Deisboeck [12] identified, this scenario can be achieved by increasing . Furthermore, our results suggest that this behavior is robust to changes in , whenever remains small (i.e., high copying fidelity ). Notice that an increase of (or, alternatively, a decrease of ) makes the population densities of and to increase and decrease respectively, while for this range of parameters the values of remain low. This highlights the relevance of targeted mutagenic therapies against tumor cells, while the most stable cells should keep replicating with a high fidelity. Such a therapy may slow down tumor growth, which could be eventually be cleared due to demographic stochasticity because the small population numbers for found for this combination of parameters (i.e., large together with low ).

We observe that there is a frontier between the different fixed points reached by the system. The pass from one fixed point to another can be given by a bifurcation. In Fig. 4 we appreciate these bifurcations (see dashed vertical lines) more clearly (in the next section we will study the bifurcations in detail). For this purpose, we have set four sections of Fig. 2 (white arrows), i.e., we have fixed different values for and . In this case these values are and .

We see the bifurcations of the fixed points of the system in Fig. 4. In and has fixed values and respectively, so bifurcations are represented depending on the replication fidelity of in both cases. We notice that for high replication fidelity only and coexist. But in the case is different: if we study this graphic in terms of the mutation rate of we can observe that for high values of it, the system stays at an equilibrium point where is greater. That means that if we make the mutation rate higher through therapy it is possible to make the tumor to achieve an equilibrium state where the most unstable cells are near to and the whole population is mainly dominated by .

In cases and we have considered fixed values and and, therefore, now bifurcations are represented in terms of the replication fidelity of . Notice that in both cases the higher is the probability of to stay stable (equivalently, the mutation rate of is low) the higher is the population of stable cells at the equilibrium point. This is crucial since they correspond to final scenarios with an important presence of genetically stable cells. Comparing these cases, we observe that, if the mutation rate of is higher, i.e., case (c) (, i.e., ), the population in the equilibrium point reached is also higher. This suggests the existence of a threshold in the unstability of bringing more stable cells into the final equilibria.

### 2 Stability analysis and bifurcations

The linear stability analysis of the fixed points characterized in the previous Section is performed by using the Jacobi matrix:

 Lμ(x∗)=⎛⎜⎝f0Q−Φ(x)−x0f0−x0f1−x0f2f0(1−Q)−x1f0f1Q′−Φ(x)−x1f1−x1f2−x2f0f1(1−Q′)−x2f1f2−Φ(x)−x2f2⎞⎟⎠. (7)

We are specially concerned with the domain, in the parameter space, where the malignant cells become dominant and the stability of such equilibrium state. Thus, taking , the Jacobi matrix has the following eigenvalues:

 λ1=f0Q−f2,xxλ2=f1Q′−f2,xandxλ3=−f2. (8)

So, is an attractor if the two inequalities, and , are satisfied. From the expression of the eigenvalues we conclude that there is a critical condition for the mutation rates and , given by: and , respectively (see also [25]). These conditions separate the domain where only remains from the other two cases. From Fig. 5(a) we can confirm that the fixed point is an attractor point until this critical condition i.e., the error threshold, (shown by a vertical dotted line) is reached. From that value, the fixed point is unstable, with local -dimensional stable invariant manifold and a -dimensional unstable one.

 λ1=f0Q−f1Q′,xxλ2=f1(f2−f1Q′)f1−f2,xandxλ3=2f1f2Q′−f21Q′−f22f1−f2. (9)

If we evaluate the stability of the fixed point , where populations and coexist, we get the following eigenvalues from the Jacobi matrix:

Notice that, as we mention in Remark 2, when , such fixed point does not exist, hence the eigenvalues of the Jacobi matrix do not provide any information. This can be appreciated in Fig. 5(b), when , thus .

Finally, we want to study the stability of the fixed point where all populations coexist. An analytical expression for the eigenvalues exists, since we have the Jacobi matrix and the analytical expression for the fixed point itself. But they have really complicated expressions and, even more, we are considering . This means that the greater , the analytical expressions for the eigenvalues are more complicated to find. This is the reason why it is interesting to compute them numerically.

To compute the eigenvalues of this -dimensional matrix we proceed in the following way. This simplified procedure turns to be quite fast in our case (dimension ) but may not be applicable for higher dimensions. It works as follows:

• We first apply the power method to compute an approximation for the eigenvalue of maximal modulus. As it is known, it is based on the idea that the sequence should behave for large values of as the direction of the eigenvector associated to the largest (in modulus) eigenvalue of . To compute it one starts from an arbitrary initial vector (in our case, for instance, , and computes . Provided the difference between the largest eigenvalue and the second largest eigenvalue (in modulus, always) is not too small, the quotients of the components of , that is , converge to . To avoid problems of overflow one normalizes , i.e. , at any step. Problems of convergence appear when the two largest eigenvalues of are close.

• We apply the same method to to obtain the smallest eigenvalue of , since is eigenvalue of iff is eigenvalue of . To compute we have used its -decomposition, which writes as the product of two matrices: an orthogonal matrix and an upper triangular matrix , i.e., .

• Provided , are accurate enough, the third eigenvalue is determined from the value of the determinant of the matrix , that we have derived from its -decomposition.

We apply this procedure to the study of possible bifurcations in the stability of the equilibrium points obtained when we do not have analytic expressions for the associated eigenvalues. To show it, we fix a value for and move . We have chosen a value for corresponding to the line of the diagram in Fig 2. Observe that, when moving at the green zone our equilibrium point is of the form while it is of the form when we move along the blue one. For each equilibrium point, whose components depend on , we compute its differential matrix and its associated eigenvalues. These eigenvalues, computed numerically as mentioned above, are plotted in Fig. 5(c). They are all three real, and thus no transient oscillations are found in the dynamics. Observe that there is an interchange of the number of positive and negative eigenvalues around the bifurcation value , but no change in their stability. They are always unstable since we have at least one positive eigenvalue.

It is interesting to highlight the change in the geometry of such equilibrium point: for values of under the bifurcation value its invariant stable manifold is -dimensional and its invariant unstable manifold -dimensional; on the contrary, after the bifurcation the associated dimensions of the invariant manifolds are exchanged. In both cases, the orbit starting at initial conditions finishes at the stable manifold of that equilibrium point. This is why, although being a unstable fixed point, the trajectory under analysis asymptotically travels towards the fixed point.

### 3 Transient times

In this section we analyze the time taken by the system to reach a given (small) distance to one of the fixed points found in section 1. Typically, the behaviour of transient times change near bifurcation threshold, and, particularly for our system, we are interested in possible changes in transients due to changes in mutation rates. These phenomena could be relevant in patient response under mutagenic therapy.

Figure 7 shows the same as Fig. 2, in the sense that we tune the same parameters of the system, but now we compute the transient times. Here we use the same method to integrate the ODE system, but taking into account only the time taken to reach a distance of the equilibrium point. Notice that if we observe the contour lines represented at the bottom of Fig. 6 we have the same as in Fig. 2. This means that the “time to equilibrium” of the system increases when we are near bifurcation points, which are represented with black lines in Fig. 2.

When slightly modifying the mutation rate of near the error threshold, we see that the time taken by the tumor to reach an equilibrium state is significantly lower. In particular, when increasing the mutation rate, not only time decreases but also we see that the tumor reaches an equilibrium point where the stable cells population is higher. This can be appreciated in Fig. 7, where red, green and blue lines represent , and respectively and the dashed, orange line represents the transient time, which is rescaled to be able to relate it with the corresponding equilibrium point. In case of possible mutagenic therapies directly targeting the most unstable cells, our results reveal that pushing mutation rates of unstable cells beyond the error threshold (where population starts increasing and population starts decreasing) the achievement of the equilibrium point would be very fast, and no further increase of mutation rates may be needed to achieve a faster response (notice that the transient time for in Fig. 6 does not significantly decrease for smaller values). In other words, when rescaling Fig. 6 (b) in terms of logarithm in base 10 (as we can see in Fig. 7) if we increase the mutation rate of , time decreases faster than an exponential function, which is interesting result since it implies an equilibrium scenario with a dominant presence of that can be reached in a short time.

## 4 Conclusions

In this article we study an ODE system modeling the behavior of a population of tumor cells with a mean-field quasispecies model introduced by Solé and Deisboeck [12]. Thus, in our simplified model we have assumed the following: the model does not consider stochasticity and does not take into account spatial correlations between cells as well. Also, we do not consider cell death (due to e.g., the immune system), otherwise we model competition in terms of replication and mutation using the quasispecies framework [14, 15].

First, we have found the fixed points of the system analytically and we have studied their stability both analytically and numerically. We have also characterized the bifurcations between the different fixed points. We conclude that, depending on the parameters, the system can reach different equilibrium states one of them involving low populations of tumor cells while keeping large populations of genetically stable cells. This scenario can be achieved by increasing the level of genetic instability conferred through the mutation rate of the mutator-phenotype population . If this mutation rate exceeds an error threshold, the replication rate of the more malignant subpopulation is reduced to a point where it has no longer competitive advantage.

Further analysis of the effects of mutation rate on the dynamics of the system have allowed us to characterize a transcritical bifurcation. Under such a bifurcation, the time taken by the system to reach the desired equilibrium state is shown to drastically decrease with slight changes in the mutation rate near the error threshold. This result can give us a clue about how medication and therapy may affect the tumor behavior in case of direct mutagenic therapies against tumor cells. In other words, it is possible to modify the mutation rate of the cells through therapy, and we have seen that if we slightly increase the mutation rate of near the error threshold we can obtain an equilibrium point with a dominant population of stable cells rapidly. Such an scenario could become relevant since populations of tumor cells could be decreased with mutagenic therapy, and such populations could be more prone to extinction due to stochastic fluctuations inherent to small populations sizes.

Further research should explicitly consider stochastic fluctuations by transforming the differential equations in stochastic differential equations e.g., Langevin equations. Also, the addition of spatial correlations would indicate if our results are still present in theoretical models for solid tumors. In this sense, it would be interesting to analyze the transient times tied to increasing mutation values considering space in an explicit manner. It is important to notice that transient times are extremely important in terms of responses velocities after therapy. Different causes and mechanisms have been suggested to be involved in transient duration, which have been especially explored for ecological systems [34]. For instance, it is known that space can involve very long transients until equilibrium points are reached. This phenomenon could be investigated extending the model we investigated by means of partial differential equations.

## Acknowledgements

J.T.L. has been partially supported by the Spanish MICIIN/FEDER grant MTM2012-31714 and by the Generalitat de Catalunya number 2014SGR-504. J.S. has been funded by the Fundación Botín.

### Footnotes

1. The codes used to obtain the results presented in this work are available upon request

### References

1. Cairns, J., 1975, Mutation selection and the natural history of cancer. Nature 255, 197-200.
2. Merlo, L.M.F., Pepper, J.W., Reid, B.J., Maley, C.C., 2006. Cancer as an evolutionary and ecological process. Nature Rev. Cancer 6, 924-935.
3. Cahill, D.P., Kinzler, K.W., Vogelstein, B., Lengauer, C., 1999. Genetic instability and darwinian selection in tumors. Trends Genet. 15, M57-M61.
4. Loeb, L.A., 2001. A Mutator Phenotype in Cancer. Cancer Res. 61, 3230-3239.
5. Oliver, A., Canton, R., Campo, P., Baquero, F., et al., 2000. High frequency of hypermutable Pseudomonas aeruginosa in cystic fibrosis lung infection. Science 288, 12514.
6. Bjedov, I., Tenaillon, O., Gérard, B., Souza, V., et al., 2003. Stress- induced mutagenesis in bacteria. Science 300, 14049. 35.
7. Sniegowski, P.D., Gerrish, P.J., Lenski, R.E., 1997. Evolution of high mutation rates in experimental populations of E. coli. Nature 387, 7035.
8. Anderson, G.R., Stoler, D.L., Brenner, B.M., 2001. Cancer as an evolutionary consequence of a destabilized genome. Bioessays 23, 103746.
9. Bielas, J.H., Loeb, K.R., Rubin, B.P., True, L.D., et al., 2006. Human cancers express a mutator phenotype. Proc. Natl. Acad. Sci. USA 103, 18238-18242.
10. Hoeijmakers, J.H., 2001. Genome maintenance mechanisms for preventing cancer. Nature 411, 366-374.
11. Benedikt, B., Siebert, R., Traulsen, A., 2014. Cancer initiation with epistatic interaction between driver and passenger mutations. J. Theor. Biol. 358, 52-60.
12. Solé R.V., 2001. Phase transitions in unstable cancer cells populations. The European Physics Journal B. 35, 117-123.
13. Brumer, Y., Michor, F., Shakhnovich, E.I., 2006. Genetic instability and the quasispecies model. J. Theor. Biol. 241, 216-222.
14. Eigen, M., 1971. Selforganization of matter and evolution of biological macromolecules. Naturwissenschaften 58, 465-523.
15. Eigen, M., Schuster, P., 1979. The hypercycle: a principle of natural selforganization (Springer-Verlag New York). Naturwissenschaften 58, 465-523.
16. Solé, R.V., Sardanyés, J., Díez, J., Mas, A., 2006. Information catastrophe in RNA viruses through replication thresholds. J. Theor. Biol., 240, 353-359.
17. Marcus, P. I., Rodriguez, L.L, Sekellick, M. J. Interferon induction as a quasispecies marker of Vesicular stomatitis virus populations. J. Virol., 72(1), 542-549.
18. Cichutek, K., Merget, H., Norlay, S., Linde, R., et al., 1992. Development of a quasispecies of human immunodeficiency virus type 1 in vivo. Proc. Natl. Acad. Sci. USA, 89, 7365-7369.
19. Domingo, E., Biebricher, C.K., Eigen, M., Holland, J.J., 2001. Quasispecies and RNA virus evolution: principles and consequences. Landes Bioscience, Georgetown Texas.
20. Bull J.J., Sanjuán R., Wilke C.O., 2007. Theory of lethal mutagenesis for viruses. J. Virol. 81, 2930-2939.
21. Wylie, C.S., Shakhnovich, E.I., 2012. Mutation Induced Extinction in Finite Populations: Lethal Mutagenesis and Lethal Isolation. PLoS Comput. Biol 8(8), e1002609.
22. Lauring, A.S., Andino, R., 2010. Quasispecies Theory and the Behavior of RNA Viruses. PLoS Pathog. 6(7), e1001005.
23. Sardanyés, J., Simó, C., Martínez, R., Solé, R. V. and Elena, S.F., 2014. Variability in mutational fitness effects prevents full lethal transitions in large quasispecies populations. Nature Scientific Reports 4, 4625.
24. Crotty, S., Cameron, C.E., Andino, R., 2001. RNA virus error catastrophe: Direct molecular test by using ribavirin. Proc. Natl. Acad. Sci. USA, 98, 6895-6900.
25. Solé R.V., Deisboeck T.S., 2004. An error catastrophe in cancer?. J. Theor. Biol. 228, 47-54.
26. Swetina, J., Schuster, P., 1982. Self-replication with errors. A model for polynucleotide replication. Biophys. Chem. 16, 329-345.
27. Jorba, Á., Zou, M., 2004. A software package for the numerical integration of ODEs by means of high-order Taylor methods. See http://www.maia.ub.es/ angel/taylor/taylor.pdf
28. Nowak, M.A., 1992. What is a quasispecies?. Trends in Ecology and Evolution 7, 118-121.
29. Usher, M.B., 1971. Developments in the Leslie model in mathematical models in ecology. Blackwell Scientific Publications.
30. Lux, T., Morales-Arias, L., Sattarhoff, C., 2011. A Markov-switching Multifractal Approach to Forecasting volatility. Kiel Working Paper, 1737.
31. Forney, G.D., 1973. The Viterbi algorithm. Proceedings of the IEEE, 61, 268-278.
32. Sergey, B., Lawrence, P., 1998. The anatomy of a large-scale hypertextual Web search engine. Computer Networks and ISDN Systems, 33, 107-117.
33. Bryan, K., Leise, T., 2006. The \$25,000,000,000 eigenvector. The linear algebra behind Google. SIAM Review, 48 (3), 569-81.
34. Hastings, A., 2004. Transients: The key long-term ecological understanding, Trends in Ecol. Evol. 19, 39-45.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters