We point out that a simple and generic strategy in order to lower the risk for extinction consists in developing a dormant stage in which the organism is unable to multiply but may die. The dormant organism is protected against the poisonous environment. The result is to increase the survival probability of the entire population by introducing a type of zero reproductive fitness. This is possible, because the reservoir of dormant individuals act as a buffer that can cushion fatal fluctuations in the number of births and deaths, which without the dormant population would have driven the entire population to extinction.
Drug resistance is a very serious problem not least for chemotherapeutic treatment of cancer, and it is very important to understand, as far as possible, the ways to countermeasure drug resistance in a number of biomedical contexts. In order to unravel the possible mechanisms responsible for this phenomena, Iwasa and co-workers (Iwasa et al. 2003, 2004) developed a prominent theoretical model of the evolutionary dynamics of escape. They base their approach on the assumption that n point mutations in some crucial parts of the genome are necessary for escape. They further assumed that the different mutants can be described by binary strings (with entries of +1 or −1) of length n. There are 2n − 1 such mutants. It is assumed that treatment reduces the proliferation ratios of sensitive mutants, R < 1; whereas resistant mutants are such that R > 1. The corresponding evolutionary dynamics is modelled in terms of Galton–Watson multi-type branching process (GWMBP; Kimmel & Axelrod 2002), where at each generation each individual of each type has a given (in general, mutant-dependent) probability of mutating and producing offspring belonging to a different type. The problem is to calculate the probability that a resistant mutant is reached within a population of size N. The model proposed by Iwasa and co-workers has been analysed in more detail by Serra (2006) and Serra & Haccou (2007). While this approach is both biologically and mathematically sound, we believe that the efficiency with which evolutionary escape allows organisms to avoid extinction, when attacked by lethal drug, can be understood in a more simple way as a consequence of a dormant phase. To assess the generic importance of dormancy as an escape strategy, we develop a very simple model framework, in which we only focus on the most essential aspects of population and evolution dynamics.
Here, we show that the existence of a dormant phase can be crucial for a population to escape extinction. We consider the following simplified scenario. We study the survival probability of a population of organisms that can exist in the form of three different types (figure 1). The types differ in their response to the presence of a drug. Types (1) and (2) have similar reproduction and death rates when no drug is present. However, the drug is supposed to be lethal to type (2) but neutral to type (1). We are interested in how the existence of a dormant mode, type (3), affects the survival probability of the entire population. The dormant type can neither reproduce, nor is it susceptible to the drug. However type (3) can die and it can undergo a transformation back to type (1). To understand the effect of the dormant type, we consider different realizations (indicated in figure 1) of the possible flow between the three different types.
Our model is relevant to a number of biological cases, in particular to entities such as cancer cells, which have been observed to evolve resistance to therapy: treatments impose a selective pressure that eliminates the lesser fit strands of the corresponding populations but also drives an evolutionary process, whereby better adapted individuals, i.e. individuals immune to the effects of the corresponding drug, eventually take over. Further rationale for our model is provided by the response of tumour cells to hypoxia (i.e. oxygen starvation). It is a well-known fact that hypoxia induces arrest of the cell cycle (Gardner et al. 2001; Alarcón et al. 2004). This means that the rate at which cells replicate under hypoxia is drastically reduced, as hypoxia downregulates the activity of the pathway regulating the progression through the cell cycle (Gardner et al. 2001). Moreover, another well-understood fact about cell survival/death regulation is that there exists cross-talk between the pathways regulating cell death and cell division (Nishioka & Welsh 1994; Padmanabhan et al. 1999; Alberts et al. 2002). This means that, in normal circumstances, the downregulation of the cell cycle machinery implies the downregulation of the apoptotic (programmed cell death) machinery. It can be argued that such cross-talk is disrupted in cancer cells, and that cancer cells are such that cell death regulatory pathways are downregulated. Therefore, under hypoxic conditions, cancer cells undergo a drastic reduction of their division and death rates.
The study of the influence of hypoxic cancer cells fits rather well within the original remit in which the dynamics of evolutionary escape was put forward, as hypoxic cells are known to play a major role in the resistance to chemotherapy and radiotherapy in tumours (Bristow & Hill 2008). To study the influence of such a subpopulation on the global dynamics of the total population, we model it as a dormant or quiescent population with negligible proliferation rate and small death rate.
Previous works have analysed some of the effects of quiescent cells on tumour growth (e.g. Alarcón et al. 2003; Brikci et al. 2008). All these previous studies hint on the role of quiescence cells in the dynamics of tumour growth and its role in resistance to tumour growth, but the issue of the evolutionary dynamics involved is not directly addressed.
We now turn to the study of the survival probability. We disentangle the interplay between drug susceptibility, type (2) and dormancy, type (3), by analysing three different versions of the population dynamics, all of which are depicted in figure 1 and formulated as the following three models.
Model A. This is the ‘normal’ situation, where no drug is present. Types (1) and (2) are essentially equivalent, except that when type (1) reproduces, it may undergo a ‘mutation’ and end up as type (3). When type (2) reproduces it may mutate and become type (1).
Model B. This is the situation in the presence of the drug. The death rate of type (2) is now significantly bigger than the death rates of types (1) and (3). Moreover, the drug makes type (2) unable to reproduce and therefore the flow from type (2) to type (1) is absent.
Below, we will find that the dormant stage increases the survival probability of the population in the presence of the drug. To emphasize that the increased ability of the population to escape extinction is in fact caused by the presence of a dormant stage, we also consider an extreme version of the dynamics in which type (1) is unable to die.
Model C. This version of the the dynamics is equivalent to model B except that type (1) now is assumed not to be able to die directly, but has to flow through either type (2) or type (3) to do so. We then demonstrate that even in this extreme situation does the availability of the dormant stage enhances the population's chance for avoiding extinction.
We notice that the different types (1), (2) and (3) can be thought in epidemiological terms as susceptible, infected and immune. The version of the dynamics defined as model B can be thought to represent an age-structured population. Type (3) is then juveniles, type (1) are mature reproduction active individuals and type (2) are individuals in the post-reproductive stage.
An economical and precise way to present the dynamics is in terms of the generator functions for the corresponding Galton–Watson multi-type branching process. We include these generator functions in the tables in the appendix A. Within the theory of multi-type branching process, the condition for eventual non-extinction is given in terms of the spectral radius, ρ, of the matrix A = (aij), whose entries are the expected values of the number of offspring of type j produced of an individual of type i. If ρ > 1, there is a finite probability of N(t → ∞) > 0. The quantities aij are calculated in terms of the corresponding generating functions (table 2): Aij = ∂jGi|x̄=1. From these standard arguments of the theory of multi-type branching process, which carry on in a straightforward manner when size-dependence is taken into account (Jagers 1999), we can establish from tables 1 and 2 in appendix A, that the condition for asymptotic survival, that is for P[N(t → ∞) > 0] > 0, is 1.1 i.e. r12 < 1/2. Otherwise, the probability of eventual extinction is 1. Now, consider the presence of a third type (type (3)) within the population, namely, quiescent individuals. These individuals are not allowed to proliferate but are resilient to the environment and can survive under hostile conditions, hence we assume that their death rate ε ≪ 1. In addition, quiescent cells are assumed to be able to revert back to type (1) at a given rate, which depends on the availability of resources. The issue we intend to analyse is whether the introduction of a quiescent subpopulation helps to escape from the whole population being extinct. More precisely, the question we aim to address is: assume r12 ≥ 1/2, is it possible for a population whose dynamics is described in model B and table 2 to elude eventual extinction? specifically, this scenario has been considered within the context of modelling of tumour growth, in particular, in the response of cancer cells to hypoxia (low levels of oxygen; Alarcón et al. 2003). Cancer cells appear to become quiescent in response to hypoxia, which is thought to give them an advantage in their competition with their normal counterparts as well as resistance to radiotherapy and chemotherapy (Brown & Wilson 2004).
We now present the results of simulations of the survival probability for the three different scenarios represented by models A, B and C.
2. Simulation results
Numerical simulations of the underlying multi-type branching process (see appendix A for details) confirm that, in fact, the presence of a quiescent population as a mechanism to escape from harsh conditions is indeed feasible. Figure 2 shows how the survival probability depends on the probability of an individual of type (1) to become quiescent, r13. We can see that as r13 decreases, the threshold for survival in terms of r31, i.e. the probability of a quiescent individual to revert to type (1), moves towards smaller values of the parameter r31. This means that, following a decrease in the flux of individuals from the type (1) population to type (3), survival is only possible by reducing the inverse flux. This observation reveals that the mechanism by which quiescence helps escape is by acting as a reservoir, where part of the population can be safely ‘stored’. If the flux from type (3) back to type (1) is too big, it effectively increases the flux from type (1) to type (2), thus increasing lethality. This behaviour, however, is sensitive to the decay rate of type (3) individuals, ε: as this parameter increases, the survival probability decreases (figure 4). This means that the death rate of quiescent individuals being small is an instrumental factor for quiescence-induced escape from harsh environments⇔.
Figure 3 shows that as r12 increases further into the regime, where equation (1.1) predicts sure extinction of the type (2) population, survival is only possible by decreasing r31, i.e. by increasing the average time individuals spent in the quiescent state.
To test this scenario further, we perform simulations in which the following situation is considered. First, we let a population whose dynamics is given by model A, evolve for some time in an environment free of the hostile agent. In this environment, individuals of both types (1) and (2) can thrive. The dynamics in an environment free of any hostile agent is described in terms of the generating functions shown in table 1. The introduction of the hostile agent is described by changing the dynamics of the population to the one described by model B (figure 1 and table 2), which is essentially the same as model A (table 1) but with type (2) individuals doomed to perish. We further assume that the hostile agent is active only for a given period of time, after which its detrimental effect on type (2) organisms ceases and the dynamics of the population reverts to model A (according to table 1). Parameter values are chosen so the dynamics of the agent-free population is super-critical (i.e. the corresponding survival probability is bigger than zero). After the hostile agent has been removed, we calculate the probability that the population survives both with and without quiescence. This scenario is highly relevant to evolutionary escape from drugs, which have a finite lifespan. This particular escape problem was considered (Iwasa et al. 2003).
The corresponding results are shown in figures 5 and 6. Figure 6 shows the survival probability, P[N > 0], as a function of r13 and, indeed, we observe that when quiescent individuals are not considered, i.e. r13 = 0, the probability of surviving after the end of the activity period of the hostile agent is null. As r13 increases, so does the corresponding survival probability, thus further confirming that quiescence is a feasible mechanism to help biological populations to escape harsh environmental conditions. Moreover, figure 5, which shows particular realizations of the process with (figure 5b), and without (2a) quiescence, also helps in understanding more about the mechanism, whereby quiescence allows escape: as can be seen in figure 5b, during the period of activity of the hostile agent most of the population is actually quiescent with a little proportion of the population in types (1) and (2). This means that quiescence mediates survival by acting as a reservoir or buffer.
3. An alternative model
Although the results discussed until here provide compelling evidence in favour of hypothesis that the presence of quiescence can induce escape without the need of an increase in the reproductive fitness of any particular type of organism, some doubt could be cast on our argument so far. One might feel inclined to argue that our model still relies on a multi-type population, where one type (type (1)) is effectively more fit than the others. To address this possibility, we consider model C in figure 1. This version is closely related to the current discussion concerning the response of cancer cells to oxygen starvation. This model is also formulated in terms of a GWMBP and characterized in terms of the generating functions given in table 3. Its rationale is as follows. Let us assume that a population of cells is divided into two types: those in which the pathways regulating cell division are activated (type (1)) and those in which the active pathways are those regulating apoptosis, i.e. cell death (type (2)). Type (1) cells can either proliferate, or stay as they are, or suffering activation of the apoptotic pathways, thus becoming type (2) cells. The flux of population between types (1) and (2) is controlled by the parameter r12. Type (2) cells are those marked for cell death, but the model provides for their staying within the population for a while (not dying immediately), but they cannot proliferate. Type (3) cells are, as in the previous model, quiescent cells. In this context, the introduction of a hostile agent (drug, removal of oxygen, etc.) corresponds simply to increasing the value of r12.
In the absence of quiescence, standard arguments reveal that extinction with probability 1 occurs when r12 ≥ 1/2. The question is once again whether quiescence can rescue the population from extinction under such conditions. Figure 7 shows that this is indeed the case, provided that the flux out of the quiescent state back into type (1) does not exceed a critical value. These results are completely analogous to those obtained for the previous model.
4. Invasion dynamics
We now discuss under which conditions a small proportion of individuals, y, which can undergo quiescence, modelled by models A and B depending on whether a drug is present or not, respectively, and with r13 ≠ 0, can take over a population of fully growing individuals, i.e. a population modelled by models A and B with r13 = 0. In particular, we will analyse the corresponding invasion probability, PF(y), by direct simulation of the population dynamics and by using the analytical approximation provided by the so-called evolutionary formalism (Arnold et al. 1994; Demetrius 1997; Alarcón & Jensen 2009; Demetrius et al. 2009), which allows a more thorough exploration of the behaviour of the invasion probability as a function of the model parameters.
The problem we address here relates to whether a small population of mutants can take over an incumbent population. Demetrius and co-workers (Demetrius 1997; Demetrius et al. 2009) have developed a formalism based on the application of ideas from ergodic theory to evolutionary problems (Arnold et al. 1994) and the diffusion approximation (Feller 1951). This formalism allows us to estimate the fixation probability of a mutant population in the presence of an incumbent species. Here, we generalize this formalism to apply it to the problem of whether a genetic inactivation generating new phenotypes can invade the incumbent population.
The starting point of this formalism is the following fundamental equation: 4.1 where r = log(Λ0) is the growth rate (or Malthusian parameter), with Λ0 being the dominant eigenvalue of A, whereas H and F are the entropy and the proliferative potential, defined as: 4.2 where is defined as: 4.3 with , and π is the stationary distribution associated to , with given by . Equations (4.1–4.3) are derived from a variational principle (Arnold et al. 1994), namely 4.4 where ν is a Markov measure, Hν = −∑j=1nν(Aj) log ν(Aj) and γ = logax0x1, where Aj, j = 1, … , n is a partition of the phase space of the system and ax0x1 is the transition probability between two states of the system, x0 and x1. The solution to this variational problem produces equations (4.1–4.3) (Billingsley 1965; Arnold et al. 1994).
According to Demetrius et al. (2009), the diffusion approximation yields the following equation for the fixation or invasion probability as a function of the initial concentration of mutants, y is given by: 4.5 where the total population N is assumed to be constant and s is defined by: 4.6 where Δr ≡ rq − r, Δσ2 = σq2 − σ2, with rq(r) is the growth rate of the quiescence (normal) population, and σq2(σ2) is the variance of the quiescent (normal) population, 〈N〉 is the average stationary population.
The demographic parameters (i.e. r, σ2, rq, and σq2) can be estimated from the evolutionary formalism (see Alarcón & Jensen 2009; Demetrius et al. 2009). The growth rates are given by r = log λ0, where λ0 is the dominant eigenvalue of the of the matrix A = (aij) = (∂jGi(x̄)|x̄=1). In Demetrius et al. (2009), it is shown that the parameter σ2 can be obtained by slightly perturbing the parameters that determine the dynamics of the system, i.e. the mean-field dynamics being given by A(δ) = (aij1+δ), and then doing an expansion for small δ. Accordingly, σ2 is given by: 4.7
Hence, in the linear approximation with H(δ) given by H(δ) = H + δHδ, we have that σ2(γ) = −Hδ, which is given by: 4.8 where πi(δ) = πi + δπi(δ) and pij(δ) = pij + δpij(δ) are the corresponding linear approximations to π and P = (pij), when we take A(δ) = (aij1+δ) δ ≪1. The details of how these quantities are actually calculated is given in appendix A of Alarcón & Jensen (2009).
In order to show that this formalism can be used to analyse the invasion of a population possessing a dormant type, we have done simulations where one population, which is sensitive to the drug, competes with two different population: one capable of undergoing quiescence and another one which is unable to enter such state. We compare the numerics with the analytical results obtained from equation (4.5) in figure 8. Both analytical and numerical results show that whereas in the former case the quiescent population takes over the sensitive population almost surely, in the latter case it is unlikely that the non-resident population takes over the resident one. It is worth remarking that the good agreement between simulations and analytical result seen in figure 8 implies that the evolutionary formalism developed by Demetrius and co-workers (Demetrius 1997; Demetrius et al. 2009) is able to make the invasion problem considered here analytically tractable.
5. Summary and discussion
We have shown that quiescence is a feasible mechanism for biological populations to escape hostile environments. This mechanism is expected to be very relevant to the important issue concerning population dynamics of cancer cells and the competition mechanisms between cancer cells and their normal counterparts. In particular, our model addresses issues related to the response of cancer cells to oxygen starvation (Brown & Wilson 2004; Bristow & Hill 2008). The mechanism involved in quiescence-dependent escape is in essence simple: it provides a buffer for the population to be safe from the hostile environment. That is, in cases where the population would go extinct owing to a fluctuation in the birth and death events, the reservoir in the quiescence buffer population can bring the population back from the brink of extinction. The efficiency of this mechanism for escape is sensitive to the entry and exit rates to and from the quiescent state, respectively.
The main difference between our mechanism for escape and the one proposed in Iwasa et al. (2003, 2004) is the following. Iwasa and co-workers propose a mechanism based on random search of the state of types of individuals for one that is fitter than the others in the presence of a given selection pressure (drug, etc.). In contrast, we have demonstrated here that the overall survival probability of a population may increase by introducing a dormant type. Since this type is unable to reproduce, its reproductive fitness is zero; nevertheless, the existence of this drug-resistant stage is able to improve the fitness of the entire population inasmuch as the population obtains a higher probability for survival.
The mechanism proposed here shares a number of features in common with the phenomenon of bacterial persistence (Balaban et al. 2004; Levin & Rozen 2006; Gardner et al. 2007; Lewis 2007). Persistence is a form of resistance to antibiotics exhibited by bacterial colonies, where resistance is not acquired by a gene mutation that allows the bacteria to grow exponentially fast in the presence of a particular drug. Instead, as schematically shown in figure 9, killing of bacteria goes through a fast phase, where the population decreases exponentially until the killing rate slows down leaving behind a remnant of cells, which differ from the sensitive ones in that they are in a dormant state but are otherwise genetically identical to their sensitive counterparts. In fact, when the drug is removed, a colony of ‘normal’ bacteria is regrown from the persistent cells. The mechanism for persistence appears to involve a phenotypic switch (Balaban et al. 2004), where cells switch from a rapidly growing, but drug-sensitive phenotype to a dormant but drug-resistant phenotype, although there has recently been the suggestion that persistence might be a social trait (Gardner et al. 2007).
According to Balaban et al. (2004), there exist two types of persister. Type I persister, according to their terminology, are only produced during the exponential growth phase and therefore their numbers are fixed at the time of inoculation and determined by the size of the inoculum. Type II persister, on the contrary, divide and grow continuously, but an order of magnitude slower than their non-persistent counterparts, and their numbers are determined by the total population numbers. The quiescence mechanism put forward here is therefore distinct from type I persister, but very similar to the behaviour exhibited by type II persister, which means that the analysis methods, in particular, the evolutionary formalism used to study the competition between subpopulations exhibiting different behaviour, can be extended to study this type of bacterial persistence.
Another issue in which our model differs from previous work on bacterial persistence is the following. We are considering different populations adopting different strategies, namely, a wild-type population that thrives in the absence of drug composed of two subpopulations: one that is capable of undergoing quiescence and another one that is not. By adopting this scenario, we can study the evolutionary dynamics of these populations and analyse which strategy is evolutionarily stable and which one is susceptible to be invaded. In that respect, our analysis goes beyond, for example, the population models proposed by Balaban et al. (2004).
Within the context of the problem of how the hypoxic subpopulation affects the dynamics of the whole tumour, our model sheds further light on this issue. It is commonly thought that hypoxia increases the probability of survival of the tumour by providing a selective pressure that favours the evolution and survival of more aggressive phenotypes (Bristow & Hill 2008). Our model shows that quiescence by itself is enough to increase the survival probability of the population without further increase in the phenotypic variety of the population. Obviously, the residual population provides a springboard for these evolutionary processes to ensue.
T.A. and H.J.J. gratefully acknowledge the EPSRC for funding under grant EP/D051223.
Appendix A: Simulations
In this appendix, we briefly describe the method we have used to produce our simulation results. Our simulations start with one single individual of type (1). Its offspring is determined by the probabilities of producing descendants as prescribed by the generating functions given in tables 1⇔–3, corresponding to models A, B and C, respectively. In general, the coefficients of the Taylor expansion of Gi(x̄) = ∑l,n,m=0 Pi;lnmx1l x2n x3m are the probabilities per generation per individual (of type i) that a type i-individual produces l descendants of type (1), n descendants of type (2) and m descendants of type (3). In subsequent generations, we go over all the individuals in existence in the last generation and the numbers and types of their descendants within the next generation are calculated in the same way. Simulations are run over 1000 realizations of 1000 generations each. The survival probability, P[N > 0], is calculated as the ratio between the number of simulations such that N(T = 1000) > 0 and the total number of realizations.
- Received March 4, 2010.
- Accepted May 19, 2010.
- © 2010 The Royal Society