## Abstract

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*.

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.

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

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).

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

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$|.

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

where |$C( t )$| is the concentration of the drug at time *t, C _{max}* 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

*C*and

_{max}*k*as pharmacokinetic or pharmacokinetic parameters.

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:

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

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$|.

## Introduction

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.

## Materials and Methods

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 *b _{0}*(

*t*) and

*d*, 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

_{0}*u*per cell division, thus generating new clones harboring specific resistance mechanisms, which are again characterized by their birth, death and mutation rates

_{i}*b*(

_{i}*t*),

*d*, and

_{i}*u*(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

_{i}*N*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).

_{i}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).

## Results

### 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.

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.

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.

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).

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 “*C*1” and “*C*2,” 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 *C*1, *C*2, 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, *C*1*, *C*2*, 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 (*C*1*, *C*2*, 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 (*C*1, *C*2, and *I*) and proposed (*C*1*, *C*2*, 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.

## Discussion

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.

## Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

## Authors' Contributions

**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

## Acknowledgments

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

## Grant Support

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.