The identification of optimal drug administration schedules to battle the emergence of resistance is a major challenge in cancer research. The existence of a multitude of resistance mechanisms necessitates administering drugs in combination, significantly complicating the endeavor of predicting the evolutionary dynamics of cancers and optimal intervention strategies. A thorough understanding of the important determinants of cancer evolution under combination therapies is therefore crucial for correctly predicting treatment outcomes. Here we developed the first computational strategy to explore pharmacokinetic and drug interaction effects in evolutionary models of cancer progression, a crucial step towards making clinically relevant predictions. We found that incorporating these phenomena into our multiscale stochastic modeling framework significantly changes the optimum drug administration schedules identified, often predicting nonintuitive strategies for combination therapies. We applied our approach to an ongoing phase Ib clinical trial (TATTON) administering AZD9291 and selumetinib to EGFR-mutant lung cancer patients. Our results suggest that the schedules used in the three trial arms have almost identical efficacies, but slight modifications in the dosing frequencies of the two drugs can significantly increase tumor cell eradication. Interestingly, we also predict that drug concentrations lower than the MTD are as efficacious, suggesting that lowering the total amount of drug administered could lower toxicities while not compromising on the effectiveness of the drugs. Our approach highlights the fact that quantitative knowledge of pharmacokinetic, drug interaction, and evolutionary processes is essential for identifying best intervention strategies. Our method is applicable to diverse cancer and treatment types and allows for a rational design of clinical trials. Cancer Res; 77(14); 3908–21. ©2017 AACR.

Major Findings

Drugs used in cancer treatment exhibit varied temporal dynamics of absorption and clearance, summarized by their pharmacokinetic parameters. Typically, these dynamics are neglected in evolutionary models of cancer progression under combination therapy. By developing a computational framework to incorporate pharmacokinetic and drug interaction effects, we showed that the complex interplay of dynamics of drugs and various cancer clones significantly affect dosing strategies. Applied to an ongoing phase Ib clinical trial, our framework predicted efficacies of different dosing schedules and suggested new strategies to significantly improve cancer cell eradication and reduce associated toxicities.

Quick Guide to Equations and Assumptions
The Expected Number of Cancer Cells in Different Clones after Treatment Initiation

We model tumor evolution as a multitype branching process (Fig. 1A). After the initiation of combination therapy, the number of cancer cells evolve according to their individual birth and death rates, which are themselves modified as a function of time by the changing drug concentrations (Fig. 1B–E; see Materials and Methods and Supplementary Data for details). As outlined in the Supplementary Data, the expected number of cells at any time after treatment initiation $\langle{n_i}( t )\rangle$ is given by

formula
formula

where $i = 0$ denotes the original cell type containing all (epi)genetic alterations necessary for conferring the cancer phenotype and $i = 1, \ldots ,N$ denotes the $N$ different resistant cell types. The cell-intrinsic birth rates ${b_i}$ are time-dependent as the drug concentration may modulate the rate of cell division. In the simplest case, we consider constant death rates ${d_i}$⁠, but the above expressions are equally valid for time-dependent death rates. The mutation rates are denoted by ${u_i}$⁠, and ${N_i}$ represents the number of preexisting cells of type $i\ $at the time of treatment initiation (Fig. 1F). The quantity $\mathop \sum \nolimits_z {u_z}$ denotes a sum over the mutation rates of all $N$-resistant types. The above expression is derived under the assumption that the waiting time between events (birth or death of cells in various clones) is distributed exponentially, thus amounting to an assumption of Markovian dynamics (18). The expected total number of cells $\mathop \sum \nolimits_{i = 0}^{i = N} \langle{n_i}( t )\rangle$ is calculated to compare the efficacies of six dosing schedules (Supplementary Fig. S1) as functions of time. The expressions above cannot be solved analytically for realistic pharmacokinetic profiles, and hence have to be solved numerically (see Supplementary Fig. S2 for comparisons with exact simulations). We also derived expressions for the expected number of cancer cells when considering cross-resistance to drugs (Supplementary Data).

The Probability of Resistance

If at the time of treatment initiation there exist no mutant cells of any type (Fig. 1G), we can estimate the probability of such a mutation arising later as a function of time. As outlined in the Supplementary Data, the probability that there exists at least one resistant cell of any type (⁠$i = 1, \ldots ,N$⁠) at time T after treatment initiation is given by

formula

where ${P_{ext,i}}( {t,T} )$ is the probability (19) that a clone generated from a single resistant cell of type $i$ produced at time $t$ goes extinct by time $T$⁠.

The Model Incorporating Drug Pharmacokinetics

To describe the time evolution of the drug concentrations in vivo, we use the following simple exponential model

formula

where $C( t )$ is the concentration of the drug at time t, Cmax is the maximum concentration reached after a single drug dose, and k is the clearance rate of the drug. This model well describes the pharmacokinetics of many drugs, especially those that are absorbed on much faster timescales as compared with their clearance times. For such drugs, the absorption kinetics can be neglected and an exponential clearance is a suitable description of the kinetics. Throughout the text, we will refer to Cmax and k as pharmacokinetic or pharmacokinetic parameters.

A Simple Model for Drug Interactions

We developed a simple model for the interaction of two drugs when there are two resistant clones (⁠$N = 2$⁠). Both drugs are considered to decrease the birth rates of the sensitive clone (type-$0$ cells), while each resistant clone is sensitive to either one of the drugs (Fig. 1A). The concentrations of the two drugs $A$ and $B$ are denoted by ${C_A}( t )$ and$\ {C_B}( t )$⁠. The total concentration of drug at any time is ${C_T}( t ) = {C_A}( t )\ + \ {C_B}( t )$⁠. We first defined a variable $Q$⁠, which is the product of the two drug concentrations: $Q( t ) = {C_A}( t ){C_B}( t )$⁠. We then defined a drug interaction parameter $\alpha $ as follows:

formula

where$\ S$⁠, the strength of the drug interactions, is a constant such that $S\ \gt - 1$⁠. As $Q$ is positive at any time $t$⁠,$\ 0\, \lt\, \alpha\, \lt \,1 \,+\, S$⁠. Finally, we coupled $\alpha $ and the drug concentrations to the growth rates of the various clones as

formula

where ${b_{00}},{b_{10}},{b_{20}}$ are the birth rates in the absence of any drug. Coupling $\alpha $ to ${C_T}$ allows for synergy among the two drugs if $\alpha\, \gt \,1$⁠, antagonism if $0\, \lt\, \alpha\, \lt\, 1$⁠, and purely additive interactions if $\alpha = 1$⁠. The strength constant $S$ controls the degree of synergy or antagonism: for $S \gg 0$⁠, the drugs are strongly synergistic while for $ - 1 \,\lt \,S \,\lt \,0$⁠, the drugs are antagonistic. When $S = 0$⁠, we have$\ \alpha = 1$⁠, describing purely additive interactions. The dependence of $\alpha $ on $Q$ as described above introduces concentration dependence to the drug interactions: for a fixed total drug concentration ${C_T}$ and fixed$\ S\, \gt\, 0$⁠, $Q$ is close to zero if either one of the two drugs is present at a very low concentration and hence$\ \alpha \sim 1$⁠. On the other hand, for the same fixed values of ${C_T}$ and$\ S$⁠, if both drugs are present at moderate concentrations, $Q$ is greater than zero and $\alpha\, \gt\, 1$⁠.

As cancer arises due to a process of somatic Darwinian evolution (1–3), researchers have to grapple with the inexorable reality of development of resistance to anticancer drugs (4–6). As heterogeneous populations of cancer cells evolve in the presence of drugs, chance mutations may arise that lead to clones with greater proliferative abilities, or initially undetectable drug-resistant clones can increase in frequency. This evolutionary process allows cancer cells that are fitter in an environment of drugs to evade drug action, eventually leading to treatment failure (5–8). For example, in non–small cell lung cancer (NSCLC) caused by activating mutations in the EGFR, continuous exposure to tyrosine kinase inhibitors like gefitinib and erlotinib invariably selects for mutant clones (9–11). The expansion of these clones leads to treatment failure and necessitates the development of third-generation inhibitors able to inhibit cells driven by specific resistance mechanisms (12, 13).

Many drugs are clinically administered at or close to the MTD using a continuous dosing strategy, with the aim of eradicating the maximum possible number of cancer cells. For example, the clinical standard of care for the drug erlotinib is at the MTD of 150 mg per day administered without any treatment breaks (14, 15). However, theoretical considerations and in vitro experiments suggest the intriguing possibility of significantly increasing the number of cancer cells killed and prolonging the time to resistance, by simply changing the drug dosing strategy (16, 17). Moving beyond single resistance mechanisms and monotherapies, it has now become evident that a multitude of resistance mechanisms may emerge in any particular cancer type (20–23). This necessitates the use of a combination of drugs, making the task of searching for an optimal drug dosing strategy even more challenging than the single drug situation. This realization has given rise to many quantitative models of cancer evolution to determine the effects of various combinations and their administration strategies (8, 17, 24–33). In particular, a number of these studies have explored the idea that general principles may exist for optimally dosing drugs in combination. For example, in a series of studies (24, 34), Goldie and Coldman proposed that an alternating non-cross-resistant chemotherapy regimen is optimal when using two drugs, which was then extended by Chen and colleagues to the case of asymmetric drugs with different efficacies (32). Building on the Goldie-Coldman framework, Day introduced a “worst-drug rule” (29) that proposed administering the less effective drug earlier or more frequently than other drugs in a combination therapy regimen. This rule was later extended to the case of cross-resistant drugs by Katouli and Komarova (30). Lavi and colleagues used a structured population approach to explore the role of heterogeneity and density dependence in multidrug resistance (35), based on previous work by Lorz and others, which treats resistance as a continuous physiologic variable (36). Bozic and colleagues developed a multitype branching process model of cancer evolution and predicted that simultaneous therapy with two drugs is a better alternative than sequential therapy (27).

Although these models for combination therapy have become more sophisticated over time and approach the challenging problem of drug resistance from various viewpoints, none of them take into account the crucial aspect of time evolution of the drugs themselves. These pharmacokinetic parameters (37) determine how fast drugs are absorbed, metabolized, and cleared from an in vivo system (Fig. 1B; Materials and Methods). We have previously shown that accounting for such processes is important for determining optimal dosing strategies in the context of single-agent therapies (2, 17, 38). As the drug concentrations in the blood plasma may significantly influence tumor response, accounting for pharmacokinetic is crucial for predictions of clinically relevant dosing schedules. In addition, when multiple drugs are involved in the treatment strategy, the drugs may interact additively, synergistically, or antagonistically (39–41). Although the importance of drug interactions has been well established in classical methods of drug development and dosing (42), very few evolutionary models have explored its role in shaping the dynamics of cancer. Exceptions exist, but largely in the field of bacterial antibiotic resistance (43).

Here we developed a comprehensive, multiscale computational framework to investigate the effects of drug pharmacokinetics and concentration-dependent drug interactions (see Fig. 1 and Materials and Methods for the basic model), as well as their variability across patients, on the identification of optimum multidrug dosing strategies. With the introduction of drug kinetics, evolutionary parameters like the birth/death rates of various cancer clones become time-dependent (Fig. 1C–E), a significant departure from prevalent modeling approaches of drug combinations. We developed analytic solutions for the expected number of cells in each clone and the probability of resistance, for arbitrarily complex temporal profiles of the drugs. Using this framework, we showed that the best dosing strategy crucially depends on patient variability in pharmacokinetic and evolutionary parameters of clonal growth as well as concentration-dependent synergistic versus antagonistic effects among drugs. Our results also highlight a number of counterintuitive scenarios of combination therapy strategies, which occur as a result of the complex interplay between the pharmacokinetic parameters, population dynamics parameters, and toxicity constraints. Finally, as a demonstration of the applicability of our method in clinical settings, we analyzed the relative efficacies of various dosing strategies currently evaluated in TATTON, an ongoing phase Ib clinical trial of EGFR-mutant lung cancer patients (44). We focused on two drugs, AZD9291 and selumetinib, and identified several easily implementable changes to the current dosing schedules that would lead to a significant increase in cancer cell killing as well as reduced toxicities. Our study highlights the need to account for drug pharmacokinetics and drug interactions, without which clinical use of drugs in combination might be far from optimal.

To quantitatively describe the clonal evolution of cancer cell populations and their responses to therapeutic intervention over time, we developed a novel multiscale computational framework (Fig. 1A; Supplementary Data). Our approach is based on a stochastic model known as a multitype branching process and considers an arbitrarily large number of cell types. The original cell type, denoted by the index i = 0, contains all (epi)genetic alterations necessary for conferring the cancer phenotype but is sensitive to all drugs considered. Cells of this type divide and die at rates b0(t) and d0, respectively; in our model, the birth rate is a function of time as it is determined by the drug concentrations, which may themselves vary over time. Death rates can also be a function of time but are kept constant here for ease of explanation. Type 0 cells may accumulate different mutations at rate ui per cell division, thus generating new clones harboring specific resistance mechanisms, which are again characterized by their birth, death and mutation rates bi(t), di, and ui (Fig. 1A). The birth and death rates can be roughly approximated as the inverse of the average time for cells to undergo the next birth and death event, respectively, from the time they are born. These quantities therefore cannot become negative. However, the difference between the two can become negative, for instance when there is a high concentration of drug in the system and the total cell number decreases with time. The various clones respond to a particular drug each, but are resistant to the other drugs. Each of these cell types can then again mutate to accumulate further alterations, giving rise to cross-resistance (Supplementary Data). This evolutionary approach is coupled with a pharmacokinetic model to form a multiscale description of drug metabolism and cancer evolution (Fig. 1B–E; Supplementary Data). Figure 1C–E depict how the birth rates of individual cell types change over time in response to a particular two-drug regimen, as shown in Fig. 1B. Our quantitative approach describes both situations in which there is preexisting resistance, such that the tumor at diagnosis consists of Ni cells of each type (Fig. 1F), or no preexisting resistance, such that new cell types arise over time according to the evolutionary model specified above (Fig. 1G). We derived analytic solutions for the expected number of cells of any type as a function of time, for an arbitrary number of cell types, mutation rates, and cross-resistance to multiple drugs (Quick Guide to Equations and Assumptions; Supplementary Data), thus allowing rapid predictions and optimizations without the need for detailed stochastic simulations of the evolutionary process. The expected number of cells at the end of a dosing period was used to compare the efficacies of different drug dosing schedules in all figures, except where specifically mentioned. We also derived analytic results for the probability that there exists at least one resistant cell of any type (⁠$i = 1, \ldots ,N$⁠) at time T after treatment initiation (Supplementary Data). For all our analyses, we compared the efficacies of six different dosing schedules of two drugs (Supplementary Fig. S1). The efficacy of a particular dosing schedule was judged either by the expected number of total tumor cells remaining or the probability of resistance development after a certain time. Finally, as these analytic expressions need to be evaluated numerically when time-dependent drug concentrations are considered, we checked our numerical procedures against exact computer simulations (details in Supplementary data; Supplementary Fig. S2).

Figure 1.

Model of cancer evolution accounting for drug pharmacokinetics. A, The sensitive cell (type 0, red) proliferates with rate ${b_0}$ and can mutate to give two different cell types, type 2 (green) and 3 (purple). The type 0 cell is sensitive to both drugs A and B while the other two types are sensitive to either drug A or drug B. B, The drug concentrations change over time, governed by the pharmacokinetic parameters of the drugs. The dosing regimen shown here corresponds to regimen (ii) (see Supplementary Data for details). C–E, The varying drug concentrations in B cause the birth rates of all the cancer cell types to be functions of time. F and G, Preexisting resistant cells (F) as well as de novo mutations (G) coming up during treatment can be handled using our framework.

Figure 1.

Model of cancer evolution accounting for drug pharmacokinetics. A, The sensitive cell (type 0, red) proliferates with rate ${b_0}$ and can mutate to give two different cell types, type 2 (green) and 3 (purple). The type 0 cell is sensitive to both drugs A and B while the other two types are sensitive to either drug A or drug B. B, The drug concentrations change over time, governed by the pharmacokinetic parameters of the drugs. The dosing regimen shown here corresponds to regimen (ii) (see Supplementary Data for details). C–E, The varying drug concentrations in B cause the birth rates of all the cancer cell types to be functions of time. F and G, Preexisting resistant cells (F) as well as de novo mutations (G) coming up during treatment can be handled using our framework.

Close modal

When several drugs are administered at the same time, they may interact (39) such that the combined inhibitory effect on their targets may be more (i.e., synergism) or less (i.e., antagonism) than the sum of their individual effects. This phenomenon has been investigated, generating different methods to quantitatively characterize additivity, synergy, and antagonism (40, 41). However, drug interactions are usually not studied in conjunction with evolutionary models of cancer progression, although some simple models have been explored in the realm of bacterial drug resistance (43). Here, we developed a simple parametric model of drug interactions where “additivity” and “synergy/antagonism” are defined in the sense of the Loewe additivity model (40, 41). This approach allowed us to incorporate drug interactions into our evolutionary framework and explore the combined effects of time-dependent concentration changes and drug interactions on the efficacy of multidrug dosing schedules.

We defined a drug interaction or synergy factor $\alpha ( {{C_A},\ {C_B},S} )$⁠, which is a function of the individual drug concentrations and a single “strength” parameter $S$ (Supplementary Data). The birth rate of sensitive cells is then coupled to $\alpha $ and the total drug concentration, such that $\alpha = 1$ defines additivity, $\alpha \ \gt \ 1$ represents synergistic drug interactions and $0\ \lt \ \alpha \ \lt \ 1$ gives rise to antagonism (SI).

Variability in pharmacokinetic parameters can significantly influence outcomes

To quantitatively explore the effects of tumor cell evolution, pharmacokinetics, and drug interactions, we explored six individual treatment schedules of two anticancer agents (Supplementary Fig. S1); these serve as illustrative examples while our method is applicable to any dosing schedule. We included simultaneous dosing schemes, high-dose/low-frequency schemes, alternating doses with the drug administration times interchanged and also a low-dose continuous (daily dose) strategy. For ease of comparing the outcomes of the different schemes, we maintained two constraints: (i) a constant total concentration of administered drug over the dosing period and (ii) a constant maximum concentration allowed in one dose.

We then investigated the effects of pharmacokinetic variability by comparing combination treatment strategies using two drugs with fast clearance times (Fig. 2A–C), two drugs with slow clearance times (Fig. 2D–F), and one of each (Fig. 2G–I). The drug concentrations over time for dosing schedule (ii) are shown in Fig. 2A, D, and G and the resulting birth rates of sensitive cells are displayed in Fig. 2B, E, and H. Although each scenario investigates the same dosing strategy, the difference in pharmacokinetic parameters leads to a big difference in the temporal trajectories of drug concentrations, which in turn affect the time dependence of the birth rates. Ultimately these differences in pharmacokinetic parameters cause changes in the dosing schedule predicted to lead to optimum efficacy among the six different strategies, as shown in Fig. 2C, F, and I.

Figure 2.

Effect of drug pharmacokinetics on the optimal schedule. The same dosing strategy (ii) shown in Supplementary Fig. S1 is used in the three columns AC, DF, and GI. Results in each column were generated using three different parametrizations of the drug pharmacokinetics, as shown in A, D, and G. The different pharmacokinetic parameters lead to very different time profiles for the birth rate of sensitive cells, shown in B, E, and H. Finally, C, F, and I show comparisons of the outcomes of six different dosing schedules. The best dose (least expected number of cancer cells at T = 20 weeks) has the highest bar, while the worst dose has zero height. Error bars were generated by varying ${b_{00}}$ by 2% of its original value of 0.185. Parameters used are as follows: ${b_{00}} = {b_{10}} = {b_{20}} = 0.185$⁠, ${d_0}$ = 0.14, ${d_1}$ = 0.2, ${d_2}$ = 0.01, ${u_1} = {u_2}$ = 0.001,$\ {N_0} = 5\ \times {10^7},\ {N_1} = {N_2} = {10^4}$⁠. For the first column, the drug clearance rates are ${k_A} = {k_B}$= 3 per week, for the second, they are ${k_A} = {k_B}$= 0.1 per week, and for the final column, they are ${k_A}$= 0.1, ${k_B}$= 3 per week. The drugs are purely additive, hence S = 0 and $\alpha = 1$⁠.

Figure 2.

Effect of drug pharmacokinetics on the optimal schedule. The same dosing strategy (ii) shown in Supplementary Fig. S1 is used in the three columns AC, DF, and GI. Results in each column were generated using three different parametrizations of the drug pharmacokinetics, as shown in A, D, and G. The different pharmacokinetic parameters lead to very different time profiles for the birth rate of sensitive cells, shown in B, E, and H. Finally, C, F, and I show comparisons of the outcomes of six different dosing schedules. The best dose (least expected number of cancer cells at T = 20 weeks) has the highest bar, while the worst dose has zero height. Error bars were generated by varying ${b_{00}}$ by 2% of its original value of 0.185. Parameters used are as follows: ${b_{00}} = {b_{10}} = {b_{20}} = 0.185$⁠, ${d_0}$ = 0.14, ${d_1}$ = 0.2, ${d_2}$ = 0.01, ${u_1} = {u_2}$ = 0.001,$\ {N_0} = 5\ \times {10^7},\ {N_1} = {N_2} = {10^4}$⁠. For the first column, the drug clearance rates are ${k_A} = {k_B}$= 3 per week, for the second, they are ${k_A} = {k_B}$= 0.1 per week, and for the final column, they are ${k_A}$= 0.1, ${k_B}$= 3 per week. The drugs are purely additive, hence S = 0 and $\alpha = 1$⁠.

Close modal

As is evident from Fig. 2C, when both drugs are metabolized rapidly compared with the frequency of dose administration, there is a negligible difference between the dosing strategies, with the best strategy (iv) better than the worst (vi) by only 0.07%. In stark contrast, if at least one of the drugs accumulates due to slow clearance (Fig. 2F and I), there is a 25%–30% difference in outcomes of different strategies. Strategies (iv) and (v), which are both high-dose alternating schedules, reduce tumor cell numbers about 25% more than the low-dose continuous strategy (vi; Fig. 2F). However, with a different set of pharmacokinetic parameters, (v) is no longer a good dosing strategy and proves to be almost 30% worse than the optimum (Fig. 2I). Type 2 cells have the largest net proliferation rate under zero drug conditions. Drug B targets type 2 and type 0 cells, and hence, when it has a fast clearance rate as in Fig. 2G, dosing it first [as in schedule (iii) and (v)] amounts to a period of time with very little drug in the system, thereby allowing all types of cells to proliferate at essentially their zero drug rates and subsequently poor behavior of schedule (v) (Fig. 2I). When the order of the two drugs is switched, as in schedule (iv), drug A builds up, leading to a reduction of type 0 and type 1 cells and good performance of this schedule (Fig. 2I). If both drugs have slow clearance rates, however, as shown in Fig. 2D, the situation is symmetrical and both schedules (iv) and (v) work equally well (Fig. 2F). Finally, the reason schedule (ii) does not work as well as (iv) even though drug A is dosed first in both schedules (Fig. 2F and I) is that the lower concentration in each administered dose in (ii) allows less accumulation of drug as compared with (iv), even though the frequency of dosing in (ii) is higher than that in (iv). For essentially the same reason, schedules (ii) and (iii) do not work as well as schedule (i) in Fig. 2F, as both drugs have slow clearance rates in this example, simultaneous dosing in (i) allows the plasma drug concentration to increase after every dose administration as compared with (ii) and (iii), thus leading to greater cancer cell eradication.

Notwithstanding these large variations in drug pharmacokinetics, Fig. 2C, F, and I shows that it may be possible that a particular schedule is robustly the best way to administer drug combinations. Our analyses uncovered that schedule (iv) is optimum in all three scenarios, which highlights the crucial role evolutionary models could play in selecting optimum dosing strategies, particularly in the presence of interpatient heterogeneity in pharmacokinetic parameters.

To obtain more general insights regarding the role of pharmacokinetics in combination therapies, we then explored other parameter combinations and model assumptions for the same set of six dosing schedules tested above (Supplementary Figs. S3–S6). We explored a number of parameter sets within two model choices: different birth rates of the various clones (Supplementary Fig. S3A and S3B), different death rates (Supplementary Fig. S3C and S3D), different clearance rates of the two drugs (Supplementary Fig. S3E and S3F) for a model in which drugs decrease birth rates (Supplementary Fig. S3A, S3C, and S3E) as well as a model where drugs increase death rates (Supplementary Fig. S3B, S3D, and S3F). We also explored scenarios with different mutation rates of the two resistant types (Supplementary Fig. S4), minimizing the probability of resistance development (as opposed to the expected tumor size) when preexisting mutant cells do not exist at time of therapy initiation (Supplementary Fig. S5) and a model that allows for cross resistance to both drugs (Supplementary Fig. S6). These calculations highlight a universal feature that holds irrespective of the details of the parameter set being used: when at least one of the two drugs has a slow clearance rate such that it accumulates over time in the various dosing schedules being tested, we can expect to see significant differences among the different schedules and hence searching for the optimal schedule is a worthwhile endeavor. The specifics of the particular cell types under investigation and the drugs being administered then determine which treatment strategy results in optimal cancer cell eradication. On the other hand, if both drugs have fast clearance rates and cannot be induced to accumulate over time in the dosing regimens, then there will be only tiny differences in the outcomes of these various dosing schemes. Our work therefore provides a general principle regarding when we can expect different dosing schemes to lead to significantly different outcomes in the two-drug combination therapy scenario.

Drug interaction effects influence treatment efficacy

We then used our quantitative description of drug interactions (Supplementary Data) within our multiscale framework to investigate the effects of additivity (Fig. 3A–E) versus synergy (Fig. 3F–J) between two drugs on the efficacy of different dosing schedules. Figure 3A and F show the three-dimensional surface of drug inhibition in the additive and synergistic scenarios, respectively. To further highlight the nature of the drug interactions, Fig. 3B and G depict the contour lines (isobolograms) of 20%, 50%, and 80% inhibition. The curvature of the contours in Fig. 3A and B are well known signatures of additivity and synergy, respectively, within the framework of the Loewe model (40, 41). Figure 3C and H show how the synergy factor $\alpha $ behaves under additive and synergistic scenarios, while Fig. 3D and I display birth rates of sensitive cells in response to treatment regimen (ii). Finally, the outcomes of the six dosing schedules under the different drug interaction scenarios are shown in Fig. 3E and J. For both figures, the initial condition used was ${N_0} = 5\ \times {10^7},\ {N_1} = {N_2} = {10^4}$⁠, representing a case of preexisting resistant clones. The number of tumor cells decreases with time to varying degrees depending on the dose and drug interaction. In Fig. 3E, the number of tumor cells remaining after 20 weeks of the best drug schedules, (iv) and (v), is approximately $9.2\ \times {10^6}$⁠, which is 25% less than the worst schedule (vi). Thus, a high-dose alternating strategy works best under this scenario of additive drugs (see next section and Fig. 4 for a detailed analysis of the underlying reason). In contrast, when the drugs are highly synergistic with S = 1,000, Fig. 3J shows that the high-dose low-frequency simultaneous scheme (i) works 31% better than the alternating strategies (iv) and (v). This result is in agreement with the intuitive expectation that strongly synergistic drugs work best when dosed simultaneously.

Figure 3.

The optimal dosing scheme depends on drug interactions. A–E show purely additive interactions while F–J show synergistic interactions. A and F show 3D plots of drug inhibition as functions of the two drug concentrations. B and G show isobolograms (contours of equal drug inhibition) generated from the 3D surfaces. The curvature of the isobolograms define additive and synergistic interactions, respectively, according to the Loewe definition of drug interactions. C, For drugs that are purely additive, $\alpha = 1$ for all concentrations of the drugs (Supplementary Data). H, For drugs that are synergistic, $\alpha \gt 1$ and are concentration dependent as described by our model (Supplementary Data). E and J show the time evolution of a tumor starting out with approximately 50 million cells under the six different dosing schemes (i)–(vi). All parameters are identical, except $S$⁠: $S = 0\ $in A–E while $S = 1,000$ in F–J. Schedule (iv) works best in E while schedule I is best in J, representing two very different ways of dosing two drugs. Parameters used are as follows: ${b_{00}} = {b_{10}} = {b_{20}} = 0.185$⁠, ${d_0}$ = 0.14, ${d_1}$ = ${d_2}$ = 0.1, ${u_1}$=0.01, ${u_2}$ = 0.0001, ${k_A} = {k_B} = \ $0.1 per week,$\ {N_0} = 5\ \times {10^7},\ {N_1} = {N_2} = {10^4}$⁠. For the first column, $S = 0$ implies purely additive drug interactions, while for the second column, S = 1,000 indicates strongly synergistic drug interactions. Dosing schedule (ii) was used to generate C, D, H, and I.

Figure 3.

The optimal dosing scheme depends on drug interactions. A–E show purely additive interactions while F–J show synergistic interactions. A and F show 3D plots of drug inhibition as functions of the two drug concentrations. B and G show isobolograms (contours of equal drug inhibition) generated from the 3D surfaces. The curvature of the isobolograms define additive and synergistic interactions, respectively, according to the Loewe definition of drug interactions. C, For drugs that are purely additive, $\alpha = 1$ for all concentrations of the drugs (Supplementary Data). H, For drugs that are synergistic, $\alpha \gt 1$ and are concentration dependent as described by our model (Supplementary Data). E and J show the time evolution of a tumor starting out with approximately 50 million cells under the six different dosing schemes (i)–(vi). All parameters are identical, except $S$⁠: $S = 0\ $in A–E while $S = 1,000$ in F–J. Schedule (iv) works best in E while schedule I is best in J, representing two very different ways of dosing two drugs. Parameters used are as follows: ${b_{00}} = {b_{10}} = {b_{20}} = 0.185$⁠, ${d_0}$ = 0.14, ${d_1}$ = ${d_2}$ = 0.1, ${u_1}$=0.01, ${u_2}$ = 0.0001, ${k_A} = {k_B} = \ $0.1 per week,$\ {N_0} = 5\ \times {10^7},\ {N_1} = {N_2} = {10^4}$⁠. For the first column, $S = 0$ implies purely additive drug interactions, while for the second column, S = 1,000 indicates strongly synergistic drug interactions. Dosing schedule (ii) was used to generate C, D, H, and I.

Close modal
Figure 4.

Detailed analysis of drug and tumor cell population dynamics under three treatment regimens from Fig. 3E — (i), (iv), and (vi). A–C, pharmacokinetics of drugs A, B, and total drug for treatment schedules (i), (iv), and (vi), respectively. As the schedules (iv) and (vi) involve simultaneous dosing of drugs A and B, their curves overlap and the blue curve representing the time course of drug A is not visible. D–F, Temporal evolution of the net proliferation rates of type 0, type 1, and type 2 cells, respectively. Blue curves represent dosing schedule (i), pink schedule (iv), and yellow schedule (vi). G–I, Number of type 0, type 1, type 2 cells remaining after T = 20 weeks, respectively. The color code is as described for D–F. Parameters used are identical to Fig. 3E.

Figure 4.

Detailed analysis of drug and tumor cell population dynamics under three treatment regimens from Fig. 3E — (i), (iv), and (vi). A–C, pharmacokinetics of drugs A, B, and total drug for treatment schedules (i), (iv), and (vi), respectively. As the schedules (iv) and (vi) involve simultaneous dosing of drugs A and B, their curves overlap and the blue curve representing the time course of drug A is not visible. D–F, Temporal evolution of the net proliferation rates of type 0, type 1, and type 2 cells, respectively. Blue curves represent dosing schedule (i), pink schedule (iv), and yellow schedule (vi). G–I, Number of type 0, type 1, type 2 cells remaining after T = 20 weeks, respectively. The color code is as described for D–F. Parameters used are identical to Fig. 3E.

Close modal

Interestingly, our analysis demonstrates that between the two simultaneous strategies, (i) and (vi), the high-dose low-frequency schedule (i) works approximately 25% better than the low-dose continuous strategy (vi). As the synergy factor $\alpha $ is larger for higher concentrations of drugs (Supplementary Data), the high-dose schedule allows $\alpha $ to build up to larger values faster (Fig. 3H) in comparison with the low-dose continuous schedule. As a result, the birth rate of the sensitive cells decreases faster in regimen (i), leading to greater tumor cell death.

Simultaneous administration of drugs may not be optimum even when multiple resistant clones preexist at the time of therapy initiation

Finally, intuition suggests that when different resistant clones exist at the time of treatment initiation, the best dosing strategy should be to use all drugs simultaneously instead of alternating them. However, as the outcome of the drug treatment is a complex function of many factors, the birth and death rates of individual clones, the pharmacokinetics of the drugs, the drug interactions as well as the toxicity constraints imposed on the dosing schedules, only detailed evolutionary analyses can predict the optimum dosing schedule. This expectation is highlighted by the surprising fact that the high-dose alternating strategy (iv) was optimum in Fig. 3E, significantly surpassing the low-dose, continuous simultaneous strategy (vi) and also the other simultaneous strategy (i). Note that the initial conditions were ${N_0} = 5\ \times {10^7},\ {N_1} = {N_2} = {10^4}$ cells, which models preexistence of both types of resistant clones at the time of treatment initiation.

We now dissect the reasons for the surprising optimality of dosing schedule (iv). Figure 4A–C depict the concentrations of the two drugs as well as the total drug over the dosing period in schedules (i), (iv), and (vi), respectively. Although the total amount of drug administered is identical in each of the three schedules over the entire period of dosing, the higher concentration pulses in schedule (iv) (pink curve in Fig. 4B) initially allow for a more rapid build-up of the total drug concentration as compared with schedule (vi) (yellow curve in Fig. 4C). As a consequence, the net proliferation rate of type 0 cells under schedule (iv) (pink curve in Fig. 4D) becomes negative, implying that the number of these cells starts to reduce. On the other hand, the net proliferation rate of type 0 cells under schedule (vi) (yellow curve in Fig. 4D) initially remains positive (although smaller than the zero drug scenario), implying that these cells still increase in number albeit at a slower rate. As time progresses and the lower dosing frequency of schedule (iv) causes greater variability in drug concentration as compared with schedule (vi), the net proliferation rate of type 0 cells under the former regimen transiently becomes larger than the net proliferation rate under the latter regimen (pink and yellow curves in Fig. 4D). Nevertheless, the net proliferation rate remains more negative under regimen (iv) for longer periods of time, as is evident from a comparison of the pink and yellow curves in Fig. 4D. This argument applies to type 1 cells as well (Fig. 4E), but not type 2 (Fig. 4F, see below). This reason ultimately leads to a significantly smaller number of type 0 and type 1 cells at the end of the dosing period under schedule (iv) as compared with schedule (vi) (Fig. 4G and H, respectively).

A similar argument applies to a comparison of schedules (iv) and (i), although in this case the two schedules are much closer to each other as can be seen from the pink and blue curves in Fig. 4D and pink and blue bars in Fig. 4G. As mentioned above, this argument does not hold for the type 2 cells: schedule (i) does marginally better than (iv) in terms of cell kill, as the latter involves starting the dosing with drug A (Supplementary Fig. S1), which does not affect type 2 cells. This results in an initial period when the net proliferation rate of type 2 cells is significantly smaller when using schedule (i) (blue curve in Fig. 4F) as compared with schedule (iv) (pink curve in Fig. 4F), leading to a marginally lower number of type-2 cells after the end of the dosing period (Fig. 4I). In all, as type 0 and type 1 cells dominate the total tumor population, the high-dose alternating strategy (iv) works best for reducing the total tumor burden. This nonintuitive result would not have been easily deducible without a thorough mathematical treatment.

An in silico clinical trial characterizes effects of interpatient pharmacokinetic variability

To explore the effects of heterogeneity among patient pharmacokinetic parameters, we then used our framework to calculate the expected number of cancer cells as a function of time for a cohort of 20 patients with drug clearance rates varying between 0.1 and 0.3 per week (6 of the 20 patient tumor trajectories are shown in Fig. 5A), for the dosing schedule (i) shown in Supplementary Fig. S1. In each trajectory, there initially are $5\ \times {10^7}$ sensitive cells and ${10^4}$cells of each resistant type, thus mimicking conditions of preexisting resistance. Note that while the different cell types have different birth/death kinetics (Fig. 5), they remain identical between different patients. Only the clearance rates of the drugs are allowed to vary from patient to patient, which therefore drives the observed interpatient heterogeneity in treatment response. As drugs A and B are administered (Fig. 1A), the tumor rapidly shrinks as sensitive cells decrease in number. However, the selection pressure exerted by the drugs leads to fitness advantages for the resistant clones, thus causing a rebound in each trajectory. Fig. 5B displays the distribution of time to progression for the cohort, highlighting the variability in patient responses as a result of heterogeneity in pharmacokinetic parameters.

Figure 5.

Interpatient pharmacokinetic variability can lead to different clinical outcomes. A, Six example tumor growth trajectories are shown, highlighting tumor rebound after initial decrease in size for different pharmacokinetic parameters. The two drugs have identical clearance rates (⁠${k_A} = {k_B}$⁠) in one trajectory and vary between 0.1 (bottom curve) to 0.3 (top curve) per week between the different trajectories. B, Distribution of time to progression, defined as the time at which the expected number of tumor cells reaches a minimum. Twenty patient trajectories were calculated, with the drug clearance rates varying between 0.1 and 0.3 per week. Dosing schedule (i) was used for A and B and parameters used are: ${b_{00}} = 0.185,\ {b_{10}} = 0.3,\ {b_{20}} = 0.2,\ {d_0} = 0.14,\ {d_1} = 0.01,\ {d_2} = 0.01,\ {u_1} = 0.01,\ {u_2} = 0.0001,\ {N_0} = 5\ \times {10^7},\ {N_1} = {N_2} = {10^4}.$C, An extended version of the model in Fig. 1A, with four resistant cell types and drugs. Type 0 cells are sensitive to all drugs while all other types are sensitive to only one drug each. D, Clonal composition for different patient pharmacokinetic parameters, in the presence of four drugs. The dosing schedule used is 0.02 units of each drug A, B, C, and D administered simultaneously on days 0, 4, 8, 12, and 16. PK-1 and PK-2 correspond to drug clearance rates of 0.1 and 0.2 per week, respectively. The different colors correspond to the different cell types as shown in C. Other parameters are: ${b_{00}} = 0.185,\ {b_{10}} =$$ {b_{20}} = {b_{30}} = {b_{40}} = 0.35,\ {d_0} = 0.14,\ {d_1} = 0.2,\ {d_2} = 0.01,\ {d_3} = 0.1,\ {d_4} = 0.14,\ {u_1} = {u_2} = {u_3} = {u_4} = 0.001,\ {N_0} = 5\ \times {10^7},\ {N_1} = {N_2} = {N_3} = {N_4} = {10^4}.$

Figure 5.

Interpatient pharmacokinetic variability can lead to different clinical outcomes. A, Six example tumor growth trajectories are shown, highlighting tumor rebound after initial decrease in size for different pharmacokinetic parameters. The two drugs have identical clearance rates (⁠${k_A} = {k_B}$⁠) in one trajectory and vary between 0.1 (bottom curve) to 0.3 (top curve) per week between the different trajectories. B, Distribution of time to progression, defined as the time at which the expected number of tumor cells reaches a minimum. Twenty patient trajectories were calculated, with the drug clearance rates varying between 0.1 and 0.3 per week. Dosing schedule (i) was used for A and B and parameters used are: ${b_{00}} = 0.185,\ {b_{10}} = 0.3,\ {b_{20}} = 0.2,\ {d_0} = 0.14,\ {d_1} = 0.01,\ {d_2} = 0.01,\ {u_1} = 0.01,\ {u_2} = 0.0001,\ {N_0} = 5\ \times {10^7},\ {N_1} = {N_2} = {10^4}.$C, An extended version of the model in Fig. 1A, with four resistant cell types and drugs. Type 0 cells are sensitive to all drugs while all other types are sensitive to only one drug each. D, Clonal composition for different patient pharmacokinetic parameters, in the presence of four drugs. The dosing schedule used is 0.02 units of each drug A, B, C, and D administered simultaneously on days 0, 4, 8, 12, and 16. PK-1 and PK-2 correspond to drug clearance rates of 0.1 and 0.2 per week, respectively. The different colors correspond to the different cell types as shown in C. Other parameters are: ${b_{00}} = 0.185,\ {b_{10}} =$$ {b_{20}} = {b_{30}} = {b_{40}} = 0.35,\ {d_0} = 0.14,\ {d_1} = 0.2,\ {d_2} = 0.01,\ {d_3} = 0.1,\ {d_4} = 0.14,\ {u_1} = {u_2} = {u_3} = {u_4} = 0.001,\ {N_0} = 5\ \times {10^7},\ {N_1} = {N_2} = {N_3} = {N_4} = {10^4}.$

Close modal

Finally, an example scenario of administering an identical schedule of four drugs, A–D (Fig. 5C and D), exemplifies the general nature of our approach. As there is increased interest in the clinic in increasing the number of drugs in combination therapies, Fig. 5C and D is an illustration of how our model can be extended to such cases. Figure 5D displays the clonal size composition of the tumors after 20 weeks of drug administration, for two pharmacokinetic scenarios: all drugs have a clearance rate of (i) 0.1 (PK-1) or (ii) 0.2 per week (PK-2). As evident in both scenarios, the resistant clones undergo rapid clonal expansion and by the end of 20 weeks reach numbers comparable with the sensitive clone, causing a tumor relapse. Although the dosing schedule is the same for both scenarios, the faster rate of drug clearance in scenario (ii) leads to a poorer outcome in terms of tumor volume.

Identification of superior dosing strategies for the ongoing TATTON clinical trial

Having established our general computational framework in the previous sections, we now demonstrate how it can be used to predict outcomes of dosing regimens using clinical and preclinical data on two drugs, AZD9291 and selumetinib, for the treatment of NSCLC. AZD9291 is a third-generation EGFR inhibitor that is effective against both cells carrying the EGFR tyrosine kinase inhibitor (TKI)-sensitizing alterations as well as cells harboring the EGFR T790M–resistant mutation (45, 46). Compared with other EGFR TKIs, AZD9291 was shown to have significantly less toxicity (47), and hence was combined with the MEK inhibitor selumetinib to decrease the chance of development of resistance. Selumetinib reduces cell proliferation by inhibiting the RAS–MAPK pathway (46), amplification of which is a frequent cause of resistance in NSCLC cell lines. Here, we combined growth kinetics data from in vitro cell culture experiments and pharmacokinetic data from clinical trials to predict outcomes of the various dosing schedules tested in TATTON, “A Multi-arm, Phase Ib, Open-Label, Multicentre Study to Assess the Safety, Tolerability, Pharmacokinetics and Preliminary Anti-tumour Activity of AZD9291 in Combination With Ascending Doses of Novel Therapeutics in Patients With EGFRm+ Advanced NSCLC Who Have Progressed Following Therapy With an EGFR TKI,” a currently ongoing phase Ib clinical trial (44). Our model suggests that simple modifications to the currently employed dosing strategies have the potential to significantly improve clinical outcomes.

We first used data from the PC9 and NCI-H1975 NSCLC cell lines (Fig. 4A and B in Eberlein and colleagues; ref. 46) to parameterize growth and death rates of cell types 0, 1, and 2 in the presence and absence of drugs. The method is given in detail in the Supplementary Data and the results are shown in Fig. 6A and B. Type 0 cells proliferate significantly faster than the other cell types when no drug is present (Fig. 6A). The net growth rates of type 0 and type 1 cells were also quantitatively consistent with the results of a previous study (17). As the earlier study had also characterized the birth and death rates of these two cell types separately, we used the same values for the death rates as estimated in (17) (⁠${d_0} = 0.12\ da{y^{ - 1}},\ {d_1} = 0.06\ da{y^{ - 1}}$⁠) to calculate birth rates of the various clones. Under the influence of 160 nmol/L AZD9291, type 0 cells die faster than they grow, and hence have a net negative growth rate (Supplementary Fig. S7). The net growth of type 1 cells slows down when treated with AZD9291, but remains positive (Supplementary Fig. S7). On the other hand, as type 2 cells are insensitive to this drug, they have a net proliferation rate that is higher than that of type 1 cells, explaining the observation that KRAS or NRAS mutations, but not T790M mutations, are observed in all parental PC9 populations (type 0 cells) treated with AZD9291 (46). Finally, the data showed that a combination of 160 nmol/L AZD9291 and 100 nmol/L selumetinib completely prevented the emergence of resistance during the period of observation (46). This finding implies that 100 nmol/L selumetinib reduces the net growth rate of type 2 cells to at least 0 and possibly to negative values (see Supplementary Data for further discussions on this point).

Figure 6.

Analysis of clinical and preclinical data on AZD9291 and selumetinib and prediction of improved dosing schedules for the TATTON phase-Ib trial. A, Linear regression of in vitro growth assay data provides net growth rates of the three cell types. B, Estimated birth rates of the three cell types as a function of drug concentration (AZD9291 for type 0 and 1; selumetinib for type 2 cells). C and D, Pharmacokinetics after multiple doses of the C2 and I dosing strategies of the TATTON trial, respectively. Both schedules involve 80 mg once a day daily AZD9291 (pink). In the case of selumetinib, the doses are 75 mg twice a day (green) treatment (C) and 75 mg twice a day (D), 4 days a week (green) treatment. E, Predicted reduction in the total number of cancer cells for the original and proposed dosing strategies C2 and C2*, respectively, at three time points after the start of treatment. F, Predicted improvement in outcomes of our proposed dosing strategies C1*, C2*, and I* compared with the original strategies. This metric of comparing the effects of original and proposed dosing schedules is invariant with time, unlike measures shown in E. Our proposed dosing strategy involves administering AZD9291 twice a day (instead of once) and selumetinib once a day (instead of twice). The days of dosing remain identical to the original schemes used in TATTON. G, Other dosing strategies tested to investigate saturating benefits of increasing AZD9291 concentration. Selumetinib was administered at 75 mg per dose in each of these schedules. Initial conditions in E–G are ${N_0} = 50\ \times {10^6},{N_1} = 1\ \times {10^6},{N_2} = {10^4}$ cells. Other parameters used (see Supplementary Data for details) are ${b_0} = 0.77 - 0.0041[ {{\rm{AZD}}9291} ],{b_1} = 0.414 - 0.0019[ {{\rm{AZD}}9291} ],{b_2} = 0.294 - 0.00234[ {{\rm{selumetinib}}} ],\ {d_0} =$$ 0.12\ {\rm{da}}{{\rm{y}}^{ - 1}},\ {d_1}\ = \ {d_2}\ = {\rm{\ }}0.06\ {\rm{da}}{{\rm{y}}^{ - 1}},\ {k_{{\rm{AZD}}9291}} = 0.302\ {\rm{da}}{{\rm{y}}^{ - 1}},\ {k_{{\rm{selumetinib}}}} = 3.538\ {\rm{da}}{{\rm{y}}^{ - 1}}.$

Figure 6.

Analysis of clinical and preclinical data on AZD9291 and selumetinib and prediction of improved dosing schedules for the TATTON phase-Ib trial. A, Linear regression of in vitro growth assay data provides net growth rates of the three cell types. B, Estimated birth rates of the three cell types as a function of drug concentration (AZD9291 for type 0 and 1; selumetinib for type 2 cells). C and D, Pharmacokinetics after multiple doses of the C2 and I dosing strategies of the TATTON trial, respectively. Both schedules involve 80 mg once a day daily AZD9291 (pink). In the case of selumetinib, the doses are 75 mg twice a day (green) treatment (C) and 75 mg twice a day (D), 4 days a week (green) treatment. E, Predicted reduction in the total number of cancer cells for the original and proposed dosing strategies C2 and C2*, respectively, at three time points after the start of treatment. F, Predicted improvement in outcomes of our proposed dosing strategies C1*, C2*, and I* compared with the original strategies. This metric of comparing the effects of original and proposed dosing schedules is invariant with time, unlike measures shown in E. Our proposed dosing strategy involves administering AZD9291 twice a day (instead of once) and selumetinib once a day (instead of twice). The days of dosing remain identical to the original schemes used in TATTON. G, Other dosing strategies tested to investigate saturating benefits of increasing AZD9291 concentration. Selumetinib was administered at 75 mg per dose in each of these schedules. Initial conditions in E–G are ${N_0} = 50\ \times {10^6},{N_1} = 1\ \times {10^6},{N_2} = {10^4}$ cells. Other parameters used (see Supplementary Data for details) are ${b_0} = 0.77 - 0.0041[ {{\rm{AZD}}9291} ],{b_1} = 0.414 - 0.0019[ {{\rm{AZD}}9291} ],{b_2} = 0.294 - 0.00234[ {{\rm{selumetinib}}} ],\ {d_0} =$$ 0.12\ {\rm{da}}{{\rm{y}}^{ - 1}},\ {d_1}\ = \ {d_2}\ = {\rm{\ }}0.06\ {\rm{da}}{{\rm{y}}^{ - 1}},\ {k_{{\rm{AZD}}9291}} = 0.302\ {\rm{da}}{{\rm{y}}^{ - 1}},\ {k_{{\rm{selumetinib}}}} = 3.538\ {\rm{da}}{{\rm{y}}^{ - 1}}.$

Close modal

The last step in parameterizing our model is to obtain the pharmacokinetic parameters of the two drugs, AZD9291 and selumetinib. As phase I clinical trials have already been performed on patient cohorts with both drugs, we obtained pharmacokinetic parameters from the literature. We used ${C_{max}}$ (the maximum observed plasma concentration) and ${t_{1/2}}$ (the elimination half life) values after a single dose for both drugs, as this allowed us to predict the pharmacokinetic time course for the entire period of drug dosing (see Supplementary Fig. S8 for details).

Three dosing strategies are being used in the AZD9291 + selumetinib arm of the TATTON trial (44). In all three strategies, 80 mg AZD9291 is administered once daily. For selumetinib, either 50 or 75 mg is administered twice daily in two continuous dosing strategies that we call “C1” and “C2,” respectively, for convenience. In the third, intermittent strategy, referred to as “I,” selumetinib is administered at 75 mg twice daily for four days a week, accompanied by three days of drug holiday. Figure 6C and D show the pharmacokinetic time course of AZD9291 and selumetinib in $C2\ $and $I\ $as predicted from the respective single-dose parameters.

As the total drug administered in $I\ $is less than that in either $C1\ $or $C2$⁠, and as the drug holidays in $I$ allow the selumetinib concentration to decrease to zero periodically (Fig. 6D), we first investigated whether the intermittent strategy would be significantly less effective in tumor shrinkage as compared with the continuous schemes. The results shown in Supplementary Fig. S9 predict that there is very little difference between the dosing schedules C1, C2, and I—drug holidays (the three days when selumetinib is not dosed) in schedule I alter the reduction of the total number of tumor cells by less than $0.5\% $⁠. The three drug holidays of selumetinib only make a small difference, in absolute terms, in the number of type 2 cells remaining at the end of the treatment period. The selumetinib concentration becomes almost zero during the drug holiday, which allows type 2 cells to start proliferating. As a result, at the end of the entire intermittent dosing period, there are approximately 26,000 type 2 cells, whereas there are approximately 2,800 type 2 cells when the continuous scheme is used (assuming there are 10,000 type-2 cells to begin with in both cases; see Supplementary Data for a discussion on initial clone sizes). Therefore, the intermittent scheme is worse than the continuous scheme in absolute terms. However, relative to the size of the entire tumor, this increase from 10,000 to 26,000 type 2 cells is negligible. As selumetinib is a highly toxic drug and drug holidays are essential for patients, our result suggests that these holidays may not make much difference to the entire tumor dynamics.

The concentrations of selumetinib achieved in vivo (Fig. 6C and D) are very high compared with what is required to reduce the birth rate of type 2 cells. As is evident from Fig. 6B, ${b_2}$ is reduced to zero in vitro by treatment with as little as 125 nmol/L selumetinib (blue line); hence higher concentrations do not decrease the net proliferation rate of type 2 cells any further. On the basis of this observation, we hypothesized that reducing the administered selumetinib dose should not significantly alter the extent of tumor reduction, but would be desirable as the drug is associated with significant toxicities—the MTD is 100 mg twice daily. On the other hand, clinical trials (47) have shown the absence of any dose-limiting toxicity even at 240 mg daily dose of AZD9291, and the in vivo concentrations of AZD9291 achievable are comparable with the active range of the drug (Fig. 6B and C). Thus, we also hypothesized that increasing the AZD9291 dose would be both feasible from a toxicity standpoint and could also increase tumor cell kill. We therefore investigated three modified dosing schedules, C1*, C2*, and I*, where AZD9291 is dosed twice a day instead of once, and selumetinib is reduced from twice to once a day. The days of administration of these new doses remain identical to the original doses. The results of our proposed schedules are shown in Fig. 6E and F. The total amount of cancer cell killing achieved by our proposed dosing strategies (C1*, C2*, and I*) is higher than for the original dosing schedules, for three arbitrary time points after the initiation of dosing (Fig. 6E). The uncertainty regarding the proliferation rate of type 2 cells in the presence of selumetinib does not significantly affect the results (see Supplementary Fig. S10).

However, the percentage reduction in cancer cells from baseline is not a good measure to characterize the difference between the original (C1, C2, and I) and proposed (C1*, C2*, and I*) dosing strategies, as the difference changes with time. A better (time-invariant) measure to judge the efficacy of our proposed dosing schemes is the percentage decrease in the total cancer cell number between the original and proposed dosing schemes (Fig. 6F). The proposed dosing schedule is about 25% better than the dosing schemes being currently used in TATTON. Varying the initial clone size distributions predicts a range of about 20%–26% improvement (Supplementary Fig. S11).

Given that selumetinib has been shown to have dose-limiting toxicities, we hypothesize that reducing selumetinib doses is preferable for the management of side effects and the adherence of patients to the recommended schedule. As the drug holiday of three days per week is predicted to not significantly affect the number of cancer cells (Fig. 6F; Supplementary Figs. S9 and S10), the clinical suggestion based on our analyses is to use schedule I* as a better alternative to the ones currently implemented in TATTON. Given that 80 mg AZD9291 and 75 mg selumetinib tablets are already available clinically, our proposed modifications are easily implementable and would not require any additional manufacturing.

Finally, we sought to identify the highest doses of AZD9291 that would significantly increase the efficacy of this treatment. We found that the improvement due to drug saturates at about 160 mg once-daily AZD9291 as shown in Fig. 6G (interestingly, 160 mg once daily is about 10% better than 80 mg twice daily)—increasing the administered concentration any further does not significantly improve the outcomes any more, indicating that at these concentrations AZD9291 accumulates to a higher concentration than is required for it to be active. This observation suggests that drugs need not necessarily always be dosed at their MTD levels—lower doses could be equally efficacious as MTD dosing strategies, while leading to lower side effects and associated toxicities.

Here we have developed a framework for exploring the effects of pharmacokinetics and concentration-dependent drug interactions on the outcome of multidrug dosing schedules within a multiscale evolutionary medicine framework. Our analyses uncovered important and nonintuitive possibilities arising in the considerations of optimum dosing schedules. The time-dependent dynamics of the drugs coupled to the population dynamics parameters of the tumor cell clones lead to a complex interplay between time scales, whose outcomes are nontrivial to predict a priori. We showed that a particular combination therapy regimen, which is optimal under one set of pharmacokinetic conditions, may be the worst strategy of drug administration when the drug pharmacokinetics are slightly different. Drug interactions add yet another layer of complexity, and our analysis highlights the need for including these effects into a systematic, rational framework for correctly identifying and analyzing optimal dosing strategies of multidrug treatments. For instance, we showed that when different resistant clones exist at the time of treatment initiation, the best dosing strategy may or may not be simultaneous dosing of all the drugs, depending on whether the drugs behave in a synergistic or additive manner. By combining such a computational framework with careful experimental parameterization of growth and death rates of tumor cell clones (17), and optimization techniques like simulated annealing to search through a high-dimensional space of treatment strategies (26), investigators will be able to design robust combination dosing strategies for heterogeneous patient cohorts.

To demonstrate how our quantitative computational framework can be applied in clinical settings, we analyzed the evolutionary dynamics of non–small cell lung cancer (NSCLC) in response to the two targeted drugs AZD9291 and selumetinib. By combining results from preclinical and previous clinical studies, we were able to suggest a simple way of significantly improving the efficacies of the treatment strategies being used in an ongoing phase Ib clinical trial of EGFR-mutant lung cancer patients. Instead of the currently used dosing strategy of once a day AZD9291 and twice a day selumetinib, we predicted that administering AZD9291 twice a day and selumetinib once a day would be significantly better in terms of cancer cell eradication as well as lowering toxicities. In addition, we showed that reducing selumetinib dosing to four days a week (instead of every day) would make no difference with regard to the rate of cancer cell eradication while reducing toxicities. Our analysis of the TATTON dosing strategies also highlighted a crucial point: that schedules administering lower drug concentrations can sometimes be equally efficacious as MTD schedules, thereby possibly reducing side effects and associated toxicities. Our work highlights the role that evolutionary modeling can play in the development of rational and more efficacious drug dosing strategies while taking complexities such as drug interactions and patient heterogeneities into account.

Our proposed strategies and general principles require validation in the clinic in the future. However, we believe that the strength of this work lies in the fact that we highlight, for the first time, the importance of accounting for pharmacokinetics and drug interactions in building predictive models of combination therapies. This approach is therefore an important step towards validating evolutionary mathematical models in the clinic. Finally, our modeling approach does not take into account effects of the immune system, the extracellular matrix and many other cell-extrinsic factors that are known to affect the dynamics of cancer cells. These are active areas of research, and as experimental model systems allow for better characterization of such effects, evolutionary models will become more accurate and translatable to the clinic.

No potential conflicts of interest were disclosed.

Conception and design: S. Chakrabarti, F. Michor

Development of methodology: S. Chakrabarti, F. Michor

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): S. Chakrabarti

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): S. Chakrabarti, F. Michor

Writing, review, and/or revision of the manuscript: S. Chakrabarti, F. Michor

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): S. Chakrabarti

Study supervision: F. Michor

We would like to thank Geoffrey Oxnard for interesting discussions on dosing strategies and an unknown reviewer for a number of valuable suggestions.

This work was supported by the Dana-Farber Cancer Institute Physical Sciences-Oncology Center (U54CA193461 and U54CA143798 to F. Michor).

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.

1.
Nowell
PC
. 
The clonal evolution of tumor cell populations
.
Science
1976
;
194
:
23
8
.
2.
Greaves
M
,
Maley
CC
. 
Clonal evolution in cancer
.
Nature
2012
;
481
:
306
13
.
3.
Foo
J
,
Michor
F
. 
Evolution of acquired resistance to anti-cancer therapy
.
J Theor Biol
2014
;
355
:
10
20
.
4.
Gottesman
MM
. 
Mechanisms of cancer drug resistance
.
Annu Rev Med
2002
;
53
:
615
27
.
5.
Garraway
LA
,
Jänne
PA
. 
Circumventing cancer drug resistance in the era of personalized medicine
.
Cancer Discov
2012
;
2
:
214
26
.
6.
Holohan
C
,
Van Schaeybroeck
S
,
Longley
DB
,
Johnston
PG
. 
Cancer drug resistance: an evolving paradigm
.
Nat Rev Cancer
2013
;
13
:
714
26
.
7.
Michor
F
,
Nowak
M
,
Iwasa
Y
. 
Evolution of resistance to cancer therapy
.
Curr Pharm Des
2006
;
12
:
261
71
.
8.
Komarova
NL
,
Wodarz
D
. 
Drug resistance in cancer: principles of emergence and prevention
.
Proc Natl Acad Sci U S A
2005
;
102
:
9714
9
.
9.
Kobayashi
S
,
Boggon
TJ
,
Dayaram
T
,
Jänne
PA
,
Kocher
O
,
Meyerson
M
, et al
EGFR mutation and resistance of non-small-cell lung cancer to gefitinib
.
N Engl J Med
2005
;
352
:
786
92
.
10.
Pao
W
,
Miller
VA
,
Politi
KA
,
Riely
GJ
,
Somwar
R
,
Zakowski
MF
, et al
Acquired resistance of lung adenocarcinomas to gefitinib or erlotinib is associated with a second mutation in the EGFR kinase domain
.
PLoS Med
2005
;
2
:
e73
.
11.
Engelman
JA
,
Zejnullahu
K
,
Mitsudomi
T
,
Song
Y
,
Hyland
C
,
Park
JO
, et al
MET amplification leads to gefitinib resistance in lung cancer by activating ERBB3 signaling
.
Science
2007
;
316
:
1039
43
.
12.
Zhou
W
,
Ercan
D
,
Chen
L
,
Yun
C-H
,
Li
D
,
Capelletti
M
, et al
Novel mutant-selective EGFR kinase inhibitors against EGFR T790M
.
Nature
2009
;
462
:
1070
4
.
13.
Walter
AO
,
Sjin
RTT
,
Haringsma
HJ
,
Ohashi
K
,
Sun
J
,
Lee
K
, et al
Discovery of a mutant-selective covalent inhibitor of EGFR that overcomes T790M-mediated resistance in NSCLC
.
Cancer Discov
2013
;
3
:
1404
15
.
14.
Hidalgo
M
,
Siu
LL
,
Nemunaitis
J
,
Rizzo
J
,
Hammond
LA
,
Takimoto
C
, et al
Phase I and pharmacologic study of OSI-774, an epidermal growth factor receptor tyrosine kinase inhibitor, in patients with advanced solid malignancies
.
J Clin Oncol
2001
;
19
:
3267
79
.
15.
Lee
SM
,
Khan
I
,
Upadhyay
S
,
Lewanski
C
,
Falk
S
,
Skailes
G
, et al
First-line erlotinib in patients with advanced non-small-cell lung cancer unsuitable for chemotherapy (TOPICAL): a double-blind, placebo-controlled, phase 3 trial
.
Lancet Oncol
2012
;
13
:
1161
70
.
16.
Foo
J
,
Michor
F
. 
Evolution of resistance to targeted anti-cancer therapies during continuous and pulsed administration strategies
.
PLoS Comput Biol
2009
;
5
:
e1000557
.
17.
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
.
Sci Transl Med
2011
;
3
:
90ra59
.
18.
Van Kampen
NG
.
Stochastic processes in physics and chemistry
.
Philadelphia
:
Elsevier
; 
2007
.
19.
Foo
J
,
Michor
F
. 
Evolution of resistance to anti-cancer therapy during general dosing schedules
.
J Theor Biol
2010
;
263
:
179
88
.
20.
Baguley
BC
. 
Multiple drug resistance mechanisms in cancer
.
Mol Biotechnol
2010
;
46
:
308
16
.
21.
Gillet
J-P
,
Gottesman
MM
. 
Mechanisms of multidrug resistance in cancer
.
Methods Mol Biol
2010
;
596
:
47
76
.
22.
Groenendijk
FH
,
Bernards
R
. 
Drug resistance to targeted therapies: Déjà vu all over again
.
Mol Oncol
2014
;
8
:
1067
83
.
23.
Szakács
G
,
Paterson
JK
,
Ludwig
JA
,
Booth-Genthe
C
,
Gottesman
MM
. 
Targeting multidrug resistance in cancer
.
Nat Rev Drug Discov
2006
;
5
:
219
34
.
24.
Coldman
AJ
,
Goldie
JH
. 
A model for the resistance of tumor cells to cancer chemotherapeutic agents
.
Math Biosci
1983
;
65
:
291
307
.
25.
Gatenby
RA
,
Silva
AS
,
Gillies
RJ
,
Frieden
BR
. 
Adaptive Therapy
.
Cancer Res
2009
;
69
:
4894
903
.
26.
Leder
K
,
Pitter
K
,
Laplant
Q
,
Hambardzumyan
D
,
Ross
BD
,
Chan
TA
, et al
Mathematical modeling of PDGF-driven glioblastoma reveals optimized radiation dosing schedules
.
Cell
2014
;
156
:
603
16
.
27.
Bozic
I
,
Reiter
JG
,
Allen
B
,
Antal
T
,
Chatterjee
K
,
Shah
P
, et al
Evolutionary dynamics of cancer in response to targeted combination therapy
.
eLife
2013
;
2
:
e00747
.
28.
Lee
MJ
,
Ye
AS
,
Gardino
AK
,
Heijink
AM
,
Sorger
PK
,
MacBeath
G
, et al
Sequential application of anticancer drugs enhances cell death by rewiring apoptotic signaling networks
.
Cell
2012
;
149
:
780
94
.
29.
Day
RS
. 
Treatment sequencing, asymmetry, and uncertainty: protocol strategies for combination chemotherapy
.
Cancer Res
1986
;
46
:
3876
85
.
30.
Katouli
AA
,
Komarova
NL
. 
The worst drug rule revisited: mathematical modeling of cyclic cancer treatments
.
Bull Math Biol
2011
;
73
:
549
84
.
31.
Mumenthaler
SM
,
Foo
J
,
Leder
K
,
Choi
NC
,
Agus
DB
,
Pao
W
, et al
Evolutionary modeling of combination treatment strategies to overcome resistance to tyrosine kinase inhibitors in non-small cell lung cancer
.
Mol Pharm
2011
;
8
:
2069
79
.
32.
Chen
J-H
,
Kuo
Y-H
,
Luh
HP
. 
Optimal policies of non-cross-resistant chemotherapy on Goldie and Coldman's cancer model
.
Math Biosci
2013
;
245
:
282
98
.
33.
Basanta
D
,
Gatenby
RA
,
Anderson
ARA
. 
Exploiting evolution to treat drug resistance: combination therapy and the double bind
.
Mol Pharm
2012
;
9
:
914
21
.
34.
Goldie
JH
,
Coldman
AJ
,
Gudauskas
GA
. 
Rationale for the use of alternating non-cross-resistant chemotherapy
.
Cancer Treat Rep
1982
;
66
:
439
49
.
35.
Lavi
O
,
Greene
JM
,
Levy
D
,
Gottesman
MM
. 
The role of cell density and intratumoral heterogeneity in multidrug resistance
.
Cancer Res
2013
;
73
:
7168
75
.
36.
Lorz
A
,
Lorenzi
T
,
Hochberg
ME
,
Clairambault
J
,
Perthame
B
. 
Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies
.
arXiv
:
1207.0923 2012
.
37.
Atkinson
AJ
 Jr
,
Kushner
W
. 
Clinical pharmacokinetics
.
Annu Rev Pharmacol Toxicol
1979
;
19
:
105
27
.
38.
Foo
J
,
Chmielecki
J
,
Pao
W
,
Michor
F
. 
Effects of pharmacokinetic processes and varied dosing schedules on the dynamics of acquired resistance to erlotinib in EGFR-mutant lung cancer
.
J Thorac Oncol
2012
;
7
:
1583
93
.
39.
Rand
KH
,
Brown
P
. 
Concentration-dependent synergy and antagonism between cefoperazone and imipenem against methicillin-resistant Staphylococcus aureus
.
Antimicrob Agents Chemother
1995
;
39
:
1173
7
.
40.
Prichard
MN
,
Shipman
C
. 
A three-dimensional model to analyze drug-drug interactions
.
Antiviral Res
1990
;
14
:
181
205
.
41.
Yeh
PJ
,
Hegreness
MJ
,
Aiden
AP
,
Kishony
R
. 
Drug interactions and the evolution of antibiotic resistance
.
Nat Rev Microbiol
2009
;
7
:
460
6
.
42.
Fitzgerald
JB
,
Schoeberl
B
,
Nielsen
UB
,
Sorger
PK
. 
Systems biology and combination therapy in the quest for clinical efficacy
.
Nat Chem Biol
2006
;
2
:
458
66
.
43.
Torella
JP
,
Chait
R
,
Kishony
R
. 
Optimal drug synergy in antimicrobial treatments
.
PLoS Comput Biol
2010
;
6
:
e1000796
.
44.
Clinical Trials.gov
AZD9291 in combination with ascending doses of novel therapeutics
.
Available from
: https://clinicaltrials.gov/ct2/show/NCT02143466.
45.
Cross
DAE
,
Ashton
SE
,
Ghiorghiu
S
,
Eberlein
C
,
Nebhan
CA
,
Spitzler
PJ
, et al
AZD9291, an irreversible EGFR TKI, overcomes T790M-mediated resistance to EGFR inhibitors in lung cancer
.
Cancer Discov
2014
;
4
:
1046
61
.
46.
Eberlein
CA
,
Stetson
D
,
Markovets
AA
,
Al-Kadhimi
KJ
,
Lai
Z
,
Fisher
PR
, et al
Acquired resistance to the mutant-selective EGFR inhibitor AZD9291 is associated with increased dependence on RAS signaling in preclinical models
.
Cancer Res
2015
;
75
:
2489
500
.
47.
Jänne
PA
,
Yang
JC-H
,
Kim
D-W
,
Planchard
D
,
Ohe
Y
,
Ramalingam
SS
, et al
AZD9291 in EGFR inhibitor-resistant non-small-cell lung cancer
.
N Engl J Med
2015
;
372
:
1689
99
.