Dose-Dependent Mutation Rates Determine Optimum Erlotinib Dosing Strategies for EGFR Mutant Non-Small Cell Lung Cancer Patients

Abstract

Background

The advent of targeted therapy for cancer treatment has brought about a paradigm shift in the clinical management of human malignancies. Agents such as erlotinib used for EGFR-mutant non-small cell lung cancer or imatinib for chronic myeloid leukemia, for instance, lead to rapid tumor responses. Unfortunately, however, resistance often emerges and renders these agents ineffective after a variable amount of time. The FDA-approved dosing schedules for these drugs were not designed to optimally prevent the emergence of resistance. To this end, we have previously utilized evolutionary mathematical modeling of treatment responses to elucidate the dosing schedules best able to prevent or delay the onset of resistance. Here we expand on our approaches by taking into account dose-dependent mutation rates at which resistant cells emerge. The relationship between the serum drug concentration and the rate at which resistance mutations arise can lead to non-intuitive results about the best dose administration strategies to prevent or delay the emergence of resistance.

Methods

We used mathematical modeling, available clinical trial data, and different considerations of the relationship between mutation rate and drug concentration to predict the effectiveness of different dosing strategies.

Results

We designed several distinct measures to interrogate the effects of different treatment dosing strategies and found that a low-dose continuous strategy coupled with high-dose pulses leads to the maximal delay until clinically observable resistance. Furthermore, the response to treatment is robust against different assumptions of the mutation rate as a function of drug concentration.

Conclusions

For new and existing targeted drugs, our methodology can be employed to compare the effectiveness of different dose administration schedules and investigate the influence of changing mutation rates on outcomes.

About the Authors:

Lin L. Liu

Affiliations Department of Biostatistics, Harvard T. H. Chan School of Public Health, Boston, MA 02115, United States of America, Department of Biostatistics and Computational Biology, Dana-Farber Cancer Institute, Boston, MA 02215, United States of America

Fei Li

Affiliation: Department of Biostatistics, Harvard T. H. Chan School of Public Health, Boston, MA 02115, United States of America

William Pao

Affiliation: Vanderbilt-Ingram Cancer Center, Nashville, TN 27232, United States of America

Franziska Michor

* E-mail: [email protected]

Affiliations Department of Biostatistics, Harvard T. H. Chan School of Public Health, Boston, MA 02115, United States of America, Department of Biostatistics and Computational Biology, Dana-Farber Cancer Institute, Boston, MA 02215, United States of America

Introduction

Recent advances have improved our understanding of the molecular alterations that drive particular cancer types and have thus enabled the development of targeted agents that specifically inhibit these lesions [1]. Examples of targeted therapies include small molecule inhibitors of the epidermal growth factor receptor (EGFR) pathway in lung cancer (e.g. erlotinib (Tarceva)) and inhibitors of the BCR-ABL tyrosine kinase in chronic myeloid leukemia (e.g. imatinib (Gleevec), dasatinib (Sprycel), and nilotinib (Tasigna)). These small molecule inhibitors are taken up into cancer cells where they interfere with abnormal signaling. Targeted therapy differs from traditional cytotoxic chemotherapy in that it not only leads to more specific effects with reduced toxicity, but also promises a future of personally tailored anti-cancer treatment [2].

The development of targeted anti-cancer therapies requires the design of optimal treatment strategies so that responses are maximized while toxicity remains tolerable [3]. Because of the combinatorial complexity of this problem, systematic and mathematical approaches have been employed in the past to identify best treatment modalities. In a seminal paper in 1977, Norton and Simon proposed a model of kinetic (i.e., non-genetic) resistance to cell-cycle specific therapy in which tumor cells followed a Gompertzian growth law [4]. This work led the authors to propose a dose-intensification strategy, supported with historical data [5] and later on Implemented as a prospective clinical trial [6]. Their model and its predictions have become known as the Norton-Simon hypothesis and inspired many subsequent studies of kinetic resistance [7–14]. In parallel, several investigations addressed the emergence of genetic resistance, i.e. resistance driven by genetic alterations in cancer cells. Coldman and co-authors were the first to introduce stochastic models of resistance against chemotherapy to guide treatment schedules [15], which led to many subsequent studies by these authors [15–17] and others [18–28]. Several other papers addressed the question about optimum dosing of targeted therapy by including the effect of quiescence on the kinetics of response to treatment [29, 30]. We recently developed a stochastic framework to optimize dosing strategies of targeted drugs [31, 32]; when applied to EGFR-mutant non-small cell lung cancer, this model allowed us to identify a treatment schedule predicted to maximally delay the onset of T790M-driven resistance [33], which is the most common mechanism of disease progression. This schedule is currently being validated in a clinical trial at Memorial Sloan-Kettering Cancer Center (NCT01967095), in which patients receive high dose erlotinib on day 1 and day 2 with daily low dose erlotinib on days 3 through 7. Dose levels for the high dose pulse are escalated starting from 600mg oral daily on day 1 and day 2 until the maximally tolerated dose is reached. First results are expected in 2016.

In this paper, we aim to extend these investigations by taking into account the ability of a targeted drug to alter the kinetics of mutation acquisition in cancer cells: the rate at which resistant cells emerge might now depend on the dose of drug administered to the patient. A possible mechanism driving this phenomenon might be the fact that DNA damage can result from the intermediates of oxygen reduction in a cell. Reactive oxygen species (ROS) are short-lived oxygen intermediates and are generated primarily by the mitochondrial respiratory chain in inflammatory cells. By interacting with free nucleotides, ROS promotes DNA damage [34], which can subsequently result in a higher frequency of mutations in the cell. Koptyra et al. [35], for instance, showed that by inhibiting ROS in leukemia cells, the mutagenesis rate was decreased, which in turn resulted in a decreased frequency of resistance against imatinib. Recent studies on the effects of traditional cytotoxic drugs have also revealed a connection between ROS and induced apoptosis in cancer cells [36]. Furthermore, traditional cytotoxic drugs were found to be associated with the generation of ROS [37–39]. The administration of anti-cancer drugs, both traditional chemotherapeutic agents and targeted drugs, may thus modulate the rate at which mutations arise in cancer cells, thereby influencing the dynamics of resistance. Here we derive a quantitative framework of dynamic mutation rates during the treatment of cancer and use this framework to identify best treatment modalities for the clinical management of human malignancies.

Methods

Consider a population of tumor cells proliferating within a tissue. We model the evolutionary dynamics of the tumor population as a multi-type non-homogeneous continuous-time birth-death stochastic process. The population consists of two types of cells: sensitive and resistant cancer cells. These cells might describe the small subset of cells capable of propagating the entire tumor cell population (i.e., “cancer stem cells”) or alternatively the entire tumor mass (i.e.“tumor bulk”). At any given time point t, the number of sensitive cancer cells is denoted as X_s(t), whereas the number of resistant cancer cells is represented by X_r(t). Sensitive cancer cells proliferate at rate λ_s(t) and die at rate μ_s(t). During each sensitive cell division, a mutation that confers resistance occurs at probability u(t); this quantity might depend on time as well. Resistant cancer cells proliferate and die at rates λ_r(t) and μ_r(t). Note that we use λ_s(t), λ_r(t), μ_s(t), μ_r(t), and u(t) as per capita rates or probabilities per time unit. The possibility of back-mutation from resistant to sensitive cells is excluded. Treatment alters growth and/or death rates of sensitive and potentially also resistant cells, thereby modulating these rates over time depending on the treatment schedule used (Fig 1).

[Figure omitted. See PDF.]

Fig 1. Schematic overview of the model.

We developed a mathematical model to investigate the effects of changing mutation rates during treatment on the evolution of resistance. When treated with targeted or traditional cytotoxic chemotherapy, a sensitive cancer cell (left) might give rise to a resistant cell at a rate that increases with the drug dose administered (case 1), is independent of the dose (case 2), or decreases with the dose (case 3). For simplicity, we show linear relationships between drug dose and mutation rate; however, more complex relationships can easily be considered using our general framework.

https://doi.org/10.1371/journal.pone.0141665.g001

We define the birth-death process as X(t) ≡ (X_s(t),X_r(t)). Let us first consider situations in which the initial cancer population merely consists of M sensitive cells, i.e. X(t) = (M, 0). This assumption will be relaxed in later sections to better describe biological situations [40, 41]. This stochastic process is defined by the following infinitesimal transition probabilities:(1)

Let us now derive some quantities of interest including the expected number of resistant cells, the variance of the number of resistant cells, and the probability of resistance. In all cases, we generalize our calculations to cases with a pre-existing resistant clone before treatment. The cases without any pre-existing resistant clones can be derived by setting X_r(0) = 0 for all following equations.

To perform exact computer simulations of this process, we employed adaptive thinning to overcome the difficulties resulting from time-dependent birth, death and mutation rates [42]. In this algorithm, exponential waiting times between events are generated by first calculating a constant rate which majorizes the true instantaneous rate at any time t. Then, for each event, an exponential waiting time is generated with this majorizing rate. This event is accepted if a uniform random variable generated on [0,1] falls below the ratio of the true rate to the majorizing rate; otherwise, the event is rejected. The simulations are used to validate the approximations developed below.

The expected number of resistant cells

Since the mutation rate is always much less than 1, u < < 1, we can approximate the birth rate of sensitive cells, λ_s(1−u), with λ_s in the derivations of the important characteristics of our model. We will later verify the validity of this approximation through the consistency of our theoretical derivations with the exact stochastic computer simulations; see the Results section. Let us first calculate the rate of transition from sensitive to resistant cells, given by(2)

Then the expected number of resistant cells arising from the sensitive cell population is given by(3)

Hence the expected number of resistant cells at time t is given by(4)where X_rr(t) is the number of resistant cells generated from resistant cells at time t, and X_sr(t) is the number of resistant cells generated from sensitive cells at time t respectively.

The probability of resistance

Another quantity of interest is the probability that there is at least one resistant cell present at time t. For an infinitesimal time interval [t_i,t_i + Δt], the rate of emergence of resistant cells is given by r(t_i)Δt. As shown previously [32], under the assumption of no back mutation, the probability of resistance arising from the sensitive cell population is given by(5)where P₀(T) denotes the probability of having no resistant cells at time T, the probability of having no resistant cells at time T that is generated from the sensitive clone, the probability that the clone arising from the initial resistant cell has become extinct by time T, and(6)

The variance of the number of resistant cells

The variance of X_r(T) is also of interest because it provides a sense of the uncertainty of our estimate of the number of resistant cells at time T. Let us first calculate the variance of resistant cells at time T that originated from sensitive cells:(7)

The variance of resistant cells originating from pre-existing resistant clones at time T is then given by(8)

Hence the total variance is given by Var(X_r(T)) = Var_s(X_r(T)) + Var_r(X_r(T)).

Estimating pharmacokinetic parameters and growth kinetics

Let us now estimate the pharmacokinetic parameters of erlotinib. We obtained data from pharmacokinetic studies from OSI/Astellas. In total, 28 subjects were given 100, 200, 400, 800, 1000, 1200, 1400, and 1600 mg of erlotinib, and the concentration of drug in the serum of these patients was measured before administration and at 2 hours, 8 hours and 24 hours after administration; see S1 Table for the data. We then employed an exponential decay function to model the concentration over time; this function is given by C(t) = C_max e^−κt, where C_max is the maximum concentration and κ is the elimination rate of the drug. Although both parameters vary across subjects, we used only one set of parameters for each dose for our primary results because of the relatively large sample size versus the small number of measured time points. The relationship between the dose amount and C_max was estimated, using least square methods, as C_max(d) = 1.393 + 0.0198d. The rate κ = 0.05 per hour was estimated as the mean decay rate for the different dosing groups. Quantitative estimates of birth and death rates of sensitive and resistant cells were experimentally determined as in [33], using a pair of isogenic PC-9 human EGFR-mutant cell lines with and without the T790M point mutation that were treated with different concentrations of erlotinib. Using the cell counts of viable and dead cells, we then determined the exponential growth and death rates of the two cell types during different concentrations of the drug for use in the stochastic model outlined above and in [33], where μ_Xs(t) ≈ 0.005 hour⁻¹, μ_Xr(t) ≈ 0.002 hour⁻¹, λ_Xs(t) ≈ exp(−4.4⋅C(t)−3.17) hour⁻¹, and λ_Xr(t) ≈ −0.001⋅C(t) + 0.03 hour⁻¹.

All numerical calculations and simulations were coded in C++ and all statistical analyses were performed using open source R [43] software.

Results

We then sought to validate our analytical approximations using exact stochastic computer simulations. As an example, let us consider the process specified in Eq (1) and define the birth and death rates for sensitive and resistant cells as follows:

These equations provide an example of a dose administration strategy leading to a drug concentration that varies like a sine function over time. Using this hypothetical treatment strategy, we then explored three different dependencies of the mutation rate on the dose: (1) a time-invariant mutation rate (independent of the drug concentration); (2) a mutation rate that increases with the drug concentration (or equivalently, decreasing with the birth rate); and (3) a mutation rate that decreases with the drug concentration (or equivalently, increasing with the birth rate). We can also interpret scenario (2) as a situation in which the sensitive cells do not have any delay in responding to the treatment effects of increasing the mutation rate, while scenario (3) represents a situation in which there is a one-period delay in the response. Thus we have

To validate the predictions of our analytical approximations, we first compared them with the output of exact stochastic computer simulations and observed good agreement (Fig 2). Our results show that the effects of the drug dose on the mutation rate (independent, increasing or decreasing mutation rate with drug concentration) influences the expected number of resistant cells (Fig 2A), the probability of resistance (Fig 2B), and the variance of resistant cells (Fig 2C). In particular, in this example, when the mutation rate changes in the opposite direction as the birth rate, as in the right panel of Fig 2A, the mutation rate increases with the drug concentration. When the mutation rate increases with the dose, the initially homogeneous tumor cell population is more prone to become heterogeneous as compared to the other scenarios. Nonetheless, the decreasing birth rate of the resistant cells impairs their growth. This fact is reflected in the lower expected number of resistant cells and probability of resistance over time.

[Figure omitted. See PDF.]

Fig 2. Evolutionary dynamics of sensitive and resistant clones during continuous therapy with different mutation dynamics.

(A) Expected number of resistant cancer cells as a function of time during continuous therapy. Blue line: mutation rate is constant during treatment; black line: mutation rate increases with drug dose; green line: mutation rate decreases with drug dose. Red circles: simulation results. Grey shaded area indicates one analytic standard deviation from the analytic mean. (B) Probability of resistance as a function of time during continuous therapy. (C) Variance of resistant cancer cells as a function of time during continuous therapy. The following parametrizations were used for both simulation and analytic approximations: A_s = 0.05, B_s = 0.1, C_s = 0.005, A_r = 0.05, B_r = 0.12, C_r = 0.002, and θ = 0.10. The values for A_u and B_u are denoted in the panels for each corresponding scenario.

https://doi.org/10.1371/journal.pone.0141665.g002

We then considered a tumor cell population containing a proportion s of pre-existing resistant cells. We investigated a pulsed therapy schedule, where drug is administered for 14 days followed by a 14-day treatment holiday (step function), as well as the sine functional form of the dosing strategy. The step functions are given by:where I(⋅) denotes the indicator function. In this case, the death rate is held as a constant as in Fig 2, while the mutation rate is assumed to be either constant, changing in the same direction, or in the opposite direction as the birth rate. We again compared our theoretical approximations to the results of the exact stochastic simulations and found good agreement (Fig 3). The same observation as in Fig 2 can be obtained when (1) there is a pre-existing resistant clone and (2) a step function of birth, death and mutation rates is chosen: when the mutation rate changes in the opposite direction as the dosing regime, then the tumor cells are more likely to develop a resistance mutation as compared to the other two scenarios, i.e. the time-invariant mutation rate and the mutation rate increasing with the drug concentration. Interestingly, we noticed that when there are pre-existing resistant clones, the impact of the dose-dependent mutation rate is not as remarkable as when the initial cell population is homogeneous.

[Figure omitted. See PDF.]

Fig 3. Evolutionary dynamics of sensitive and resistant clones during continuous therapy with pre-existing resistant cells.

(A), (B) and (C) are examples for sine wave functional forms of birth, death and mutation rates. (A) Expected number of resistant cancer cells as a function of time during continuous therapy. Blue line: mutation rate is constant during treatment; black line: mutation rate increases with drug dose; green line: mutation rate decreases with drug dose. Red circles: simulation results for no pre-existing resistant clones; orange circle: simulation result for 3% proportion of pre-existing resistant clones; purple circle: simulation result for 5% proportion of pre-existing clones. (B) Probability of resistance as a function of time during continuous therapy. (C) Variance of resistant cancer cells as a function of time during continuous therapy. (D), (E), and (F) are examples for piecewise functional forms of birth, death and mutation rates. (D) Expected number of resistant cancer cells as a function of time during continuous therapy. (E) Probability of resistance as a function of time during continuous therapy. (F) Variance of resistant cancer cells as a function of time during continuous therapy. The death rates in Fig 2 were used here as the death rates for both upper and lower panels. The birth rates in the continuous therapy in Fig 2 were also used in this figure. The birth rates in the piecewise strategy were λ_Xs(t) = 0.15⋅I(t / 14 mod 2 = 0) + 0.05⋅I(t / 14 mod 2 ≠ 0), λ_Xr(t) = 0.17⋅I(t / 14 mod 2 = 0) + 0.07⋅I(t / 14 mod 2 ≠ 0).

https://doi.org/10.1371/journal.pone.0141665.g003

Next, we investigated the pharmacokinetic effects of drug accumulation in vivo. From this point on, all results were obtained from analytic approximation, which has been proved to be consistent with the exact simulation in the previous sections. We modeled the drug elimination in the human system after one dose as an exponential decay with rate κ. When the dosing amount is given as D with corresponding C_max using the equation developed in the Methods section, then the drug concentration in the body is given by Concentration(t) = C_max e^−κt/T, where T is the dosing interval. Note that the drug concentration reaches a steady state around C_max/[1−e^−κt/T] after a sufficiently long time.

A loading dose is sometimes used to cause the drug concentration to reach its steady level more quickly. We thus compared the effects of different mutation rate scenarios (independent, increasing and decreasing rates with dose) for drug administration schedules with and without loading doses. The birth and death rates of non-small cell lung cancer cells as a function of the erlotinib concentration were estimated as in [33]; the birth rate is a linear function of the drug concentration and the death rate is given by an exponential decay of the drug concentration. All coefficients were obtained through least square estimation. The dosing regime in terms of drug concentration in vivo and the related birth and death rates are displayed in Fig 4A and 4B, respectively. We considered a scenario in which there are 10⁶ tumor (“stem”) cells capable of propagating the entire tumor cell population, with no pre-existing resistant cells. For each dosing strategy, we again considered three different functional forms for the mutation rate: (a) the mutation rate is independent of the drug concentration and hence is fixed throughout treatment, u(t) = u₀; (b) the mutation rate linearly increases with the drug concentration, u(t) = β ⋅ Concentration(t) + u₀; and (c) the mutation rate linearly decreases with the drug concentration: u(t) = −β ⋅ Concentration(t) + u₀, where β > 0 denotes the effect of the drug concentration on the mutation rate in vivo, while u₀ is the baseline mutation rate, chosen as u₀ = 10⁻⁸ (Fig 4C). Here we controlled the independent mutation rate as the time-average of the mutation rate increasing (or decreasing) with the drug concentration. The probability of resistance over time for these scenarios is shown in Fig 4D. In all cases, a dosing strategy with a loading dose outperforms the strategy without a loading dose. Based on the current parametrization, the lowest probability of resistance was obtained in the case in which the mutation rate monotonically increases with the drug concentration, suggesting a different treatment performance under different mutation rate assumptions. To test the robustness of our results, we performed similar analyses with β = 10⁻⁹ (S1 Fig) and found robust results. Together, these findings demonstrate that the magnitude of the mutation rate does not change the relative order of the performance of the three mutation regimes (a)–(c): the mutation rate linearly increasing with dose delivers the lowest probability of emerging resistance.

[Figure omitted. See PDF.]

Fig 4. Pharmacokinetic effects influence the evolutionary dynamics of resistance.

Here we consider a baseline mutation rate per cell division of 10⁻⁸ and the effect of drug concentration on the mutation rate β = 10⁻⁹. (A) The drug concentration in vivo based on the pharmacokinetic model over time for a 100mg per day continuous dosing regime. Dotted line: loading dose; solid line: no loading dose. (B) The birth rates as a function of time t. Red line: the birth rate of the resistant cells; black line: the birth rate of the sensitive cells. (C) The mutation rate of sensitive cells as a function of time t. Blue line: constant mutation rate; green line: mutation rate monotonically decreases with the drug concentration; black line: mutation rate monotonically increases with the drug concentration. (D) The probability of resistance as a function of time t. Values for birth and death rates: μ_Xs(t) ≈ 0.005 hour⁻¹, μ_Xr(t) ≈ 0.002 hour⁻¹, λ_Xs(t) ≈ exp(−4.4 ⋅ C(t)−3.17) hour⁻¹, and λ_Xr(t) ≈ −0.001 ⋅ C(t) + 0.03 hour⁻¹.

https://doi.org/10.1371/journal.pone.0141665.g004

We then compared the performance of different dosing regimes under the three different types of dynamic mutation rate assumptions as described above. These dosing regimes include: 100mg/day, 1600mg/week, and 1600mg/week combined with 100mg/day, together with loading and no loading doses. The dosing regimes as a function of time t are shown in Fig 5A and 5B. The mutation rates over time under different assumptions are displayed in Fig 5C. Overall, when comparing across different mutation rate assumptions, without loading doses, we found that 150mg/day and the combined strategy of 1600mg/week + 100mg/day lead to the lowest probability of resistance without pre-existing resistant clones (Fig 5D–5F). When a loading dose is given, the relative performance of different strategies is altered, especially between the 150mg/day and combined strategies, where 150mg/day outperformed the combined strategy under all mutation rates assumptions. When there are pre-existing resistant clones, however, the presence of a loading dose did not significantly alter the dynamics (Fig 5G–5I). Finally, the different mutation rates assumptions also did not have an effect of the relative performance of different dosing schedules. Fig 5 displays results where β = 10⁻¹⁰, and similar results are shown in S2 Fig for β = 10⁻⁹ to corroborate our conclusion: within the scenarios tested, different mutation rates do not alter the relative performance of different dosing schedules.

[Figure omitted. See PDF.]

Fig 5. Effects of varying dosing regimes on the evolution of resistance.

Here we consider β = 10⁻¹⁰. (A) The dosing regimes with no loading dose for 100mg/day, 150mg/day, 1600mg/week, 1600mg/week combined with 100mg/day during the week, and 1600mg/week combined with 75mg/day during the week. (B) The dosing regimes with loading dose for 100mg/day, 150mg/day, 1600mg/week, 1600mg/week combined with 100mg/day during the week, and 1600mg/week combined with 75mg/day during the week. (C) Mutation rate as a function of treatment concentration under different assumptions: blue: independent with treatment concentration; black: increasing with treatment concentration; green: decreasing with treatment concentration. (D)-(F) Without pre-existing resistance, the probability of resistance monitored up to one month under (D) constant mutation rate, (E) mutation rate increasing with the drug concentration, and (F) mutation rate decreasing with the drug concentration. (G)-(I) With pre-existing resistance, the expected number of resistant cells monitored up to one month under (G) constant mutation rate, (H) mutation rate increasing with the drug concentration, and (I) mutation rate decreasing with the drug concentration. Dotted line: with loading dose; solid line: without loading dose. Values for birth and death rates: μ_Xs(t) ≈ 0.005 hour⁻¹, μ_Xr(t) ≈ 0.002 hour⁻¹, λ_Xs(t) ≈ exp(−4.4 ⋅ C(t)−3.17) hour⁻¹, and λ_Xr(t) ≈ −0.001 ⋅ C(t) + 0.03 hour⁻¹

https://doi.org/10.1371/journal.pone.0141665.g005

Next, we took pre-existing resistance into consideration, assuming that there are 100 pre-existing resistant cells in the tumor. We then used the expected number of resistant cells in the population as a proxy for evaluating the performance of each dosing strategy. In this situation, different mutation rate assumptions did not lead to significantly different performance results of different dosing regimes (Fig 5F–5H): the combined therapy leads to the best performance under all three different assumptions.

Finally, we considered the optimum dosing strategy for the three different forms of mutation rate functions discussed above under the constraint of side effects and toxicity. We used a similar constraint as derived in [32], where ToleratedDose(Freq) = 50exp(−0.4 × Freq). For each dosing frequency corresponding to the dosing amount (100mg, 800mg, 1600mg) explored in our study, we rounded them to the unit of half day. In addition, the FDA-approved 150mg/day dosing regime and two combined strategies 1600mg/week + 75mg/day and 1600mg/week + 100mg/day were also included in the analysis. Only schedules with no loading dose were explored. We then monitored the following measures as the utilities of each dosing scheme: (i) the area under the curve of the trajectory of probability of resistance, (ii) the area under the curve of the expected number of resistant cell over time, (iii) the probability of resistance and (iv) the expected number of resistant cells after one and three months of treatment (Fig 6 and S3 Fig). In all cases, different assumptions of the mutation rates did not change the ranking of the dosing schemes based on either of the different measures.

[Figure omitted. See PDF.]

Fig 6. The optimum dosing strategy varies with dose-dependent mutation rates.

Here we consider β = 10⁻¹⁰. (A)–(C): (A) The area under the curve (AUC) of the probability of resistance over three months of treatment, (B) the probability of resistance after one month of treatment, and (C) the probability of resistance after three months of treatment when there are no pre-existing resistant cells for different treatment schedules (separated by different colors) and different mutation rate assumptions (independent; increasing with drug concentration; decreasing with drug concentration) indicated by the x-axis label. (D)–(F): (D) The AUC of the expected number of resistant cells over three months of treatment, (E) the expected number of resistant cells after one month of treatment, and (F) the expected number of resistant cells after three months of treatment when there are pre-existing resistant clones for different treatment schedules (separated by different colors) and different mutation rate assumptions (independent; increasing with drug concentration; decreasing with drug concentration) indicated by the x-axis label. Values for birth and death rates: μ_Xs(t) ≈ 0.005 hour⁻¹, μ_Xr(t) ≈ 0.002 hour⁻¹, λ_Xs(t) ≈ exp(−4.4 ⋅ C(t)−3.17) hour⁻¹, and λ_Xr(t) ≈ −0.001 ⋅ C(t) + 0.03 hour⁻¹.

https://doi.org/10.1371/journal.pone.0141665.g006

Discussion

In this study, we have built a mathematical framework by using a non-homogeneous multi-type birth-and-death process model with dosing-dependent mutation rates to investigate the evolution of resistant clones to targeted cancer therapies. The consideration of dose-dependent mutation rates can lead to different stochastic dynamics of the sensitive and resistant cancer cell populations as compared to the dose-independent case [44]. Thus, its inclusion might lead to different dosing strategies that are best able to prevent or delay the emergence of the resistant cell population. Throughout this paper, we have also assumed that the drug concentration inside the tumor is uniform since there is insufficient data available to parameterize other functional forms. Should the drug concentration be more uneven across the tumor, then regions with lower concentration might lead to different rates of generating resistant cells. Such conditions will be analyzed when more data becomes available. We applied our mathematical framework to non-small cell lung cancer treated in particular with the the targeted agent erlotinib to identify best dosing strategies for this cancer type; we found that in this application, different considerations about the relationship between the mutation rate and the drug concentration do not change the finding that resistance is optimally delayed using a low dose continuous plus high dose pulsed approach [33]. This strategy is currently being tested as a prospective clinical trial at Memorial Sloan-Kettering Cancer Center (NCT01967095). While this clinical trial of our original optimized dosing strategy (see Chmielecki et al, Science Translational Medicine 2011 [33]) is ongoing, the current work is meant to provide a more in-depth analysis of situations in which the drug concentration might influence the kinetics of mutation acquisition, in order to evaluate whether such scenarios might change the optimum dosing strategies. No direct clinical validation is currently possible.

Currently, the emergence of resistant cells is attributed to clonal expansion of pre-existing populations of cells that harbor such resistance mutations [45–50]. For example, EGFR T790M can be detected at low levels in patients with EGFR-mutant NSCLC prior to treatment [50]. However, an alternative explanation is that kinase inhibitor treatment indirectly induces resistance through DNA damage incurred as a result of increased cell killing and subsequent release of ROS [34, 35]. To account for this possible scenario, our mathematical framework incorporates non-homogeneous multi-type birth-and-death processes with a dosing-dependent mutation rate to investigate the evolution of resistant clones to targeted cancer therapies under general dosing regimens. We assume that a dose-dependent mutation rate may produce different stochastic dynamics of sensitive and resistant cancer cell population as compared to the dose-independent case, and thus may suggest different dosing strategies to prevent or at least slow down the emergence of the resistant cell population. In addition to the T790M mutation, cells might accumulate other genetic and/or epigenetic changes that might change the fitness of the cell. To address this situation, in our model, we implicitly assume that each cell population has a distribution of growth rates and death rates, and the kinetically dominant subpopulation of each mutation type (sensitive those cells that do not harbor the T790M resistance mutation but potentially other alterations, and resistant those cells that do harbor the T790M resistance mutation and also potentially other changes) drives the outgrowth of such clones. Specifically, this new model uses detailed pharmacokinetic data to examine three possible clinical scenarios involving EGFR mutant lung cancer treated with erlotinib: mutation rates are independent of treatment, increase with treatment, or decrease with treatment. In all situations considered, including pre-existing or no pre-existing resistant clones, and loading or no loading dose in the beginning of the treatment schedules, the relative performance of different dosing strategies is robust against assumptions of mutation rates. These findings suggest that considering constant mutation rate is sufficient for modeling the effects of different dosing strategies in EGFR-mutant lung cancer treated with erlotinib. Currently, little data is available for cell kinetics in response to DNA damaging agents and therefore we have only focused on the treatment response to erlotinib in this paper. However, our modeling framework can be applied to DNA damaging agents as well as long as the relevant growth and death rates are available.

For more precise parameterization of our model, time-series data of the growth, death and mutation rates during different drug combinations and doses are required. Next-generation sequencing has been widely used to determine the mutation frequencies under different treatments or cancer subtypes [51], but it will require a very large sequencing depth and information on the total number of cell divisions to obtain estimates for mutation rates per cell division. Future studies will attempt to confirm our models based upon actual mutation rate data induced by inhibitor treatment. Our model can also be extended to investigate combination drug therapy, which has been suggested as a possible solution for overcoming resistance to cancer therapy [52].

Supporting Information

[Figure omitted. See PDF.]

S1 Fig. Pharmacokinetic effects influence the evolution of resistance.

Here we consider the baseline mutation rate per cell division to be 10⁻⁸ and the effect of drug dose on the mutation rate to be β = 10⁻¹⁰. (A) The drug concentration in vivo based on the pharmacokinetic model over time for 100mg per day dosing regime. Dotted line: loading dose; solid line: no loading dose. (B) The birth rates as a function of time t. Red line: the birth rate of the resistant cells; black line: the birth rate of the sensitive cells. (C) The mutation rate of sensitive cells as a function of time t. Blue line: constant mutation rate; green line: mutation rate monotonically decreases with the drug concentration; black line: mutation rate monotonically increases with the drug concentration. (D) The probability of resistance as a function of time t. Values for birth and death rates: μ_Xs(t) ≈ 0.005 hour⁻¹, μ_Xr(t) ≈ 0.002 hour⁻¹, λ_Xs(t) ≈ exp(−4.4 ⋅ C(t)−3.17) hour⁻¹, and λ_Xr(t) ≈ −0.001 ⋅ C(t) + 0.03 hour⁻¹.

https://doi.org/10.1371/journal.pone.0141665.s001

(TIFF)

S2 Fig. Effects of varying dosing regimes on the evolution of resistance.

Here we consider β = 10⁻⁹. (A) The dosing regimes with no loading dose for 100mg/day, 150mg/day, 1600mg/week, 1600mg/week combined with 100mg/day during the week, and 1600mg/week combined with 75mg/day during the week. (B) The dosing regimes with loading dose for 100mg/day, 150mg/day, 1600mg/week, 1600mg/week combined with 100mg/day during the week, and 1600mg/week combined with 75mg/day during the week. (C) Mutation rate as a function of treatment concentration under different assumptions: blue: independent with treatment concentration; black: increasing with treatment concentration; green: decreasing with treatment concentration. (D)-(F) Without pre-existing resistance, the probability of resistance monitored up to one month under (D) constant mutation rate, (E) mutation rate increasing with the drug concentration, and (F) mutation rate decreasing with the drug concentration. (G)-(I) With pre-existing resistance, the expected number of resistant cells monitored up to one month under (G) constant mutation rate, (H) mutation rate increasing with the drug concentration, and (I) mutation rate decreasing with the drug concentration. Dotted line: with loading dose; solid line: without loading dose. Values for birth and death rates: μ_Xs(t) ≈ 0.005 hour⁻¹, μ_Xr(t) ≈ 0.002 hour⁻¹, λ_Xs(t) ≈ exp(−4.4 ⋅ C(t)−3.17) hour⁻¹, and λ_Xr(t) ≈ −0.001 ⋅ C(t) + 0.03 hour⁻¹.

https://doi.org/10.1371/journal.pone.0141665.s002

(TIFF)

S3 Fig. The optimum dosing strategy varies with dose-dependent mutation rates.

Here we assume β = 10⁻⁹. (A)–(C): (A) The AUC of probability of resistance over three months of treatment, (B) the probability of resistance after one month of treatment, and (C) the probability of resistance after three months of treatment when there are no preexisting resistant clones for different treatment schedules (separated by different colors) and different mutation rate assumptions (independent; increasing with drug concentration; decreasing with drug concentration) indicated by the x-axis label. (D)–(F): (D) The AUC of expected number of resistant cells over three months of treatment, (E) the expected number of resistant cells after one month of treatment, and (F) the expected number of resistant cells after three months of treatment when there are preexisting resistant clones for different treatment schedules (separated by different colors) and different mutation rate assumptions (independent; increasing with drug concentration; decreasing with drug concentration) indicated by the x-axis label. Values for birth and death rates: μ_Xs(t) ≈ 0.005 hour⁻¹, μ_Xr(t) ≈ 0.002 hour⁻¹, λ_Xs(t) ≈ exp(−4.4 ⋅ C(t)−3.17) hour⁻¹, and λ_Xr(t) ≈ −0.001 ⋅ C(t) + 0.03 hour⁻¹.

https://doi.org/10.1371/journal.pone.0141665.s003

(TIFF)

S1 Table. Pharmacokinetic data for the concentration of erlotinib in the serum of 28 patients from OSI/Astellas study.

T_max denotes the average maximum of time to reach C_max, where C_max is the maximum concentration measured in the serum. The patients are followed at 2 hours, 8 hours, and 24 hours after erlotinib delivery.

https://doi.org/10.1371/journal.pone.0141665.s004

(XLSX)

Acknowledgments

The authors would like to thank Jasmine Foo, Rui Zhao, Abhishek Chakrabortty, Emmanuel Dimont, Donghoon Shin, and Michor Lab for insightful discussions. The authors also acknowledge support from the Dana-Farber Cancer Institute Physical Sciences-Oncology Center (NCI U54CA143798), from US National Cancer Institute grants R01CA121210 and P01CA129243, and OSI/Astellas for sharing the pharmacokinetic data.

Author Contributions

Conceived and designed the experiments: LLL FL WP FM. Performed the experiments: LLL FL. Analyzed the data: LLL FL FM. Contributed reagents/materials/analysis tools: LLL FL. Wrote the paper: LLL FL WP FM.

Citation: Liu LL, Li F, Pao W, Michor F (2015) Dose-Dependent Mutation Rates Determine Optimum Erlotinib Dosing Strategies for EGFR Mutant Non-Small Cell Lung Cancer Patients. PLoS ONE 10(11): e0141665. https://doi.org/10.1371/journal.pone.0141665

References

1. Allgayer H, Fulda S. An introduction to molecular targeted therapy of cancer. Adv Med Sci. 2008;53(2):130–138. pmid:18650147

2. Gerber DE. Targeted therapies: a new generation of cancer treatments. Am Fam Physician. 2008 Feb;77(3):311–319. pmid:18297955

3. Hutson TE, Figlin RA, Kuhn JG, Motzer RJ. Targeted therapies for metastatic renal cell carcinoma: an overview of toxicity and dosing strategies. Oncologist. 2008 Oct;13(10):1084–1096. pmid:18838439

4. Norton L, Simon R. Tumor size, sensitivity to therapy, and design of treatment schedules. Cancer Treat Rep. 1977;61:1307–1317. pmid:589597

5. Norton L, Simon R. The Norton-Simon hypothesis revisited. Cancer Treat Rep. 1986;70:163–169. pmid:3510732

6. Citron ML, et al. Randomized trial of dose-dense versus conventionally scheduled and sequential versus concurrent combination chemotherapy as postoperative adjuvant treatment of node-positive primary breast cancer: first report of Intergroup Trial C9741/Cancer and Leukemia Group B Trial 9741. J Clin Oncol. 2003;21(11):1431–1439. pmid:12668651

7. Dibrov B, Zhabotinsky A, Neyfakh Y, Orlova M, Churikova L. Optimal Scheduling for Cell Synchronization by Cycle-Phase-Specific Blockers. Math Biosci. 1983;66:167–185.

8. Agur Z. The effect of drug schedule on responsiveness to chemotherapy. Ann Acad New York Sci. 1986;504:274–277.

9. Panetta JC. A mathematical model of periodically pulsed chemotherapy: Tumor recurrence and metastasis in a competitive environment. Bull Math Biol. 1996;58:425–447. pmid:8688836

10. Costa MI, Boldrini JL, Bassanezi RC. Drug kinetics and drug resistance in optimal chemotherapy. Math Biosci. 1995;125:191–209. pmid:7881194

11. Hahnfeldt P, Folkman J, Hlatky L. Minimizing Long-Term Tumor Burden: The Logic for Metronomic Chemotherapeutic Dosing and its Antiangiogenic Basis. J Theor Biol. 2003;220:545–554. pmid:12623285

12. Gardner S, Fernandes M. New tools for cancer chemotherapy: computational assistance for tailoring treatments. Mol Cancer Therapeutics. 2003;2:1079–1084.

13. Swan GW. Role of optimal control in chemotherapy. Math Biosci. 1990;101(2):237–284. pmid:2134485

14. Martin R, Teo KL. Optimal control of drug administration in cancer chemotherapy. World Scientific Publishin; 1994.

15. Coldman A, Goldie J. A stochastic model for the origin and treatment of tumors containing drug-resistant cells. Bull Math Biol. 1986;48(3/4):279–292. pmid:3828558

16. Coldman A, Murray J. Optimal control for a stochastic model of cancer chemotherapy. Math Biosciences. 2000;168:187–200.

17. Day RS. Treatment sequencing, asymmetry, and uncertainty: Protocol strategies for combination chemotherapy. Cancer Research. 1986;46:3876–3885. pmid:3731062

18. Iwasa Y, Michor F, Nowak MA. Evolutionary dynamics of escape from biomedical intervention. Proc Roy Soc Lond B. 2003;270:2572–2578.

19. Komarova NL, Wodarz D. Drug resistance in cancer: Principles of emergence and prevention. PNAS. 2005;102:9714–9719. pmid:15980154

20. Komarova NL. Stochastic modeling of drug resistance in cancer. J Theor Biol. 2006;239:351–366. pmid:16194548

21. Iwasa Y, Nowak MA, Michor F. Evolution of resistance during clonal expansion. Genetics. 2006;172:2557–2566. pmid:16636113

22. Haeno H, Y YI, Michor F. The Evolution of Two Mutations During Clonal Expansion. Genetics. 2007;177:2209–2221. pmid:18073428

23. Comen E, Morris PG, Norton L. Translating Mathematical Modeling of Tumor Growth Patterns into Novel Therapeutic Approaches for Breast Cancer. Journal of mammary gland biology and neoplasia. 2012;17(3-4):241–249. pmid:23011603

24. De Pillis LG, Radunskaya A. A mathematical tumor model with immune resistance and drug therapy: an optimal control approach. Computational and Mathematical Methods in Medicine. 2001;3(2):79–100.

25. de Pillis LG, Radunskaya AE, Wiseman CL. A validated mathematical model of cell-mediated immune response to tumor growth. Cancer Research. 2005;65(17):7950–7958. pmid:16140967

26. Costa M, Boldrini J, Bassanezi R. Drug kinetics and drug resistance in optimal chemotherapy. Mathematical biosciences. 1995;125(2):191–209. pmid:7881194

27. Kohandel M, Sivaloganathan S, Oza A. Mathematical modeling of ovarian cancer treatments: sequencing of surgery and chemotherapy. Journal of theoretical biology. 2006;242(1):62–68. pmid:16580022

28. Gardner SN. Modeling multi-drug chemotherapy: tailoring treatment to individuals. Journal of theoretical biology. 2002;214(2):181–207. pmid:11812172

29. Komarova NL, Wodarz D. Effect of cellular quiescence on the success of targeted CML therapy. PLoS One. 2007;2(10):e990. pmid:17912367

30. Komarova NL, Wodarz D. Stochastic modeling of cellular colonies with quiescence: an application to drug resistance in cancer. Theoretical population biology. 2007;72(4):523–538. pmid:17915274

31. Foo J, Michor F. Evolution of resistance to targeted anti-cancer therapies during continuous and pulsed administration strategies. PLoS Comput Biol. 2009 Nov;5(11):e1000557. pmid:19893626

32. Foo J, Michor F. Evolution of resistance to anti-cancer therapy during general dosing schedules. J Theor Biol. 2010 Mar;263(2):179–88. pmid:20004211

33. Chmielecki J, Foo J, Oxnard GR, Hutchinson K, Ohashi K, Somwar R, et al. Optimization of dosing for EGFR-mutant non—small cell lung cancer with evolutionary cancer modeling. Science translational medicine. 2011;3(90):90ra59. pmid:21734175

34. Simon HU, Haj-Yehia A, Levi-Schaffer F. Role of reactive oxygen species (ROS) in apoptosis induction. Apoptosis. 2000;5(5):415–418. pmid:11256882

35. Koptyra M, Falinski R, Nowicki MO, Stoklosa T, Majsterek I, Nieborowska-Skorska M, et al. BCR/ABL kinase induces self-mutagenesis via reactive oxygen species to encode imatinib resistance. Blood. 2006;108(1):319–327. pmid:16527898

36. Davis W, Ronai Z, Tew KD. Cellular thiols and reactive oxygen species in drug-induced apoptosis. Journal of Pharmacology and Experimental Therapeutics. 2001;296(1):1–6. pmid:11123355

37. Ghibelli L, Coppola S, Rotilio G, Lafavia E, Maresca V, Ciriolo M. Non-oxidative loss of glutathione in apoptosis via GSH extrusion. Biochemical and biophysical research communications. 1995;216(1):313–320. pmid:7488106

38. Hug H, Strand S, Grambihler A, Galle J, Hack V, Stremmel W, et al. Reactive oxygen intermediates are involved in the induction of CD95 ligand mRNA expression by cytostatic drugs in hepatoma cells. Journal of Biological Chemistry. 1997;272(45):28191–28193. pmid:9353266

39. England K, ’Driscoll CO, Cotter TG. ROS and protein oxidation in early stages of cytotoxic drug induced apoptosis. Free Radical Research. 2006;40(11):1124–1137. pmid:17050166

40. Tortora G. Mechanisms of resistance to HER2 target therapy. J Natl Cancer Inst Monogr. 2011;2011(43):95–8. pmid:22043051

41. Diaz LA Jr, Williams RT, Wu J, Kinde I, Hecht JR, Berlin J, et al. The molecular evolution of acquired resistance to targeted EGFR blockade in colorectal cancers. Nature. 2012 Jun;486(7404):537–40. pmid:22722843

42. Lewis PA, Shedler GS. Simulation of nonhomogeneous Poisson processes by thinning. Naval Research Logistics Quarterly. 1979;26(3):403–413.

43. R Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria; 2013. Available from: http://www.R-project.org/

44. Dempsey JL, Seshadri RS, Morley AA. Increased mutation frequency following treatment with cancer chemotherapy. Cancer research. 1985;45(6):2873–2877. pmid:3986814

45. Hofmann WK, Komor M, Wassmann B, Jones LC, Gschaidmeier H, Hoelzer D, et al. Presence of the BCR-ABL mutation Glu255Lys prior to STI571 (imatinib) treatment in patients with Ph+ acute lymphoblastic leukemia. Blood. 2003;102(2):659–661. pmid:12663457

46. Shah NP, Nicoll JM, Nagar B, Gorre ME, Paquette RL, Kuriyan J, et al. Multiple BCR-ABL kinase domain mutations confer polyclonal resistance to the tyrosine kinase inhibitor imatinib (STI571) in chronic phase and blast crisis chronic myeloid leukemia. Cancer cell. 2002;2(2):117–125. pmid:12204532

47. Roche-Lestienne C, Soenen-Cornu V, Grardel-Duflos N, Laï JL, Philippe N, Facon T, et al. Several types of mutations of the Abl gene can be found in chronic myeloid leukemia patients resistant to STI571, and they can pre-exist to the onset of treatment. Blood. 2002;100(3):1014–1018. pmid:12130516

48. Roche-Lestienne C, Laï JL, Darré S, Facon T, Preudhomme C. A mutation conferring resistance to imatinib at the time of diagnosis of chronic myelogenous leukemia. New England Journal of Medicine. 2003;348(22):2265–2266. pmid:12773665

49. Turke AB, Zejnullahu K, Wu YL, Song Y, Dias-Santagata D, Lifshits E, et al. Preexistence and Clonal Selection of MET Amplification in EGFR Mutant NSCLC. Cancer cell. 2010;17(1):77–88. pmid:20129249

50. Maheswaran S, Sequist LV, Nagrath S, Ulkus L, Brannigan B, Collura CV, et al. Detection of mutations in EGFR in circulating lung-cancer cells. New England Journal of Medicine. 2008;359(4):366–377. pmid:18596266

51. Jia P, Jin H, Catherine B, Xia J, Ohashi K, Liu L, et al. Next-generation sequencing of paired tyrosine kinase inhibitor-sensitive and-resistant EGFR mutant lung cancer cell lines identifies spectrum of DNA changes associated with drug resistance. Genome research. 2013;. pmid:23733853

52. DeVita VT. Single agent versus combination chemotherapy. CA: A Cancer Journal for Clinicians. 1975;25(3):152–158.