Europe PMC
Nothing Special   »   [go: up one dir, main page]

Europe PMC requires Javascript to function effectively.

Either your web browser doesn't support Javascript or it is currently turned off. In the latter case, please turn on Javascript support in your web browser and reload this page.

This website requires cookies, and the limited processing of your personal data in order to function. By using the site you are agreeing to this as outlined in our privacy notice and cookie policy.

Novel computational method for predicting polytherapy switching strategies to overcome tumor heterogeneity and evolution.

Author information

Affiliations

  • 1. Department of Computing and Mathematical Sciences, California Institute of Technology, Pasadena, CA, USA.
      Authors
    • Jonsson VD 1
    • Matni N 1
    (2 authors)
  • 2. Department of Medicine, University of California, San Francisco, CA, USA.
      Authors
    • Blakely CM 2, 3
    • Lin L 2, 3
    • Asthana S 2, 3
    • Olivas V 2, 3
    • Pazarentzos E 2, 3
    • Gubens MA 2, 3
    • Bivona TG 2, 3
    (7 authors)
  • 3. Helen Diller Family Comprehensive Cancer Center, University of California San Francisco, San Francisco, CA, USA.
      Authors
    • Blakely CM 2, 3
    • Lin L 2, 3
    • Asthana S 2, 3
    • Olivas V 2, 3
    • Pazarentzos E 2, 3
    • Gubens MA 2, 3
    • Bastian BC 3
    • Bivona TG 2, 3
    (8 authors)
  • 4. Human Oncology and Pathogenesis Program, Memorial Sloan Kettering Cancer Center, USA.
      Authors
    • Taylor BS 4
    (1 author)
  • 5. Department of Control and Dynamical Systems, California Institute of Technology, Pasadena, CA, USA.
      Authors
    • Doyle JC 5
    (1 author)

ORCIDs linked to this article

Scientific Reports, 13 Mar 2017, 7:44206
https://doi.org/10.1038/srep44206 PMID: 28287179 PMCID: PMC5347024

This article is in the Europe PMC Open access subset. Refer to the copyright information in the article for licensing details.
Free full text in Europe PMC

This article is based on a previously available preprint.

Abstract 

The success of targeted cancer therapy is limited by drug resistance that can result from tumor genetic heterogeneity. The current approach to address resistance typically involves initiating a new treatment after clinical/radiographic disease progression, ultimately resulting in futility in most patients. Towards a potential alternative solution, we developed a novel computational framework that uses human cancer profiling data to systematically identify dynamic, pre-emptive, and sometimes non-intuitive treatment strategies that can better control tumors in real-time. By studying lung adenocarcinoma clinical specimens and preclinical models, our computational analyses revealed that the best anti-cancer strategies addressed existing resistant subpopulations as they emerged dynamically during treatment. In some cases, the best computed treatment strategy used unconventional therapy switching while the bulk tumor was responding, a prediction we confirmed in vitro. The new framework presented here could guide the principled implementation of dynamic molecular monitoring and treatment strategies to improve cancer control.

Free full text 

Logo of scirepAboutEditorial BoardFor AuthorsScientific Reports
Sci Rep. 2017; 7: 44206.
Published online 2017 Mar 13. https://doi.org/10.1038/srep44206
PMCID: PMC5347024
PMID: 28287179

Novel computational method for predicting polytherapy switching strategies to overcome tumor heterogeneity and evolution

Vanessa D. Jonsson,1,* Collin M. Blakely,2,3,* Luping Lin,2,3 Saurabh Asthana,2,3 Nikolai Matni,1 Victor Olivas,2,3 Evangelos Pazarentzos,2,3 Matthew A. Gubens,2,3 Boris C. Bastian,3,4 Barry S. Taylor,5,6,7 John C. Doyle,8,9,10 and Trever G. Bivonaa,2,3,11
Author information Article notes Copyright and License information Disclaimer

Associated Data

Supplementary Materials

Abstract

The success of targeted cancer therapy is limited by drug resistance that can result from tumor genetic heterogeneity. The current approach to address resistance typically involves initiating a new treatment after clinical/radiographic disease progression, ultimately resulting in futility in most patients. Towards a potential alternative solution, we developed a novel computational framework that uses human cancer profiling data to systematically identify dynamic, pre-emptive, and sometimes non-intuitive treatment strategies that can better control tumors in real-time. By studying lung adenocarcinoma clinical specimens and preclinical models, our computational analyses revealed that the best anti-cancer strategies addressed existing resistant subpopulations as they emerged dynamically during treatment. In some cases, the best computed treatment strategy used unconventional therapy switching while the bulk tumor was responding, a prediction we confirmed in vitro. The new framework presented here could guide the principled implementation of dynamic molecular monitoring and treatment strategies to improve cancer control.

Targeted cancer therapies are effective for the treatment of certain oncogene-driven solid tumors, including non-small cell lung cancers (NSCLCs) with activating genetic alterations in EGFR (epidermal growth factor receptor), ALK (anaplastic lymphoma kinase), BRAF, and ROS1 kinases1,2,3. However, inevitably resistance to current targeted therapies emerges, typically within months of initiating treatment and remains an obstacle to long-term patient survival1,2,3,4. The presence and evolution of tumor genetic heterogeneity potentially underlies resistance and also limits the response to successive therapeutic regimens that are used clinically in an attempt to overcome resistance in the tumor after it has emerged4,5,6,7. Indeed, while a targeted therapy may be effective in suppressing one genomic subclone within the tumor, other clones may be less sensitive to the effects of the drug. Thus, through selective pressures, resistant populations can emerge and promote tumor progression. Moreover, the current paradigm of solid tumor treatment is largely based on designing fixed (static) treatment regimens that are deployed sequentially as either initial therapy or after the clear emergence of drug-resistant disease, detected by clinical and radiographic measures of tumor progression. In contrast, designing dynamic treatment strategies that switch between targeted agents (or combinations thereof) in real time in order to suppress the outgrowth of rare or emergent drug-resistant subclones may be a more effective strategy to continually suppress tumor growth and extend the duration of clinical response. Thus, there is a need to identify principled approaches for the predictive design of effective combination (poly)therapy strategies to pre-empt the growth of multiple tumor subclones actively during treatment.

Mathematical modeling, analysis and computational simulations of tumor growth, heterogeneity and inhibition by various therapeutic modalities has long been employed as a method to provide insight into evolutionary outcomes and effective treatment strategies. Such modeling may include the use of stochastic8,9,10 or deterministic differential equation implementations11,12 to propose static or sequential treatment strategies that delay resistance in various cancer models. Recent studies by Zhao et al.13,14 incorporate the use of mathematical optimization, a fundamental subject in engineering design to predict static combination therapies that effectively address heterogeneity in a lymphoma model. Complementary engineering techniques from optimal control theory provide an additional theoretical framework to design dynamic drug scheduling regimens in the context of dynamical systems models of cancer heterogeneity and evolution. The application of optimal control theory to treatment design has a history dating back to the 1970s15,16 with more recent examples including that of scheduling angiogenic and chemotherapeutic agents17 or immuno- and chemotherapy combinations18. While mathematical modeling and engineering methods have been used extensively to inform treatment strategy design, a significant drawback to prior work in the field is that the underlying computational framework(s) have not conjointly accomplished the following important aims: (1) allowing for the systematic principled design of dynamic treatment strategies using experimentally identified models of tumor dynamic behaviors; and (2) developing quantitative methods that allow for the exploration of the robustness of predicted treatment strategies with respect to multiple common challenges in real-world patients, such as tumor heterogeneity and fluctuations in drug concentrations.

Here, we present a novel approach that combines a mathematical model of the evolution of tumor cell populations with parameters identified from our experimental data and an engineering framework for the systematic design of polytherapy scheduling directed at the following unresolved issues in the field: (1) how tumor genetic composition and drug dose constraints affect the long term efficacy of combination strategies, (2) how optimal scheduling of combination small molecule inhibitors can help to overcome heterogeneity, genomic evolution and drug dose fluctuations, and (3) how serial tumor biopsy or blood-based tumor profiling scheduling in patients can be timed appropriately. To tackle these questions, we developed an integrated experimental and computational framework that solves for candidate combination treatment strategies and their scheduling given an initial polyclonal tumor and allows the exploration of treatment design trade offs such as dosage constraints and robustness to small fluctuations in drug concentrations. This methodology is rooted in optimal control theory and incorporates an experimentally derived mathematical model of evolutionary dynamics of cancer growth, mutation and small molecule inhibitor pharmacodynamics to solve for optimal drug scheduling strategies that address tumor heterogeneity and constrain drug-resistant tumor evolution. Our key new insights include (1) heterogeneous tumor cell populations are better controlled with switching strategies; indeed, static two-drug strategies are unable to effectively control all tumor subpopulations in our study; (2) constant combination drug strategies are less robust to perturbations in drug concentrations for heterogeneous tumor cell populations, and hence more likely to lead to tumor progression; (3) countering the outgrowth of subclonal tumor populations by switching polytherapies even during a bulk tumor response can offer better tumor cell population control, offering a non-intuitive clinical strategy that pro-actively addresses molecular progression before evidence of clinical or radiographic progression appears.

Results

The presence and evolution of intratumoral genetic heterogeneity in a patient with EGFR-mutant lung adenocarcinoma

To explore the utility of our approach, we focused on EGFR-mutant lung adenocarcinoma. Many mechanisms of resistance to EGFR-targeted therapies in lung adenocarcinoma are well characterized19. Furthermore, tumor heterogeneity and multiple resistance mechanisms arising in a single patient can occur2,19. Thus, overcoming polygenic resistance is of paramount importance in this disease and will likely require a non-standard approach. To illustrate this point, we investigated the molecular basis of targeted therapy resistance in a 41-year old male never-smoker with advanced EGFR-mutant (L858R) lung adenocarcinoma. This patient responded to first-line treatment with erlotinib but progressed on this therapy within only four months after initial treatment, instead of the typical 9–12 month progression free survival observed in EGFR-mutant lung adenocarcinoma patients. We reasoned that genomic analysis of this patient’s outlier clinical phenotype could reveal the molecular pathogenesis of suboptimal erlotinib response. Using a custom-capture assay20,21, we deeply sequenced the coding exons and selected introns of 389 cancer-relevant genes in both the pre-treatment and the erlotinib-resistant tumor specimen and matched normal blood to identify somatic alterations that could mediate resistance (Materials and Methods). Exome sequencing of the pre-treatment specimen confirmed the presence of the EGFRL858R mutant allele that was identified through prior clinical PCR-based sequencing of this EGFRL858R specimen (data not shown), and additionally revealed mutant allele-specific focal amplification of the EGFR coding locus that resulted in a high allelic frequency (95% variant frequency) (Fig. 1B,C). We discovered a rare concurrent subclone in the treatment-naïve tumor with a BRAFV600E mutation (6% variant frequency; Fig. 1B). This observation is consistent with a recent report of a BRAFV600E mutation in an erlotinib-resistant lung adenocarcinoma specimen22 and recent data indicating that EGFR-mutant lung adenocarcinoma cells can often develop EGFR TKI resistance through RAF-MEK-ERK pathway activation23. The frequency of the subclonal BRAFV600E mutation increased approximately 10-fold upon acquired erlotinib resistance, from 6% to 60% in the primary and recurrent tumor, respectively (Fig. 1C). This increase in the BRAFV600E allelic fraction was likely due to the expansion of the BRAFV600E subclone, given that we found no evidence that this increased frequency occurred as a result of focal BRAF amplification in the resistant tumor (Fig. 1C and Fig. S1). Beyond the outgrowth of mutant BRAF, we identified two additional genetic alterations in the resistant tumor that could contribute to EGFR TKI resistance: focal amplification of 7q31.2 encoding MET in the resistant tumor cells (Fig. 1B and Fig. S1), a low frequency EGFRT790M mutation (14% variant frequency) (Fig. 1B,C). All candidate somatic mutations and focal copy number amplifications conferring resistance to erlotinib therapy (in EGFR, BRAF, and MET) were confirmed by independent, validated DNA sequencing and FISH assays (data not shown and Fig. S1). By using mutant specific antibodies recognizing EGFRL858R and BRAFV600E, we found co-occurrence of these mutant proteins in individual tumor cells (Fig. S1). FISH analysis also suggested that concurrent genetic alteration of EGFR and MET occurred in individual tumor cells (Fig. S1). Thus, erlotinib therapy acted as a selective pressure for the evolution of multiple concurrent clonal and subclonal genetic alterations that could cooperate to drive rapid drug-resistant disease progression in EGFR-mutant lung adenocarcinoma.

An external file that holds a picture, illustration, etc. Object name is srep44206-f1.jpg
Concurrent genetic alterations drive rapid resistance to EGFR TKI treatment in EGFR-mutant lung adenocarcinoma.

(A) Computed tomography indicates the clinical course and timeline of disease in the patient with rapid progression on EGFR TKI therapy and shows the EGFR-mutant lung adenocarcinoma (red arrows) analyzed both prior to erlotinib treatment and upon resistance at 4 months. (B) Key somatic mutations identified by exon-capture and deep sequencing of the pre- and post-treatment tumor in (A) demonstrating concurrent alterations in EGFR and BRAF and the frequency of each mutation in pre- and post- treatment tumor samples. P-values indicated as determined by a two-tailed Fischer’s exact test. (C) DNA copy number alterations inferred from exon-capture and sequencing data indicate the focal amplification of the EGFRL858R-mutant allele was lost upon acquired resistance while the patient’s resistant tumor gained a focal amplification of MET, with no change in BRAF (relative positions indicated, chromosome 7).

Analysis of clonal concurrence and resistance

While BRAFV600E, MET activation, and EGFRT790M can individually promote EGFR TKI resistance22,24,25, the therapeutic impact of the concurrence of these alterations we uncovered has not been characterized. Therefore, we studied the effects of BRAFV600E, MET activation, and EGFRT790M, alone or in combination, on growth and therapeutic response in human EGFR-mutant lung adenocarcinoma cellular models. First, we found that expression of V600E but not wild-type (WT) BRAF promoted resistance to erlotinib in 11–18 cells that endogenously express EGFRL858R (Fig. S2). This erlotinib resistance in BRAFV600E-expressing EGFR-mutant 11–18 cells was overcome by concurrent treatment with erlotinib and selective inhibitors of either BRAF or MEK (vemurafenib26 and trametinib27 respectively (Figs S2 and S3). We next used the 11–18 system to test the effects of MET activation by hepatocyte growth factor (HGF), which phenocopies the effects of MET amplification in EGFR TKI resistance25,28 on therapeutic sensitivity. We found that MET activation not only promoted erlotinib resistance in parental 11–18 cells but also enhanced the effects of BRAFV600E expression on erlotinib resistance in these cells (Fig. S3). This resistance induced by MET activation in 11–18 parental and BRAFV600E-expressing cells was accompanied by increased phosphorylation of MEK, ERK, and AKT (Fig. S4). Treatment with the MEK inhibitor trametinib, but not the BRAF inhibitor vemurafenib or the MET inhibitor crizotinib, overcame erlotinib resistance and inhibited phospho-ERK in MET-activated BRAFV600E-expressing 11–18 cells (Fig. S4), providing a rationale for polytherapy against EGFR and MEK in EGFR-mutant tumors with activating co-alterations in MET and BRAF.

Given that we found a rare EGFRT790M subclone in the polyclonal resistant tumor, we next explored whether BRAFV600E expression could promote resistance to EGFR TKI treatment in H1975 human lung adenocarcinoma cells that endogenously express EGFRT790M and EGFRL858R. We observed that BRAFV600E modestly decreased sensitivity to afatinib, an approved irreversible EGFR kinase inhibitor effective against EGFRT790M29, and that this effect of BRAFV600E on afatinib sensitivity was blunted by vemurafenib (Fig. S5). Together, our data indicate that erlotinib therapy induced the evolution of multiple concurrent events that re-shaped the polyclonal tumor genetic landscape during the onset of resistance; resistance could be overcome by polytherapy against both EGFR and MAPK signaling in preclinical models.

Polytherapy Provides Temporary Response in Heterogeneous or MET Activated Tumors

While we conducted a finite set of experiments to test various rational drug combinations that could address the heterogeneous basis of resistance in this patient’s disease, this approach is not easily scaled; further, it is not readily feasible to explore all possible drug combinations and drug doses over a continuous range or anticipate the effects of the myriad of possible tumor subcompositions on tumor control under treatment using cell-based assays alone. Therefore, we sought to provide a more general and scalable framework for understanding the impact of each genetically-informed targeted therapy strategy on the temporal evolution of the multiple concurrent EGFR-mutant tumor cell subclones present in this patient, as a potentially more generalizable platform. We developed an ordinary differential equation (ODE) model of tumor growth, mutation and selection by small molecule inhibitors with parameters identified from experimental data (Fig. S11A,B and Equation S1) and interrogated it to uncover the limitations of the targeted treatments in the context of tumor heterogeneity and evolution. We first confirmed that our model was able to capture the essential tumor population dynamics by showing a qualitative equivalence between the patient’s clinical course and our model simulation of similar tumor subpopulations consisting of 94% EGFRL858R, 6% BRAF V600E and assuming the existence of a very low initial frequency of 0.01% MET amplification of EGFRL858R, BRAFV600E and EGFRT790M in the presence of 1 μM erlotinib (Fig. 2A,B).

An external file that holds a picture, illustration, etc. Object name is srep44206-f2.jpg
Mathematical simulation qualitatively captures the patient’s evolution on erlotinib.

(A) A simulation of the mathematical model of lung adenocarcinoma evolution (SI, Equation (S1)) in the presence of 1 μM erlotinib, given the patient-derived pretreatment initial tumor cell subpopulations (94% EGFRL858R, 6% BRAF V600E, 0.01% MET amplification of EGFRL858R, BRAFV600E and EGFRT790M). Parameters used in the simulation were derived from growth and viability assays of parental 11–18 EGFRL858R-positive lung adenocarcinoma cells or those cells engineered to express mutations listed above and treated with 0 or 50 ng/ml HGF, in the presence of varying concentrations of erlotinib and fit according to Equations S8, S9 and S11. (B) Tumor cell populations present at day 0, 6 and 17 of the simulation in (A), including the total HGF+ cell population at day 17 (gray). The model qualitatively captures a possible evolutionary trajectory and results in a similar final tumor cell composition as that of the patient, (B) day 17 vs. Fig. 1(B and C).

To systematically explore the utility of many different drug combination regimens to overcome polygenic resistance, we used our computational model to calculate the efficacy of clinically relevant doses of erlotinib and afatinib in combination with either crizotinib, trametinib or vemurafenib on the growth of parental 11–18 and H1975 cells EGFR mutant cell lines. We found that most polytherapies could address only certain subpopulations (Fig. 3A). For example, the afatnib/trametinib combination elicited a complete response for a representative heterogeneous MET-negative tumor cell population comprised of (89% EGFRL858R, 10% EGFRL858R BRAFV600E, 1% EGFRL858R, T790M) compared to rapid progression for its MET activated analog (Fig. 3C). Moreover, we computed the concentrations of erlotinib or afatinib in combination that could guarantee a progression-free response for both MET activated or MET neutral tumor cell populations (SI, Mathematical Methods) and found that in many cases, the concentrations were considerably higher than clinically feasible (due to either known pharmacokinetic limitations or dose limiting toxicities) (Fig. 3B).

An external file that holds a picture, illustration, etc. Object name is srep44206-f3.jpg
Modeling pharmacodyamic effects of concurrent BRAFV600E expression and MET activation in EGFR-mutant lung adenocarcinoma cells and their implication on progression.

(A) Drug efficacy as measured by the effect of 1.5 μM erlotinib or 0.5 μM afatinib in combination with either 0.5 μM MET inhibitor crizotinib, 0.5 μM MEK inhibitor trametinib or 5 μM BRAF inhibitor vemurafenib on cell growth (SI, Equation S1) of parental 11–18 EGFRL858R-positive lung adenocarcinoma cells or those cells engineered to express mutations listed above and treated with 0 or 50 ng/ml HGF. (B) Concentrations of EGFR TKIs afatinib and erlotinib in combination with either 0.5 μM crizotinib, 0.5 μM trametinib or 5 μM vemurafenib that guarantee progression free tumor reduction for any HGF− or HGF+ initial tumor subpopulations according to the model, measured by the minimum concentration of erlotinib or afatinib that results in exponential stability of the evolutionary dynamics model (SI, Section 3.2). (C) Simulations of the lung adenocarcinoma model for combinations of 0.5M afatinib + 0.5 μM trametinib and 1.5 μM erlotinib + 0.5 μM μcrizotinib for the HGF− and HGF+ tumors specified. (D) (Left) Simulations of the evolutionary dynamics of different HGF− lung adenocarcinoma initial tumor subpopulations with a constant treatment of 0.7 μM, 0.5, 0.3 or 0.1 μM afatinib in combination with 0.5 μM of trametinib (red) and of different HGF+ lung adenocarcinoma initial tumor subpopulations with a constant treatment of 8.32 μM, 3.2 μM, 1.5 μM or 0.75 μM erlotinib in combination with 0.5 μM crizotinib (blue). (Right) Maximum eigenvalue decompositions (SI, Section 3.2) classify which subpopulations can lead to progression at different concentrations of EGFR TKI for the afatinib + trametinib combination and the erlotinib + crizotinib combination.

To better understand the efficacy of the combination therapy over time, we sought to classify which initial tumor cell subpopulations could eventually lead to therapeutic failure when treated with different concentrations of EGFR TKIs in combination with crizotinib, trametinib or vemurafenib. We defined the evolutionary stability of a subpopulation as the worst-case evolutionary outcome, in each case where the particular subpopulation is present upon treatment initiation. More precisely, the evolutionary stability is the maximum eigenvalue of each evolutionary branch downstream of the subclone (SI, Section 3.2). This approach provides an assessment of which subclones present in the initial tumor cell population are likely to lead to overall progression (a positive evolutionary stability) versus those that lead to response (a negative evolutionary stability) when treated with a particular combination therapy. Our analysis confirms that progression-free response on combination therapies is sensitive to both EGFR TKI concentration and dependent on whether pre-existent subpopulations are effectively targeted at these concentrations (Fig. 3D and Figs S6–S8). Overall, this analysis revealed that combinations of two signal transduction inhibitors had limited effectiveness in durably controlling resistance over a longer time horizon.

Engineering Drug Scheduling to Control Tumor Evolution

We next explored how the rational design of combination drug scheduling strategies could address this issue. Experimental studies have recently proposed drug pulsing30 or drug switching10 as a strategy to delay the growth of certain cancers. To this end, we proposed a novel methodology rooted in engineering principles to design drug scheduling strategies that best control the growth and evolution of tumor cell populations. In particular, we apply concepts from optimal and receding horizon control theory to our experimentally integrated model of lung adenocarcinoma evolution to compute treatment strategies that minimize tumor cell populations over time. Our algorithm allows for the specification of treatment design constraints such as maximum dose, the time horizon over which the treatment strategy is applied and the switching horizon, that is the minimum time over which one particular treatment can be applied. This algorithm can be extended to include other drug related characteristics and treatment design constraints. In addition, the framework allows for the analysis of tradeoffs between these aspects of the design space as well as others, such as how robust the predicted treatment strategies are with respect to uncertainties in the model or perturbations in drug dosages.

For a predetermined time and minimum switching horizon, we define an optimal control problem (SI, Algorithm 1) and solve for the drug combination that best minimizes the existing tumor cell subpopulations for every receding switching horizon. Given that any one polytherapy is unlikely to be simultaneously effective against all subpopulations, the resulting optimal strategy, which maximizes the response of the tumor cells present at every time horizon (SI, Mathematical Methods), is potentially one that switches between drug combinations, at defined time points during the treatment course.

As proof-of-principle, we determined which drug scheduling regimens could maximally reduce different initial tumor cell populations by solving our control problem for different allowable switching horizons over a thirty day period (Fig. 4). The afatinib/trametinib combination was the optimal constant strategy for tumor cell populations harboring the EGFRL858R, T790M mutation, and although this strategy invoked progression free response in HGF− tumor cell populations, most EGFRL858R HGF+ tumor cell populations progressed on the therapy over thirty days (Fig. 4A vs 4C and Fig. S7AB). For the HGF− tumor population comprised of 89% EGFRL858R, 10% EGFRL858RBRAFV600E and 1% EGFRL858R, T790M, the optimal constant strategy provided overall response leaving a dominant EGFRL858RBRAFV600E tumor subpopulation present whereas the optimal ten day switching strategy provided an enhanced response over the constant strategy by alternately targeting EGFRL858R and EGFRL858R, T790M tumor cell subpopulations (Fig. 4B). In the case of the HGF treated tumor cell distribution consisting of 90% EGFRL858R and 10% EGFRL858R, T790M, a constant combination of afatinib/trametinib was effective against the EGFRL858R, T790M, HGF+ subpopulation despite overall progression due to the outgrowth of the EGFRL858R, HGF+ tumor cell population, whereas a five day switching regimen between afatinib/trametinib and erlotinib/crizotinib combinations alternately targeted the HGF+ EGFRL858R, T790M and the EGFRL858R subpopulations (Figs 3A and ​and4B)4B) leading to overall response.

An external file that holds a picture, illustration, etc. Object name is srep44206-f4.jpg

Optimal drug scheduling strategies solved by Algorithm 1 (SI, Section 2.2) for representative initial tumor cell distributions (A),(C), for a 30 day timeframe and 30, 15, 10, 5, 3 and 1 day minimum switching horizons, give one EGFR TKI, either 1.5 μM erlotinib (ERL) or 0.5 μM afatinib (AFA) in combination with either 5 μM vemurafenib (VEM), 0.5 μM trametinib (TRA) or 0.5 μM crizotinib (CRI) and corresponding simulations (B,D) of the lung adenocarcinoma evolutionary dynamics for a subset of optimal drug scheduling strategies.

More generally, the optimal constant strategies determined by our algorithm are combinations that best minimize existing tumor cell subpopulations at every switching horizon. In particular, a greater reduction in tumor cells can be achieved by switching between therapies that alternately target different subpopulations, even while there is overall response in the tumor (Fig. 4A). This finding suggests a non-intuitive approach to the clinical management of solid tumors that would represent a departure from the current standard clinical practice. Our model suggests an advantage to switching treatments pro-actively even during a bulk tumor response, while the current paradigm in the field is to switch from the initial treatment to a new drug(s) only after there is clear evidence of radiographic or clinical progression on the initial treatment.

To understand the potential benefits of switching strategies in tumors with different initial genetic heterogeneity, we computed the optimal switching strategies for a subset of tumor cell distributions and compared them to their corresponding computed optimal constant strategies. We found that the larger the number of subclones present in the initial tumor, the more beneficial even a small number of switches could be for overall tumor cell population control (Fig. 5A and Fig. S9A). For a highly heterogeneous tumor cell population comprised of HGF treated 89% EGFRL858R, 10% EGFRL858RBRAFV600E, 1% EGFRL858R, T790M mutations, the predicted fifteen day switching therapy (afatinib/trametib followed by erlotinib/crizotinib) provides an immediate benefit versus the predicted constant treatment strategy (afatinib/trametinib), yielding a 10-fold decrease in final tumor population. By contrast, for a more homogeneous tumor consisting of 90% EGFRL858R, 10% EGFRL858R, T790M, the optimal predicted 30, 15 and 10 day switching strategies are indistinguishable from the constant therapy strategy for population control. Our predictions indicate that a similar 10-fold reduction in final population (similar to that achieved in the heterogeneous tumor instance analyzed above) is achieved only with a more rapid, five day switching strategy for this more homogeneous tumor population (afatinib/trametinib, then alternating between erlotinib/trametinib and afatinib/vemurafenib). These findings emphasize our results that while polytherapy may a provide response in some subsets of tumor cell populations, it provides only a temporary or no response in heterogeneous or MET activated tumors; in these cases, even minimal therapy switching can provide an immediate and more substantial benefit for overall tumor population control.

An external file that holds a picture, illustration, etc. Object name is srep44206-f5.jpg
Exploring the robustness of treatment strategies through model simulation.

(A) Switching strategies are more beneficial to tumor cell populations with more initial heterogeneity. (Left) Fold change in final lung adenocarcinoma tumor cell populations at day 30 versus day 0 over the course of the optimal 30, 15, 10, 5, 3, and 1 day treatment strategies solved by algorithm 1 (SI, Section 2.2) and normalized by fold change in final tumor cell population for the constant 30 day treatment strategy for an initial tumor cell population comprised of (90% EGFRL858R, 10% H1975 EGFRL858R, T790M) and another comprised of (89% EGFRL858R, 10% BRAFV600E, 1% EGFRL858R, T790M) subclones. (Right) Sum of fold change for the final lung adenocarcinoma populations (SI, Equation S5) for select initial tumor cell distributions (SI, Table 1) and their corresponding optimal 30, 15, 10, 5, 3, and 1 day treatment strategies, categorized by the number of subclones in the initial tumor cell population. Smaller fold change sums indicate that more switching is beneficial to reduce final populations, whereas larger fold changes indicate that more switching does not necessarily help in reducing the final tumor populations. (B) EGFR TKI dose perturbations. (Left) Fold change in number of lung adenocarcinoma cells between day 30 and day 0, as a function of percent EGFR TKI dose reduction for the optimal 30, 15, 10, 5 and 1 day strategies solved by algorithm 1 (SI, Section 2.2) for tumor cell populations indicated above. The shaded areas indicate the regions of the perturbation space where the treatment strategy reduces the initial tumor cell population by more than 30% (response, light blue), increases the size of the original tumor population size by more than 20% (progression, red), or maintains the original tumor population size between the two (stability, white). (Right) Bar graphs indicate the maximum reduction in EGFR TKI dose supported by the optimal strategy such that there is still reduction in tumor size at day 30 with respect to day 0 for the V600E and the pretreatment MET tumor. (C) The average maximum percent EGFR TKI dose reduction supported before progression for lung adenocarcinoma tumors with different number of initial tumor cell subpopulations and for predicted optimal 30, 15, 10, 5, and 1 day switching strategies.

Robustness Analysis of Switching Strategies

Motivated by studies indicating that tissue to plasma ratios for certain drugs such as erlotinib can be low31, we sought to computationally explore how dose reductions of TKI combinations could affect the evolution of tumor cell populations. This is a particularly relevant clinical issue, as many drugs when used in combination often require a reduction in the recommended monotherapy dose due to toxicity of the dual drug therapy in patients. To examine this question, we simulated the optimal switching strategies corresponding to 30, 15, 10, 5 and 1 day switching horizons subject to EGFR TKI dose reductions for a set of initial tumor cell populations and studied the effects on the final and average tumor populations over the course of the treatment (SI, Mathematical Methods).

For a tumor with a smaller number of initial subclones, such as one comprised of 90% EGFRL858R and 10% EGFRL858R, T790M, all switching strategies induced a response for EGFR TKI dose reductions of up to 50% (Fig. 5B). In contrast, with the more complex HGF treated tumor cell population comprised of 89% EGFRL858R, 10% EGFRL858RBRAFV600E, 1% EGFRL858R, T790M, only combination strategies with switching horizons of 10 day or shorter induced a response (Fig. 5B). Notably, we observed that the shorter the switching horizon, the higher dose reduction that could be supported while still maintaining a progression free response (Fig. 5B and Fig. S9B). We observed this phenomenon more generally when we simulated different tumor cell initial distributions (Fig. 5C). Thus, we find that the greater number of subclones present in the initial tumor, the greater the benefit there is in increasing switching frequency in terms of the achieving robustness to perturbations in EGFR TKI drug concentration.

Switching Strategies Control or Delay Progression in vitro

Motivated by the results of our treatment strategy algorithm, we tested drug scheduling strategies on select tumor subpopulations in an in vitro model of EGFR mutant lung adenocarcinoma. Specifically, we synthesized the optimal treatment strategy for a heterogeneous HGF treated tumor cell population consisting of 89% EGFRL858R, 10% EGFRL858RBRAFV600E, 1% EGFRL858R, T790M, and imposed a constraint that at most one switch could occur, as a starting point to simulate what might be most clinically feasible. The resulting optimal treatment strategy predicted by our modeling, consisting of the erlotinib/crizotinib (days 0–5) followed by the afatinib/trametinib (days 5–30) combination, was shown to elicit the best response in vitro, validating our predictive model (Fig. 6B).

An external file that holds a picture, illustration, etc. Object name is srep44206-f6.jpg
Engineering optimal treatment strategies for concurrent, clonal genetic alterations in EGFR-mutant lung adenocarcinoma and predicting their therapeutic impact.

(A) Simulations of the optimal treatment strategy predicted by algorithm 1 (SI, Section 2.2) consisting of 1.5 μM erlotinib + 0.5 μM crizotinib for days (0–5) followed by 0.5 μM afatinib + 0.5 μM trametinib for days (5–30); the same strategy but with the switch occurring at day 10 and, constant strategies of 0.5 μM afatinib + 0.5 μM trametinib or 1.5 μM erlotinib + 0.5 μM crizotinib for 30 days, for an initial tumor cell population of 89% EGFRL858R, 10% EGFRL858RBRAFV600E, 1% EGFRL858R, T790M, HGF treated. (B) Evolution experiment shows that the predicted strategy for an initial tumor cell population of 89% EGFRL858R, 10% EGFRL858RBRAFV600E, 1% EGFRL858R, T790M, treated with 50 ng/ml HGF, is optimal. Overlaid numbers indicate the relative cell density of each well at day 30 compared to the erlotinib + crizotinib well (magenta). Computational simulations in (A) show that the predicted optimal strategy has the greatest reduction in tumor cells in vitro (B, red) compared to the same strategy with a 10 day switch (yellow). A simulation of the model predicts that a constant treatment of afatinib + trametinib produces little change in number of tumor cells (B, blue) and that a constant treatment of erlotinib + crizotinib predicts the exponential outgrowth of the initial EGFRL858R, T790M MET amplified subpopulation, experimentally validated in (B, magenta).

To show how a delay in the switching time might affect response to therapy, we tested equivalent initial tumor cell populations but changed the treatment strategy to start the afatinib/trametinib combination at day 10 instead of at day 5. This resulted in worse overall response than the 5 day switching regimen (Fig. 6B). The corresponding model simulation highlights that although the erlotinib/crizotinib combination effectively targeted the HGF treated EGFRL858R mutation during the first 10 days, it allowed the HGF treated EGFRL858R, T790M subclone to dominate for a longer period of time, thereby impeding overall response.

Discussion

One of the fundamental challenges in the principled design of combination therapies is the pre-existence and temporal expansion of intratumor genetic heterogeneity that can often lead to rapid resistance with first-line targeted therapies. To address this problem, we sought to develop a new modeling framework to systematically design principled tumor monitoring and therapeutic strategies. We applied a receding horizon optimal control approach to an evolutionary dynamics and drug response model of lung adenocarcinoma that was identified from experimental and clinical data. Based on the clinical and experimental data, our computational method generated optimal drug scheduling strategies for a comprehensive set of initial tumor cell subpopulation distributions.

Our initial insight was that constant drug combination strategies that guarantee progression free response for tumor cell populations with considerable heterogeneity and/or MET activation, required EGFR TKI concentrations that were considerably higher than are typically clinically feasible. At clinically relevant doses, these constant combination strategies were not effective against all tumor cell subpopulations and inevitably, those subpopulations with even slight evolutionary advantages could undergo clonal expansion and cause resistance. To overcome this issue, we used our algorithm to generate optimal drug scheduling strategies that could preempt the outgrowth of these subpopulations over fixed switching periods, and showed that these strategies outperformed constant combination strategies for most tumor cell subpopulation distributions. Notably, our computational analysis showed there was more benefit in applying switching strategies in the context of increasing pre-existing genetic heterogeneity and these switching strategies provided more robustness guarantees in the presence of perturbations in drug concentrations that can occur in patients. We demonstrated successful in vitro validation of our optimal control approach for selected tumor subpopulation distributions. In particular, for an in vitro analog of our clinical case, a non-intuitive combination therapy switching strategy offered better tumor control than constant treatment strategies.

We found that the most effective drug scheduling strategies were ones that addressed existing subpopulations as they emerged during the course of the treatment, even during a bulk tumor response. In contrast, current standard of care clinical practice is generally to delay switching to second-line therapy until after there is clear evidence of radiographic or clinical progression. Our approach suggests a paradigm shift that would require regular monitoring of an individual patient’s tumor mutational status, for instance by mutational analysis of plasma cell-free circulating tumor DNA, so-called “liquid biopsies”32,33,34,35,36. Our modeling strategy could potentially synthesize this genetic information to yield both the design and prioritization of specific drug regimens and the optimal time for clinical deployment, informed by the molecular findings in a particular patient. Such treatments may need to be applied (non-intuitively) during the initial tumor response, instead of later during therapy or after drug resistance is readily apparent by standard clinical measures in some cases. We envision that our approach could help contribute to the shift from a reactive to pro-active, dynamic management paradigm in solid tumor patients in the molecular era. Drug scheduling strategies synthesized by the algorithm for the initial tumor cell population could be adapted to account for genetic alterations that are detected by the analysis of serial liquid (or tumor) biopsies, leading to a dynamic learning model through iterative refinements; as such, the model could suggest more effective strategies with time. Additional validation and refinement of our model is required to address several important open questions, that include defining the role and impact, if any, of subclones of unclear significance, and the degree and implications of mutational and clonal concordance and discordance in liquid versus tumor biopsies. If a dose reduction is required due to drug toxicity, this framework allows for the computational exploration of alternative strategies that can best control tumor progression for a particular patient’s dosage constraints. In particular, we can quantify the effectiveness of these strategies and begin to understand the tradeoffs between dose reductions, tumor heterogeneity and frequency of switching. Additional considerations such as pharmacokinetics, the tumor microenvironment and metastatic processes37,38 could extend this model to add more clinical relevance. Our approach, if further validated in prospective clinical trials, could guide the optimal timing of serial clinical specimen sampling (plasma, tumor) and radiographic analysis to streamline clinical management. Overall, the combination of techniques stemming from mathematical optimization and control theory combined with more clinically applicable tumor dynamics models is a promising approach to aid the rational design, clinical testing, and clinical adoption of dynamic molecular monitoring and drug scheduling strategies to better control complex solid cancers such as lung cancer in real-time and improve clinical outcomes.

Materials and Methods

Computational Methods

The details of mathematical models and experimental methods may be found in SI Mathematical Methods. The mathematical model of lung adenocarcinoma growth mutation and selection by small molecule inhibitors was formulated as system of ordinary differential equations (ODEs). The treatment strategy algorithm was formulated as a receding horizon optimal control problem with the objective of minimizing lung adenocarcinoma populations at every horizon and implemented using python version 3.4.3, scipy version 1.11.0.

Experimental Methods

Patient sample preparation and sequence capture

Formalin fixed paraffin embedded (FFPE) NSCLC fine needle aspirate biopsy specimens and a normal blood sample were obtained from the patient under institutional informed consent both prior to erlotinib treatment and upon erlotinib resistance. Lung tumor biopsy specimens contained >75% tumor cells upon histopathological analysis by a board-certified pathologist. Barcoded sequence libraries were generated using genomic DNA from FFPE tumor material and matched normal blood using the NuGEN Ovation ultralow library systems and according to manufacturer’s instructions (NuGEN, San Carlos, CA). These libraries were among an equimolar pool of 16 barcoded libraries generated and subjected to solution-phase hybrid capture with biotinylated oligonucleotides targeting the coding exons of 389 cancer-associated genes using Nimblegen SeqV.D.J.Cap EZ (Roche NimbleGen, Inc, Madison, WI). Each hybrid capture pool was sequenced in a single lane of Illumina HiSeq2000 instrumentation producing 100 bp paired-end reads (UCSF Next Generation Sequencing Service). Sequencing data was demultiplexed to match all high-quality barcoded reads with the corresponding samples.

Sequencing Analysis

Paired-end sequence reads from normal blood, pre-treatment tumor, and erlotinib-resistant tumor samples were aligned against build hg19 of the reference genome with BWA39. Duplicate reads were marked, alignment and hybridization metrics calculated, multiple sequence realignment around candidate indels performed, and base quality scores recalibrated across all samples with the Picard suite (http://picard.sourceforge.net/) and the Genome Analysis Toolkit (GATK)40. Somatic point mutations were detected in the treatment-naïve and resistant tumors using MuTect41, while small insertions and deletions (indels) were identified with GATK. Given the depth of sequencing achieved and the presence of low-frequency oncogenic mutations in the normal sample likely due to circulating tumor DNA, mutations were excluded as germline if they exceeded a frequency of 10% in the normal sample. Non-synonymous mutations were annotated for their sequence context, effect, and frequency in lung adenocarcinoma and squamous cell tumors from The Cancer Genome Atlas (TCGA) project and Imielinski et al.42. All previously characterized oncogenic alleles in NSCLC or mutations previously linked to erlotinib resistance were also manual inspected in both treatment-naïve and resistant tumors. This analysis revealed a single sequencing read bearing the T790M mutation in the primary tumor (total coverage at this locus: 1300x). This was insufficient evidence from sequencing data to formally call the mutation, but we cannot exclude the possibility that EGFRT790M exists pre-treatment in a very rare clone (<0.08%) for which our target depth of coverage limited our sensitivity. DNA copy number alterations where inferred from the mean sequence coverage for each target region in each sample corrected for overall library size. Amplifications and deletions were determined from ratios of coverage levels between the pre- and post-treatment tumors and the matched normal blood sample. Due to the elevated signal to noise from targeted capture and sequencing of FFPE material from lower input amounts, overt genomic amplifications and deletions were required to affect multiple target regions (exons) of a given gene before being called as detected. The EGFRT790M and BRAFV600E variants were confirmed by a standard Clinical Laboratory Improvement Amendments (CLIA)-approved PCR-based shifted termination assay (data not shown). The changes in EGFR and MET copy number were validated using established fluorescence in situ hybridization clinical assays.

Cell Lines and Reagents

Human lung cancer cell lines were acquired as previously described43,44. Cells were grown in RPMI 1640 supplemented with 10% (high serum) or 0.5% (low serum) fetal bovine serum (FBS), penicillin G (100 U/ml) and streptomycin SO4 (100 U/ml). Erlotinib, afatinib, vemurafenib, crizotinib, and trametinib were purchased from Selleck Chemicals (Houston, TX). Drugs were resuspended in DMSO at a concentration of 10 mM and stored at −20 °C. Erlotinib and afatinib were used at working concentrations ranging from 0.010–1.5 μM. Vemurafenib was used at a working concentration of 5.0 μM, and trametinib and crizotinib were used at a 0.5 μM. HGF was purchased from Peprotech (Rocky Hill, NJ) and resuspended at 50 g/ml in sterile PBS + 0.1% BSA. Cells were treated with HGF at 50 ng/ml.

Generation of stable cell lines

293-GPG viral packaging cells were transfected with pBABE (empty vector), pBABE-mCherry-BRAF-WT and pBABE-mCherry-BRAFV600E constructs (kindly provided by Dr. Eric Collision, UCSF, San Francisco, CA) using Lipofectamine-2000 (Life Technologies, Pleasanton, CA) per manufacturer’s instructions. Virus containing media was harvested three days post transfection and used to infect 11–18 and H1975 lung cancer cell lines. Cells were incubated with virus containing media supplemented with 6 μg/ml of polybrene for 24 hours. Media was changed to standard cell growth media (RPMI-1640 + 10% fetal bovine serum and 100 U/ml penicillin G and 100 U/ml streptomycin SO4) and cells were expanded for 48 hours, at which point puromycin (2 μg/ml) was added to the media and cells were allowed to grow for an additional 4 days. Cells that survived puromycin selection (stable cell lines) were used in all subsequent experiments.

Cell Viability and Growth Assays

Assays were performed as previously described43,44. Briefly, cells were seeded overnight at a density of 5,000 cells per well in 96-well plates in RPMI containing 10% FBS and treated with indicated reagents for 72 hours. Viable cell numbers were determined using the CellTiterGLO assay according to manufacturer’s instructions (Promega). Each assay consisted of six replicate wells and was repeated at least twice in independent experiments. Cell viability is presented as the mean (±s.e.m.) erlotinib or afatinib inhibitory concentration 50 (IC50). Statistical significance between treatment groups was determined by the Bonferroni’s multiple comparisons ANOVA statistical test.

Immunoblot analysis

Cells were harvested 24 h after initiation of treatment with reagents. Cells were scraped and lysed in lysis buffer (50 mM Tris·HCl pH 8.0, 150 mM sodium chloride, 0.1% SDS, 0.5% sodium deoxycholate, 1% Triton X 100, 5 mM EDTA containing protease and phosphatase inhibitors (Roche Diagnostics. Indianapolis, IN). After quantitation by Pierce BCA assays (Thermo Scientific, Rockford, IL), 25 μg of each sample was separated by gel electrophoresis on 4–15% Criterion TGX precast gels (BioRad, Hercules, CA) and transferred to nitrocellulose membrane. For immunoblots, the following antibodies were used: anti-total EGFR (1:1000 dilution, Bethyl Laboratories, Inc., Montgomery TX), anti-pEGFR, anti-total Met, anti-pMet, anti-total Mek, anti-pMek, anti-total Akt, anti-pAkt, anti-total Erk, anti-pErk (1:1000, Cell Signaling Technology Inc., Danvers, MA), BRAFV600E Monoclonal Antibody (Clone VE1, 1:1000, Spring Bioscience, Pleasonton, CA), BRAF WT (1:1000, Santa Cruz Biotech, Santa Cruz, CA) and anti-actin (1:5000 dilution, Sigma-Aldrich, Saint Loius, MO), HRP-conjugated anti-rabbit Ig (used at a 1:3000 dilution, Cell Signaling), and HRP-conjugated anti-mouse IgG (used at a 1:3000 dilution, Cell Signaling). Specific proteins were detected by using either ECL Prime (Amersham Biosciences, Sunnyvale, CA) or the Odyssey Li-Cor (Lincoln, NE) with the infrared dye (IR Dye 800, IR Dye 680)-conjugated secondary antibodies (1:20,000, Li-Cor).

Immunofluorescence

5 micron FFPE sections were adhered to glass microscope slides (Fischer), and dewaxed 2× with Histoclear II (National Diagnostics). Tissues were rehydradated through an ethanol series (100%, 95%, 70% 2 min each), and then washed in with phosphate buffered saline (PBS). Antigen retrieval was performed with 0.1 M Citrate pH 6.0 with the 2100 Retriever (Electron Microscopy Sciences). After antigen retrieval, slides were washed with PBS, followed by blocking for 30 minutes with PBS + 20% normal bovine serum, 4% bovine serum albumin (Sigma Life Sciences). Sections were then incubated with EGFRL858R rabbit (43B2, Cell signaling) and BRAFV600E mouse (VE1, Spring Biosciences) primary antibodies at 1:50 in blocking buffer overnight at 4 °C in a humidified chamber. A negative control in which the primary antibodies were omitted was included. Slides were washed 3x with PBS and then incubated with secondary antibodies Alexa-fluor 488 goat anti-rabbit IgG and Alexa-fluor 594 goat anti-mouse IgG secondary antibodies (Life Technologies) (1:200 in blocking buffer) for 60 minutes at room temperature. Slides were then washed 3x with PBS and coverslipped with Prolong Gold Antifade with DAPI (Invitrogen). Images were obtained using the Eclipse C1si spectral confocal laser scanning microsope (Nikon) with identical laser settings for pre-treatment and resistant tumor sections.

Fluorescence in situ hybridization (FISH)

FISH analysis was performed using probes targeting MET, EGFR, BRAF and CEP 6 as previously described45,46. Slides were analyzed with an epi-fluorescence microscope equipped with DAPI single band-pass, Aqua single band-pass, Gold single band-pass and Red single band-pass filter sets (Abbott Molecular Inc.). Tumor-bearing areas were identified using the DAPI filter at low magnification 10x lens) with the aid of hematoxylin and eosin stained sections. All areas involving tumor cells were scanned and within this area the cells that did not overlap and had bright signals for the probes were used for enumeration. For identification of focal amplification, the signals targeting BRAF, MET and EGFR were compared to Centromere 6 (CEP6) signal as in Supplemental Fig. 1 with a ratio greater than 2 considered evidence of amplification47.

All methods and experimental protocols were performed in accordance with UCSF laboratory biosafety regulations under a protocol approved by the UCSF Biosafety Committee. Informed consent was obtained for the analysis of the human clinical specimens in accordance with UCSF institutional regulations under an IRB-approved research protocol approved by the UCSF Committee on Human Research.

Additional Information

How to cite this article: Jonsson, V. D. et al. Novel computational method for predicting polytherapy switching strategies to overcome tumor heterogeneity and evolution. Sci. Rep. 7, 44206; 10.1038/srep44206 (2017).

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Material

Supplementary Information:

Acknowledgments

The authors acknowledge funding support from NIH, the Pew-Stewart Charitable Trust, the Kinship-Searle Foundation, and the Van Auken Foundation (to TGB). C.M.B. was supported by grants from The Lung Cancer Research Foundation (P0060520/A123243) and AACR (14-40-18-BLAK). J.C.D. provided funding for V.D.J. and N.M.

Footnotes

T.G.B. is a consultant to Astrazeneca, Novartis, Array, Revolution Medicines and has received research funding from Ignyta. C.M.B. has received funding from Clovis Oncology, Ignyta, and MedImmune.

Author Contributions T.G.B., V.D.J. and C.M.B. conceived and designed the study. V.D.J. conceived and developed the math model and control theoretic algorithm. J.C.D. advised on control theoretic problem approach. V.D.J. and N.M. implemented the algorithm. C.M.B. and V.D.J. designed experiments, performed experiments and analyzed data. V.O., L.L. and E.P. performed experiments. L.L., and B.C.B. performed sequencing. S.A., B.C.B., and B.S.T. analyzed sequencing data. M.A.G. provided tumors for analysis. V.D.J., C.M.B. and T.G.B. wrote the manuscript, with input from all authors. We thank Nikoletta Sidiropoulos for pathology assessment and independent confirmation of the BRAF V600E mutation. We thank Tyrrell Nelson for assistance in sequencing library preparation. We thank Swapna Vemula for assistance with FISH analysis. We thank Russell Johnson for assistance with figure formatting.

References

  • Swanton C. & Govindan R. Clinical implications of genomic discoveries in lung cancer. New England Journal of Medicine 374, 1864–1873 (2016). [Abstract] [Google Scholar]
  • Bivona T. G. & Doebele R. C. A framework for understanding and targeting residual disease in oncogene-driven solid cancers. Nature medicine 22, 472–478 (2016). [Abstract] [Google Scholar]
  • Sawyers C. L. The 2011 Gordon Wilson lecture: overcoming resistance to targeted cancer drugs. Transactions of the American Clinical and Climatological Association 123, 114 (2012). [Europe PMC free article] [Abstract] [Google Scholar]
  • Venkatesan S. & Swanton C. Tumor evolutionary principles: How intratumor heterogeneity influences cancer treatment and outcome. In American Society of Clinical Oncology educational book/ASCO. American Society of Clinical Oncology. Meeting vol. 35, e141–9 (2015). [Abstract] [Google Scholar]
  • de Bruin E. C. et al.. Spatial and temporal diversity in genomic instability processes defines lung cancer evolution. Science 346, 251–256 (2014). [Europe PMC free article] [Abstract] [Google Scholar]
  • Zhang J. et al.. Intratumor heterogeneity in localized lung adenocarcinomas delineated by multiregion sequencing. Science 346, 256–259 (2014). [Europe PMC free article] [Abstract] [Google Scholar]
  • Alizadeh A. A. et al.. Toward understanding and exploiting tumor heterogeneity. Nature medicine 21, 846–853 (2015). [Europe PMC free article] [Abstract] [Google Scholar]
  • Chmielecki J. et al.. Optimization of dosing for EGFR-mutant non–small cell lung cancer with evolutionary cancer modeling. Science translational medicine 3, 90ra59–90ra59 (2011). [Europe PMC free article] [Abstract] [Google Scholar]
  • Bozic I. et al.. Evolutionary dynamics of cancer in response to targeted combination therapy. Elife 2, e00747 (2013). [Europe PMC free article] [Abstract] [Google Scholar]
  • Zhao B. et al.. Exploiting temporal collateral sensitivity in tumor clonal evolution. Cell 165, 234–246 (2016). [Europe PMC free article] [Abstract] [Google Scholar]
  • Gatenby R. A., Silva A. S., Gillies R. J. & Frieden B. R. Adaptive therapy. Cancer research 69, 4894–4903 (2009). [Europe PMC free article] [Abstract] [Google Scholar]
  • de Pillis L. G., Gu W. & Radunskaya A. E. Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations. Journal of theoretical biology 238, 841–862 (2006). [Abstract] [Google Scholar]
  • Zhao B., Pritchard J. R., Lauffenburger D. A. & Hemann M. T. Addressing genetic tumor heterogeneity through computationally predictive combination therapy. Cancer discovery 4, 166–174 (2014). [Europe PMC free article] [Abstract] [Google Scholar]
  • Zhao B., Hemann M. T. & Lauffenburger D. A. Intratumor heterogeneity alters most effective drugs in designed combinations. Proceedings of the National Academy of Sciences 111, 10773–10778 (2014). [Europe PMC free article] [Abstract] [Google Scholar]
  • Bahrami K. & Kim M. Optimal control of multiplicative control systems arising from cancer therapy. Automatic Control, IEEE Transactions on 20, 537–542 (1975). [Google Scholar]
  • Swan G. W. & Vincent T. L. Optimal control analysis in the chemotherapy of igg multiple myeloma. Bulletin of mathematical biology 39, 317–337 (1977). [Abstract] [Google Scholar]
  • dOnofrio A., Ledzewicz U., Maurer H. & Schättler H. On optimal delivery of combination therapy for tumors. Mathematical biosciences 222, 13–26 (2009). [Abstract] [Google Scholar]
  • Chareyron S. & Alamir M. Mixed immunotherapy and chemotherapy of tumors: Feedback design and model updating schemes. Journal of theoretical biology 258, 444–454 (2009). [Abstract] [Google Scholar]
  • Helena A. Y., Riely G. J. & Lovly C. M. Therapeutic strategies utilized in the setting of acquired resistance to EGFR tyrosine kinase inhibitors. Clinical Cancer Research 20, 5898–5907 (2014). [Europe PMC free article] [Abstract] [Google Scholar]
  • Gao S. P. et al.. Mutations in the EGFR kinase domain mediate STAT3 activation via IL-6 production in human lung adenocarcinomas. The Journal of clinical investigation 117, 3846–3856 (2007). [Europe PMC free article] [Abstract] [Google Scholar]
  • Wagle N. et al.. High-throughput detection of actionable genomic alterations in clinical tumor samples by targeted, massively parallel sequencing. Cancer discovery 2, 82–93 (2012). [Europe PMC free article] [Abstract] [Google Scholar]
  • Ohashi K. et al.. Lung cancers with acquired resistance to EGFR inhibitors occasionally harbor braf gene mutations but lack mutations in KRAS, NRAS, or MEK1. Proceedings of the National Academy of Sciences 109, E2127–E2133 (2012). [Europe PMC free article] [Abstract] [Google Scholar]
  • Tricker E. M. et al.. Combined EGFR/MEK inhibition prevents the emergence of resistance in EGFR-mutant lung cancer. Cancer discovery 5, 960–971 (2015). [Europe PMC free article] [Abstract] [Google Scholar]
  • Pao W. et al.. Acquired resistance of lung adenocarcinomas to gefitinib or erlotinib is associated with a second mutation in the EGFR kinase domain. PLoS Med 2, e73 (2005). [Abstract] [Google Scholar]
  • Turke A. B. et al.. Preexistence and clonal selection of MET amplification in EGFR mutant NSCLC. Cancer cell 17, 77–88 (2010). [Europe PMC free article] [Abstract] [Google Scholar]
  • Gazdar A. F. & Minna J. D. Deregulated EGFR signaling during lung cancer progression: mutations, amplicons, and autocrine loops. Cancer Prevention Research 1, 156–160 (2008). [Europe PMC free article] [Abstract] [Google Scholar]
  • Bonavia R. et al.. Tumor heterogeneity is an active process maintained by a mutant EGFR-induced cytokine circuit in glioblastoma. Genes and development 24, 1731–1745 (2010). [Europe PMC free article] [Abstract] [Google Scholar]
  • Wilson T. R. et al.. Widespread potential for growth-factor-driven resistance to anticancer kinase inhibitors. Nature 487, 505–509 (2012). [Europe PMC free article] [Abstract] [Google Scholar]
  • Solca F. et al.. Target binding properties and cellular activity of afatinib (BIBW 2992), an irreversible ErbB family blocker. Journal of Pharmacology and Experimental Therapeutics 343, 342–350 (2012). [Abstract] [Google Scholar]
  • Thakur M. D. et al.. Modelling vemurafenib resistance in melanoma reveals a strategy to forestall drug resistance. Nature 494, 251–255 (2013). [Europe PMC free article] [Abstract] [Google Scholar]
  • Lankheet N. A. et al.. Concentrations of erlotinib in tumor tissue and plasma in non–small-cell lung cancer patients after neoadjuvant therapy. Clinical lung cancer 16, 320–324 (2015). [Abstract] [Google Scholar]
  • Newman A. M. et al.. An ultrasensitive method for quantitating circulating tumor dna with broad patient coverage. Nature medicine 20, 548 (2014). [Europe PMC free article] [Abstract] [Google Scholar]
  • Oxnard G. R. et al.. Association between plasma genotyping and outcomes of treatment with osimertinib (azd9291) in advanced non–small-cell lung cancer. Journal of Clinical Oncology, JCO667162 (2016). [Europe PMC free article] [Abstract] [Google Scholar]
  • Chabon J. et al.. Circulating tumour dna profiling reveals heterogeneity of EGFR inhibitor resistance mechanisms in lung cancer patients. Nature communications 7, 11815 (2016). [Europe PMC free article] [Abstract] [Google Scholar]
  • Newman A. M. et al.. Integrated digital error suppression for improved detection of circulating tumor dna. Nature biotechnology (2016). [Europe PMC free article] [Abstract] [Google Scholar]
  • Lanman R. B. et al.. Analytical and clinical validation of a digital sequencing panel for quantitative, highly accurate evaluation of cell-free circulating tumor dna. PloS one 10, e0140712 (2015). [Europe PMC free article] [Abstract] [Google Scholar]
  • Anderson A. R., Weaver A. M., Cummings P. T. & Quaranta V. Tumor morphology and phenotypic evolution driven by selective pressure from the microenvironment. Cell 127, 905–915 (2006). [Abstract] [Google Scholar]
  • Scott J. G. et al.. Mathematical modeling of the metastatic process. In Experimental Metastasis: Modeling and Analysis 189–208 (Springer, 2013). [Google Scholar]
  • Li H. & Durbin R. Fast and accurate short read alignment with burrows–wheeler transform. Bioinformatics 25, 1754–1760 (2009). [Europe PMC free article] [Abstract] [Google Scholar]
  • McKenna A. et al.. The genome analysis toolkit: a mapreduce framework for analyzing next-generation dna sequencing data. Genome research 20, 1297–1303 (2010). [Europe PMC free article] [Abstract] [Google Scholar]
  • Cibulskis K. et al.. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples. Nature biotechnology 31, 213–219 (2013). [Europe PMC free article] [Abstract] [Google Scholar]
  • Imielinski M. et al.. Mapping the hallmarks of lung adenocarcinoma with massively parallel sequencing. Cell 150, 1107–1120 (2012). [Europe PMC free article] [Abstract] [Google Scholar]
  • Bivona T. G. et al.. FAS and NF-[kgr]B signalling modulate dependence of lung cancers on mutant EGFR. Nature 471, 523–526 (2011). [Europe PMC free article] [Abstract] [Google Scholar]
  • Zhang Z. et al.. Activation of the AXL kinase causes resistance to EGFR-targeted therapy in lung cancer. Nature genetics 44, 852–860 (2012). [Europe PMC free article] [Abstract] [Google Scholar]
  • Gerami P. et al.. Fluorescence in situ hybridization (FISH) as an ancillary diagnostic tool in the diagnosis of melanoma. The American journal of surgical pathology 34, 1146–1156 (2009). [Abstract] [Google Scholar]
  • Dalton M. AJ et al.. Use of fluorescence in-situ hybridization (FISH) to distinguish intranodal nevus from METastatic melanoma. The American journal of surgical pathology 34 (2010). [Europe PMC free article] [Abstract] [Google Scholar]
  • Helena A. Y. et al.. Analysis of tumor specimens at the time of acquired resistance to EGFR-TKI therapy in 155 patients with EGFR-mutant lung cancers. Clinical Cancer Research 19, 2240–2247 (2013). [Europe PMC free article] [Abstract] [Google Scholar]
Articles from Scientific Reports are provided here courtesy of Nature Publishing Group

Full text links 

Read article at publisher's site: https://doi.org/10.1038/srep44206

Read article for free, from open access legal sources, via Unpaywall: https://www.nature.com/articles/srep44206.pdfUnpaywall

Citations & impact 

Impact metrics

19
Citations
Jump to Citations

Citations of article over time

Alternative metrics

Altmetric item for https://www.altmetric.com/details/17305319
Altmetric
Discover the attention surrounding your research
https://www.altmetric.com/details/17305319

Article citations

Go to all (19) article citations

Data 

Data behind the article

This data has been text mined from the article, or deposited into data resources.

BioStudies: supplemental material and supporting data

Similar Articles 

To arrive at the top five similar articles we use a word-weighted algorithm to compare words from the Title and Abstract of each citation.

Funding 

Funders who supported this work.

NCI NIH HHS (1)