## Abstract

In recent experiments on isogenetic cancer cell lines, it was observed that exposure to high doses of anticancer drugs can induce the emergence of a subpopulation of weakly proliferative and drug-tolerant cells, which display markers associated with stem cell–like cancer cells. After a period of time, some of the surviving cells were observed to change their phenotype to resume normal proliferation and eventually repopulate the sample. Furthermore, the drug-tolerant cells could be drug resensitized following drug washout. Here, we propose a theoretical mechanism for the transient emergence of such drug tolerance. In this framework, we formulate an individual-based model and an integro-differential equation model of reversible phenotypic evolution in a cell population exposed to cytotoxic drugs. The outcomes of both models suggest that nongenetic instability, stress-induced adaptation, selection, and the interplay between these mechanisms can push an actively proliferating cell population to transition into a weakly proliferative and drug-tolerant state. Hence, the cell population experiences much less stress in the presence of the drugs and, in the long run, reacquires a proliferative phenotype, due to phenotypic fluctuations and selection pressure. These mechanisms can also reverse epigenetic drug tolerance following drug washout. Our study highlights how the transient appearance of the weakly proliferative and drug-tolerant cells is related to the use of high-dose therapy. Furthermore, we show how stem-like characteristics can act to stabilize the transient, weakly proliferative, and drug-tolerant subpopulation for a longer time window. Finally, using our models as *in silico* laboratories, we propose new testable hypotheses that could help uncover general principles underlying the emergence of cancer drug tolerance. *Cancer Res; 75(6); 930–9. ©2015 AACR*.

This study provides a new perspective on the inherent risks of interventional chemotherapy in patients with cancer by showing how the adaption of even nongenetically unstable cell populations exposed to antiproliferative drugs can be acted upon by selective forces, which drive the outgrowth of rapidly proliferative drug-resistant cell populations.

We describe the evolution of phenotype in a well-mixed PC9 cancer cell population exposed to cytotoxic drugs [whose concentration at time *t* is *c*(*t*)] using both an I-B and an IDE formalism. In general, we consider the case of constant infusion of cytotoxic drugs (i.e., *c*(·):=*c* > 0) and characterize the state of each PC9 cell by its expression levels of two phenotypic traits, survival potential and proliferation potential. In this framework, we identify a PC9 cell as having a low value of survival potential and a high value of proliferation potential, a DTP cell as having a high value of survival potential and a low value of proliferation potential, and a DTEP cell as having a high value of survival and proliferation potential (schematized as in Fig. 1B).

We model the key biologic mechanisms of phenotype evolution in the PC9 cancer cells, namely selection, stress-induced adaptation, and nongenetic phenotype instability, by using three separate mathematical strategies, which are tailored to fit the I-B and the IDE formalisms. Selection is modeled through a proliferation probability *p*(·,·,·) and a death probability *d*(·,·) (see Eqs. 1 and 3), which depend on the levels of proliferative and survival potentials of the cells, as well as on the cell microenvironment at time *t*. Importantly, because DTPs and DTEPs exhibit vastly different proliferation rates, we assume that *d*(·,·) does not vary with proliferation potential. On the other hand, we assume that maintaining a high survival potential in a drug-free environment is costly to a cell and will act to reduce *p*(.,.,.) (refs. 13, 32, 33). Finally, we assume that *p*(·,·,·) is dependent on the total population size to represent competition between cells for space and nutrients. This dependence could also be introduced in *d*(·,·) without altering the qualitative results of the model.

Stress-induced adaptation of cell proliferation level is modeled by an advection term that leads to a decrease in the level of proliferation. The speed of adaptation *v*(·,·;·) depends on the cell microenvironment at time *t*, the level of survival potential, and the average sensitivity of the cell proliferative potential to stress-inducing agents |$\bar \upsilon$|. We assume that adaptation can only occur in the presence of the drugs. Hence, we specify that *v*(·,·;·) = 0 when *c* = 0.

Finally, we assume that small (large) epimutations correspond to small (large) changes in cell properties and also that small epimutations occur at a much higher frequency than large epimutations (34). Therefore, nongenetic phenotype instability is modeled as Brownian motion, or diffusion, in the levels of cell proliferation and survival potentials.

In the I-B formalism, we consider each cell as an individual agent and label it by an index 0 ≤ *i* ≤ *N*(*t*), where *N*(*t*) ≥ 0 is the size of the population at time *t* ∈ [0, *T*], and *T* is the end time of the simulation. The normalized expression levels of the survival-potential and proliferation-potential traits in each cell *i* at time *t* are modeled, respectively, by the random variables *X _{i}*(

*t*):[0, ∞] → [0, 1] and

*Y*(

_{i}*t*):[0, ∞] → [0, 1]. We simulate the evolution of the cell population in discrete time, according to the algorithm pictured in Supplementary Fig. S1A and S1B. Over the time interval between two successive time instants

*t*and

*t*+ Δ

*t*, we first allow each cell

*i*either to proliferate, undergo apoptosis, or remain in a quiescent state according to the respective probabilities

If a cell proliferates, we assume both daughter cells inherit the parent's trait values. After all cells had undergone one iteration of the proliferation and death process and returned to the quiescent state, we then let each cell update its trait values according to the following system of discretized stochastic differential equations:

where |$W_i^1 (t)$| and |$W_i^2 (t)$| are standard normal random variables for all 0 ≤ *i* ≤ *N*(*t*) and *t* ≥ 0, and *D* is the average rate of phenotypic fluctuations.

In the IDE formalism, we consider the cell population to be structured by two continuous, real variables *x* ∈ [0, 1] and *y* ∈ [0, 1], which represent, respectively, the normalized expression levels of survival potential and proliferation potential. The population density of cancer cells is modeled by the function *n*(*x*, *y*, *t*) ≥ 0, where the global population density at time *t* ∈ [0, *T*] is computed as |$\varrho(t)\, = \,\int_0^1 {\int_0^1 {n(x,\,y,\,t)} \,{\rm d}x{\rm d}y}$|, and the time evolution of *n* is governed by the following equation:

We set |$\beta = D^2/2$|, so that a link between the solution of the IDE model and the outcome of the I-B model can be established (at least formally) through the Feynman–Kac formula (35).

In the extended I-B model, we keep the birth and death process unchanged from the original I-B model. However, instead of allowing the DTPs to update at each time step, we consider them to be stem-like so that they do not change their phenotype, unless they proliferate. When a DTP undergoes proliferation, we allow its offspring to update according to *p*_{1}, *p*_{2}, and *p*_{3}, which represent, respectively, the probabilities for symmetric self-renewal, asymmetric self-renewal, and symmetric differentiation. The update mechanisms of the PC9s and DTEPs remain unchanged. This new algorithm is pictured in Supplementary Fig. S1C and S1D.

The assumptions on the functions *p*, *d*, and *v* are summarized in Supplementary Tables S1 and S2. All definitions of the parameter functions are provided in the Supplementary Material.

## Introduction

The interplay between intratumor heterogeneity (1, 2) and perturbations in the tumor microenvironment (3–6) plays a key role in the emergence of drug resistance in cancer cell populations. Intratumor heterogeneity results mainly from genetic modifications (e.g., mutations), leading, through a Darwinian-like process, to the selection of tumor cells expressing phenotypes adapted to the local environment (2, 7). However, intratumor heterogeneity can also emerge from nongenetic processes mediated by either stochastic events or epigenetic modifications (8–10). In fact, perturbations of the tumor microenvironment induced by cytotoxic agents may “instruct” cells to enter a more stress-resistant phenotype (11).

In recent experiments, performed on a number of genetically homogeneous populations of cancer cells (including lung cancer, melanoma, colorectal cancer, breast cancer, and gastric cancer cells), Sharma and colleagues (12) showed that epigenetically regulated changes in phenotype can play an important role in the development of reversible drug tolerance. During these experiments, a small subpopulation of drug-tolerant cells was consistently detected, that could maintain viability in the presence of high-dose drug therapy. These drug-tolerant persisters (DTP) were shown to be nonproliferative and display markers specific to stem cell–like cancer cells. After a period of time, approximately 20% of DTPs changed their phenotype to resume normal proliferation and lost stem-like markers, still in the constant presence of the drugs. The resulting cells were labeled the drug-tolerant expanded persisters (DTEP; Fig. 1A). Interestingly, both DTPs and DTEPs could be drug resensitized by drug-free passaging.

Therefore, the three distinct subpopulations—the parental cancer cells (PC9), DTPs, and DTEPs—that compose the whole cancer cell population at various times during drug treatment, although genetically identical, possess different functional phenotypes. Most notably, they can be characterized by their respective levels of survival potential (i.e., the level of robustness toward life-threatening events in extreme conditions, which in this case can be identified as the level of drug tolerance) and proliferation potential (i.e., the rate of cell proliferation; ref. 13). But what is driving the evolution of phenotypes observed in the PC9 cancer cells? Is it simply a case of selection, in which cells with certain properties survive and proliferate better in a given environment (14)? Are individual cells changing their properties in response to environmental cues (11, 15)?

During drug therapy, the average level of survival potential in the PC9 cell population is nondecreasing over time. We argue that this can be explained by the interaction between two mechanisms. The first mechanism is nongenetic instability (16), which results in phenotypic fluctuations without alteration of the genotype due to processes that may be intrinsic or extrinsic to the cell (17), and may include noise related to gene expression (18, 19), protein production (20), and DNA methylation (9, 10). Importantly, nongenetic fluctuations in phenotype can occur without a cell dividing (unlike changes in phenotype due to genetic mutations) and can be passed on to subsequent generations (20). Fluctuations in phenotype give rise to phenotype heterogeneity and can result in the transition of drug-sensitive cells into a robust phenotype (4, 10, 21). These cells are able to survive drug exposure, and hence the second mechanism, selection, can cause the most robust cells to become the majority in the population (14).

The dynamics of the average level of proliferation potential is even more interesting. Initially, as the majority phenotype transitions from PC9 to DTP, the average proliferation level drops. However, after a period of time, as DTEPs emerge from DTPs, the average level of proliferation increases back to normal. The interplay between nongenetic instability and selection could explain the emergence of the DTEPs from the DTPs. Nevertheless, these mechanisms alone cannot account for the initial drop in the average level of proliferation of the surviving cell population.

A natural explanation is that there are a small number of DTPs present in the initial population of PC9 cells that survive the initial drug therapy, and are subsequently selected for in the presence of the drugs (14). This is exactly the scenario suggested by the authors in ref. 12, despite the lack of definitive evidence that there are DTPs present in the initial population of PC9 cells. Therefore, it is possible that there are no DTPs present in the initial population of PC9 cells. If this is the case, then we propose that the DTP phenotype must be stress-induced (11, 15, 22).

In this adaptive scenario, PC9 cells are induced to lower their proliferative potential in response to stress caused by exposure to the drugs. This idea is supported by several observations. First, in response to stress caused by exposure to chemotherapeutic agents, damage to DNA is a common event and can activate the DNA damage response—a process that has been associated with cell-cycle arrest in lung cancer cells (22, 23)—until either the DNA can be repaired, or apoptosis is triggered. Second, one of the main chemotherapy drugs used in the experiments of ref. 12, gefitinib, has been shown to induce cell-cycle arrest in G_{1} phase in lung cancer cell lines, including PC9 cells (24). Therefore, we hypothesize that in the presence of the drugs, cells with a low survival potential will experience significant stress, and respond by lowering their proliferation rate (23, 25). On the other hand, robust cells, with a high survival potential, experience much less stress in the presence of the drugs and will not be induced to lower their proliferation rate.

Motivated by these considerations, here, we propose an individual-based (I-B) computational model (11, 26–28) and an integro-differential equation (IDE) model (29–31) of the phenotype evolution observed in ref. 12. Such models can be used as *in silico* laboratories to highlight stylized facts, and uncover mechanisms that underlie emergent features of cancer cell populations. The I-B computational model allows an intuitive and flexible description of the system at hand, whereas the IDE model makes it possible to study the system in terms of qualitative and asymptotic analysis, and is computationally less expensive.

Our models rely on the assumption that PC9 cells possess significant phenotypic plasticity. Furthermore, we assume that, during drug treatment, the surviving cell population can consist of cells residing within a spectrum of intermediate phenotypic states, ranging from PC9 through to DTP, and to DTEP. We also assume that the PC9s can gradually acquire drug tolerance at the same time as they gradually change their proliferation ability.

Unlike previous models of resistance in cancer cell populations that focus on the evolution of just one phenotypic trait (11, 27–31), here, we introduce a novel strategy to model the effects of stress-induced adaptation, and we focus on the evolution of two phenotypic traits that show substantial variability during drug treatment and after drug washout: a cell's survival potential and proliferation potential. Throughout, we consider the proliferation potential and survival potential of cells separately. However, mathematically they are not strictly independent traits, due to the strategy we use to model the effects of stress-induced adaption.

In the presence of high-dose drug therapy, we illustrate the ability of the models to recover the evolution of the PC9 cells into DTPs and finally into DTEPs, and also the drug resensitization of DTPs and DTEPs. Showing this, we can then argue that the models have validity for suggesting plausible evolutionary mechanisms underlying drug tolerance. Therefore, we use our I-B computational and IDE models in parallel to address the following questions:

**Q1.**Is nongenetic instability necessary in the development of drug tolerance?**Q2.**How do the evolutionary dynamics of the cell population compare between the regimes of low and high drug dose?**Q3.**Is it possible to determine whether therapy induces the DTP phenotype?**Q4.**What is the role of stem-like characteristics of the DTP cells in the development of drug tolerance?

Our study highlights the important role of nongenetic instability in the emergence of DTEPs, and also leads us to propose experiments that may help to assess whether there are DTPs present *ab initio*, and whether PC9 cells are undergoing stress-induced adaptation of their phenotype. Finally, our models predict that the transient dominance of the DTP subpopulation is regulated by their stem-like characteristics and is a direct result of the high doses of chemotherapeutic drugs used in ref. 12.

## Results

For each scenario presented below, we consider the concentration of cytotoxic drugs (i.e., the parameter *c*) to be high (low) if it is greater (less) than, or equal to, the LC_{90} (LC_{50}) dose, which is defined as the drug dose required to kill 90% (50%) of the total cell population in the initial stage of drug therapy, before the population starts to recover. A sample of parameter sets and their corresponding LC_{90} and LC_{50} doses is provided in Table 1. The end-time of the simulations is represented by the parameter *T*, and satisfies |${\frac {\left| N^\infty - N(t\, = \,T)\right|} {N^\infty}} \lt \,0.05$| or |${\frac{\left| {\varrho^\infty - \varrho(t\, = \,T)} \right|} {\varrho^\infty }}\, \lt \,0.05$|, where *N*^{∞} and |$\varrho^\infty$| are the respective asymptotic values of the total cell number and global population density, which are computed for a given set of parameters values.

(a_{1}, a_{2}, a_{3})
. | (b_{1}, b_{2}, b_{3})
. | β
. | |$\bar \upsilon$| . | x*
. | K
. | λ
. | LC_{90}
. | LC_{50}
. |
---|---|---|---|---|---|---|---|---|

(0.03, 0.25, 0.05) | (0.15, 1.3, 0) | 0.5 × 10^{−5} | 0.012 | 0.8 | 1.0 × 10^{5} | 0 | 0.138 | 0.090 |

(0.03, 0.25, 0.05) | (0.15, 1.3, 0) | 0.5 × 10^{−5} | 0.005 | 0.8 | 1.0 × 10^{5} | 0 | 0.189 | 0.114 |

(0.03, 0.25, 0.05) | (0.15, 1.3, 0) | 0.5 × 10^{−5} | 0 | N/A | 1.5 × 10^{5} | 0.05 | 0.213 | 0.121 |

(0.03, 0.25, 0.05) | (0.15, 1.3, 0) | 0.5 × 10^{−5} | 0 | N/A | 1.0 × 10^{5} | 0.02 | 0.212 | 0.121 |

(0.03, 0.25, 0.05) | (0.15, 1.3, 0) | 2.5 × 10^{−5} | 0 | N/A | 1.0 × 10^{5} | 0.02 | 0.220 | 0.122 |

(0.03, 0.25, 0.05) | (0.15, 1.3, 0.02) | 1.0 × 10^{−4} | 0 | N/A | 0.5 × 10^{5} | 0.05 | 0.215 | 0.115 |

(0.03, 0.25, 0.05) | (0.15, 1.3, 0.02) | 1.0 × 10^{−4} | 0.070 | 0.8 | 0.5 × 10^{5} | 0 | 0.100 | 0.067 |

(a_{1}, a_{2}, a_{3})
. | (b_{1}, b_{2}, b_{3})
. | β
. | |$\bar \upsilon$| . | x*
. | K
. | λ
. | LC_{90}
. | LC_{50}
. |
---|---|---|---|---|---|---|---|---|

(0.03, 0.25, 0.05) | (0.15, 1.3, 0) | 0.5 × 10^{−5} | 0.012 | 0.8 | 1.0 × 10^{5} | 0 | 0.138 | 0.090 |

(0.03, 0.25, 0.05) | (0.15, 1.3, 0) | 0.5 × 10^{−5} | 0.005 | 0.8 | 1.0 × 10^{5} | 0 | 0.189 | 0.114 |

(0.03, 0.25, 0.05) | (0.15, 1.3, 0) | 0.5 × 10^{−5} | 0 | N/A | 1.5 × 10^{5} | 0.05 | 0.213 | 0.121 |

(0.03, 0.25, 0.05) | (0.15, 1.3, 0) | 0.5 × 10^{−5} | 0 | N/A | 1.0 × 10^{5} | 0.02 | 0.212 | 0.121 |

(0.03, 0.25, 0.05) | (0.15, 1.3, 0) | 2.5 × 10^{−5} | 0 | N/A | 1.0 × 10^{5} | 0.02 | 0.220 | 0.122 |

(0.03, 0.25, 0.05) | (0.15, 1.3, 0.02) | 1.0 × 10^{−4} | 0 | N/A | 0.5 × 10^{5} | 0.05 | 0.215 | 0.115 |

(0.03, 0.25, 0.05) | (0.15, 1.3, 0.02) | 1.0 × 10^{−4} | 0.070 | 0.8 | 0.5 × 10^{5} | 0 | 0.100 | 0.067 |

NOTE: A sample of sets of parameter values used in the simulations, with their corresponding LC_{90} and LC_{50} doses, which were calculated numerically. *x** is the parameter that specifies the value of *x*, above which, stress-induced adaptation is switched off (details in the Supplementary Material).

### Mathematical models recover the population level dynamics of the PC9 cells

To determine whether our models can reproduce the experimental observations of the PC9 cancer cell lines during drug therapy, as reported in ref. 12, we numerically solve the IDE model and simulate the I-B model for a high dose of the cytotoxic drugs. In Fig. 2, we present numerical solutions from the I-B and IDE models for the case when there are a few DTPs present in the initial population (i.e., 5% of the global population density), and when there is no stress-induced adaptation of the proliferation potential. Figure 2A displays the trajectories of the mean phenotypic expression of the population. It is clear that, as time progresses, the mean phenotypic expression changes from being PC9-like, to being DTP-like and finally to being DTEP-like. Such dynamics are also shown in Fig. 2C–J, which displays the phenotypic distribution at different times during drug therapy (see also Supplementary Movie S1, which shows the phenotypic distribution evolving in time for one simulation of the I-B model). In Fig. 2B, the plot of the global population density as a function of time reveals that, shortly after drug exposure, there is a large reduction in cell number, followed by a time of relatively constant population level, before the population recovers back to normal levels. Similar results can be obtained when there are no DTPs present in the initial population, and when we include stress-induced adaptation of proliferation levels (we present numerical solutions for this case, from the IDE model in Fig. 3A and B and the top rows of Supplementary Figs. S3 and S4, whereas Supplementary Movie S2 shows the phenotypic distribution evolving in time for one simulation of the I-B model in this case).

However, only the I-B model can shed light on the individual cell-level dynamics, which, at times, is quite different from the population-level dynamics. A comparison between the trajectories of the phenotypic expression of six individual cells and the population mean trajectory from the I-B model is provided in Supplementary Fig. S5A. In this case, there are no DTPs present in the initial population, and we include stress-induced adaptation of the proliferation potential. Supplementary Fig. S5A shows one cell deviating from the average population-level behavior, by rapidly evolving a DTEP phenotype directly from the PC9 phenotype. Therefore, although the majority of surviving cells are predicted to follow the phenotypic trajectory from PC9 to DTP and through to DTEP (see Supplementary Movies S1–S3), due to stochasticity, this will not always be the case.

We present numerical solutions from the I-B and IDE models when *c* = 0, to see what happens to the PC9 cells if they are left to evolve in the absence of drugs in Supplementary Fig. S5C and S5D. It is clear that, although selection forces the majority of cells to remain in the PC9 phenotype, the effect of the diffusion term is to make the population more heterogeneous over time. Therefore, there are a small number of cells that are more slowly proliferating and/or with higher resistance, despite the fact that the fast-proliferating and less-resistant cells have a fitness advantage in the drug-free conditions.

Figure 2 also shows that we can obtain a good match for the dynamics predicted by the I-B model with that predicted by the IDE model. We also achieve the same qualitative match for the other cases considered in this section. Therefore, with the exception of our study on the stem cell characteristics of DTPs (where we use the extended I-B model), from now on we present numerical results obtained from the IDE model only.

In the experiments reported in ref. 12, both the resistant DTPs and DTEPs can be drug resensitized by drug-free passaging. Therefore, we numerically solve the IDE model and simulate the I-B model for a period of time with high doses of cytotoxic drugs, so that either DTPs or DTEPs are left as the dominant subpopulation. We subsequently switch the drug concentration off to simulate drug washout, and allow the system to continue to evolve. We find that the models can successfully reproduce drug resensitization of the DTPs and DTEPs (see Fig. 4 and Supplementary Fig. S6).

Finally, in ref. 12, the authors suggest that reversible drug tolerance in the PC9 cancer cell lines is due to an epigenetic mechanism, rather than resistance-conferring genetic mutations. To determine whether a genetic mechanism could generate the same phenotype dynamics observed in ref. 12, we define an IDE model of phenotype evolution, in which only selection and genetic mutations are driving cell dynamics (described in the Supplementary Material). Our simulations (see Supplementary Fig. S2) show that a detectable subpopulation of DTPs cannot emerge during treatment if there are no DTPs present in the system initially. In this case, either the PC9 population will move toward extinction, or the DTEPs will evolve directly from the PC9s. Furthermore, if there are some DTPs present in the system initially, and selection allows them to transiently emerge as the dominant subpopulation during drug treatment, then the time it takes for the DTEPs to subsequently emerge from the DTPs is significantly longer, than if phenotypic changes were the result of nongenetic mechanisms. These results support the idea that epigenetic mechanisms, rather than genetic mutations, are responsible for the evolution of phenotype observed in the cell lines studied in ref. 12.

### Nongenetic instability is crucial for the emergence of DTEPs

To investigate the role of nongenetic instability in the evolution of drug tolerance in PC9 cells, we analyze the phenotype evolution in the PC9 population for different values of the average rate of phenotypic fluctuations (i.e., the parameter *β*) for a high dose of cytotoxic drugs. We consider two hypothetical scenarios:

Scenario 1, when the initial population consists of only PC9s and stress-induced adaptation of the proliferation level is present.

Scenario 2, when a small number of DTPs are present initially and stress-induced adaptation of the proliferation level is absent.

For the first scenario, our model predicts that if the average rate of phenotypic fluctuations is too large, then the PC9s will evolve directly into the DTEPs. Furthermore, if this rate is too small, then the PC9s will become extinct before they are able to express the more robust DTP phenotype. This is illustrated in Fig. 3A and B, in which we plot the trajectories of the mean phenotypic expression of the population and the corresponding total cell density as a function of time, for three values of *β* (see also Supplementary Fig. S3A and S3B).

For the second scenario, the model predicts that in the absence of nongenetic instability, the DTEPs cannot emerge from the DTPs. Therefore, the DTP phenotype is stabilized inside the surviving population. Consequently, compared with the case when nongenetic instability is present, the long-term total global population density can be much smaller. In fact, if the concentration of the drugs is high enough, then the population can go extinct when nongenetic instability is absent. These predictions are illustrated in Fig. 3C and D (see also Supplementary Fig. S3C and S3D).

### A high dose of cytotoxic drugs is necessary for the transient dominance of DTPs

We use our IDE model to investigate the evolutionary dynamics of the PC9 cells in a low drug-dose regime. For scenario 1, described above, our model predicts that, *ceteris paribus*, decreasing the drug concentration can prevent the emergence of DTPs and speed up population recovery. These dynamics are illustrated in Fig. 5A and B, in which we plot the trajectories of the mean phenotypic expression of the population and the corresponding global population density as a function of time, for both a high dose and a low dose of the cytotoxic drugs. This is also true for scenario 2 (see Fig. 5C and D). Therefore, we can conclude that a higher dose of the cytotoxic drugs is a key ingredient to see the transient dominance of DTPs. A more complete set of comparisons between the high and low drug-dose regimes, for a range of parameter values, is provided in Supplementary Figs. S3 and S4.

### Low cytotoxic drug-dose experiments could establish if the DTP phenotype is not stress induced

Another interesting prediction from our model in the low-dose regime is that, if there are a few DTPs present in the initial PC9 population, and stress-induced adaptation is absent, then it is possible to see an increase in the average proliferative potential in the surviving population, followed by a decrease before a second increase (see Fig. 6). This behavior is not predicted in the high-dose regime, nor is it predicted if there are no DTPs present initially, and is driven by the concurrent emergence of the DTEPs from both of the initial subpopulations (PC9s and DTPs). In the low drug-dose regime and without stress-induced adaptation of the proliferative potential, some of the initial PC9 cells and their progeny are able to survive the drug exposure and gradually adopt a phenotype with a higher survival potential. Hence, some DTEPs begin to emerge directly from the PC9s. At the same time, the small DTP initial subpopulation, which has a significant survival advantage in the presence of the drug compared with the PC9s, is transitioning into a more proliferative (DTEP-like) phenotype and, after a time-lag, begins to outcompete the PC9s. Therefore, the majority population switches from the PC9s (which are transitioning into DTEPs) to the DTPs (also transitioning into DTEPs), thus forcing the average proliferative potential to decrease for a period of time before increasing again.

Therefore, one way to determine whether there are preexisting DTPs within the initial PC9 population, and if stress-induced adaptation is not occurring, would be to repeat the experiments of ref. 12 with a low dose of the cytotoxic drugs, and observe the average mitotic rate of the surviving population. Our models predict that, if there is an initial increase in the average mitotic rate, followed by a decrease and a further increase, then it is likely that DTPs are present in the initial population, and stress-induced adaptation is not occurring.

### Stem-like characteristics temporarily stabilize the subpopulation of DTPs

Sharma and colleagues (12) showed that the DTPs display stem cell–like surface markers, which are absent in the DTEPs. Therefore, it is reasonable to assume that the DTPs express some stem-like characteristics. In particular, stem cells have the potential to generate more stem cells, as well as differentiated daughter cells (36)—they are capable of symmetric cell divisions, in which a dividing stem cell can produce either only stem cell daughters (self-renewal) or only differentiated daughters, or asymmetric cell divisions, in which a dividing stem cell produces one stem cell daughter and one differentiated daughter.

Motivated by these considerations, we use the extended I-B model to investigate the role of self-renewal and asymmetric cell divisions in the development of drug tolerance. We analyze the dynamics of the cells for a range of different values of the symmetric self-renewal probability *p*_{1}, when the probability of asymmetric self-renewal *p*_{2} equates to 1 − *p*_{1}, so that the probability of symmetric differentiation *p*_{3} is zero (see Supplementary Fig. S1C and S1D).

Again, we consider the two hypothetical scenarios outlined above. In both scenarios, if DTPs do transiently dominate the surviving population, then a higher probability of self-renewal increases the time taken for the emergence of the DTEPs from the DTPs. Furthermore, after the DTEPs have emerged, a higher probability of self-renewal corresponds to a larger and more stable subpopulation of DTPs in the surviving population. These predictions are illustrated in Fig. 7, in which, for scenario 1, we plot the trajectories of the mean phenotypic expression of the population and the corresponding numbers of PC9s, DTPs, and DTEPs as a function of time for three different values of *p*_{1} (see also Supplementary Fig. S7C and S7D for scenario 2).

On the other hand, if the DTPs do not transiently dominate the surviving population, either because the population becomes extinct before it can become drug tolerant (see Supplementary Fig. S7G and S7H), or because the DTEPs emerge directly from the PC9s (see Supplementary Fig. S7A and S7B), then, not surprisingly, the probability of self-renewal *p*_{1} does not play any role in the dynamics of the population.

## Discussion

The PC9 cell line experiments reported in ref. 12 are performed in an isolated and relatively homogeneous environment, and involve only a few constituents—a genetically identical cell population, culture media, and combinations of drugs. Furthermore, each experiment has clear observables, namely the percentage of surviving cells and their phenotypic distributions, which make the PC9 cancer cell lines an ideal system to study from a mathematical perspective.

Here, we have presented an I-B computational model and an IDE model of the evolution of phenotype observed in ref. 12, which rely on the assumption that this evolutionary process is driven by the interplay between nongenetic phenotype instability, stress-induced adaptation, and selection. Distinct mathematical strategies were used to model each evolutionary mechanism, to enable us to better understand the effect of each mechanism on phenotypic evolution. The models reproduce the main experimental observations detailed in ref. 12, and support the idea that epigenetic mechanisms, rather than genetic mutations, are responsible for the evolution of cell phenotype. Our analysis highlights the important role of nongenetic fluctuations in phenotype in the emergence of drug tolerance in PC9 cancer cell lines. In particular, we suggest that the absence of nongenetic instability can result in the stabilization of the DTP phenotype in the surviving population, so that DTEPs do not emerge, or even in extinction. This is a key result because it supports the idea that epigenetic therapy may be a promising therapeutic strategy in the war against cancer (37–39).

Another important prediction of our models is that the transient dominance of DTPs is strictly related to the use of high doses of cytotoxic drugs. If experimentalists apply a lower dose of cytotoxic agents to the PC9 cell population during drug therapy, we propose that it would be highly unlikely to observe DTPs. Rather, we would expect the DTEPs to emerge directly from the PC9 population. Note that this is the usual way to yield stable drug-tolerant lineages (40).

Our analyses also suggest that stem-like, self-renewing cell divisions of DTPs can act to increase the number of DTPs in the system, and is consistent with other mathematical models of stem cell dynamics (41). In particular, we expect that DTP self-renewal can act to stabilize the DTP subpopulation for a longer period of time during drug therapy, compared with the case without self-renewal, as well as to delay the emergence of DTEPs.

On the basis of our models, we can conclude that if there are no DTPs present in the initial population of PC9 cells, then it is likely that a proper interplay between nongenetic phenotype instability, stress-induced adaptation, and selection is mandatory for the transient appearance of the DTP phenotype during high-dose drug therapy. On the other hand, if there are some DTPs present in the initial population, then nongenetic fluctuations in phenotype and selection are enough to explain the experimental observations reported in ref. 12. Therefore, the next biologically meaningful question is: are DTPs present in the initial population? Our analysis enabled us to propose a low cytotoxic drug-dose experiment, which could answer this question. Another possible way to determine whether DTPs are present *ab initio* would be to use the experimental technique described in ref. 42 to identify and isolate slow-cycling cells in the parental PC9 cell population, and subsequently compare the drug sensitivity of these cells with that of the more rapidly dividing PC9 cells.

The low-dose experiment we proposed could also determine that stress does not induce the DTP phenotype. However, another way to assess whether PC9s are induced by stress to acquire the DTP phenotype would be to perform a flow-cytometry analysis on the PC9s with respect to the cluster of differentiation 24 (CD24) surface marker, at a time soon after drug delivery. If the surviving population has an intermediate level of expression soon after drug exposure, relative to the expression of CD24 in PC9s and DTPs, then this would suggest that the DTP phenotype is a result of nongenetic fluctuations in phenotype, selection, and stress-induced adaptation.

To keep the models as simple as possible, we chose to include only mechanisms that were necessary to reproduce the experimental results observed in ref. 12. However, if appropriate future experimental evidence emerges for this system, such as time series data for proliferation rates or cell density, then this could be used to determine whether an additional mechanism, like stress-induced adaptation of the survival potential, should be added to the basic framework established here.

A natural way to extend our model would be to study the emergence of epigenetic drug tolerance in a primary tumor. In this framework, we could assume cells to be organized in a nonvascularized and radially symmetric microspheroid (43), and introduce an additional structuring variable that stands for the normalized linear distance of cells from the center of the spheroid. We could also introduce an additional evolution equation for the local concentration of cytotoxic drugs. In this case, the interplay between diffusion and consumption of cytotoxic drugs could lead to the creation of distinct niches differentiated by the local environment, which would provide ecologic opportunities for diversification.

In our study, we have considered a two-dimensional structuring phenotype. It is then natural to wonder about what could be a biologic quantitative variable underlying such a continuous phenotype. We propose, as a reasonable candidate, the degree of epigenetic modifications of the DNA (methylation and histone acetylation) in relevant parts of the genome. Such modeling has been previously proposed in ref. 44, using a one-dimensional, phenotype-like structure variable. By using a multidimensional phenotype, we consider differential effects of epimutations on genes responsible for the classical intracellular pathways of proliferation and survival.

## Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

## Authors' Contributions

**Conception and design:** R.H. Chisholm, T. Lorenzi, A. Lorz, L.N. de Almeida, A. Escargueil, J. Clairambault

**Development of methodology:** R.H. Chisholm, T. Lorenzi

**Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.):** R.H. Chisholm, T. Lorenzi

**Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis):** R.H. Chisholm, T. Lorenzi, A. Lorz, A.K. Larsen, L.N. de Almeida, A. Escargueil, J. Clairambault

**Writing, review, and/or revision of the manuscript:** R.H. Chisholm, T. Lorenzi, A. Lorz, A.K. Larsen, L.N. de Almeida, A. Escargueil, J. Clairambault

**Study supervision:** A. Lorz, L.N. de Almeida, A. Escargueil, J. Clairambault

## Acknowledgments

The authors thank both reviewers for their most helpful and insightful comments on the article.

## Grant Support

This research was funded in part by the French “ANR blanche” project Kibord: ANR-13-BS01-0004. T. Lorenzi was supported by the Fondation Sciences Mathematiques de Paris (FSMP), by a public grant overseen by the French National Research Agency (ANR) as part of the “Investissements d'Avenir” program (reference: ANR-10-LABX-0098), and by the labex LMH through the grant no. ANR-11-LABX-0056-LMH in the “Programme des Investissements d'Avenir”.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked *advertisement* in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

## References

^{+}cells in lung adenocarcinoma by activating notch signaling