Pharmacokinetic Profiles Determine Optimal Combination Treatment Schedules in Computational Models of Drug Resistance

A major obstacle in cancer research is the development of resistance to anticancer drugs (1, 2). Although, this phenomenon has long been studied, it has gained even more attention after the introduction of targeted therapies and technological advances such as next-generation DNA/RNA sequencing (3–5). 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 treatment. Unfortunately, the initial clinical response to targeted therapies is almost always temporary, as acquired resistance to these drugs frequently develops (2, 3).

Clonal evolution represents a key element for understanding and predicting the dynamics of resistance (5, 6). In this context, cancer cells evolve specific (epi)genetic alterations during tumor growth that lead to the outgrowth of different subclones within a tumor cell population. This intratumor heterogeneity has important negative implications for targeted therapies because somatic mutations continue to occur even in the presence of drugs, leading to increased aggressiveness of the cancer as new drug-resistant cells with improved malignant potential emerge (1). Moreover, resistance might also be due to therapy-induced selection of a small resistant subpopulation of cells that was already present in the original tumor (7).

This realization has given rise to many stochastic mathematical models of genetic resistance, which consider the emergence of these alterations as random events during cell division (8–13). Goldie and Coldman presented some of the earliest work in the 1980s (14, 15), proposing the use of branching process models to study preexisting or acquired resistance to chemotherapy in tumor cell populations. Branching process models are a class of stochastic models that describe the growth and composition of populations made up of stochastically reproducing individuals (16, 17). A subclass of these models is referred to as birth–death processes because each cell experiences a random birth or death event. Multi-type branching processes are convenient for modeling clonal evolution of tumor cells, as during the division process the emergence of random (epi)genetic alterations can give rise to tumor subclones with different fitness (i.e., proliferative capacity or resistance to apoptosis) than their ancestors (18). This approach has inspired other groups to mathematically characterize drug resistance dynamics and to investigate potential administration schedules that could delay the emergence of resistance (19–21). For instance, Hata and colleagues (7) used a branching process to model gefitinib-tolerant and gefitinib-resistant/persistent clones in EGFR-mutant non–small cell lung cancer and showed that acquired resistance can occur either by de novo acquisition of the EGFR-T790M mutation or by expansion of preexisting EGFR-T790M–positive subclones under the selective pressure of gefitinib therapy. Branching process models have also been used to study multi-drug resistance (11, 13) and the dynamics of metastasis formation (12, 22).

We previously developed a comprehensive computational strategy to explore the evolutionary dynamics of heterogeneous tumor cell populations while taking pharmacokinetic and drug interaction effects into account (23–25). This approach consists of a cell-level description of the changes in sensitive and resistant cells over time and in response to treatment, using a birth-death-mutation process to model population dynamics. Our framework represents a departure from previously developed models where the birth and death rates of sensitive and resistant cells were not directly influenced by drug concentrations and were assumed to be constant over time. In addition, the pharmacokinetic processes included in our previous work (20, 24–26) were exponential declines of drug concentrations, a modeling choice that might not hold for all treatments and cell types. Furthermore, in vitro growth inhibition response of cells to a range of different drug concentrations is generally nonlinear, which needs to be incorporated in a flexible framework (27, 28).

In this article, we developed a generalized modeling methodology that accounts for complex pharmacokinetic models and nonlinear effects of drugs on the growth and death kinetics of tumor cells. We integrated this framework in an R package called ACESO (A Cancer Evolution Schedule Optimizer), providing users with an accessible tool to rationally identify optimum single-agent and combination treatment administration strategies for oncogene-driven cancers. We demonstrate the use of ACESO to explore optimum dosing strategies using publicly available data from the Harvard Medical School (HMS) Library of Integrated Network-based Cellular Signatures (LINCS) database (29). In addition, we demonstrate how the new features and flexibility of ACESO can identify optimum administration schedules for erlotinib in non–small cell lung cancer (NSCLC; refs. 20, 21) by extending the previous model to account for central nervous system (CNS) metastases. Within these models, calculating the probability of developing resistance is a standard feature in ACESO. Our multi-scale framework represents a crucial step toward making clinically relevant predictions as it considers the most important aspects governing treatment response and cancer evolution.

The underlying model in ACESO is a multi-type branching process with time-dependent rates, which consists of multiple different types of cells; see the Quick Guide to Equations and Assumptions and Supplementary Data for a detailed description of the model and mathematical equations, respectively. The tumor-initiating cells with oncogene-activating mutations, denoted as type-0 cells, are referred to as sensitive cells because they contain all (epi)genetic alterations necessary for conferring the tumor phenotype, but are sensitive to anticancer therapies. These alterations can be a point mutation, amplification, deletion, inversion, methylation or other epigenetic change, or any other alteration that arises at a certain rate per cell division. Type-0 cells proliferate and die with rates |{b_0}( t )$| and |{d_0}( t )$|⁠, respectively. These rates can be modulated by the cytostatic and/or cytotoxic effect exerted by the drug(s) as a function of their concentrations over time.

Type-0 cells accumulate (epi)genetic alterations at probability |{u_i}$| per cell division, hence at rate |{u_i}{b_{\it 0}}( t )$|⁠, to generate new clones harboring specific resistance mechanisms. This parameter can be a constant value or be modified by the effects of drugs. Newly emerging resistant cell types (denoted as type-i cells, with |i\ = \ 1,2, \ldots $|⁠) are again characterized by their birth and death rates, |{b_i}( t )$| and |{d_i}( t )$|⁠, respectively, which may also exhibit a concentration-dependent profile. Each of these resistant cell types can then accumulate further (epi)genetic alterations and become resistant to additional drugs. Figure 1A shows a schematic representation of a four-type branching process with cross-resistance. The current implementation of ACESO is limited to combinations of two drugs due to the difficulty of discerning the effect of each drug in combinations of more than two agents as well as interaction kinetics between them.

This quantitative approach considers situations in which there is preexisting resistance at the time of tumor diagnosis and treatment initiation, as well as de novo resistance emerging during treatment. In previous work, we derived analytic solutions to calculate the expected number of cells of any type as a function of time and the probability of developing resistance after treatment initiation (24, 25). These results are summarized in the Supplementary Data and allow for faster predictions without the need for stochastic simulations of the evolutionary process. This cell-level description of how heterogeneous cell populations evolve over time is then coupled to pharmacokinetic models, which describe how the drug concentrations change over time (Fig. 1B), and with information quantifying drug–drug interactions (Fig. 1C).

Data from in vitro cell viability and death assays were required to determine the growth and death rate parameters for each cell type under varying drug concentrations. Here, for each treatment condition (i.e., concentration level of one or more drugs) the corresponding rate constants of birth and death were estimated. One example of a viability assay that can be used to characterize cell proliferation is the MTS assay (30), where the number of viable cells is measured over time. This data is then used to identify a net growth rate constant assuming an exponentially growing cell population. The dynamics of tumor cell death can be characterized using apoptosis assays such as Annexin V/propidium iodide (PI) FACS assays (31), where cells with positive Annexin V staining are considered dead. Our framework uses the fraction of dead cells resulting from these assays to determine the baseline death rate, as well as the death rates in the presence of different drug concentrations (see Supplementary Data for a detailed description). Finally, the net birth rate is calculated as the sum of the estimated cell growth and death rates for the different drug concentrations. Other assays for proliferation (e.g., KI67 protein assays; ref. 32) and death (e.g., Caspase activation assays; ref. 33) can also be used for rate estimation. The structure of an ideal database for ACESO is described in the Supplementary Data.

To estimate drug effects on birth and death rates, we perform the following steps: (i) in case of exposure to a single drug, ACESO includes several nonlinear models (see Supplementary Data; ref. 28) that are fitted to the concentration–response data, which provide the values of the drug effect parameters. The selection of the best model is based on the lowest Akaike information criterion (AIC) value (34). (ii) The relationship between estimates of b or d and drug concentrations obtained from drug combination studies is established using nonparametric models (35, 36). These methods do not assume any specific model structure a priori to fit the data, and thus, the birth/death rate versus concentration profiles are described without assuming any mechanism of interaction, providing a large degree of flexibility.

To implement more complex pharmacokinetic models and to simulate different routes of administration and dosing schedules, we integrated the mrgsolve R package (37) in ACESO, which allows for rapid simulations of ordinary differential equation–based pharmacokinetic models. When patient-specific pharmacokinetic data are available, user-defined or compartmental models available in ACESO using the mrgsolve model specification can be used to select the best model describing the drug concentrations over time. If data from different individuals are available, parameter estimates describing the median tendency of the data are reported, as estimation of interindividual variability is not currently supported in ACESO. Selection between models is based on the minimum AIC value and inspection of the individual fits.

ACESO helps to identify whether the combined action of two drugs is synergistic or antagonistic or whether there is no interaction. To quantify the degree of synergy/antagonism between two compounds, the most common approach is to compare the combination effect to a null reference model of no interaction. If the combination response is greater than what is expected from the reference model, the combination is classified as synergistic, while antagonism is defined as the combination producing less than the expected effect. There are several conventional approaches that define different null models to assess drug synergy/antagonism. The Loewe additivity model (38) quantifies a zero-interactive state for the combination of two drugs and defines synergy/antagonism as a combined inhibitory effect that is greater/lower than the sum of the individual effects of the drugs. Highest single agent (HSA), also known as Gaddum noninteraction model (39), is another popular model that defines “independent action” of drugs when the predicted effect of a combination is that of the most effective drug alone. According to this model, any combined effect stronger than the effect of a single drug is called synergism and a weaker effect is antagonism. Additional information related to these models is provided in the Supplementary Data. These two models are used in ACESO to define the pairwise drug combination effects. The Loewe/HSA scores reflect the difference between the measurement and the surface obtained under the no-interaction models, such that values less than zero represent antagonism and values greater than zero represent synergism.

Raw cell viability data used in this analysis can be retrieved from the HMS LINCS database (http://lincs.hms.harvard.edu/db/datasets/20259), where nine small-molecule kinase inhibitors given in combination were chosen for analysis of combinatorial drug sensitivities in the BT-20 triple-negative breast cancer (TNBC) cell line (obtained from ATCC). Briefly, in these experiments, cells were exposed during 72 hours to constant concentration levels of the drugs either as single agents or in combinations. Because of the lack of information regarding drug-resistant cell lines in LINCS, we assumed that these cells had decreased birth and death rate parameters compared with their sensitive counterparts, and hence the proliferation and death rate parameters obtained for the parental BT-20 cell line were divided by different values depending on the drugs being analyzed and tested for sensitivity.

We used cell viability (MTS assay) and apoptosis (Annexin V/PI FACS assay) data from isogenic erlotinib-sensitive and erlotinib-resistant pairs of PC-9 NSCLC cell lines from a previous publication (40). In these experiments, the total numbers of viable and dead cells in culture were counted at 48, 60, and 72 hours after treatment initiation for different erlotinib concentrations and in the absence of the drug. In contrast to the TNBC data above, for the NSCLC data we used had response kinetics to a single drug only and also did not have information on the growth patterns of resistant cell types.

Using the building blocks outlined above, we created a multi-scale computational modeling platform that incorporates tumor evolution, pharmacokinetics, and drug interaction dynamics (Fig. 1D). This computational platform is publicly available as an R package called ACESO. This method can be applied to investigate best dose administration schedules for any combination treatment.

We analyzed the dynamics of heterogeneous TNBC cell lines during combination treatment. We obtained cell viability data resulting from the exposure of the BT-20 TNBC cell line to different concentrations of nine small-molecule kinase inhibitors used as single agents or in combination (see Data availability section, above). We focus our Results section on the discussion of the combination of alpelisib and neratinib; results of analyses using all drug combinations in the database are shown in the Supplementary Data. For the purpose of the model, BT-20 cells sensitive to both drugs are referred to as type-0 cells, whereas cells resistant to one of the drugs are defined as type-1 or type-2 cells.

We first estimated tumor cell proliferation and death rates. To this end, we used an exponential model to characterize cell growth and estimated a net growth rate parameter for each treatment condition. Figure 2A shows the relationship between viable cells measured 72 hours after exposure and the different concentration levels of the targeted therapies alpelisib and neratinib, as well as the net growth rate estimates obtained from these measurements. Because no information regarding apoptosis assays was available to estimate the death rate of type-0 cells (d₀) and assuming that these targeted therapies do not induce cell death per se, a different d₀ value was defined for each concentration of drugs, which was equal to the minimum value needed to ensure that all birth rates of type-0 cells (b₀) were positive regardless of the treatment condition. The values of b₀ for the cells treated with alpelisib and neratinib are shown in Fig. 2A (right). These parameters were fitted using a generalized additive model (GAM) to generate the predicted surface required for the simulation of the evolutionary process (Fig. 2B). GAMs flexibly model the relationship between the b₀ values and each drug concentration by using a number of basis functions that can have a wide variety of shapes (36). The resulting isobolograms highlight points of equal drug contribution to the decrease of sensitive cell birth rates; Supplementary Fig. S1 shows the result of fitting GAMs to the b₀ values obtained under additional drug combinations in our database.

Once drug-resistant cells have emerged from drug-sensitive cells, they follow their own branching process model characterized by their corresponding birth and death rates. We considered situations in which drug concentrations did not have any influence on drug-resistant cell proliferation, and therefore, these cells were only affected by the varying concentrations of the other drug in the combination. The birth rate–drug concentration profiles were then characterized using the nonlinear models for single agents included in ACESO. Figure 2B shows the fitting of neratinib-resistant (type-1) and alpelisib-resistant (type-2) cell birth rates, where an exponential decay model and a five parameter log-logistic function, respectively, best fit the data as indicated by the lowest AIC values (Table 1). The nonlinear models used to characterize the other cell types in our database are also summarized in Table 1 along with the corresponding parameter estimates and SEs. The best-fitting models based on AIC values were exponential, log-logistic, and Weibull functions.

In addition, a value for the probability of emergence of (epi)genetic resistance mechanisms of 10⁻⁷ was defined for cells resistant to one drug, and of 10⁻⁸ for cross-resistance. These values were chosen because the baseline point mutation rate for genetically stable cells is around 10⁻¹⁰–10⁻¹¹ per base per cell division (41) and we considered 1,000 loci in the genome giving rise to resistance; similar rates may apply for the emergence of other resistance mechanisms such as amplifications, deletions, and methylation changes.

We next sought to understand whether drug interactions played an important role in the combination studies. Specifically, we aimed to determine whether the interaction of the drugs causes increased inhibition of the proliferation of drug-sensitive cells compared with an additive effect of the drugs (Loewe additivity principle) or the effect of each drug alone (HSA). The matrices in Fig. 3A and B show the values obtained by computing the differences between the measurements from the combination of alpelisib and neratinib and the surface obtained with the Loewe or HSA models, respectively, such that negative values represent antagonism and positive values synergism. In this case, both models revealed synergistic interactions between the drugs. Interaction matrices for the remaining drug combinations in this study are shown in Supplementary Figs. S2 and S3. The analysis using the Loewe model (Supplementary Fig. S2) reported that the drugs with the greatest synergistic effect were the combination of alpelisib and neratinib, alpelisib with trametinib, and dactolisib with trametinib; whereas the highest antagonism values were found for the dactolisib and lapatinib combination. The matrices from the HSA model (Supplementary Fig. S3), on the other hand, presented additional synergistic effects that were not found with the Loewe model (e.g. alpelisib with saracatinib or selumetinib) and less antagonistic interactions (although the highest antagonism was still found for the dactolisib and lapatinib combination). The discrepancy between the two figures is not surprising as each model has a different definition of synergism/antagonism.

We obtained relevant pharmacokinetic models and parameter estimates from previous studies (42, 43) and used these estimates for the simulations of drug concentration profiles (Fig. 3C; Supplementary Table S1). In the case of neratinib and dactolisib, pharmacokinetic parameters were estimated from the median concentration–time data scanned from the original publication using ACESO (Fig. 3D). For the examples outlined here, we repeated the pharmacokinetic dynamics for each cycle, but more complex models can also be simulated using ACESO.

Once we defined the pharmacokinetic models, the birth/death/mutation rates of the different cell types, and the effect of different drugs concentrations and combinations on these rates, we used ACESO to identify optimum treatment combination schedules. The quality of a schedule was judged by its impact on the total number of different tumor cell types. For these simulations, we considered a population of 1 × 10⁶ sensitive cells and 1,000 drug-resistant cells prior to treatment initiation. Simulations were implemented for each combination of drugs; the findings for alpelisib and neratinib combination are outlined below and the results for other drug combinations are summarized in the Supplementary Data and Supplementary Figs. S4 and S5.

We tested different dosing schedules for alpelisib and neratinib combination therapy (Table 2). For these schedules, we then determined the dynamics of resistant cell populations over time (Fig. 4A) and compared the expected number of tumor cells 30 days after treatment obtained from each regimen (Fig. 4B and E). In the case of alpelisib, any modification from the MTD administration of 300 mg once a day resulted in an increased number of alpelisib-sensitive cells after 30 days. This effect arises due to the large birth rate of alpelisib-sensitive cells in the absence of drug and the rapid clearance of alpelisib, both of which allow tumor cells to grow quickly during treatment breaks (Fig. 4A; regimen 8). In the case of neratinib, administering high-dose pulses with low daily maintenance dosing gave similar results to the MTD schemas, because low concentrations of the drug had a considerable effect on the birth rate of type-2 cells (see Fig. 2) and therefore the schedule involving a high-dose pulse of 720 mg/week with an additional daily dose of 160 mg for the remainder of the week was comparable with the approved schedule of 240 mg once a day (<1% difference in Fig. 4E). The synergistic effect between these two drugs, however, led to the best schedule of 300 mg once a day alpelisib with 240 mg once a day neratinib as it is the one that further decreases the number of type-0 cells and consequently, the total number of cells (43.1% improvement from the worst schedule in Fig. 4B), although the differences with respect to other regimens were small.

We then sought to investigate optimum erlotinib dosing strategies for the treatment of NSCLC using the evolutionary cancer model from (20, 40). NSCLCs that harbor mutations within EGFR gene are sensitive to the tyrosine kinase inhibitor (TKI) erlotinib (44). Unfortunately, all patients treated with this drug acquire resistance, most commonly as a result of a secondary mutation within EGFR (T790M mutation). In addition, approximately one-third of the patients develop CNS metastases after initial response to EGFR TKIs (45).

Here, we extended the original two-type branching process to account for CNS metastases and to show the flexibility of ACESO to handle different scenarios of clonal evolution. As it has been previously suggested, the metastasized cells retain EGFR TKI sensitivity if sufficient drug concentration can be achieved in brain parenchyma or the cerebrospinal fluid (CSF; ref. 45). We thus considered three cell types (Fig. 5A): primary tumor drug-sensitive cells, primary tumor drug-resistant cells, and drug-sensitive metastasized cells—drug-sensitive cells that have disseminated to a new site and have the same proliferation and death rates as the original cells.

We used publicly available data of erlotinib-sensitive and erlotinib-resistant cell birth and death rates (see Data availability section; ref. 40). We then fit the concentration–response curves for the birth rates of each cell type using ACESO. The selected model for erlotinib-sensitive cells was a five parameter logistic function with a cell intrinsic birth rate of 0.04 h⁻¹, an E_max of 0.003 h⁻¹, an EC₅₀ value of 0.76 μmol/L, a Hill coefficient h of 1.5, and an asymmetry factor f of 8.26. A model was also fitted for the influence of erlotinib on the inhibition of drug-resistant cell growth. Although, it was not the model with lower AIC, a linear equation with an intercept of 0.034 h⁻¹ and a slope parameter of −0.0005 (μmol/L·h)⁻¹ was used to model the birth rates of the resistant PC-9 cell line as in the original article. The death rates obtained from the apoptosis assay counts for both cell types were almost identical for the different drug concentrations used in the experiments, and hence, we defined a constant death rate of 0.12 day⁻¹ and 0.06 day⁻¹ for the drug-sensitive and -resistant cell populations, respectively.

We then generated the erlotinib plasma drug concentration versus time profiles using the model. The corresponding parameters were obtained from a pharmacokinetic analysis of data from 1,047 patients with solid tumors (46) instead of using simpler equations for the pharmacokinetics obtained from scarce data points from the literature as in (20). Simulations represented in Fig. 5B were performed considering the administration of a 1,600 mg dose of erlotinib once a week over 1 month of treatment. The simulations demonstrated differences in drug exposure obtained with each model, which had a large impact on the risk of developing resistance (Fig. 5C); depending on the drug clearance, drug concentrations drop faster and thereby cause cells to resume proliferation more quickly. This fact highlights the importance of properly characterizing the pharmacokinetics of drugs as they can have a large influence in selecting optimum dosing strategies. On the other hand, previous studies demonstrated that the concentration of erlotinib measured in the CSF during standard daily dosing of 150 mg was too low (5% of that in plasma according to ref. 47) to prevent EGFR-mutant NSCLC cell proliferation. Thus, we assumed a linear pharmacokinetic process with a partition coefficient of 0.05 to simulate erlotinib concentrations in the metastatic site.

This analysis explores whether a higher dose of erlotinib provides an increased concentration in CNS that may prove to be more efficient in metastatic disease while preventing the development of resistance at the primary site. On the basis of the above considerations, the standard oral daily dosing schedule of 150 mg was compared with a high-dose pulse of 1,600 mg once a week and four different combinations of low-dose continuous and high-dose pulsed strategies (Fig. 5D–F) assuming 1 month of treatment. Additional simulations increasing the daily dosing for erlotinib were not investigated, as daily doses above 200 mg are associated with unacceptable toxicity (48). In contrast, weekly high-dose pulses up to 2,000 mg have been tolerated by patients with NSCLC despite persistent nausea reported in patients receiving higher doses than 1,200 mg (49).

Simulations were performed for cases of (i) an absence of preexisting tumor cells at the metastatic site at the time of tumor diagnosis, and (ii) established CNS metastases at the time of the primary tumor diagnosis, assuming no preexisting resistance in both cases. We used a mutation probability of 10⁻⁸ per erlotinib-sensitive cell division and a dissemination probability to CSF of 10⁻⁹.

The resulting numbers of the different cancer cell populations are shown in Fig. 5D and E. Variation in the number of sensitive cells in the primary tumor site from one schedule to the next was minimal. However, for erlotinib-resistant cells and especially for metastasized erlotinib-sensitive cells, using high-dose pulses considerably improved the results. We conclude that a combined low-dose continuous and high-dose pulsed erlotinib schedule was the most successful dosing strategy at preventing progression in patients with CNS metastases, but did not show significantly delayed emergence of resistance due to the T790M EGFR mutation; however, once resistant cells arise, these schedules were also demonstrated to be more effective at killing the resistant cell population.

In this work, we have built a multi-scale computational framework in R, called ACESO, which incorporates a continuous time multi-type branching process model to investigate the evolution of resistant clones. This model incorporates concentration-dependent birth and death rates, as well as the probability per cell division that a resistance mechanism arises due to an (epi)genetic alteration. The use of mathematical models to search for optimized treatment schedules that delay the appearance of resistance in cancer cells is part of a growing effort to improve clinical trial design for patients with cancer. Here, we provide an open-source tool to contribute to this research by exploring the evolution of heterogeneous tumor cell populations while taking pharmacokinetics and drug interaction effects into account.

There are several examples in the literature introducing different models to describe drug pharmacokinetics (Supplementary Fig. S6), tumor progression, and drug effects (38, 40, 42, 46). However, the implementation of our current goal previously required the use of different platforms for parameter estimation, model simulation or numerical/graphical diagnostics, or programming skills, delaying the application of integrative quantitative analysis of cancer data. The objective of this work was to integrate the tools required to characterize all necessary processes into a single computational framework. The application of our method to individual cancer types requires the availability of data on cell growth and death during different clinically tolerated concentrations of the drugs to parameterize the behavior of sensitive and resistant cancer cell populations. ACESO uses in vitro data because of the difficulty of obtaining such concentration–response curves from in vivo settings. This fact complicates the predictions of clinically relevant time scales for the emergence of resistance because growth kinetics in vitro occur on a time scale orders of magnitude faster than in vivo; however, as long as the relative differences in growth rates are constant between in vivo and in vitro situations or an in vitro–in vivo correlation model is included in the framework, the ranking of schedules based on in vitro rates is applicable to in vivo situations.

In this work, we first applied ACESO to data from the HMS LINCS database. Because of the lack of data regarding apoptosis assays and drug-resistant cell lines, our results serve as an example needing validation before clinical implementation. We also applied ACESO to an extension of a previously used model of EGFR-mutant cancer cells (20, 40). The previous effort identified a low-dose/high-dose pulse erlotinib administration schedule as the predicted best intervention strategy for stage IV NSCLC, which was then tested in a prospective clinical trial (21) at the Memorial Sloan-Kettering Cancer Center (New York, NY; http://clinicaltrials.gov/show/NCT01967095). This trial established the MTD as 1,200 mg for 2 days a week and 50 mg for the remainder of the week. Although, the study was not powered to demonstrate a significant delay of the emergence of T790M-driven resistance due to larger than expected variability in patient pharmacokinetics, a significant reduction in the rate of progression due to CNS disease was observed (0% vs. up to 33% given historic control). Using ACESO, we now expanded the original model (20) to include a population of drug-sensitive metastasized cells that are able to migrate to the CNS; we also used a population pharmacokinetics model (46). ACESO again predicted pulse dosing to be optimal because this strategy allows for increased drug concentrations in the CNS while the continued daily dosing controlled the progression of disease. This regimen is now being studied in a cohort of patients with EGFR-mutant lung cancer with untreated CNS metastases.

The models represented in the current version of ACESO assume exponential cell growth, no spatial effects, and no competing interactions between cells. Therefore, microenvironmental phenomena such as tumor–immune cell interactions or a changing architecture of blood vessels are currently not considered but could, in the future, be incorporated by extending the methodology to more complex models. The model is able to describe any resistance mechanism that arises at a certain rate per cell division; such changes can include (epi)genetic alterations such as point mutations, amplifications, deletions, inversions, methylation of histone modification changes, or any other alterations. To parameterize such models, the kinetics of different cell types and interaction kinetics among them need to be measured in appropriate model systems. Apart from describing a growing population of tumor cells, our approach can describe any kind of cell population or bacterial infection where (epi)genetic changes cause drug resistance. Using ACESO, different strategies can also be tested for robustness due to variability in pharmacokinetic parameters among patients, variable growth and death rates of sensitive and resistant cells, as well as different compositions of the tumor at the start of therapy. We propose use of our methodology for other cancer types and for therapeutic agents at the forefront of clinical development.

I.F. Trocóniz reports grants from MSD, Menarini, and Eli & Lilly outside the submitted work. No potential conflicts of interest were disclosed by the other authors.

I. Irurzun-Arana: Data curation, software, formal analysis, visualization, methodology, writing-original draft. T.O. McDonald: Conceptualization, supervision, methodology, writing-review and editing. I.F. Trocóniz: Supervision, funding acquisition, writing-review and editing. F. Michor: Conceptualization, resources, supervision, funding acquisition, writing-review and editing.

This work was supported by the Dana-Farber Cancer Institute Physical Science-Oncology Center; grant number NIH U54CA193461.

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