A model for the intrinsic limit of cancer therapy: Duality of treatment-induced cell death and treatment-induced stemness.

Intratumor cellular heterogeneity and non-genetic cell plasticity in tumors pose a recently recognized challenge to cancer treatment. Because of the dispersion of initial cell states within a clonal tumor cell population, a perturbation imparted by a cytocidal drug only kills a fraction of cells. Due to dynamic instability of cellular states the cells not killed are pushed by the treatment into a variety of functional states, including a "stem-like state" that confers resistance to treatment and regenerative capacity. This immanent stress-induced stemness competes against cell death in response to the same perturbation and may explain the near-inevitable recurrence after any treatment. This double-edged-sword mechanism of treatment complements the selection of preexisting resistant cells in explaining post-treatment progression. Unlike selection, the induction of a resistant state has not been systematically analyzed as an immanent cause of relapse. Here, we present a generic elementary model and analytical examination of this intrinsic limitation to therapy. We show how the relative proclivity towards cell death versus transition into a stem-like state, as a function of drug dose, establishes either a window of opportunity for containing tumors or the inevitability of progression following therapy. The model considers measurable cell behaviors independent of specific molecular pathways and provides a new theoretical framework for optimizing therapy dosing and scheduling as cancer treatment paradigms move from "maximal tolerated dose," which may promote therapy induced-stemness, to repeated "minimally effective doses" (as in adaptive therapies), which contain the tumor and avoid therapy-induced progression.

Introduction

The single major cause of treatment failure in cancer therapy is the emergence of a treatment-resistant tumor that drives recurrence. Other than in the case of some early-stage tumors, it is tacitly accepted that relapse is inevitable during the course of drug treatment. This assumption has translated into the unquestioned practice of quantifying the efficacy of all treatments by how long one extends the time period between diagnosis and either progression, in the form of a clinical relapse, or death [1]. The former metric defines a progression-free survival (PFS) time, which quantifies not the prevention of progression, but a delay, as evident in the proportional hazard model [2]. Treatment success is therefore measured as a “shift to the right” of the decaying Kaplan-Meier curve, which represents the fraction of progression-free surviving patients in the treated cohort compared to that of the control group. The hence derived extension of median PFS time (or loosely equivalent, of the median time to progression, or TTP) has become a de facto measure for success of a therapy [3].

In theoretical and in vitro experimental models of post-treatment regrowth of a tumor cell population, the “recovery time” in which the surviving tumor cells regrow to reach the population threshold present at the onset of treatment is a biological characteristic of the tissue. It can be considered the equivalent of the clinical TTP [4].

Recurrence, or tumor cell regrowth after treatment, is generally thought to be the result of selection in a process of Darwinian somatic evolution: Given sufficient genetic variability in a sufficiently large initial (pretreatment) cell population, it is considered statistically certain (possibly as a result of increased mutation rate in cancer cells) that the population contains cells carrying genomic mutations that confer drug-resistance and stem-like traits [5-8]. A single cell with such a mutation will survive the treatment and clonally expand, thus driving the tumor regrowth under treatment.

This genetic explanation implicitly acknowledges de facto inevitability of relapse for a certain set of parameters including mutation rate, cell population size and selection pressure [5, 6, 8, 9]. In addition, elimination of drug-sensitive cells by treatment releases drug-resistant cells from spatial and nutritional constraints and facilitates their proliferation, thereby creating an apparent causal link between treatment and recurrence [4, 10–13].

In recent years, non-genetic cell phenotype plasticity and the resulting cell population heterogeneity has been recognized as a source of the resistant cell phenotype, which could underlie recurrence without implication of genetic mutations [8, 14–19]. Extensive phenotypic heterogeneity within a population of cells is generated by the variability of an individual cell’s biochemical state. Such non-genetic variability emanates, in part, from the ability of the cell’s gene regulatory network to produce a diversity of stable gene expression patterns (attractors), resulting in multistability within a single clonal, isogenic population. Gene expression noise stochastically disperses individual cells in gene expression space, allowing them to occupy a range of these stable cell states. Because such stochasticity of gene expression causes continuous phenotype switching and equilibration of phenotypes, this type of heterogeneity is not subject to the rule of extinction of a phenotype, as is the case with genetic mutation.

The resulting phenotypic diversity of the isogenic cell population is, while also probabilistic, more prevalent than that generated by genetic mutations and it produces distinct, robust, recurrent, and biologically relevant phenotypic sub-states in clonal cell populations [20-22]. Such sub-states include mesenchymal, persister, or stem-like states, etc., as single-cell RNAseq now amply reveals [19, 23–30]. Some of these states can confer resistance and are sufficiently robust to be inherited across several cell generations [22, 31]. In this manner, non-genetic probabilistic variation can also drive Darwinian selection of resistant cells, at least for a number of cell generations.

While both genetic and non-genetic variation arise in a probabilistic manner, there is a key difference. Because non-genetic variant attractor states are the result of regulatory mechanisms, they can also be directly induced by environmental signals. Such instruction to change gene expression programs in a directed manner by an external input via a deterministic (or strongly biased probabilistic) control, as opposed to selection of an undirected probabilistic internal change, plays a dominant role in a tumor’s acquisition of stem-like drug-resistant cells at short time scales (days) compared to the clonal expansion of rare mutant clones [9]. Such regulated change of cell state may be part of a robust, evolved cellular defense mechanism against cellular toxins [32].

A growing body of evidence suggests that emergence of stem-like and therapy-resistant cells along with the associated changes in gene expression are induced (as opposed to selected) by the cytotoxic stress of treatment [29, 33, 34]. In other words, there is a causal biological mechanism linking treatment to stress to the stem-like phenotype. The recurring appearance of stem-cell-like gene expression programs, or “stemness,” the speed of response and involvement of canonical signaling pathways that confer multidrug resistance (such as Wnt signaling-mediated upregulation of the ABC membrane pumps) collectively support the idea of stress-induced activation of preexisting xenobiotic resistance programs in cells by treatment [35-37].

Therefore, any cytocidal treatment may be a double-edge sword: while a one fraction of cells in the heterogeneous population is killed, another fraction of cells is induced by treatment stress to enter a stem-like or more resistant persister state—planting the seed for recurrence [8, 38]. Drug-induced resistance thus poses an intrinsic limit to curability of tumors under any treatment that involves cell stress, as is the case with chemotherapy, target-selected therapy or radiation. The role of somatic Darwinian evolution in recurrence relative to that of non-genetic cell state transitions has been analyzed using mathematical models in order to minimize selection pressure during treatment [4, 5, 39–41]. However, the intrinsic limit that induced resistance places on therapy has only recently been systematically evaluated [39].

Here we analyze a most elementary, generic, mathematical model for the conflicting processes that are inherent to cancer therapy: treatment-induced cell death and treatment-induced transition from a drug-sensitive phenotype to a drug-resistant (stem-like) one. Under this formulation, recurrence is “wired-into” the population dynamics and, depending on quantitative details, can become manifest. Over a relevant parameter space of an ordinary differential equation (ODE) model, we analytically evaluate the activity profiles of a drug in inducing cell death vs. transition to the resistant state. We quantify how these features of treatment relate to the intrinsic inevitability of recurrence, measured as TTP. We thus provide a formal survey of the consequence of the core process of non-genetic induction of resistance by treatment, irrespective of the ensuing selection and micro-environmental influences.

We find that depending on the relative susceptibilities for cell death versus induction of the resistant state there can be a non-monotonic dependence of TTP on drug dose, such that beyond a narrow optimal dose, the induction of resistance dominates and increasing treatment intensity will drastically shorten TTP. Thus, knowledge of the measurable propensities of cells to die or activate resistance mechanisms as a function of dose is critical information for optimizing therapy. Our focus on treatment-induced non-genetic tumor cell state change complements the evolutionary framework, fills a conceptual gap to help explain why it is so hard to cure advanced cancer and can be used for modeling scheduling to avoid treatment-associated progression as sought by adaptive therapy [11].

Materials and methods

Dynamical model of tumor growth

We consider an ODE model that describes the population dynamics of cancer cells that interconvert between two distinct cell states: drug-sensitive and drug-resistant. Let x(t) = [x1(t), x2(t)] denote the population state vector, where x1(t) is the number of sensitive cells, and x2(t) is the number of resistant cells. Following the formulation given by Zhou et al., the sizes of these two populations are governed by the following linear ODE (Fig 1) [42]:

Fig 1

Levels of parameterization in the dynamical model of tumor growth.

At the “highest” (i.e., most general) level, there are the rate constants that govern the growth and state transition dynamics of the cancer cell population. One level down, we introduce drug treatment into the model by assuming that the rate constants d, d, and k are logistic functions of drug dose m. The parameters that determine the shape of these dose-response curves—the pharmacodynamic constants, e.g., EC50—form the “lowest” (i.e., most specific) level of the model parameters.

We analyze this model of ongoing treatment with the assumption that treatment acts by raising the per capita death rate of cancer cells. Herein, the parameter b (b) denotes the fixed division rate of sensitive (resistant) cells, whereas d (d) denotes the total death of sensitive (resistant) cells undergoing drug treatment. The parameter k denotes the transition rate at which sensitive cells become resistant, and k is the transition rate at which resistant cells become sensitive. All of these parameters are strictly positive and depend in a particular way from the drug dose as discussed below.

We use TTP—the time it takes for the tumor to surpass its pretreatment population size—as a quantitative measure of recurrence under drug treatment. In this model, we denote TTP by t, which is defined as where N(t) = x1(t) + x2(t) is the total number of cells in the tumor (Fig 1). The goal of this model is to understand how t depends on the model parameters.

Stability analysis

Aside from the degenerate case in which one eigenvalue of A is zero, the only fixed point of the ODE in Eq 1 is the point x = 0: extinction of all cancer cells. The eigenvalues λ1,2 of A determine the local stability of this fixed point, whereas the eigenvectors v(1,2) of A give the primary directions of growth and/or decay in state space. For this model, we can derive general conditions on the growth and transition rates under which the origin is an unstable node, a stable node, or a saddle point (Table A in S1 Text). Under appropriate initial conditions, these three cases correspond to unchecked tumor growth, tumor extinction, and tumor regression followed by regrowth. The last scenario provides the simplest mathematical conception for the relapse phenomenon.

The saddle point nature of the extinction state suggests that the dynamics are ultimately difficult to control and contain. It also points to a marked “turning point” for each trajectory: if the initial tumor population is “close enough” to the stable manifold, the trajectory will first move towards the origin before being repelled away along the unstable manifold, indicating the difficulty to eradicate the tumor (S1–S4 Figs). We can think of this behavior as temporary control of tumor size before the inevitable progression. S3 Eq in S1 Text tells us that it is possible to control tumor regrowth wherever the relative fitness of the sensitive phenotype is sufficiently large.

The extinction state is a saddle point whenever the ratio (b − d)/(b − d) is greater than the ratio (S3 Eq in S1 Text). That is, temporary control of the tumor size is possible when the relative fitness of the sensitive phenotype is sufficiently large. The constant b − d − k is the net flux of the drug-sensitive population per unit density of drug-sensitive cells, and the constant k is the flux into the drug-sensitive population per unit density of drug-resistant cells. The latter parameter, sometimes referred to as the “backflow rate,” is useful in characterizing how the pool of drug-resistant cells allows the drug-sensitive population to avoid extinction during chemotherapy [42].

Time to progression

The behavior of TTP as a function of drug dose is closely related to the saddle point dynamics of the system. For one, to have t > 0, the total population must decrease initially before recovering to its initial value, which is typically only observed when the origin is a saddle point. Moreover, under the saddle point dynamics, t = t is the unique non-zero time point for which N(t) = N(0).

Even with a closed expression for N(t), solving the expression N(t) = N(0) for t is not generally possible as it involves the sum of distinct exponential terms. Instead, we can approximate t by decomposing the saddle point dynamics into the distinct stages of population decrease (remission) and increase (regrowth). By definition of TTP, the total population N(t) is undergoing regrowth at time t = t, i.e., N′(t) > 0. The theory of linear dynamical systems tells us that the exact solution N(t) is given by a linear combination of the exponential terms and , where λ1,2 are the eigenvalues of A (S5 Eq in S1 Text). In the case of a saddle point, where λ1 > 0 and λ2 < 0, the tumor population regrowth is necessarily driven by the exponential term corresponding to the positive eigenvalue λ1. Therefore, at time t, we may assume that the total population is well approximated by this exponential term: . Under this assumption, we can approximate TTP as follows (S8 Eq in S1 Text):

Pharmacodynamic model of continuous therapy

To consider drug treatment, we assume that drug dose is constant throughout the course of therapy (i.e., continuous therapy) and that the rate constants d, d, and k depend on the amount of drug m present in the system, i.e., that chemotherapy reduces tumor burden by increasing the death rate of cancer cells. At the same time, it increases the rate at which sensitive cells become resistant, the basis for the “double-edged sword” effect of chemotherapy. In mathematical terms, we introduce a secondary parameter m that denotes the drug dose, and we assume that the primary parameters d, d, and k that capture tumor cell population dynamics are increasing functions of m (Fig 1). Scaling m as percentage of the maximum tolerated dose (MTD), i.e., 0 ⩽ m ⩽ 100, we next discuss the pharmacodynamics of cellular drug response, i.e., functional forms for the dependence of the three primary parameters from m.

Using the commonly observed sigmoid shape of biological response curves, which reflect the cumulative probabilistic response of individual cells in a heterogeneous population, we use logistic functions to describe the rate constants d(m), d(m), and k(m) (Fig 1): We assume that k, the rate constant for the re-sensitization of resistant cells, does not change with drug dose.

In the case of d(m), each of the four parameters δ, E, r and P determines the shape of the logistic curve and describes a behavior affected by the drug. For this reason, we refer to these parameters as the pharmacodynamic parameters associated with a given rate constant for cell response. In particular, the parameters E, P, and r respectively determine the maximum response (efficacy), EC50 (which is inversely related to potency), and saturation rates for each of the above responses of the drug. For example, a drug with high efficacy and high potency to kill sensitive cells is characterized by a high value for E and a low value for P.

Taken together, our model incorporates three levels of parameterization: first, the parameters of the general linear population dynamics model are the rate constants b, b, etc. (Fig 1). Second, some of these rate constants (e.g., k) depend on drug dose, which is represented by the parameter m (Fig 1). And third, the way each of these rate constants depend on drug dose is determined by their respective pharmacodynamic constants, as given in Eq 4 (Fig 1).

In drug development implicit pharmacodynamic parameters are empirically tuned to optimize treatment outcome, which, in the case of our model, is t. Teasing apart how the drug dose m affects TTP as a function of the pharmacodynamic parameters is not a straightforward problem because we do not have an explicit expression for t. An exhaustive search of the entire parameter space is impractical. Because drug-induced transition to the resistant state is a hitherto unaccounted-for factor affecting recurrence, we start with a cursory analysis of how t changes when we vary the dose-response relationship of the sensitive-to-resistant transition rate k(m).

Results

Qualitative sensitivity analysis

Representative samples of the parameter space demonstrate the possible qualitatively distinct treatment outcomes—t as a function of drug dose—that depend on the pharmacodynamics (i.e., rate of behavior as function of drug dose m) of drug-induced resistance relative to killing (Fig 2). We define “drug resistance” by taking E ≪ E: the drug’s cell-killing efficacy is lower for the drug-resistant phenotype than it is for the drug-sensitive phenotype. Although drug-resistance can be manifest in dampening the killing effect either by decreasing the drug potency (or equivalently, increasing the EC50) or efficacy to kill cells (or a combination of both), we do not vary this aspect of resistance. Instead, we model resistance as a lowering the resistant cell death rate at MTD, which reflects a more profound effect on the cell state. To incorporate the observed fitness cost of resistance in the absence of drug, we further assume that b > b and δ = δ [4, 10, 43, 44].

Fig 2

A coarse-grained sensitivity analysis of the effect of drug-induced resistance on treatment outcome.

The parameters given in S1 Table are held fixed, while the EC50 and efficacy of resistance induction are varied in the above four cases. Efficacy increases as the parameter E increases, whereas potency decreases as the EC50 value (P) increases. A. Case A: high efficacy, low potency (E = 0.2, P = 50). B. Case B: high efficacy, high potency (E = 0.2, P = 10). C. Case C: low efficacy, low potency (E = 0.02, P = 50). D. Case D: low efficacy, high potency (E = 0.02, P = 10).

Thus, we anchor the pharmacodynamics for killing of sensitive and resistant cells d(m) and d(m) and vary the parameters that determine drug-induced transition to the resistant state with respect to either EC50 (P) or efficacy (E) for the transition rate k(m) (S1 Table, Fig 2). We consider these two parameters specifically because of the different effects that they exert on the dose-response curve of a given drug. The efficacy of a drug, or the drug response at MTD, is an intrinsic property of the drug and cannot be compensated by alteringg the drug dose. On the other hand, the potency of a drug, which corresponse to the inverse of the EC50, is a property that can be compensated by altering drug dose: a low-potency drug can achieve the same response as a high-potency drug, but at a higher drug dose.

The unknown then is how the relation between rate of induction of the resistant phenotype, k(m), and the “kill curves,” d(m) and d(m), shape t as a function of drug dose m. The coarse-grained, but comprehensive, sensitivity analysis of our model model is achieved by altering the values of P and E, relative to the two fixed kill curves pharmacodynamics curves (Fig 2). The following four canonical cases represent qualitatively distinct, plausible pharmacodynamical relationships between killing and induction of resistance due to treatment.

Case A: High efficacy, low potency of resistance induction

We first consider a drug with a high efficacy and low potency (i.e., high E and P) for inducing resistance (Fig 2A). For such a drug, t is not monotone increasing in drug dose. Instead, it first increases before decreasing to a plateau (Fig 3A). Thus, increasing drug dose past a certain point significantly worsens the treatment outcome.

Fig 3

Summary of model behavior for Cases A, B, C, and D as a function of drug dose m.

The net growth rates of sensitive and resistant cells are defined as g ≔ b − d and g ≔ b − d, respectively. First row: Total growth and transition rates, g + g and k + k, respectively. Second row: Eigenvalues of the matrix A in Eq 1. The inset in cases A and C highlight the non-monotonic behavior of the eigenvalue λ1. Third row: Time to progression t plotted alongside its asymptotic approximation .

We observe that t increases with drug dose when an increase in dose corresponds to a significant decrease in the total growth rate relative to the total switching rate (Fig 3A). On the other hand, t decreases with drug dose when an increase in drug dose corresponds to a significant increase in the total switching rate relative to the total growth rate (Fig 3A). The fact that the transition dynamics dominate the growth dynamics at high drug doses indicates that under cytotoxic stress, tumor recurrence is not driven by cell growth (Fig 3). Instead, it is driven by the ability of the tumor cells to evade extinction by transitioning to a stem-like, drug-resistant state.

Returning to our analysis of the saddle point dynamics, the positive eigenvalue λ1 is also non-monotonic in drug dose, first decreasing to a minimum before increasing to a plateau (Fig 3). This inverse relationship between λ1 and t agrees with the above observation that λ1 roughly corresponds to the rate of tumor growth during recurrence.

Case B: High efficacy, high potency of resistance induction

We now consider a drug with the same efficacy in inducing the stem-like state as in Case A but with higher potency (Fig 2B). Unlike Case A, t is now monotone increasing in drug dose (Fig 3B). Therefore, if the drug dose is sufficiently high, the treatment outcome in Case B is not sensitive to fluctuations about a dosage. This robustness is favorable to the sensitivity observed under Case A, in which a small fluctuation in drug dose can significantly worsen the t. On the other hand, the maximum t in Case B is low compared to that in Case A (Fig 3A and 3B). A treatment that is robust to variation in drug dose is not necessarily a good one if progression is only, at most, slightly delayed.

As before, the behavior of t is well summarized by the positive eigenvalue λ1 and the difference between the total growth and switching rates (Fig 3B). In particular, the latter shows how increased drug potency to induce stemness affects t. As EC50 for k is closer to that for d and d than it is in Case A, the total growth and state-switching rates saturate at roughly the same drug dose (Fig 3B). Therefore, the difference between the two rates only increases over a small range of doses, in which in t increases monotonically (Fig 3B). Cases A and B reveal that the relative potency of a given drug to induce resistance versus kill cells determines whether t is monotone increasing in drug dose or not.

Case C: Low efficacy, low potency of resistance induction

For a drug with the same low potency as in Case A, but with lower efficacy, t is again, as in Case A, a non-monotonic function of drug dose (Figs 2C and 3C). However, the difference between the maximum possible t and t at MTD is much smaller for the drug than in Case A, since the maximum difference between the total growth and switching rates is smaller than in Case A (Fig 3A and 3C). In terms of treatment outcomes, we can think of this drug as an improvement from Cases A and B: t is overall less sensitive to variation in drug dose than in Case A, and the maximum possible t is greater than that in both Cases A and B.

The overall increase in t is also reflected in saddle point dynamics of the system. In Case C, the magnitude of the negative eigenvalue λ2 is overall lower than in Case A (Fig 3C). This decrease in magnitude agrees with the increase in t, as λ2 roughly corresponds to the rate of tumor remission in the early stages of treatment. That is, the remission stage is prolonged under treatment by a low-efficacy drug, as compared to a high-efficacy drug (S1–S4 Figs).

Case D: Low efficacy, high potency of resistance induction

Finally, we consider a drug that is a “combination” of Cases B and C; that is, the drug has low efficacy but high potency in inducing resistance to cell killing (Fig 2D). Under this drug, t is a similar “combination” of Cases B and C: it is monotone increasing in drug dose, and its value at MTD is the same as in Case C (Fig 3D). Therefore, a drug with low efficacy and high potency for drug-induced resistance relative to killing produces a treatment scheme that effectively delays tumor progression for a wide range of drug doses. Compared to case C, the higher potency in inducing stemness at low dose abrogates the optimal dose window, i.e., the counterintuitive peak in t at lower dose.

Parameter search

In order to determine how well Cases A, B, C, and D capture the dependence of t on drug potency (P) and efficacy (E) of resistance induction, we perform a broad parameter search. Specifically, we compute t over a wide range of P and E beyond the four representative cases from the previous section (Fig 2). Instead of plotting t as a function of drug dose m, we plot three quantities of interest as a function of P and E: the drug dose at which t attains its maximum, the maximum value of t, and the value of t at MTD (Fig 4). Doing so, we find that t attains its maximum at MTD (i.e., is monotonic in drug dose) when P is low regardless of the value of E (Fig 4A). Therefore, t is indeed a non-monotonic function of drug dose only when the potency to induce resistance is low. Furthermore, the maximum value of t shows little sensitivity to either E or P, and instead appears to depend on whether t is monotone increasing in drug dose or not (Fig 4B). We also find that the value of t at MTD decreases as E increases independent of the value of P (Fig 4C).

Fig 4

Analysis of time to progression (t) over a wide range of values for potency (P) and efficacy (E) of inducing resistance.

Plotted values are indicated by the colorbar on each plot. All other parameters are fixed at the values given in S1 Table. A. Drug dose m at which t attains its maximum. B. Maximum value of t over all drug doses. C. Value of t at MTD.

The results of this finer-grained parameter search indicate that the four cases given in the main text indeed characterize the dependence of t on the parameters P and E: the non-monotonic dependence of t on drug dose is governed by P, whereas the value of t at MTD is governed by E. This parameter search, however, is still limited to the two parameters P and E. Taken together, a key finding is that at low potency (high EC50) for inducing resistance, where higher drug doses m are needed to trigger this cell state transition, increasing the dose past a certain point reduces t. A further evaluation of the parameter space, paired with rigorous sensitivity analysis, is required to characterize how t depends on the complete set of model parameters.

Virtual cohort simulations

In statistical analysis of clinical studies, individual patient measures, such as PFS time and TTP, are typically not displayed directly; instead, the data are often presented in the form of Kaplan-Meier curves, which show the fraction of surviving and progression-free patients as a function of time [2, 3]. The process of calibrating mathematical models of cancer dynamics to clinical data typically involves applying regression, or some other parameter-fitting method, to a Kaplan-Meier curve from observed cancer patient cohorts [45, 46]. Doing so requires generating a “virtual patient cohort” from the model, usually by assuming some statistical distribution of a set of the model parameters, and sampling individual “patients” from this distribution [45, 46]. As a step towards grounding our model in clinical data, so as to make meaningful predictions about treatment courses, we generate virtual patient cohorts and present an analysis of the resulting Kaplan-Meier curves.

To generate virtual patient cohorts, we assume that the basal rates of cell birth, death, and state transition, as well as the fraction of resistant cells at tumor detection, vary from patient to patient (S2 Table). For simplicity, we assume that each of these parameters are uniform random variables. One possible way of updating this assumed prior with a more appropriate posterior distribution based on clinical data would be to use an expectation-maximization approach [46]. We assume that the remaining parameters, tumor size at diagnosis and the pharmacodynamic constants, remain fixed across all patients, as they are intrinsic properties of standard clinical practice and a given drug, respectively (S1 Table).

Keeping with our previous analysis, we generate n = 100 virtual patient cohorts of n = 103 patients each, and compute the TTP curve t(m) for each patient at values of E and P given by the four previous cases A, B, C and D (Fig 2). Once we have t(m) for each patient in a given cohort, we can compute the progression-free fraction for a range of possible drug doses, where we treat t as a progression event for each patient. We then take the mean progression-free fraction across all cohorts, giving us a set of Kaplan-Meier curves (Fig 5). Our aim is to see if the non-monotonic dependence of treatment outcome on drug dose m is manifest at the ensemble level of the patient cohort. Whereas for an individual patient, an improved treatment outcome is indicated by an increased TTP, an improved outcome for an entire cohort is indicated by an upward and/or rightward shift in the Kaplan-Meier curve.

Fig 5

Kaplan-Meier curves averaged over n = 100 virtual cohort simulations of n = 103 patients each, plotted for fixed drug doses m.

Error bars denote plus and minus one standard deviation. The x-axis denotes time t, and the y-axis denotes the fraction of patients in a cohort who are progression free by time t (i.e., t > t). Each curve is colored according to the fixed drug dose m applied to each cohort, indicated by the color bar in each plot. Panels A, B, C, and D correspond to varying the parameters P and E according to the four cases shown in Fig 2. Unless specified as being drawn from the uniform distributions given in S2 Table, all other parameters are held fixed at the values given in S1 Table.

Comparing the resulting Kaplan-Meier curves for cases A-D, we find that the mean behavior of the survival fraction mirrors that of t. As before, in cases A and C, where potency to induce resistance is low, the overall treatment outcome is non-monotonic in drug dose: as m increases from 0 to MTD, the Kaplan-Meier curve shift upwards, then downwards (Fig 5). Also in line with our TTP analysis is the observation that in cases C and D, where efficacy to induce resistance is low, the overall treatment outcome is better than that for cases A or B: the upwards shift of the Kaplan-Meier curves going from low to high drug doses is much larger in the former cases than in the latter (Fig 5). Thus, our conclusions about the qualitative behavior of the model at the level of the individual “patient’s” tumor cell population scale up to the level of the statistical ensemble (i.e., the cohort), replicating both the “expected” as well as the “counterintuive” effects of increasing drug dose. This result is significant as it demonstrates that our analysis is in some sense robust to noise. The next logical step would be to see if our analysis still holds, and is therefore clinically meaningful, after fitting our parameters to actual patient data.

Discussion and conclusion

In advanced tumors, post-treatment recurrence is almost an intrinsic feature in the course of tumor progression. It is increasingly acknowledged that even if an untreated tumor would result in rapid progression and death, recurrence after treatment is causatively or mechanistically linked to the act of treatment. The traditional explanation for recurrence after initial remission invokes selection of preexisting mutant cells in which genetic mutations confer the resistance phenotype.

More recently, the “competitive release” of the resistant cells, when these cells expand into the niche freed by the killing of sensitive cells by treatment, has been proposed as mechanism of tumor recurrence, adding the non-intuitive twist that “more killing” is not better [4, 10–13, 39]. Nonlinear models of competition between sensitive and resistant subpopulations, such as that presented by Kozlowska et al., lead to similar saddle-point dynamics of tumor depletion and progression as our model, albeit at a longer time scale that encompasses multiple cell generations [46].

Another interesting fundamental similarity between genetic and non-genetic mechanisms, at least in terms of a possible formal generalization, pertains to an additional layer of non-linearity that we have not considered here: a multi-step process that gives rise to a progressive increase (‘gradual evolution’) in resilience. Increasing genomic instability with tumor progression not only decreases the threshold for cell death exploited by many therapies, but also accelerates mutational exploration of new phenotypes in an evolutionary process. A well-studied case of self-propelling increase in resistance is the amplification of the DHFR gene repeats that confers resistance to methotrexate. Once amplified, the locus is even more prone to undergo genomic recombinations and to further amplify, including in response to treatment stress that promotes chromosomal breakage [47]. The non-genetic equivalent is that with increasing malignancy (“stemness”), the barrier for cell type transitions is reduced, evident in the increasing phenotypic plasticity of advanced tumors, and hence, the chance of the tumor cell to enter an even more dedifferentiated, resilient stem-like cells in response to treatment stress [47].

The oft observed re-sensitization of recurrent tumors, rapid rate for appearance of resistance markers, and ubiquitous cell phenotype plasticity, however, suggest a role for non-genetic, reversible phenotype switching in tumor recurrence [14–19, 48–50]. Most recently, the induction of cells to transition into a stem- or mesenchymal-like resilient state by the cell stress imparted by treatment has received increasing acceptance and has been confirmed by single-cell resolution measurements [19-30]. Population-wide resistance to treatment induced by the same perturbation intended to kill and eradicate tumor cells thus poses a conflicting situation: a double-edged sword that complicates treatment response.

The mechanisms that allow a cell to be either fated to death or to enter a resistant state following the same perturbation may depend on the initial microstate due to stochastic gene expression [38]. This uncertainty is here modeled by the sigmoidal shape of the pharmacodynamical functions, which represents the cumulative probability of the dispersed propensity of cells to respond in either way to treatment. We used a minimal model to characterize how the “relative strength” (potency and efficacy) of a drug to either kill tumor cells or convert them into resistance cells affects the population dynamics, as measured by the time it takes for the cell population to recover and grow to its pre-treatment size (time to progression t). We focused on four qualitatively distinct scenarios (A, B, C, and D) that correspond to all possible combinations of high versus low potency and high versus low efficacy of the treatment to induce the resistant state relative to killing.

Despite the elementary form of the model, interesting behaviors emerge: the four scenarios produced robust, prototypic behaviors for the dependence of t on drug dose (Fig 4). The two cases (A and C) in which the potency of the drug (dose for half-maximal effect, EC50) to induce the resistant phenotype was substantially lower (high EC values) than the potency to kill cells produced a non-monotonic t-to-drug dose relationship. In such a scenario, there is an optimal window of dose for maximal t: drug doses higher than the optimal dose will have lower benefits in terms of t.

The efficacy of the drug (the amplitude of the dose-response curve) also affected t by determining its plateau value at high drug doses. High efficacy of the drug to induce the resistant state (Cases A, B) resulted, as expected, in a low value of t at high doses, independent of whether the dose-dependence is monotonic (Case B) or non-monotonic (Case A). On the other hand, a drug with low efficacy to induce treatment resistance (Cases C, D) resulted in a high value of t at high doses. Sensitivity of t to changes in drug efficacy means that drug dose optimization is paramount in the case where the potency of the drug to induce resistance is low relative to the potency of cell killing—parameters that could be determined in preclinical studies.

Unfortunately, fine-grained dose-escalation clinical trials (e.g. with at least three doses) are generally not conducted that would expose the non-monotonic effect. To establish the connection between intrinsic biological properties of drugs in triggering state transition to the resistant state and the clinical consequences as observed in drug trials, we performed virtual patient cohort simulations. We found that the qualitative conclusions of the model about the effect of induced resistance on treatment success are robust to variation in the other model parameters and result in corresponding non-monotonic dependence of the progression free survival in the simulated patient drug trial cohorts.

Adaptive therapy, in which treatment is stopped upon regression and re-started upon regrowth, or metronomic therapy, which applies a low dose at regular intervals, may be worthwhile treatment schemes in cases A and C because there is an optimal drug dose below MTD in these cases. The current rationale behind dose-minimization treatment strategies is to avoid fixation of mutant resistant clones due to competitive release and selection [10–13, 40]. However, if resistance is inducible and reversible, the intended “containment” of the tumor (as opposed to the harder-to-achieve “eradication”) is even more readily achieved than when guided by the concept of competitive release of mutant resistant clones. Thus, if we consider the new biological rationale of non-genetic, reversible dynamics of treatment-induced resistance, a strategy such as adaptive or metronomic therapy may be further optimized to be more effective. However, in order to make any meaningful predictions about optimal treatment courses we must first ground our model in clinical data, either by using empirical estimates of parameters from the literature, or by using a statistical learning framework to fit parameter distributions from the data [45, 46, 51, 52].

Table A. Summary of all possible cases for the stability of the origin under the ODE in Eq 1.

(ZIP)

Click here for additional data file.

Dynamics of tumor population recovery for Case A.

A. Left: Time-evolution of the total population N(t) plotted on a logarithmic scale for drug doses m = 23, 35, 45, 85. The horizontal dashed line indicates initial population size N(0). Right: Turning point t, at which N(t) reaches a minimum, as a function of drug dose m. Red markers indicate reference points of m = 23, 35, 45, 85. B. Parametric curve (t(m), t(m)) relating the turning point t and the time to progression t. Dashed lines of slope 2 (blue) and 3 (red) are given for reference.

(EPS)

Click here for additional data file.

Dynamics of tumor population recovery for Case B.

A. Left: Time-evolution of the total population N(t) plotted on a logarithmic scale for drug doses m = 23, 35, 45, 85. The horizontal dashed line indicates initial population size N(0). Right: Turning point t, at which N(t) reaches a minimum, as a function of drug dose m. Red markers indicate reference points of m = 23, 35, 45, 85. B. Parametric curve (t(m), t(m)) relating the turning point t and the time to progression t. Dashed lines of slope 2 (blue) and 3 (red) are given for reference.

(EPS)

Click here for additional data file.

Dynamics of tumor population recovery for Case C.

A. Left: Time-evolution of the total population N(t) plotted on a logarithmic scale for drug doses m = 23, 35, 45, 85. The horizontal dashed line indicates initial population size N(0). Right: Turning point t, at which N(t) reaches a minimum, as a function of drug dose m. Red markers indicate reference points of m = 23, 35, 45, 85. B. Parametric curve (t(m), t(m)) relating the turning point t and the time to progression t. Dashed lines of slope 2 (blue) and 3 (red) are given for reference.

(EPS)

Click here for additional data file.

Dynamics of tumor population recovery for Case D.

A. Left: Time-evolution of the total population N(t) plotted on a logarithmic scale for drug doses m = 23, 35, 45, 85. The horizontal dashed line indicates initial population size N(0). Right: Turning point t, at which N(t) reaches a minimum, as a function of drug dose m. Red markers indicate reference points of m = 23, 35, 45, 85. B. Parametric curve (t(m), t(m)) relating the turning point t and the time to progression t. Dashed lines of slope 2 (blue) and 3 (red) are given for reference.

(EPS)

Click here for additional data file.

Initial conditions, rate constants and pharmacodynamic parameters with fixed values (units not specified).

(XLSX)

Click here for additional data file.

Parameter distributions used to general virtual patient cohorts.

(XLSX)

Click here for additional data file.

20 Jan 2022

Dear Prof. Huang,

Thank you very much for submitting your manuscript "A model for the intrinsic limit of cancer therapy: duality of treatment-induced cell death and treatment-induced stemness" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

Several reviewers expressed serious concerns about the paper. In particular, Reviewer #2 indicates that the model dynamics needs to be compared with dynamics of similar simple models, which involve just competition between resistant and sensitive cells, without invoking induction of stemness and resistance.

Reviewer #3 expresses very serious concern about ability of using minimal semi-mechanistic model for providing meaningful predictions without calibration/validation against an experimental dataset.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Mark Alber, Ph.D.

Deputy Editor

PLOS Computational Biology

Mark Alber

Deputy Editor

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The authors introduce a novel mathematical model of treatment that pushes cells either into death or an increased stem-like state. The simple model is analyzed mathematically to provide approximations for time to progression (TTP), and an analysis of the important pharmacological parameters which determine when an optimal dose exists below the max tolerable dose (MTD). In general, the manuscript is interesting and well-conceived, and PLOS Computational Biology is a suitable journal for this type of manuscript.

The main result from the model is that there is non-monotonicity of eigenvalues (and therefore, there exists an optimal drug dose) only when potency to induce drug resistance is low.

Minor concerns:

1. Authors note that figure 3A & 3C are non-monotonic in lambda 1. Is this true? It’s difficult to tell visually from the figures, as it appears that lambda 1 plateaus to a lower bound. In this case, it might be helpful to pick an exaggerated example.

2. Can the authors please describe the reasoning behind displaying the summation “kSR + kRS” in figure 3. This represents the rate of flow between Resistant/Sensitive, added to the rate of flow between Sensitive/Resistant. Is this correct, and why is the addition important quantity (and not, for example, the net total rate of Sensitive: gS - kRS).

If I understand correctly, kRS is a constant (not a function of drug dose), as it does not appear in figure 2.

3. I think it would be helpful to place at least one of the supplementary figures (S1 through S4) inside the main text, to show example trajectories under treatment.

4. Figure 4 is a nice result, although I wonder if the result can be simplified. Perhaps this simulation could be simplified for plotting the dose at which max value of TTP occurs, across the range of values for PSR and ESR. If the corresponding dose is equivalent to MTD (100), then it’s implied to be monotonic.

5. The point in the discussion about adaptive therapy is intriguing. While it is certainly true that the MTD dose for cases A & C result in inferior TTP, it does not necessarily follow that adaptive therapy would outperform MTD here, because adaptive approaches often administer MTD dosing, but with periodic holidays. Perhaps it would be good to note the worthwhile treatment options of both adaptive therapy and metronomic, or other dose minimization strategies, given that there appears to be an optimal dose, m, in cases A & C.

6. Please place a box around the legend in Figures S1-4 (right side of panel A), because it makes it appear that the legend marker is a data point which doesn’t fall on the black line.

Reviewer #2: This is an interesting paper, from a prolific and highly-powered research team, which uses an extremely simple linear ODE model to depict the limits of chemotherapy, under assumption that chemotherapy leads to induction of resistance and does not just select out pre-existing resistant cells. Analysis of the model leads to the conclusion that dynamics involves a saddle point and that 4 qualitatively different regimes of treatment can be distinguished, leading to different life outcomes for the patient. If proven true, and substantiated, this may be useful as an advice to physicians and cancer counselors, and of course patients.

Presentation of the model and results is preceded by an interesting introductory review explaining the ideal of "stemness", which is a reversal of differentiation leading potentially to evolution of resistance. These are important and interesting concepts, even if many of these indicate that stemness is a feature usually present as a continuum of types (cf. the papers from Herbert Levine's group related to EMT), and not 0/1 as it is effectively assumed in the model.

In this reviewer's opinion, the paper is interesting but a thorough treatment of several important issues is missing.

1) How the dynamics of the model compares to similarly simple models, which involve just competition between resistant and sensitive cells, without invoking induction of stemness and resistance. It will be interesting and important to show that this model results in dynamics different from the almost equally simple nonlinear ODE model of lung cancer treatment by Kozlowska et al. (PLoS CB 2020), which involves only selection of pre-existing resistant cells. Will such not lead to trajectories similar as in the "saddle-point" set-up?

(2) Enough is know by now about progression dynamics of many types of cancer under treatment, to endow the model with parameters more thoroughly grounded in data. Kozlowska et al. paper on lung cancer has such estimates, and an earlier paper in Cancer Research with the same first author features estimates for ovarian cancers. The same is the case regarding the another Cancer Research paper from Bozic lab concerning colon cancer. Which of the 4 regimes will prevail under data from these papers?

(3) In the older literature, gradual evolution of resistance was considered, such as resistance to methotrexate induced by DHFR gene amplification. How will performance of the current model change under gradual evolution of ? A thorough study would require far-reaching generalization of the model, but it is possible that relevant mathematicsl results exist in the literature

It seems fair to request at least a discussion of these issues.

The paper is generally very well-written, although the term "baked-in" seems a bit idiosyncratic compared to the more common "wired-in".

Reviewer #3: In this work the authors study a very important and timely question: how treatment response and recurrence are influenced by the dual impact of cytocidal treatments: the killing impact of the drug vs the propensity for the drug to induce the transition to a resistant phenotype. By introducing a generic elementary model, the authors undertake a quantitative analysis of how the propensity of the drug to promote death versus a transition to resistance influences the often-used clinical measurement of median time to progression (TTP). The authors perform a very reasonable set of analyses on their model, centered around the question of how a drug’s capacity to induce resistance influences TTP. While I find the quantification of the results to be interesting, I am left wondering how much impact the paper will have given the simplistic nature of the model, the lack of calibration to data, and the fairly intuitive nature of many of the results. I expand upon my concerns below.

Major comments

• I always appreciate working with a minimal semi-mechanistic model. What I wonder here is: is this model too simple to provide meaningful predictions? I would be less concerned about the simplicity of the model if there was some calibration/validation against an experimental dataset. But, in the absence of this, I’m left wondering: 1) is the assumption of exponential growth sufficient, and 2) are the functional forms in eqn. (4) reasonable choices to model the impact drug dose has on the rates of cell-kill and transition-to-resistance?

• The TTP analysis is very nicely presented. I guess I’m just stuck wondering: what has been gained from the analysis? I absolutely think the message that designing a dosing strategy cannot be optimally done without understanding the propensity of the drug to induce resistance is an important one. But, if one accepts the evidence of drug-induced resistance, is it not self-evident that optimal dosing strategies would be influenced by this drug-induced phenotypic transition? If the model were grounded in data, that would be different, as then directly actionable suggestions could be made. But without that, the result really is: induced resistance can impact TTP. No doubt I agree, but I agreed with this point just from reading the intro to the paper!

• Many of the qualitative conclusions made seem fairly intuitive, given the assumptions built into the model. While more quantitative aspects of the conclusions are certainly not intuitive, those are also a lot more dependent on functional form and the somewhat arbitrary values parameters were set to. In particular, the following did seem intuitive without an analysis of the model: 1) In the case of high efficacy of induction (Cases A, B), there is diminishing TTP gains as done increases. 2) Saturation level of TTP decreases as efficacy of resistance induction increases. The one result I found less intuitive is that non-monotonicity of TTP is only observed when the potency to induce drug resistance is low. Though, I do think it's more accurate to say "low relative to the potency of tumor killing".

Minor comments

• Lines 78-79: the authors state that “However, the intrinsic limit that induced resistance places on therapy has not been systematically evaluated.” While not as systematic analysis as done herein, this was studied in a model of induced resistance in “Mathematical Approach to Differentiate Spontaneous and Induced Evolution to Drug Resistance During Cancer Treatment” by Greene and colleagues.

• Figure 1: why isn’t the drug dose m shown influencing the death of resistant cells?

• At first I was having trouble interpreting “high potency” and “low potency”. Further (and mor careful) reading made clear that high potency means that a low concentration of drug is needed to cause the halfway response between baseline induction and maximal induction. While it makes sense to call this “high potency”, I found it confusing and thought the language could be clarified earlier in the paper.

• Line 226: “We first consider a drug with a high efficacy and low potency (i.e., high EC50) for inducing resistance to cell killing”. At first the language threw me off, because there are also terms in the model that describe the efficacy and potency of cell killing. It would be clearer to me if this just said “We first consider a drug with a high efficacy and low potency (i.e., high EC50) for inducing resistance.”

• I was wondering if the authors could explain, in Figure 3 (bottom row), why all four plots seem to first hit their maximum value around the same dose of m = 40. I would have thought that the potency parameter P_SR would have influenced this time, but maybe it’s more a function of the fixed parameters in the model?

• When I downloaded the supplemental figures, they were of low image quality, even when viewing them in a small window.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Alexander R. A. Anderson

Reviewer #2: Yes: Marek Kimmel

Reviewer #3: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

10 May 2022

Submitted filename: Angelini_PLOS_CB_response.pdf

Click here for additional data file.

20 Jun 2022

Dear Ms. Angelini,

We are pleased to inform you that your manuscript 'A model for the intrinsic limit of cancer therapy: duality of treatment-induced cell death and treatment-induced stemness' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Mark Alber, Ph.D.

Deputy Editor

PLOS Computational Biology

Mark Alber

Deputy Editor

PLOS Computational Biology

***********************************************************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: RE: #1 – thank you; the insets much improve the manuscript.

RE: #2 – thank you; the clarification added to the manuscript regarding k_RS is helpful.

All other concerns were adequately addressed. The manuscript is much improved after this round of review.

Reviewer #3: The authors' responses and edits have given me a better appreciation of the work done. I especially like the addition of the section on virtual cohort simulations. I am left with just a few minor questions:

- Line 224: the line "we anchor the pharmacodynamics for killing of sensitive and resistant cell..." is unclear to me. Am I correct that this means that the PD parameters in ds(m) and dR(m) are fixed?

- Line 293: I now understand how to read this, but my brain still freezes when I read the statement "remission is prolonged under treatment by a low efficacy drug". For me, it is clearer if it read "remission is prolonged under treatment with low efficacy for inducing resistance, as compared to treatment with a drug with a high efficacy for inducing resistance." Or, something like that!

- A typical step in the generation of a virtual population cohort, when sampling from a uniform distribution, has become to optimize the parameters so that all growth curves fall within some pre-determined number of standard deviations from the expected growth curve (see: Efficient Generation and Selection of Virtual Populations in Quantitative Systems Pharmacology Models by Allen et al, 2016). Otherwise, virtual patients generated from randomly sampling uniform distributions can have dynamics very unlikely to occur in a real patient. Though, I grant, the ranges considered for most parameters is quite small, so maybe this step wasn't needed. Curious to know the authors' thoughts.

A few typos noticed:

- Line 220: "a lowering the resistant cell death"

- Line 230 "alteringg"

- Line 460: missing period

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Alexander Anderson

Reviewer #3: No

19 Jul 2022

PCOMPBIOL-D-21-02030R1

A model for the intrinsic limit of cancer therapy: duality of treatment-induced cell death and treatment-induced stemness

Dear Dr Huang,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Zsofia Freund

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol