Evolutionary modeling of combination treatment strategies to overcome resistance to tyrosine kinase inhibitors in non-small cell lung cancer.

Many initially successful anticancer therapies lose effectiveness over time, and eventually, cancer cells acquire resistance to the therapy. Acquired resistance remains a major obstacle to improving remission rates and achieving prolonged disease-free survival. Consequently, novel approaches to overcome or prevent resistance are of significant clinical importance. There has been considerable interest in treating non-small cell lung cancer (NSCLC) with combinations of EGFR-targeted therapeutics (e.g., erlotinib) and cytotoxic therapeutics (e.g., paclitaxel); however, acquired resistance to erlotinib, driven by a variety of mechanisms, remains an obstacle to treatment success. In about 50% of cases, resistance is due to a T790M point mutation in EGFR, and T790M-containing cells ultimately dominate the tumor composition and lead to tumor regrowth. We employed a combined experimental and mathematical modeling-based approach to identify treatment strategies that impede the outgrowth of primary T790M-mediated resistance in NSCLC populations. Our mathematical model predicts the population dynamics of mixtures of sensitive and resistant cells, thereby describing how the tumor composition, initial fraction of resistant cells, and degree of selective pressure influence the time until progression of disease. Model development relied upon quantitative experimental measurements of cell proliferation and death using a novel microscopy approach. Using this approach, we systematically explored the space of combination treatment strategies and demonstrated that optimally timed sequential strategies yielded large improvements in survival outcome relative to monotherapies at the same concentrations. Our investigations revealed regions of the treatment space in which low-dose sequential combination strategies, after preclinical validation, may lead to a tumor reduction and improved survival outcome for patients with T790M-mediated resistance.

Introduction

The advent of pathway-targeted therapies has dramatically influenced cancer research and medical oncology over the past decade. Many of these therapies are capable of inducing an initial response; however, in many cases, the tumor evolves more aggressive and resistant phenotypes over time and eventually becomes insensitive to therapy.[1,2]De novo and acquired resistance to targeted therapies represent a major clinical problem that continues to challenge efforts to delay progression of disease and improve overall survival rates.[3−5] Gaining a better understanding of the evolution of resistance and identifying treatment strategies that alter the penetrance of resistance throughout a tumor are imperative for improving patient outcomes. One powerful approach to address this problem is to use mathematical modeling of the evolutionary dynamics of therapeutic resistance.[6−9] Mathematical models enable systematic exploration of the infinite-dimensional space of potential treatment strategies through variation of parameters such as drug dose, treatment timing, and combination options. Mathematical modeling can also be used to predict optimized treatment schedules based on a variety of biological end points (e.g., maximal time to progression of disease, maximal rate of tumor reduction, minimal probability of resistance, minimal tumor size, or minimal resistant cell frequency) as well as an assessment of the robustness of these biological end points to changes in the schedule and dosing. As such, mathematical modeling narrows down an infinite space of possible treatment strategies to a subset of strategies with the greatest potential that can then be validated in preclinical models before being introduced to patient care.[6,8]

In this study we focus on lung cancer, the leading cause of cancer-related deaths in the United States.(10) Non-small cell lung cancer (NSCLC) accounts for 80% of all lung cancers and consists of three main types: adenocarcinomas, squamous cell carcinomas, and large cell carcinomas. Standard first-line therapy for advanced NSCLC consists of platinum-based chemotherapy and has a modest effect on overall patient survival. Approximately 10–15% of NSCLCs in North America and 30% in Asia harbor mutations in the EGFR kinase domain that trigger activated signaling of the EGFR pathway and frequently result in responses to the EGFR tyrosine kinase inhibitors (TKI) erlotinib and gefitinib.[11−13] The majority of EGFR mutant patients exhibit tumor regression upon EGFR TKI treatment; however, of the 70% that initially respond, all relapse within about one year after initiation of therapy.[14,15] Several mechanisms of acquired resistance to TKIs are responsible for this relapse; in about 50% of cases, the T790M “gatekeeper” mutation in EGFR causes resistance.[16−18] Some data suggest that the T790M mutation may pre-exist the start of therapy in many patients.(19)

Four large phase III trials (TRIBUTE, INTACT 1, INTACT 2, and TALENT) were initiated to evaluate whether concurrent treatment of EGFR TKIs with standard chemotherapy enhances overall survival for advanced NSCLC patients. The results from these trials led to the conclusion that this combination strategy was unable to significantly improve patient survival.[20−22] At the time of these trials, there were no obvious indicators to suggest that combining these therapies would not lead to improved outcomes for patients. After all, previous clinical trials demonstrated that chemotherapy as a single agent prolongs survival of NSCLC patients when compared to placebo, and that those patients who failed first-line chemotherapy and were then administered erlotinib had improved survival relative to those not treated with erlotinib.[23−25] Due to failures of these combination trials and the results of multiple preclinical studies, a strategy for combining erlotinib with standard chemotherapy (i.e., carboplatin and paclitaxel) with sequential pulsing of the two agents was proposed.(26) Recent clinical studies have shown that intermittent dosing of EGFR TKIs with chemotherapy is superior to concurrent dosing.[27−29] This finding suggests that by simply altering the dose and schedule of currently used drugs, the efficacy of combination therapies can be improved. Therefore, quick and cost-effective methods are needed for assessing the potential of a given treatment regime before administering it to patients and before initiating large, expensive, and time-consuming clinical trials.

In this paper, we utilized mathematical modeling in conjunction with quantitative in vitro experiments to identify optimal combination treatment strategies using erlotinib and paclitaxel to prevent or delay resistance to treatment in NSCLC cells. We used an integrative approach to investigate the evolutionary dynamics of a tumor, which are determined by the composition of, and interactions between, sensitive and resistant cells in the presence of treatment. Although several mechanisms of resistance to EGFR TKIs (e.g., MET amplification) have been identified, we focused on addressing the penetrance of preexisting T790M-harboring resistant cells in NSCLC since this mechanism of resistance remains a major clinical problem. Our approach, however, can be generalized to address other mechanisms of resistance and will therefore be of clinical interest for several scenarios of resistance as well as for other cancer types in which resistance to targeted agents arises due to known mechanisms.

Our biological model of sensitive and resistant cells consisted of two NSCLC cell lines, HCC827 and H1975, with different sensitivities to the EGFR inhibitor erlotinib.(30) The HCC827 cell line harbors a mutation in the EGFR tyrosine kinase domain (E746–A750) that renders these cells sensitive to EGFR TKI therapies, whereas the H1975 cells harbor both the L858R and T790M mutations and are resistant to EGFR TKIs. We experimentally identified the parameters of our mathematical model and validated the model predictions with independent experiments. In particular, we generated predictions describing how the tumor composition evolves over time and how it depends on the initial preexisting fraction of resistant cells under various concentrations of the drug. We extended these predictions to study the evolution of resistance under sequential combination strategies using a range of paclitaxel and erlotinib concentrations. We then developed an integrated quantitative framework to identify regions of a combination treatment strategy space that provided optimal outcomes (i.e., eventual elimination of the NSCLC population or maximal delay of resistance in our system). These results serve as examples for the utility of mathematical modeling for defining treatment schedules with favorable outcomes in cell line systems. After appropriate preclinical validation, these suggestions are expected to lead to improved outcomes in the clinic.

Experimental Section

Mathematical Methods

We used a binary nonhomogeneous two-type birth and death process model to represent the dynamics of the TKI-sensitive and TKI-resistant cell populations over time under varying treatment schedules.(31) Similar models have been previously used to study the dynamics of resistance of cancer cells in a variety of settings.[7,8,32−35] Here, we used a variation of these models to investigate the dynamics of penetrance of resistance under sequential pulsed combination therapy with two drugs, erlotinib and paclitaxel. Penetrance, in contrast to emergence, was defined as the outgrowth of preexisting resistant clones. We applied this mathematical model to study combination therapy administered to mixtures of HCC827 and H1975 cells. Under this framework, sensitive and resistant cells proliferate and die at rates that depend on the current therapy being used and its concentration. The net proliferative rates (birth minus death rates) for each cell type were obtained through analysis of the experimental cell line data. The key assumptions of this model are that (1) cell populations are governed by an exponential growth model during each treatment and break phase and, (2) once a cell is born, its lineage is independent of other cells in the population. The second assumption is discussed in the model validation section of the Results. Cell populations were a mix of sensitive and resistant cells at varying ratios, and we investigated the effects of varying the initial composition on the outcome of treatment. We did not consider the mutation of sensitive cells into resistant cells in this model since resistance generated through this mechanism has been studied previously in other settings.[7,8,32−35] We also neglected the possibility of back-mutation from resistant to sensitive cells since this is an event of exceedingly small probability. This assumption of lack of back-mutation is a well-accepted model in the population genetics literature where it is known as the infinite-sites model.(36)

Calculation of Growth Rates

The sole unknown parameter used in our model was the net proliferative rate of each cell type under various drug concentrations. To obtain net proliferative rates, we fit the experimental counts of live cells at various time points under each drug concentration to the exponential growth model described above; see later sections for details of the experimental approach. A linear regression of the log-transformed data was performed to obtain fitted rates at each drug concentration for both cell lines. The experimental data confirmed a good fit with the exponential growth model (data not shown). A nonlinear least-squares regression was performed to yield the relationship between concentration and net proliferative rate for each cell line and drug pair. The net proliferative rates of the sensitive and resistant cells as a function of erlotinib concentration are denoted by the fitted functions λS(x) and λR(x), respectively. Similarly, the net proliferative rates of sensitive and resistant cells as a function of paclitaxel concentration are represented by μS(x) and μR(x), respectively.

Mean Population Dynamics over Time

We considered treatment strategies in which fixed concentrations of erlotinib and paclitaxel were sequentially pulsed in time. In the treatment strategies, t1 and t2 represent the respective lengths of the erlotinib and paclitaxel pulses and c1 and c2 represent the respective concentrations of the two drugs; these concentrations were constant for each pulse. Each treatment cycle consisted of a pair of erlotinib and paclitaxel pulses. The mean size of the sensitive cell population after m therapy cycles is given by S(m(t1 + t2)) = s0 exp[m(t1λS(c1) + t2μS(c2))], and the mean of the resistant cell population is given by R(m(t1 + t2)) = r0 exp[m(t1λR(c1) + t2μR(c2))], where m(t1 + t2) is the time at the end of cycle m. The average population size of sensitive cells at time t between the completion of cycles m and m + 1 is given by S(t) = s0 exp[m(t1λS(c1) + t2μS(c2)) + (t – m(t1 + t2))λS(c1)] if t – m(t1 + t2) ≤ t1 or by S(t) = s0 exp[(m + 1)t1λS(c1) + mt2μS(c2) + ( t – (m + 1)t1 + mt2))μS(c2)] if t – m(t1 + t2) > t1. Analogous formulas hold for the resistant cell population. This formulation was then used to predict the evolution of cell admixtures over time.

Possibility of Elimination of the NSCLC Cell Population

One biological end point for evaluating a treatment schedule is whether the schedule is capable of driving the tumor cell population to extinction. In the context of our model, tumor cell extinction occurred if and only if [S((m + 1)(t1 + t2))]/[S(m(t1 + t2))] – exp[(t1λS(c1) + t2μS(c2))] was less than one or, equivalently, if and only if (t1λS(c1) + t2μS(c2)) was negative.

Choice of Drug Concentrations

We searched a range of low-dose paclitaxel concentrations (0–30 nM) and a low to moderate range of erlotinib concentrations (0–10 μM). These dose ranges were chosen for their clinical relevance.[37,38] Phase I pharmacokinetic studies have shown that the maximum tolerated daily oral dose of erlotinib (150 mg) results in a maximum plasma concentration of 3 μM, and in a phase I/II trial, high doses up to 2000 mg weekly were generally well tolerated. These doses result in a maximum plasma concentration of around 20 μM.(39) For paclitaxel, the maximum plasma concentration achievable after a 24 h infusion of 180 mg/m2 is around 1 μM, and shorter 3 h infusions at the same dose achieve a Cmax of approximately 5 μM.(40) Paclitaxel is given once every 21 days in a standard chemotherapy regimen. However, since we considered schedules using daily paclitaxel at concentrations one to two orders of magnitude lower than these concentrations and alternated with erlotinib at low to moderate doses, we hypothesized that many of the treatment strategies in our schedule space could be well-tolerated in the clinic. This hypothesis is further supported by the results of several combination clinical trials, which combined paclitaxel at 200 mg/m2 and carboplatin with erlotinib 150 mg/day, and have not reported dose-limiting effects.[29,41]

Time to Progression of Disease

The time to progression of disease (POD) is defined as the first time after treatment initiation at which the total population size increases. For some dose combinations of paclitaxel and erlotinib, there did not exist a schedule that led to elimination of the total population. However, it was possible to identify a combination dosing strategy that maximized the time to POD. Under some treatment schedules, the total cell number was initially reduced and resistant outgrowth caused POD; in other cases, the total cell number was never reduced, and in those cases, the time of POD was defined to be zero.

Due to the cyclic nature of a sequentially pulsed dosing strategy, the determination of the time to POD requires comparisons of the population sizes, S(t) + R(t), at the same point of each cycle. Thus we defined the time to POD as inf{m ≥0: S((m + 1)(t1 + t2)) + R((m + 1)(t1 + t2)) ≥ S(m(t1 + t2)) + R(m(t1 + t2))}, which represents a critical point in the total population size. One disadvantage of this discrete formulation is that longer schedules will always be favored, since the time until discovery of POD at the end of a cycle will be delayed if longer treatment cycles are used. Therefore, we introduced a continuous-time approximation of the sensitive and resistant cell populations: if X represents the fraction of time on erlotinib, the total population size at time t is approximated by f(t) = s0 exp[t(λS(c1)X + μS(c2)(1 – X))] + r0 exp[t(λR(c1)X + μR(c2)(1 – X))]. Note that the value of f(t) matches the exact mean population size at the end of each treatment cycle. The use of this approximation provided two benefits: (1) a continuous representation of the dynamics in time enabled the calculation of explicit derivatives and thus an analytical identification of the critical points, and (2) this formulation characterized the time to POD without the aforementioned sampling effects. A consequence of this formulation is that the outcome of a treatment schedule depends explicitly on the fraction of time spent on each drug, i.e., modifying pulse lengths by a common multiplicative constant does not affect the outcome. This finding is consistent with the original model as long as time scales spanning multiple treatment cycles are considered.

The critical point of f(t) provides an estimate of the time to POD and is given as follows:where αS (X,c1,c2) = XλS(c1) + (1 –X)μS(c2) and αR(X,c1,c2) = XλR(c1) + (1 – X)μR(c2) . Values of t* that are negative or possess a nonzero imaginary part were encountered when a treatment schedule did not elicit a decrease in tumor cell numbers; in these cases, the time to POD was set to zero.

Experimental Methods

Materials

The following stains were purchased from Invitrogen: Hoechst 33342 (#H21492), propidium iodide (PI, #P1304MP), CellTracker orange CMTMR (#C2927), and CellTracker green CMFDA (#C7025).

Cell Culture

HCC827 and H1975 cells were purchased from ATCC and cultured in RPMI media supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin solution using standard growth conditions of 37 °C and 5% CO2.

Growth Assays

Approximately 5,000 cells per well were seeded in 96-well black or clear bottom plates (Corning Inc., #3904) under standard culture conditions. The following day, cells were treated with erlotinib (0, 0.1, 1, 10 μM) or paclitaxel (0, 1, 10, 100 nM). Live and dead cell counts were determined using the Cellomics ArrayScan VTI HCS Reader, a high-throughput quantitative imaging system. Briefly, cells were stained with 5 μg/mL Hoechst 33342 and 5 μg/mL PI for 45 min prior to analysis. Average intensity of Hoechst 33342 and PI was determined for each individual cell using the target activation bioapplication of the Cellomics Arrayscan. A final readout of total cell number and percentage of PI positive cells was quantitatively measured within a given well for a given treatment and for a given time. Each condition was performed in replicates of four. For mixture experiments, cells were mixed at respective ratios prior to plating in the 96-well plates with the initial seeding density (5,000 cells/well) kept constant. All data points used in the analysis were taken before any confluence effects were apparent.

Fluorescence Imaging

Cells were labeled using CellTracker probes following manufacturer’s instructions (HCC827 labeled with CellTracker orange CMTMR and H1975 labeled with CellTracker green CMFDA). Cells were seeded at 75,000 cells/well in a 24-well plate either unmixed or mixed at a ratio of 1:1. The following day, 1 μM erlotinib or 10 nM paclitaxel was added and the cells were imaged 48 h post treatment.

Results

Characterization of Growth Kinetics in Unmixed Populations

As a first step toward building our mathematical model, we measured the growth rates of the individual cell populations during drug treatment. Quantitative measurements of the numbers of live and dead HCC827 and H1975 cells during treatment with erlotinib or paclitaxel were determined using the Cellomics Arrayscan (Figures 1 and S1 in the Supporting Information). These data were then utilized to calculate the growth rates of the HCC827 and H1975 cells. H1975 cells grew more rapidly than the HCC827 cells at all concentrations of erlotinib and slightly faster in the absence of treatment. A decrease in the growth rate of HCC827 cells was observed during erlotinib treatment, whereas H1975 cells only exhibited a response to erlotinib at the highest dose of 10 μM (Figure 2A). During paclitaxel treatment, the degree of differential response between HCC827 and H1975 cells was not as dramatic as compared to erlotinib treatment. Both cell lines appeared to be sensitive to the drug; H1975 cells grew slightly slower than the HCC827 cells in the presence of paclitaxel (Figure 2B). These growth rates then served as calibrants for our mathematical model.

Figure 1

Quantification of cell numbers during drug treatments. (A, B) Quantitative measurements of the number of live and dead HCC827 and H1975 cells at 0, 24, and 48 h post erlotinib and paclitaxel treatment were determined from cell images (see Figure S1 in the Supporting Information). Each treatment condition was performed in replicates of four, and the data was displayed as an average ± 2 SD.

Figure 2

Determination of growth rates during drug treatments. HCC827 and H1975 growth rates were calculated using the live cell counts determined in Figure 1 for a time period of 0–48 h at each concentration of (A) erlotinib (0.1, 1, 10 μM) and (B) paclitaxel (1, 10, 100 nM).

Mathematical Modeling Predicted Growth Kinetics of Mixed Populations of Cells

The behavior and possible interactions between sensitive and resistant cells in response to drug-induced perturbations were important parameters to consider as inputs into the model. Using the unmixed growth rates of the HCC827 and H1975 cells, we derived a mathematical model to predict the population dynamics of mixtures of HCC827 and H1975 cells at various ratios (1:1, 1:9, 1:4, 4:1) during administration of different concentrations of erlotinib or paclitaxel. The predictions were then validated experimentally by quantifying the number of live cells in admixed populations with the Cellomics Arrayscan. We found an average relative error between predictions and experimental data of 7.8% and 9.57% for erlotinib and paclitaxel treatments, respectively (Figure 3A,B). In addition, the quantitative growth rate measurements of the mixed populations in response to drugs were corroborated using fluorescent imaging (Figure 3C). HCC827 cells were more sensitive to erlotinib treatment than H1975 cells irrespective of their admixing with the resistant population. These results suggested that the interactions between the HCC827 and H1975 cells mixed at various ratios did not significantly influence their respective growth rates.

Figure 3

Growth kinetics of admixed populations. (A, B) Mathematical model predictions of the growth kinetics of admixtures of HCC827 and H1975 cells at various ratios 1:1, 1:9, 1:4, and 4:1 (1:1 ratio shown here) under different drug treatments of erlotinib or paclitaxel were made using the individual growth rates determined in Figure 2 (denoted as open circles). These predictions were then experimentally validated by determining the total number of cells of the admixed population in the presence of erlotinib or paclitaxel treatment using the Cellomics Arrayscan (denoted as straight lines). The average relative error rate between the model predictions and experimental validation was 7.8% and 9.57% for erlotinib and paclitaxel treatment, respectively. Each treatment condition was performed in replicates of four and the data displayed as an average ± 2 SD. (C) HCC827 cells labeled with CellTracker orange and H1975 labeled with CellTracker green were admixed at a 1:1 ratio and treated with 1 μM erlotinib for 48 h. Fluorescent images were captured at 10× magnification.

Model Predictions of Treatment Outcomes

We next used the validated mathematical model to predict the growth rates of resistant and sensitive cells under various dosing strategies using the fitted growth rate curves from experiments as model parameters (Figure 2). Here we considered population dynamics in the in vitro system; in order to apply these findings to in vivo situations, more detailed kinetics of the pharmacokinetic (PK) effects as well as interactions with the microenvironment, immune system, endothelial cells, etc. are necessary. The formulation of mathematical models of such situations and the experimental determination of quantitative interaction, PK, and growth kinetics are the topic of ongoing work.

Population Dynamics during Constant-Concentration Monotherapy

Figure 4 shows the predicted cell population size and composition as a function of time under constant monotherapy with erlotinib at 0.1, 1, and 10 μM or paclitaxel at 1, 10, and 100 nM. The initial population was 106 cells with a 0.1% fraction of preexisting resistant cells. At 10 μM erlotinib, the overall population size was predicted to initially decline; however, this treatment strategy selected for an expanding resistant subpopulation that led to progression of disease after this initial decline. This behavior is consistent with clinical observations in which 100% of patients who initially respond to erlotinib develop acquired resistance, often mediated by a T790M mutation. We observed that, at a 100 nM dose, paclitaxel elicited a decline in both the HCC827 and H1975 cell numbers. Interestingly, at a low dose of 10 nM, paclitaxel treatment resulted in a net decrease in the H1975 population but a slight increase in the HCC827 population, a profile opposite that of erlotinib.

Figure 4

Mathematical model predictions of cell composition during constant-concentration monotherapy. (A) Total (black line), HCC827 (blue line), and H1975 (red line) population sizes over time under constant monotherapy with erlotinib (0.1, 1, 10 μM). (B) Total, HCC827, and H1975 population size over time under constant monotherapy with paclitaxel (1, 10, 100 nM). Initial population in both panels was considered to be 106 cells with 0.1% preexisting resistant frequency.

The Initial Resistant Fraction Affects the Time until Progression of Disease

We then evaluated the effect of varying the initial resistant fraction on the population dynamics under a sequential alternating schedule. The treatment schedule consisted of 3 μM erlotinib pulsed with 100 nM paclitaxel in a 1:20 ratio. Each pulse pair (erlotinib–paclitaxel) represented one treatment cycle. Given the assumptions of the mathematical model, only the relative fraction of time on each drug determines the outcome; thus, modifying pulse lengths by a common multiplicative constant does not alter the long-term outcome under the assumption that pulse lengths were short enough such that multiple treatment cycles were achieved. Figure 5 displays the predicted tumor composition and size over time for populations starting with 0.01%, 0.1%, and 10% resistant cells. A positive correlation between the initial resistance frequency and time until POD was observed.

Figure 5

Time to progression of disease for varying initial resistant fractions. Model predictions of total (solid) and H1975 (dashed) population size over time for populations starting with 0.01% (black), 0.1% (red), and 10% (blue) initial resistance frequencies. All predictions were conducted under a sequential alternating treatment schedule consisting of 3 μM erlotinib pulsed with 100 nM paclitaxel in a 1:20 ratio.

Optimal Sequential Combination Schedules Delay Progression of Disease

Using our mathematical modeling framework, we explored a large range of sequential combination strategies to identify optimal treatment schedules. We studied a low-dose (0–30 nM) paclitaxel treatment in sequential combination with low to moderate doses of erlotinib (0–10 μM). For each erlotinib and paclitaxel dose pair in this range, we considered the continuum of all possible sequential dosing schedules, where each schedule was identified by the fraction of time during which erlotinib was administered. For each dose pair, we then investigated whether there existed a dosing strategy that resulted in the eventual elimination of both HCC827 and H1975 cells. If no such strategy existed, this dose pair was placed in the “no elimination” region of the treatment space. We found that there were regimens that resulted in the disappearance of the tumor cell population. Figure 6A shows a map of the treatment space in which blue regions identify the region of elimination of the NSCLC cell population and red regions identify the region where NSCLC cell elimination is not possible. Note that our treatment space also included the subset of constant-concentration monotherapies at each dose in the range, since the pulse length of either therapy can be set to zero. For paclitaxel, monotherapy at concentrations over 20 nM led to elimination. If erlotinib was added to the treatment, elimination was achieved at lower doses of paclitaxel.

Figure 6

Optimal sequential combination schedules delay progression of disease. (A) Map of sequential combination treatment space shows regions in which elimination of the NSCLC cell population is possible (blue) and impossible (red). A range of doses for paclitaxel (0–30 nM) and erlotinib (0–10 μM) were considered; each point on the map represents the possibility of elimination of the tumor cell population for a specific pair of dose concentrations. Points were colored blue if there existed a treatment timing schedule with the specified dose pair that achieved overall elimination; if no suitable schedule existed, points were colored red. (B) For each point in the no-elimination (red) region of (A), the best possible outcome (defined as maximal time to progression of disease) was determined. Color in this map represents the maximal time to disease progression. (C) For each point in the no-elimination (red) region of (A), the schedule achieving the optimal outcome was identified. Color in this map represents the fraction of time spent on erlotinib in the optimal schedule. The initial population size in these maps was 106 cells with 0.1% initial resistance frequency.

Although we investigated a large range of dose pairs and a range of alternating schedules, not all are possible to administer in the clinic due to toxicity constraints. In current practice, paclitaxel is dosed once every 21 days at a level resulting in a plasma concentration in the blood of up to 10 μM. The drug is eliminated from the body with a half-life of less than 15 h.(40) Although a concentration of 20 nM, as investigated in Figure 6A, is predicted to result in a plasma concentration well below the acceptable 10 μM, no toxicity data is available to confirm whether a sustained plasma concentration of 20 nM is tolerated in patients. Therefore, it is important to investigate the optimal strategies leading to both elimination of the NSCLC cell population (when tolerated) and the optimal outcomes when elimination is not possible.

In cases when elimination was not possible, we defined the time to POD as the time at which the total tumor size ceased to decrease and began to increase. If a specific treatment strategy failed to reduce the total cell number initially, then the time of POD was defined as zero. Details on how this time was calculated are provided in the Experimental Section. For each dose pair in our treatment space, the best possible outcome in the no-elimination region was identified by the maximal time to POD. In other words, we searched the space of treatment schedules to find the schedule that resulted in maximal time to POD. The map of best possible outcomes (maximal times to POD) is shown in Figure 6B, and the strategy corresponding to this optimal outcome is displayed in Figure 6C. In the latter, treatment strategies were identified by the fraction of time during which erlotinib was administered. We observed that, with certain combinations of low-dose sequential therapies, we achieved a longer time to POD than with higher-dose therapies once every 21 days (Figure 6). In addition, many of these optimal strategies involved very low doses of paclitaxel (around 10–15 nM) and erlotinib (2–10 μM) for approximately equal lengths of time.

Identification of Sequential Paclitaxel–Erlotinib Schedules That Can Achieve Elimination of the NSCLC Cell Population

Finally, we further considered the treatment strategies predicted to result in NSCLC elimination. For each dose pair in the elimination region of Figure 6A, we identified the range of pulse-timing schedules that achieve tumor cell elimination. These ranges are shown in Figure 7A,B, illustrated by the minimum and maximum fraction of time during which erlotinib was administered. For example, consider the point X in Figure 7A which represents the dose pair 6 μM erlotinib and 15 nM paclitaxel. Figure 7A and Figure 7B demonstrate that schedules with a fraction of time on erlotinib between 0.28 and 0.37 led to eradication of the tumor.

Figure 7

Range of sequential paclitaxel–erlotinib schedules that can achieve elimination of the NSCLC cell population. (A, B) Minimal and maximal erlotinib fraction for sequential treatment strategies that achieve cell elimination for each dose-pair in the elimination (blue) region of Figure 6A. The initial population size in these maps was 106 cells with 0.1% initial resistance frequency.

Next we investigated the dynamics of the tumor composition under strategies within this optimal range, starting with an initial population of 0.1%. Figure 8A displays the predicted tumor composition and size over time when 6 μM erlotinib was administered for 28% of the time and 15 nM paclitaxel was given for 72% of the time (the “minimal erlotinib strategy”). Figure 8B shows the analogous dynamics under the strategy when 6 μM erlotinib was given 37% of the time and 15 nM paclitaxel was given 63% of the time. Under the latter strategy, the overall tumor composition decreased more quickly than under the former strategy, even though the minimal erlotinib strategy reduced the resistant population more quickly. Both dosing strategies resulted in the elimination of the sensitive and resistant cell populations.

Figure 8

Tumor evolution under optimal sequential combination schedules vs constant-concentration monotherapy. (A) Dynamics under the optimal sequential combination strategy identified in Figure 7A, using a minimal fraction of time on erlotinib (28% time on erlotinib at 6 μM, 72% of time on paclitaxel at 15 nM). (B) Dynamics under optimal sequential combination strategy identified in Figure 7B, using a maximal fraction of time on erlotinib (37% time on erlotinib at 6 μM, 63% of time on paclitaxel at 15 nM). (C) Total (black line) and H1975 (red line) population sizes over time under constant monotherapy with erlotinib (6 μM). (D) Total (black line) and H1975 (red line) population sizes over time under constant monotherapy with paclitaxel (15 nM).

For comparison, we also investigated the dynamics of tumor composition during administration of each drug alone at these same concentrations. At 6 μM erlotinib, growth of the HCC827 cell population was strongly inhibited with little effect on the H1975 population (Figure 2). At 15 nM paclitaxel, there was a slight decrease of the number of H1975 cells while the HCC827 cell population continued to increase. Figure 8C and Figure 8D display the dynamics of the population under continuous 6 μM erlotinib or continuous 15 nM paclitaxel. Either of these drugs alone at these concentrations resulted in POD in a relatively short amount of time. Notably, continuous erlotinib elicited an initial reduction in tumor size that was much more rapid than that achieved in any of the other strategies evaluated; however, we simultaneously observed a selection for a resistant subpopulation and subsequent rapid POD. This finding reflects the fact that targeting the sensitive population with a strong differential selection pressure may not be the best strategy in the long term. The alternating combination therapy using drugs with opposite weak differential selection profiles such as low-dose paclitaxel with erlotinib ultimately resulted in the elimination of the NSCLC cell population in our model.

Discussion

In this paper, we have presented a combined mathematical modeling and experimental approach to investigate the effects of combination treatment strategies and schedules on the evolution of acquired resistance in non-small cell lung cancer. Our results suggest that optimally timed combination strategies may achieve dramatic improvements in overall outcome over monotherapy with the same drugs and concentration. In addition, we found that applying high doses to achieve the fastest possible tumor reduction rate initially was not always the best strategy in the long term, as this additionally led to a maximal selective pressure, which was rapidly evaded by the acquisition of resistance mutations.

The apparent strength of targeted therapies such as erlotinib is the ability to strongly inhibit populations harboring the specific molecular target, leading to a dramatic tumor size reduction if these sensitive populations comprise the majority of the initial population. However, such strategies may lead to progression of disease (POD) due to the outgrowth of a resistant subpopulation. Using our mathematical model to search the space of schedules, we identified regions on the dose-combination map in which theoretical elimination of the NSCLC cell population was possible. We observed that, as a direct consequence of the growth rate response curve, continuous paclitaxel schedules above 18 nM resulted in eventual elimination of the tumor. However, information is not available about whether a sustained plasma concentration of the drug at this level can be tolerated in patients; thus we proceeded to investigate all regions of the dose-combination map where elimination of the cell population was possible. For example, we found that a sequential schedule combining 6 μM erlotinib for 28–37% of the time with 15 nM paclitaxel for the remaining time was predicted to lead to an eventual elimination of the tumor cells. At these doses, monotherapy with erlotinib resulted in an initial tumor reduction followed by POD, consistent with clinical observations, and monotherapy with paclitaxel resulted in a lack of response.[1,16,18] In contrast, the alternating strategy was predicted to lead to elimination of the cancer cell population. When the erlotinib concentration was increased beyond 6 μM, favorable schedules were found at even lower paclitaxel concentrations and shorter treatment times. For treatment schedules for which NSCLC cell elimination was not possible, we investigated a wide range of scheduling strategies and identified both the optimal outcome, defined as the maximal time until POD, and a schedule that achieved this outcome. We found that the strategies that maximally delayed POD were closer to the center of the dose-combination map with paclitaxel doses between 10 and 15 nM and erlotinib at doses between 1 and 10 μM approximately symmetrically pulsed (∼50% of time on each drug).

These results suggest that sequentially pairing a targeted inhibitor with a cytotoxic drug inducing a weak differential selection pressure on the TKI-sensitive and TKI-resistant populations may lead to a better overall outcome. The inhibitory effect on the resistant cell population by the cytotoxic agent was sufficient to enable the design of an alternating pulsed strategy that controlled and eventually eliminated both cell populations. Although the initial rate of tumor reduction was not as dramatic using these sequential combination strategies as was monotherapy of a molecularly targeted drug, we predicted that sequential combination strategies led to slow tumor elimination rather than POD due to resistance. Thus, the moral of the story from the tortoise and the hare, “slow and steady wins the race”, also seems to apply when designing treatment strategies. We hypothesize that there are many existing cytotoxic therapies that, at low or moderate doses, induce a slight inhibition of the growth of cells with resistance to targeted therapies. At these low doses, these drugs are good candidates for sequential combination trials. We predict that sequential combination therapy will provide a better, less costly alternative to the development of second generation molecularly targeted inhibitors for the resistant cell populations, which are themselves intrinsically vulnerable to additional resistance mechanisms.

Supporting these findings, several human and mouse trials of NSCLC suggested that sequential therapy using a cytotoxic agent and either erlotinib or gefitinib was more effective than monotherapy with either drug or with combination concurrent dosing.[26,27,29,42] It was postulated that doses or timing of such sequential therapy would greatly influence the outcome. Our mathematical model predicted that sequential therapy does indeed provide a better outcome than either therapy alone at the same doses. The overall outcome was sensitive to timing, dose, and initial ratio of sensitive to resistant cells, and we were able to identify the correct balance of pulses to overcome TKI resistance. We predicted, and validated experimentally, that the initial ratio of resistant to sensitive cells influenced the overall time until POD. Identification of noninvasive methods for monitoring the molecular genotypes of tumors throughout a course of treatment would provide valuable real time data that could be used to dynamically update the model to help further guide treatment schedules. More specifically, this information could be used to define the threshold of resistant cells at which a particular treatment strategy would succeed or fail. Several groups are working toward developing such platforms that would routinely analyze circulating tumor cells or tumor DNA in plasma to provide quantitative molecular characterization of tumors.(19)

We recognize that, as with most models, our framework has limitations. First, the two NSCLC lines used in our biological model are not isogenic. However, they were chosen to train our model for several reasons: (1) they carry known EGFR mutations that are observed clinically and confer sensitivity and resistance, respectively, to EGFR TKI therapies; (2) there exists a significant differential growth rate between these two lines during treatment with EGFR TKIs; and (3) we are interested in modeling the penetrance of resistance and therefore are not considering the rate of new mutations that would convert sensitive cells to resistant cells. Another limitation is that the data used as input to the model was derived from an in vitro system. Thus, pharmacokinetic processes present in vivo such as absorption and elimination of the drug were neglected as well as potential drug interactions. Furthermore, interactions with endothelial, mesenchymal, and immune system cells were not considered in our model. However, even with these idealized strategies and rates measured in vitro, we were able to recapitulate key clinical findings. As long as the relative relationships between growth rates in vitro are similar to those in vivo, relative benefits of various scheduling paradigms can be evaluated in our system, and our predictions will provide starting points for preclinical and clinical evaluation of sequential combination therapies. Thus, although the model has limitations in terms of precise clinical predictions, two main conclusions can be drawn: the timing of drug scheduling in combination therapies can have a striking impact on the overall outcome of therapy, and mathematical modeling provides a useful and efficient method to investigate and optimize over the multidimensional space of scheduling strategies. These realizations serve as a starting point for future investigations that will address more complex scenarios arising in vivo as well as additional resistance mechanisms and drugs. Our approach will also be useful for investigating the dynamics of resistance against targeted therapies for other tumor types.