Addressing Genetic Tumor Heterogeneity through Computationally Predictive Combination Therapy

Recent high-throughput sequencing and genomic hybridization studies have revealed substantial intratumoral heterogeneity in patients with cancer (1, 2). Sections of single biopsies as well as biopsies taken from primary and metastatic regions revealed a highly complex nonlinear branching clonal evolutionary model as the basis for cancer progression and intratumoral diversity (3–11). Matched diagnostic and relapsed patient samples have also revealed the dynamic nature of a heterogeneous tumor, with the dominant subpopulation at relapse often originating from a minor preexisting subclone (12–16). As such, spatial and temporal intratumoral heterogeneity presents a fundamental challenge for the rational design of combination chemotherapeutic regimens, which remain the primary treatment for most systemic malignancies.

There have been numerous attempts to study tumor heterogeneity and optimal therapeutic strategies. Many have focused on theoretical approaches to examine drug scheduling with generic single drugs or drug combinations on a heterogeneous population containing a sensitive and a resistant subpopulation (17–20). Some of these scheduling strategies have been followed by experimental validations in vitro (21, 22). However, a rational approach to design drug combinations to minimize the effects of intratumoral heterogeneity in a tractable model with experimental validation in vitro and in vivo in preclinical models has been lacking.

Here, we apply a computational optimization algorithm constructed on the experimental foundation of known single-drug efficacies for genetically variant cell subpopulations to predict how drug combinations will affect heterogeneous tumors. The most crucial assumption inherent in our algorithm has been demonstrated in our recent work: that commonly used combinations of chemotherapeutics act as linear averages of each component drug against homogenous tumors (23). For example, a two-drug combination of drug A and drug B creates a combination selective pressure that resembles a simple weighted sum of A + B. Integration of this experimental foundation with our mathematical framework offers an advance in predictive understanding of how drug combinations influence the fate of heterogeneous tumors, which we successfully validated both in vitro and in vivo.

Our experimental system derives from the conceptual premise that tumors undergo branched clonal evolution with deregulated oncogene expression and/or tumor suppressor loss, providing a basis for transformation and additional genetic changes accumulating during tumor progression. These additional mutations underlie the development of a heterogeneous tumor population (Fig. 1A, as a simplified example). An attractive approach to model this heterogeneity in a tractable system amenable to systematic study is the use of an RNA interference (RNAi)–based approach (Fig. 1B). Knockdown of specific genes of interest can approximate the loss-of-function that occurs in tumors, with the combination of multiple short hairpin RNA (shRNA)–expressing subpopulations modeling the diversity present in a heterogeneous tumor. We have previously used an RNAi-based approach to elucidate mechanisms of single and combination drug action (23, 24) using Eμ-Myc; p19^Arf−/− lymphoma cells. Given this dataset of known therapeutic effects of cytotoxic and targeted therapies (Supplementary Table S1) on individual shRNA-expressing subpopulations, we hypothesized that we could computationally predict superior versus inferior treatment strategies for minimizing subpopulations within a heterogeneous tumor (Fig. 1C).

Our mathematical algorithm is based on integer programming, which identifies the set of drugs that should be present to accomplish an aspired goal, or “objective function.” In the first manifestation here, the defined goal is to minimize outgrowth of specific tumor subpopulations within a known heterogeneous population (see Methods and Supplementary Data for mathematical formulations). The rationale for this goal as a first test of our hypothesis is that current clinical approaches tend to focus on targeting the predominant tumor subpopulation, and the outgrowth of originally minor subpopulations representing resistance to that treatment is at present a daunting clinical problem.

We used this algorithm, in combination with a previously published dataset on drug–genotype interactions (23, 24), to identify the most effective two-drug combinations sampled from a large number of three-component heterogeneous tumor populations (tumors that are composed of two RNAi-produced subpopulations plus the parental subpopulation; Supplementary Fig. S1). We then systemically simulated all possible combinations of three-component populations and determined the most and least effective two-drug combinations (Supplementary Table S2). In many cases, the optimal therapy was the same independent of whether we examined the entire heterogeneous population or just the predominant subpopulation. However, there was also a subset of population compositions where the solutions differ (Supplementary Fig. S2). This suggests that for some heterogeneous tumor populations, we cannot derive the optimal therapy based on solely the predominant subpopulation alone. One such example was a three-component population consisting of the parental Eμ-Myc; p19^Arf−/− lymphoma (no shRNA) and subpopulations expressing either a Chk2 (a DNA damage checkpoint regulator) or a Bok (a Bcl2-family cell death mediator) shRNA. The optimal treatment for this tumor was a vincristine (Vin) plus vorinostat (SAHA) combination, whereas the worst was an irinotecan (IRT) plus chlorambucil (CBL) combination. Interestingly, if only shChk2 or shBok populations were examined alone, the predicted optimal combination was neither Vin/SAHA nor IRT/CBL (Fig. 1D). Furthermore, SAHA was not the best single agent for either shChk2 or shBokalone, but becomes part of an optimal drug combination in a heterogeneous population containing both shChk2 and shBok. This was also the case in examining lymphomas containing an alternative population composition of parental/shChk2/shBim. Thus, consideration of a heterogeneous tumor in its entirety can result in nonintuitive optimal drug combinations, sometimes containing drugs that are not the best single agent for any subpopulations.

To experimentally validate these predictions, we used an in vitro fluorescence-based competition assay (Fig. 2A and B and Supplementary Fig. S3), in which a mixture of parental lymphoma cells and shRNA-expressing subpopulations were exposed to combinations of drugs at controlled doses. Specifically, the drug combinations were dosed such that each drug in the combination contributed equally to a cumulative LD80-90 combination cell killing. GFP- or Tomato-labeled subpopulations enabled us to track the enrichment or depletion of individual subpopulations in the population mixture. Here, we observed that Vin/SAHA effectively maximized the therapeutic response in shBok-containing tumor cells (i.e., enhanced the depletion of shBok-infected cells) while minimizing the selective outgrowth of populations of tumor cells expressing shChk2 (Fig. 2C and Supplementary Fig. S4). In contrast, IRT/CBL strongly selected for the resistant shChk2subpopulation, and the parental subpopulation remained sensitive to both combinations. Combination treatment with actinomycin D/erlotinib or Vin/CBL, which were predicted to be optimal treatments if shChk2 or shBok subpopulations were considered individually, were less effective than Vin/SAHA. These predictions were also validated in a different tumor with a distinct population composition (Fig. 2D and Supplementary Fig. S4). Taken together, these data suggest that combination therapies can be tailored to effectively minimize the effects of a heterogeneous tumor population, given some knowledge of the composite subpopulations and their responses to single drugs.

As tumors continuously evolve with the acquisition of new mutations, forming more complex hierarchical structures, we wondered whether we could extend our simple model to account for another layer of complexity. Here, we performed an additional knockdown in our existing subpopulations to form a heterogeneous population consisting of parental, shChk2 (alone), shBok (alone), and shChk2 plus shBok. We then examined the response of these cells to the same four combination therapies and observed that the results were consistent with our predictions, with the optimal therapy still Vin/SAHA (Fig. 2E and F and Supplementary Fig. S4). We can also apply principal component analysis (PCA) as a method for dimensionality reduction and visualization of the resulting population, as we have done previously for single populations (25, 26). We observed that the optimal drug combination effectively impacts the population composition toward the goal of eliminating the shRNA subpopulations (Supplementary Fig. S5). Thus, for this limited set of tested cases, we have demonstrated successful prediction of effective therapeutic strategy design by computational integration of the individual drug effects across known target subpopulations, which can produce nonintuitive outcomes.

We next sought to validate our predictions in vivo. The murine Eμ-Myc lymphoma model (27, 28) allows us to perform ex vivo transduction of tumor cells followed by transplantation into syngeneic immunocompetent recipient mice (Fig. 3A). Lymph node, thymus, and spleen are among several sites of tumor dissemination in both donor and recipient tumor-bearing mice. Using this model, we observed substantial intratumoral heterogeneity by ex vivo whole mouse fluorescence imaging (Fig. 3B and Supplementary Fig. S6). However, when we then analyzed individual lymph nodes, thymus, and spleen of untreated mice using flow cytometry, we observed heterogeneity at the level of individual tumors (Supplementary Fig. S7). To validate the therapeutic effects in vivo, we determined the optimal dose for each of the individual drugs to ensure that there was comparable therapeutic effect on parental tumors from each component drug in the combination (Supplementary Fig. S8). Using these doses, combination treatment with Vin/SAHA successfully minimized the emergence of any tumor subpopulation, while IRT/CBL enriched significantly for the intrinsically resistant shChk2 subpopulation (Fig. 3C and D). Because our combination therapies were optimized to minimize specific tumor subpopulations relative to the parental lymphoma cells, we also examined the tumor-free survival of mice with heterogeneous lymphoma tumors normalized to that of mice with homogeneous lymphoma tumors. Upon treatment, Vin/SAHA improved relative tumor-free survival compared with IRT/CBL (Fig. 3E). Taken together with our in vitro results, these data suggest that we can apply a mathematical optimization approach to predict optimal therapeutic strategies for minimizing the emergence of genetically defined tumor component populations and affect the tumor-free survival of mice when presented with a heterogeneous instead of homogeneous tumor.

We further explored the comparison between Vin/SAHA and IRT/CBL on relative tumor-free survival in shChk2/shBok/parental tumors across wide ranges of subpopulation proportions. With a three-component population, all tumor compositions can be represented on a ternary plot, with each corner corresponding to a homogeneous population of a single component and any point within the plot corresponding to some specific composition mixture of the three components. We experimentally determined the tumor-free survival of mice with homogeneous tumors containing one of the three subpopulations for treatment with either of the combination regimens (Supplementary Figs. S9 and S10). Using this information, we generated a descriptive model showing the efficacy difference between Vin/SAHA and IRT/CBL across all possible subpopulation proportions. We observed that by introducing heterogeneity (shChk2 and/or shBok subpopulations) into the parental tumor, Vin/SAHA increasingly improves relative tumor-free survival (Fig. 4A). We can also derive the comparison in terms of absolute tumor-free survival, to inform us of tumor compositions for which Vin/SAHA dominates over IRT/CBL (Supplementary Fig. S11). In both cases, an example heterogeneous tumor composition (used for the competition assays above), which was not used in generating the model, approximated well with the prediction described by the model.

Because our optimization model predicts the effect of drug combinations on the evolution of tumor composition, and our descriptive model approximates the efficacy difference between treatments at specified tumor compositions, we can combine these models to examine drug efficacy at multiple time points upon single or consecutive drug treatments. Although both drug combinations have minimal effects on changes in tumor composition with empty vector controls (Fig. 4B), the drug treatments predictably affected the heterogeneous tumor (Fig. 4C)—landing at tumor compositions at relapse for which we can also now estimate the therapeutic efficacy comparisons for the next round of treatment. As preexisting subclones before treatment may exist at extremely low frequencies in clinically observed cases, we also examined theoretically the dynamics of tumor heterogeneity with multiple rounds of combination therapy at different initial tumor compositions, including extremes such as 0.1% preexisting shChk2 and shBok. Sequential IRT/CBL strongly selected for the shChk2 subpopulation, whereas Vin/SAHA had less selective pressure on the subpopulations (Supplementary Fig. S12). Thus, based solely on tumor heterogeneity dynamics (assuming no new induced resistance mechanisms), Vin/SAHA remains as a superior effective drug combination, whereas after several rounds of IRT/CBL, the effectiveness of this drug combination can be dramatically reduced because of the outgrowth of shChk2 in becoming the dominant subpopulation.

Taken together, knowledge of the in vitro efficacy of single drugs on individual subpopulations has allowed us to optimize for a drug combination that minimizes the effects of heterogeneity and improves relative tumor-free survival. We combined additional in vivo efficacy data of the chosen drug combinations on individual subpopulations to approximate the in vivo therapeutic responses over wide ranges of tumor subpopulation proportions in a heterogeneous tumor. These two models enable the tracking of tumor trajectories upon single- and multi-course treatments and examination of the effectiveness of subsequent treatment strategies at relapse(s).

Starting with Peter Nowell's seminal 1976 article on the clonal evolution of cancer (29), we have gradually gained the ability to monitor and deconvolute the step-wise alterations in cancer evolution that underlie intratumoral heterogeneity. Not only are distinct tumor subclones found to coexist within the same tumor regions (30, 31), but analyses of biopsies in primary and metastatic tumor sites further suggest that metastatic subclones originate from a nonmetastatic parental clone in the primary tumor (3–11). Additional changes at the posttranscriptional and epigenetic level can potentially further diversify a tumor population with functional variations (32, 33). This heterogeneous tumor population is also dynamic, as has been shown in the responses to standard combination chemotherapeutic regimens, with preexisting minor subclones expanding to dominate at relapse (12–16). As such, current combination regimens can have unpredictable and/or unintended consequences on the resulting tumor diversity. The original rationale for the use and choice of combination therapies has been primarily to increase the effective and tolerable drug dose, while minimizing resistance (34–38). In theory, this leads to more cancer cell killing and decreases the likelihood that a resistant clone will compromise therapeutic success. However, in light of recent studies showing the extent of intratumoral variation and its clinical implications, it is paramount to incorporate tumor diversity and the expected evolutionary trajectories into rational drug combination design to achieve predictable tumor response, reduce chances of relapse, and prolong patient survival. Our joint theoretical and experimental approach suggests that we can derive drug combinations that can minimize outgrowth of specific subpopulations in a given heterogeneous tumor while enhancing tumor-free survival in mice—provided we have some knowledge of the tumor composition and the response of component subpopulations to single drugs.

Here, we used the Eμ-Myc lymphoma mouse model as the basis for our computational modeling and experimental validations. The presence of p19^Arf−/− in this model may remove selective pressure on specific mutations that arise during tumor progression and in response to treatment. Thus, the broad applicability of this approach to other cancer models remains to be seen. Nevertheless, the Eμ-Myc lymphoma model is an extremely well-characterized preclinical mouse model and thus a good initial system to tractably address the rational design of drug regimens in the context of genetically diverse tumors. Furthermore, Eμ-Myc models a defining translocation observed in Burkitt lymphoma. Previously targeted molecular analyses on different sites of the same patient and matched sample at diagnosis and relapse revealed shared MYC breakpoints and BCl6 mutations, but differential TP53 and MYC mutations across the different samples, suggesting clonal evolution spanning across different tumor sites and in response to treatment (39).

The complexity of heterogeneous tumors found in patients undoubtedly complicates applications of this approach. However, our studies highlight basic principles toward rational drug combination design. For example, as was the case for Vin/SAHA in this tumor model, a therapeutic strategy may involve finding a set of drugs that, when combined, maximizes the therapeutic synthetic lethality of certain subpopulation(s), while minimizing the selective pressure for outgrowth of cells bearing undesirable alteration(s). Here, the maintenance of tumor heterogeneity—representing the persistence of a “naïve” pretreatment state—may be preferable to the undesirable outgrowth of specific subpopulation(s). On the other hand, if we can potentially drive toward the enrichment of therapeutic-sensitive subpopulations while preserving overall tumor sensitivity, a sequential treatment approach may be successful as part of a multi-course regimen. Thus, the ability to apply optimization approaches to control the trajectories of tumor composition offers opportunities to maximize the effects of concurrent or successive drug combination regimens.

Murine Eμ-Myc; p19^Arf−/− B-cell lymphomas were cultured in B-cell medium [45% Dulbecco's Modified Eagle Medium (DMEM), 45% Iscove's Modified Dulbecco's Medium (IMDM), 10% FBS, supplemented with 2 mmol/L l-glutamine, and 5 μmol/L β-mercaptoethanol]. The cell line was tested and shown to be free of Mycoplasma using both PCR-based (American Type Culture Collection) and biochemical-based (Lonza) methods. All drugs were obtained from LC Laboratories, Sigma-Aldrich, Calbiochem, or Tocris Biosciences. For in vivo studies, irinotecan, chlorambucil, and vorinostat were dissolved in dimethyl sulfoxide (DMSO) as stock and further diluted in 0.9% NaCl solution before injection. Vincristine was dissolved in 0.9% NaCl solution.

All shRNAs were expressed in either MSCV–LTR–MIR30–SV40–GFP (MLS; ref. 40) or MSCV–LTR–MIR30–PGK–Tomato (MLT) retroviral vectors and were previously validated for knockdown (24). Multiple shRNAs targeting each gene were used for single-drug responses to rule out off-target effects (24). Transfection and transduction were performed as previously described (23).

Single and combination drug treatments in vitro were performed as previously described (refs. 23, 24; see Supplementary Fig. S3). Briefly, lymphoma cells were infected to approximately 20% of the total population with the indicated retroviruses and were dosed with single- or two-drug combinations at concentrations needed to achieve 80% to 90% cell death (LD80-90) of the empty vector control (MLS or MLT) cells (assessed at 48 hours). For combination treatments, each drug was dosed to ensure equal contribution of overall killing in the combination. Cell death was assessed with propidium iodide incorporation. Every 24 hours, untreated and treated cells were diluted at 1:4 and 1:2, respectively. At 72 hours, the percentage of GFP⁺ and/or Tomato⁺ cells were analyzed on a cell analyzer FACSScan, FACSCalibur, or LSRFortessa (BD Biosciences). To assess enrichment/depletion of the subpopulations, a log-transformed Resistance Index (RI; refs. 23, 24) was calculated: log₂ RI = log₂[(L_t − L_tL_u)/(L_u − L_tL_u)], where L_t and L_u represents percentage of labeled cells (GFP⁺ and/or Tomato⁺) out of the total live population with and without treatment, respectively. See Supplementary Data for derivations of log₂ RI. Experiments were performed with indicated shRNAs in both MLS and MLT vectors to rule out any vector-specific effects.

Two million infected Eμ-Myc; p19^Arf−/− lymphoma cells were tail-vein injected into 7 to 9 weeks female syngeneic recipient C57/BL6 mice. Upon palpable tumor formation (∼12–13 days after injection), mice were treated with single or combination drug regimens. Dosing for each drug was determined to ensure comparable efficacy on control uninfected lymphoma cells. IRT was administered at 120 mg/kg intraperitoneally (i.p.) once on day 1, Vin at 1 mg/kg i.p. once on day 1, CBL at 10 mg/kg i.p. once on day 1, and SAHA at 300 mg/kg i.p. on days 1 and 3 and 150 mg/kg on days 2 and 4. Vehicle control consisted of maximum equivalent doses of DMSO and 0.9% NaCl solution. For competition assay, lymphoma cells from thymus, lymph node, and spleen were harvested at sizable palpable tumor after relapse and analyzed on a cell analyzer. For survival studies, tumor-free survival was monitored by daily palpation following treatment.

Fluorescence imaging of mice was acquired using the IVIS imaging system and Living Image software v4.3.1 (PerkinElmer). Two excitations (535 and 465 nm) with nine emission filters (for a pairwise combination of 10) were used. GFP, Tomato, and background signals were spectral unmixed and the resulting composite image further enhanced with brightness and contrast.

The optimization of combination therapies for minimizing heterogeneous subpopulations was mathematically formulated as a binary integer programming problem. All analyses were performed in MATLAB R2012b (MathWorks). Binary integer programming was solved using the bintprog function in the MATLAB Optimization Toolbox. In determining the cases for experimental validation, we sorted the simulation results based on difference in treatment effectiveness (evaluated by log₂ RI) between therapy optimized on the basis of consideration of entire heterogeneity and that of just the predominant subpopulation. We validated cases with large efficacy differences, including nonintuitive cases (e.g., optimal drug combinations for heterogeneity that do not contain component drugs best for predominant subpopulation) and cases with different heterogeneous tumor compositions. Refer to Supplementary Data for more details on the mathematical model for optimization and the descriptive model on the efficacy comparison between Vin/SAHA and IRT/CBL.

Rotated PCA was used to visualize the impact of different drug combinations on tumor composition. The data consist of a single matrix with each row representing a different treatment condition (or a pure subpopulation) and each column representing an averaged log₂ RI for each drug. The averaged log₂ RI was derived from matrix multiplication of the enrichment/depletion dataset matrix, R (see mathematical models section in Supplementary Data) by population composition vector, p, determined experimentally following treatment. Thus, the derived data represent a pharmacologic profile for each resulting population following a specific treatment condition. The data were mean centered and unit variance scaled. PCA was performed using SIMCA-P v11.5 (Umetrics). For ease of interpretation, a linear transformation of the principal component space was performed using MATLAB by rotating the parental population projection completely into the second principal component.

Statistical analyses were performed using Prism v5 (GraphPad). Two-tailed Student t tests and one-way ANOVA (with Bonferroni post hoc test) were used, as indicated. Error bars represent mean ± SEM. Log-rank (Mantel–Cox) test was used for comparison of survival curves.

J.R. Pritchard has ownership interest in a patent. D.A. Lauffenburger is a consultant/advisory board member of Pfizer, Merrimack Pharmaceuticals, and Complete Genomics. No potential conflicts of interest were disclosed by the other authors.

Conception and design: B. Zhao, J.R. Pritchard, D.A. Lauffenburger, M.T. Hemann

Development of methodology: B. Zhao, J.R. Pritchard, D.A. Lauffenburger, M.T. Hemann

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): B. Zhao, J.R. Pritchard

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): B. Zhao, J.R. Pritchard, M.T. Hemann

Writing, review, and/or revision of the manuscript: B. Zhao, J.R. Pritchard, D.A. Lauffenburger, M.T. Hemann

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): J.R. Pritchard

Study supervision: J.R. Pritchard, D.A. Lauffenburger, M.T. Hemann

The authors thank Luke Gilbert and Peter Bruno for advice on experimental procedures, members of the Lauffenburger and Hemann laboratory for helpful discussions, and Peter Bruno, Eleanor Cameron, Miles Miller, and Aaron Meyer for comments on the article.

Funding for this study was provided by Integrative Cancer Biology Program Grant U54-CA112967-06 (to M.T. Hemann and D.A. Lauffenburger). B. Zhao is supported by the NIH/National Institute of General Medical Sciences (NIGMS) Interdepartmental Biotechnology Training Program 5T32GM008334. This work was also supported in part by the Koch Institute Support (core) Grant P30-CA14051 from the National Cancer Institute.