Drug resistance is the single most important driver of cancer treatment failure for modern targeted therapies, and the dialog between tumor and stroma has been shown to modulate the response to molecularly targeted therapies through proliferative and survival signaling. In this work, we investigate interactions between a growing tumor and its surrounding stroma and their role in facilitating the emergence of drug resistance. We used mathematical modeling as a theoretical framework to bridge between experimental models and scales, with the aim of separating intrinsic and extrinsic components of resistance in BRAF-mutated melanoma; the model describes tumor–stroma dynamics both with and without treatment. Integration of experimental data into our model revealed significant variation in either the intensity of stromal promotion or intrinsic tissue carrying capacity across animal replicates. Cancer Res; 77(19); 5409–18. ©2017 AACR.

Major Findings

Through the integration of a simple mathematical model with in vitro and in vivo experimental growth dynamics of melanoma cell lines (both with and without drug), we were able to dissect the relative contributions of intrinsic versus environmental resistance. Our study revealed significant heterogeneity in vivo, indicating that there is a diversity of either stromal promotion or tumor-carrying capacity under targeted therapy. We believe this variation may be one possible explanation for the heterogeneity observed across patients and within individual patients with multiple metastases. Therefore, quantifying this variation both within in vivo model systems and in individual patients could have a significant impact on the design of future treatment strategies that target both tumor and stroma. Furthermore, we present guidelines for building more effective and long-lasting therapeutic strategies utilizing our experimentally calibrated model. These strategies explicitly consider the protective nature of the stroma and utilize inhibitors that modulate it.

Quick Guide to Equations and Assumptions

The tumor is classified into two subpopulations, with respect to their sensitivity to the targeted inhibitor. S and R are, respectively, drug-sensitive and drug-tolerant populations. The stroma is divided into normal cells F (i.e., fibroblasts) and reactive cells A (i.e., CAFs). The latter compartment represents fibroblasts in a transformed, secretory phenotype that promotes survival and tumor growth under drug treatment. We assume that S grows in the absence of treatment with growth rate ρs, R grows under targeted treatment at rate ρR. They share a carrying capacity K, representing the maximum packing capacity of the tissue where the tumor is growing. Targeted therapy [BRAF inhibitor (BRAFi)] induces the stroma to switch to its reactive form at a rate θ. In turn, reactive stroma (A) will promote tumor growth by an additional growth rate η. Upon removal of the targeted inhibitor, stromal renormalization occurs at rate ϕ and cancer cells are resensitized at rate ξ. The stromal-targeted inhibitor (FAKi) is assumed to reduce the stromal promotion by rate α. These interactions occur dynamically in time (t) as defined by the following system of ordinary differential equations:

formula

In addition we use the initial conditions: |$S( 0 ) = {S_0},R( 0 ) = {R_0},F( 0 ) = {F_0},A( 0 ) = {A_0}.$| Note that, g(t), h(t), and f(t) are binary functions of time that allow for specific terms in the equations to be switched on and off, depending on treatment scheduling. Given a protocol calling for targeted therapy for the time interval |$[ {t_i^{TT},t_f^{TT}} ]$| and FAKi for |$[ {t_i^{FF},t_f^{FF}} ]$|⁠, the binary functions are defined as follows.

formula

A useful measure of tumor burden control over a window of time [tA, tB] is the inverse of the AUC, defined as follows:

formula

In the past decade, many molecular targets of oncogenic drivers have been developed and approved for the treatment of pathway-specific cancers, in the hope that they could accompany or even replace highly toxic chemotherapeutic drugs (1–4). Unfortunately, this strategy turned out to be only partially successful, with strong initial responses often followed by relapse (3). In an attempt to improve these poor long-term responses, combinations of multiple inhibitors (including immunotherapies) have been attempted (5–7). Despite successes in concurrent inhibition of several pathways in preclinical models (8, 9), it would seem that in the clinical setting, combination of targeted inhibitors does not offer cure, but can at best delay inevitable disease progression caused by the onset of drug resistance (2, 10, 11).

In an effort to understand why these treatment strategies fail, and how we might redesign better and more successful treatments, we must embrace the reality that cancer is a complex evolving system. Because cancer is an evolutionary disease, it can evolve strategies to override or circumvent the action of a given inhibitor. These strategies include producing secondary mutations (11) or exploiting preexisting genetic heterogeneity. However, mutations alone are not sufficient to explain the often rapid time scale over which cancer stops responding to therapy (12). Recent evidence suggests that cancer is able to coopt the surrounding stroma to create an environment that can facilitate treatment escape (13, 14). This phenomenon is termed environment-mediated drug resistance (EMDR; ref. 12), and includes several processes ranging from cell adhesion–mediated drug resistance (15–17) to therapy-induced secretomes, such as IGF, HGF, TGFβ (8, 18), and fibronectin (19). The mechanisms of context-driven resistance we consider here are shared across a variety of solid tumors characterized by aberrations in growth-control signaling and a high level of interaction with the surrounding tumor microenvironment. Our primary focus here is on BRAF-mutated melanoma. A particular instance of EMDR in melanoma is represented by the action of cancer-associated fibroblasts (CAF) that create a habitat favorable for drug tolerance and tumor growth. The environmental remodeling includes deposition of extracellular matrix components, upregulation of growth factor production, intensification of paracrine signaling between the stroma and the tumor cells, and rewiring of the tumor cells' proliferative and survival signaling via integrin binding (12). The effect of this transformed habitat on the cancer and stromal cells is transiently induced by application of the targeted drug and is mostly reversible (20). Given the transient nature of EMDR, there may be an opportunity to modulate it through treatment holidays by allowing renormalization of the stroma to occur, potentially facilitating a better overall treatment outcome. In addition, preliminary investigations have shown benefits in inhibiting stromal-derived processes, such as elevated FAK signaling (13). Dual targeting of tumor and stromal processes represents a promising strategy for better management of BRAF-mutated melanoma.

Understanding this complex interplay between tumor and host cells undergoing treatment is ideally suited for mathematical and computational models. Recently, several theoretical studies have addressed the role of the environment in facilitating drug resistance. Mumenthaler and colleagues have studied how gradients of nutrients and drug concentration modulate the fitness of drug-sensitive and drug-resistant cell lines, and eventually determine recurrence (21). Sun and colleagues modeled the environmental adaptation to drug treatment via drug-induced resistance factors that modulate the growth dynamics of metastatic disease (22). Silva and colleagues and, more recently, Robertson-Tessi and colleagues modeled microenvironmental heterogeneity, specifically the regulation of metabolism, to understand the evolutionary dynamics driving treatment response and leading to resistance (23, 24). A significant literature already exists for mathematical models of intrinsic resistance in cancer progression and response to treatment. Lavi and colleagues offer a comprehensive review of models of cancer resistance (25). However, the focus of the majority of these models is limited to intrinsic chemotherapy resistance (26). Models that integrate the role of the stroma, which is key in the emergence of resistance to targeted therapeutics, are less well developed, but are beginning to emerge. Many studies analyze the dynamics emerging from tumor-immune interactions (27–30). Fewer mathematical models specifically describe interactions between cancer and stromal fibroblasts and their role in drug resistance (31–34). To our knowledge, the problem of separating intrinsic resistance from EMDR, through the dynamics of response to targeted therapy, has not yet been addressed.

Here, we present a first minimal model of tumor–stroma interactions, which aims to bridge the growth dynamics of cancer from in vitro and in vivo experimental models. Specifically, by using exponential growth dynamics from cells growing in vitro, we can calibrate baseline unconstrained growth dynamics. Then, using the same cancer cell line in vivo (mouse allograft), we can capture the saturation dynamics. Using these two experimental model systems, under treatment, we can then quantify the relative contribution of the environment to tumor growth.

Our calibrated model describes the baseline growth dynamics and the relevant tumor–stroma interactions determining growth and response to treatment. This, in turn, allows a fuller exploration of the role of stroma in the promotion of drug resistance, which we propose is critical for the design of optimal treatment strategies. To this end, we will explore treatment schedules that exploit tumor–stroma interactions to limit and/or delay the emergence of EMDR. Our study gives preliminary guidelines for building more effective and long-lasting therapeutic strategies, including dose fractionation and timing.

A common paradigm for the treatment of advanced stage BRAF-mutated melanoma includes targeted therapy in the form of a BRAF inhibitor (BRAFi), such as vemurafenib, recently approved for patients carrying the V600E mutation (35). Kinase inhibitors, such as vemurafenib, specifically block a molecular pathway that the cancer cells are strongly dependent on, resulting in reduced toxicity for the whole body and increased specificity for the tumor. Although this treatment can keep the cancer in check for many months, the disease will eventually recur. Having identified the environment as a key factor in therapy failure (12), alternative blockades of stromal-derived processes are actively being investigated. Here, we specifically model FAK inhibition (FAKi) that has proved effective in the preclinical setting (13).

We propose a model of EMDR for molecularly targeted cancers. Figure 1 shows a schematic of the interactions between key players in our system: cancer cells classified as either sensitive or tolerant to the targeted drug (S and R, respectively), and stroma cells in normal or reactive form (F and A, respectively). The R compartment accounts for an initial intrinsic resistant cancer population as well as cells that are transiently drug-tolerant through the action of EMDR. It is worth noting that this “catch all” compartment does not correspond to a single biological phenotype or genotype; however, it allows us to analyze growth regimes with and without targeted treatment, and most importantly to quantify the relative contribution of the environment to tumor growth dynamics under treatment. Significant bacterial literature indicates the existence of persister phenotypes that are tolerant to a number of antibiotic agents and yet do not appear to be driven by genetic changes (36). Very recently, such populations of “cancer persister cells” have been discovered in an EGFR+ lung cancer cell line (37, 38). However, in the absence of more detailed data, we develop a simplified model with an initial R population that includes cells derived from any of these mechanisms, and allow all cells to return to sensitivity, irrespective of resistance mechanism.

Figure 1.

Interactions hypothesized in the compartmental model. Positive interactions are represented by arrows, negative ones by flat ends. The BRAFi (targeted to the tumor) inhibits growth in the drug-sensitive portion of the tumor (S) and induces activation of normal stroma (F). In turn, reactive stroma (A) promotes growth in the drug-tolerant portion of the tumor (R). The stroma-targeted inhibitor FAKi dampens the effect of stromal-induced growth promotion. Upon removal of BRAFi, the tumor reacquires sensitivity to the drug and the stroma renormalizes (gray arrows).

Figure 1.

Interactions hypothesized in the compartmental model. Positive interactions are represented by arrows, negative ones by flat ends. The BRAFi (targeted to the tumor) inhibits growth in the drug-sensitive portion of the tumor (S) and induces activation of normal stroma (F). In turn, reactive stroma (A) promotes growth in the drug-tolerant portion of the tumor (R). The stroma-targeted inhibitor FAKi dampens the effect of stromal-induced growth promotion. Upon removal of BRAFi, the tumor reacquires sensitivity to the drug and the stroma renormalizes (gray arrows).

Close modal

The interactions between the cell compartments, modulated by the two drugs (BRAFi and FAKi), are defined by a set of ordinary differential equations (ODE; Eq. A) discussed in more detail in the guide to equations. A key advantage of this simple model is that it can incorporate data from both in vitro and in vivo experimental models.

Figure 2 shows the experimental data for BRAF-mutated melanoma cell lines 5555 and 4434. These cells were both cultured in vitro (Fig. 2A) and injected in vivo (Fig. 2B). Growth was observed over time, both in the absence of drug and under treatment with PLX4720, a BRAFi. We can adapt the model (Eq. A) to represent each one of these experimental conditions. Table 1 shows a summary of the experimental conditions and corresponding models. Starting from the in vitro experimental setup, corresponding to a simplified system of equations with fewer unknown parameters, we obtain parameter estimates by data fitting and consequently use these values for the data fitting of the in vivo experimental setup. In doing so, we significantly reduce the number of unknown parameters for each fit, as well as the risk of overfitting.

Figure 2.

Unpacking the relative contributions of intrinsic resistance and extrinsic environment conferred tolerance (EMDR). A,In vitro data and fit. For each condition, we obtained one estimate that best fits the three replicates at the same time. B,In vivo data and fit. Data consist of six and four untreated 5555 and 4434 mice, respectively, and six and five BRAFi-treated 5555 and 4434 mice, respectively. Only a few representative mice are shown. For each condition, the model is fitted individually to each replicate (mouse). Solid and dashed lines correspond to untreated and BRAFi-treated tumor, respectively. Note different y-axis scale for the two cell lines. Data from ref. 13.

Figure 2.

Unpacking the relative contributions of intrinsic resistance and extrinsic environment conferred tolerance (EMDR). A,In vitro data and fit. For each condition, we obtained one estimate that best fits the three replicates at the same time. B,In vivo data and fit. Data consist of six and four untreated 5555 and 4434 mice, respectively, and six and five BRAFi-treated 5555 and 4434 mice, respectively. Only a few representative mice are shown. For each condition, the model is fitted individually to each replicate (mouse). Solid and dashed lines correspond to untreated and BRAFi-treated tumor, respectively. Note different y-axis scale for the two cell lines. Data from ref. 13.

Close modal
Table 1.

Summary of experimental conditions, corresponding model, and set of parameters to be estimated

Summary of experimental conditions, corresponding model, and set of parameters to be estimated
Summary of experimental conditions, corresponding model, and set of parameters to be estimated

The in vitro setup (with a time scale on the order of a few days, Fig. 2A) can be represented by an exponential growth regime, and lacks the stromal component. This corresponds to reducing system (Eq. A) for small time t with F0 = 0, obtaining:

formula

where the only unknown parameters are: R0, S0, ρS for the untreated case (g = 0), and R0, S0, ρR for the treated case (g = 1). Parameter estimation for these triplets is carried out by approximate Bayesian computation, which builds a discrete approximation of the posterior distribution. Data are fitted to the analytic solution of Eq. D. Analytic solutions are reported in Table 1, and a detailed description of the estimation method is reported in the Supplementary Material. Figure 3A shows the marginal distributions for the growth rates of each cell line. Comparing the estimates for control and treated conditions, we see a reduction in growth rate for the treated cancer. The deficit in growth rate reveals that under drug treatment, the R population, irrespective of the mechanism of resistance, exhibits slower growth compared with the S population in untreated conditions, consistent with the previous literature (20).

Figure 3.

Approximated posterior distribution of estimated parameters. Violin plots show probability density functions (x-axis) of parameter estimates (y-axis). A, Estimates for cancer growth rates: ρS (untreated cancer), ρR (treated cancer). The bar highlights the fitness cost of intrinsic resistance. B, Estimates for K reveal heterogeneity of carrying capacity across mice. ρS from previous estimate (Table 2). C, Estimates for |$\tilde{\eta }$| reveal heterogeneity of stromal-derived protection found in vivo. |$\theta = 0.03\ 1/{\rm day}.\,\,{\rho _R}$| from previous estimates (Table 2). D, Estimates for |$\tilde{\alpha }$| for ten 5555 mice treated with BRAFi and FAKi combination. |$\theta = 0.03\ 1/{\rm day}.\,\,\rho _R = 0.49539\ 1/{\rm day} \cdot \tilde{\eta } = 12.67\ 1/{\rm day}$| (average of previous estimates; Table 2).

Figure 3.

Approximated posterior distribution of estimated parameters. Violin plots show probability density functions (x-axis) of parameter estimates (y-axis). A, Estimates for cancer growth rates: ρS (untreated cancer), ρR (treated cancer). The bar highlights the fitness cost of intrinsic resistance. B, Estimates for K reveal heterogeneity of carrying capacity across mice. ρS from previous estimate (Table 2). C, Estimates for |$\tilde{\eta }$| reveal heterogeneity of stromal-derived protection found in vivo. |$\theta = 0.03\ 1/{\rm day}.\,\,{\rho _R}$| from previous estimates (Table 2). D, Estimates for |$\tilde{\alpha }$| for ten 5555 mice treated with BRAFi and FAKi combination. |$\theta = 0.03\ 1/{\rm day}.\,\,\rho _R = 0.49539\ 1/{\rm day} \cdot \tilde{\eta } = 12.67\ 1/{\rm day}$| (average of previous estimates; Table 2).

Close modal

We assume that the growth dynamics of the cancer cells treated in vitro can be solely attributed to preexisting drug-tolerant subpopulations. On the other hand, to quantify the role of the environment on the dynamics of resistance, we turn to the mouse allografts. When the same cell lines are injected in mice, growth is significantly constrained and experiments cover a longer time scale (Fig. 2B). The observed dynamics are more accurately captured with a logistic growth regime, as described by:

formula

By assuming that cells from the same cell line grow at the same exponential rate in an unconstrained environment, we are able to use the growth rates estimated from the in vitro data (i.e., and ρS and ρR) to help calibrate the parameter estimates for the in vivo model.

By fitting the model to the untreated mice data, we obtain estimates for the parameters R0, S0, and K. Note that in the absence of treatment (g = 0), the equations for the tumor and stromal populations are decoupled; therefore, the estimate of tissue-carrying capacity (K) is independent of the quantification of interacting stromal cells. However, K is intrinsically dependent on nutrient constraints as well as the packing capacity of the tissue. Indeed, variations of this quantity are captured in the range of estimated values (see Fig. 3B and Table 2).Posterior distributions are wider in mice with higher values of carrying capacity (e.g., mice I, II cell line 5555). For low K, logistic curves reach carrying capacity within the time window of the in vivo experiments. Curves with higher K, however, have a later inflection point, and their characteristic shape is not captured in the same time window, resulting in more uncertain estimates. It is worth noting that for some mice, the data do not capture the saturating dynamics, as the experiment had to be interrupted due to animal welfare (for details on the original experiments, see ref. 13).

Table 2.

Estimated parameter values for 5555 and 4434 cell lines

ParameterCell lineMouseEstimateRange of parameter values
|{{\rm{\rho }}_{\rm{S}}}{\rm{\ }}( {1/{\rm{day}}} )$| 5555  0.66325 0.6264–0.7002 
 4434  0.4652 0.4336–0.5011 
|{{\rm{\rho }}_{{\rm{R\ }}}}{\rm{\ }}( {1/{\rm{day}}} )$| 5555  0.49543 0.4170–0.5922 
 4434  0.23942 0.2008–0.2942 
|{\rm{K\ }}( {{\rm{m}}{{\rm{m}}^3}} )$| 5555 9,318.8984 8344.3248–10435.095 
  II 8,352.5141 7587.9093–9264.7683 
  III 5,222.2471 5086.7338–5336.5553 
  IV 2,600.883 2540.8448–2672.0427 
  2,345.1481 2316.3892–2371.1414 
  VI 1,475.6726 1449.7419–1504.5806 
  Average 4,818.62  
 4434 2,290.201 2232.7464–2334.0511 
  II 1,673.7862 1612.0517–1749.9593 
  III 215.6403 206.5229–225.9717 
  IV 1,989.7505 1919.6823–2067.3783 
  Average 1,543.6147  
|{\tilde{\eta }\ }( {1/{\rm{day}}} )$| 5555 VII 0.12569 0.000906–0.41951 
  VIII 42.6592 39.73–45.02 
  IX 26.876 24.9689–28.4394 
  0.07974 0.001035–0.22109 
  XI 1.9371 1.2953–2.63 
  XII 4.3748 4.0746–4.6022 
 4434 0.52429 0.36648–0.6682 
  VI 0.34886 0.25731–0.44138 
  VII 7.54 6.1132–11.9923 
  VIII 0.8568 0.65997–1.3352 
  IX 0.10473 0.003147–0.20759 
|\tilde{\alpha }\ \ }( {1/{\rm{day}}} )$| 5555 XIII 14.6163 14.2953–14.9547 
  XIV 15.0018 14.441–15.5339 
  XV 14.3818 14.182–14.5668 
  XVI 14.7552 14.2099–15.4362 
  XVII 14.0455 13.6516–14.3181 
  XVIII 13.9436 13.6503–14.1428 
  XIX 14.2911 13.7735–14.859 
  XX 14.1929 13.6946–14.5971 
  XXI 14.4799 14.2231–14.729 
  XXII 14.5452 14.082–15.1437 
  Average 14.42533  
ParameterCell lineMouseEstimateRange of parameter values
|{{\rm{\rho }}_{\rm{S}}}{\rm{\ }}( {1/{\rm{day}}} )$| 5555  0.66325 0.6264–0.7002 
 4434  0.4652 0.4336–0.5011 
|{{\rm{\rho }}_{{\rm{R\ }}}}{\rm{\ }}( {1/{\rm{day}}} )$| 5555  0.49543 0.4170–0.5922 
 4434  0.23942 0.2008–0.2942 
|{\rm{K\ }}( {{\rm{m}}{{\rm{m}}^3}} )$| 5555 9,318.8984 8344.3248–10435.095 
  II 8,352.5141 7587.9093–9264.7683 
  III 5,222.2471 5086.7338–5336.5553 
  IV 2,600.883 2540.8448–2672.0427 
  2,345.1481 2316.3892–2371.1414 
  VI 1,475.6726 1449.7419–1504.5806 
  Average 4,818.62  
 4434 2,290.201 2232.7464–2334.0511 
  II 1,673.7862 1612.0517–1749.9593 
  III 215.6403 206.5229–225.9717 
  IV 1,989.7505 1919.6823–2067.3783 
  Average 1,543.6147  
|{\tilde{\eta }\ }( {1/{\rm{day}}} )$| 5555 VII 0.12569 0.000906–0.41951 
  VIII 42.6592 39.73–45.02 
  IX 26.876 24.9689–28.4394 
  0.07974 0.001035–0.22109 
  XI 1.9371 1.2953–2.63 
  XII 4.3748 4.0746–4.6022 
 4434 0.52429 0.36648–0.6682 
  VI 0.34886 0.25731–0.44138 
  VII 7.54 6.1132–11.9923 
  VIII 0.8568 0.65997–1.3352 
  IX 0.10473 0.003147–0.20759 
|\tilde{\alpha }\ \ }( {1/{\rm{day}}} )$| 5555 XIII 14.6163 14.2953–14.9547 
  XIV 15.0018 14.441–15.5339 
  XV 14.3818 14.182–14.5668 
  XVI 14.7552 14.2099–15.4362 
  XVII 14.0455 13.6516–14.3181 
  XVIII 13.9436 13.6503–14.1428 
  XIX 14.2911 13.7735–14.859 
  XX 14.1929 13.6946–14.5971 
  XXI 14.4799 14.2231–14.729 
  XXII 14.5452 14.082–15.1437 
  Average 14.42533  

NOTE: For each in vitro condition, we obtained one growth rate estimate that best fits the three replicates at the same time. For each in vivo condition, the model was fitted to each replicate (mouse) to obtain individual estimates of carrying capacity, stromal promotion, and stromal inhibition. Range of values of approximated posterior distributions are also reported.

For the treated mice setup (g = 1), the equations are coupled. We can solve the last two equations of Eq. E analytically, to write A as a function of F0 and θ. Defining |$\tilde{\eta } = \eta {F_0},$| we reduce the parameter number in the analytic solution of Eq. E. At this stage, experimental quantification of the rate of stromal activation is not available; therefore, estimates for the parameters |${R_0},\,{S_0},\,\,\tilde{\eta }$| will be carried out with a range of θ values. We observed high sensitivity of the estimates of |$\tilde{\eta }$| to variations in this experimentally undefined parameter θ (Supplementary Fig. S1). Figure 3C shows the estimated values for |$\tilde{\eta }$| for each mouse. This reveals considerable variation in the stromal support across mice, hinting at an underlying heterogeneity in stromal habitats and activation. As estimates of K and |$\tilde{\eta }$| are dependent on the previously estimated growth rates (ρS and ρR, respectively), the ABC estimation was run for values of growth rates within the range captured by the fit to the in vitro data (see Table 2). The resulting posterior distributions varying in relation to the growth rates are shown in Supplementary Figs. S2 and S3, respectively. The variation in response to BRAF inhibition across replicates could be the result of underlying heterogeneity either in tissue-carrying capacity or in stromal support, or both. Our estimation protocol for the stromal promotion parameter |$\tilde{\eta }$| makes use of an average carrying capacity K previously estimated. However, the variation across replicates could also be explained by variation in carrying capacity. Therefore, we further investigated the BRAFi-treated mice data, to infer the posterior distribution of |$\tilde{\eta }$| as K is varied and vice versa. Supplementary Figure S4 shows the resulting posterior distributions in the |$( {K,\tilde{\eta }} )$| space for a sample mouse. The posterior distribution of K is highly sensitive to the variation of |$\tilde{\eta }$| (see yellow violin plots), and vice versa. However, the best overall fits of |$\tilde{\eta }$| and K are located in the same region of the space. This means that for a given mouse, we can unequivocally identify a combination of values for the carrying capacity and stromal support that best explains the data.

Finally, we can quantify the inhibiting action of the stromal-targeted drug in the form of |$\tilde{\alpha } = \alpha {F_0},$| fitting data from mice treated with both BRAFi (PLX4720) and FAKi (PF562271) to the following version of the model:

formula

Despite the variability of responses across replicates (see data and fits in Supplementary Fig. S5), the resulting estimates for |$\tilde{\alpha }$| show little variation (Fig. 3D; Table 2). This implies that the variability in treatment response may be attributed to the heterogeneous stromal composition of the tissue (highlighted in Fig. 3C), as opposed to the efficacy of the stromal inhibition.

Calibrating our model across in vitro and in vivo data allows us to gain insight into the dynamics of the system that a qualitative analysis of these experiments cannot capture. Figure 3A shows the marginal posterior distribution for growth rates ρS and ρR, with a reduction of the latter quantifying the impact that drug tolerance has on proliferative capacity.

Comparing in vitro and in vivo dynamics allows us to assess the relative contribution of the environment to drug resistance. This analysis revealed significant heterogeneity across replicates (mice), both in terms of tissue-carrying capacity, and stromal protection (Fig. 3B and C). This heterogeneity translates to a high variability of response to treatments that target both the tumor and the stroma, despite the apparent more homogeneous inhibitory effects of the stromal-targeted drug (Supplementary Figs. S5 and Fig. 3D).

Analysis of the ODE model with the combination treatment of BRAFi and FAKi (Eq. F) gives insight into the dynamics of the system as a function of stromal promotion and tumor growth rate. Specifically, we can discriminate two distinct cases:

  • (i) |${\rm If}\ ( {\alpha\, \lt\, \eta } )\ {\rm or }\ ( {\alpha \, \gt\, \eta \ {\rm and}\, {\frac{{{\rho _R}}}{{{{\rm F}_0}( {\alpha - \eta } )}}}\, \gt\, 1} )\ {\rm{ then }}\ {\frac{{{\rm d}( {S + R} )}}{{{\rm dt}}}} \ge 0\forall\ t \ge 0,$|

  • (ii) |${\rm If}\ ( {\alpha\, \gt\, \eta \ {{\rm and\ 0}}\, \lt\, {\frac{{{\rho _R}}}{{{{\rm F}_0}( {\alpha - \eta } )}}}\, \lt\, 1} )\ {\rm{ then }}\ {\frac{{d( {S + R} )}}{{dt}}} \ge 0\forall \ 0 \le t \le {t^{*}},$|

where t* = − |${\frac{1}{\theta}} {\rm log} \, \Big(1 - {\frac{\rho _R} {F_{0}( {\alpha - \eta } )}}\Big)$|⁠.

In the first case, either the stromal promotion is too strong to be compensated by the FAKi, or the stromal promotion is weak, but the tumor growth rate is elevated. Then, the overall tumor burden is monotonically increasing, although bounded by the carrying capacity, and therapy is ineffective. In the second case, when stromal promotion is weak and the tumor growth rate is reduced, then the therapy is effective provided that it is administered for a sufficiently large period of time.

As an example, consider the cohort of 5555 BRAFi-treated mice (VII through XII) and using the parameterized model (Eq. A), with α taken as the average of the previous estimates (see Table 2), we can subclassify the mice. According to our estimates, mice VII, X, XI, and XII fall into case (ii), meaning that with a combination of BRAFi and FAKi, it is possible to achieve control as long as we treat past time t*. On the other hand, mice VIII and IX fall into case (i); hence, the tumor is always growing under treatment, eventually reaching carrying capacity. Supplementary Figures S6 and S7 show a simulated treatment combination of BRAFi + FAKi calibrated on two representative mice, case (i) and case (ii), respectively.

For a tumor–stroma system falling into case (i), recurrence is inevitable, but may be delayed with alternative scheduling strategies. Given that the phenotypic changes underlying EMDR are transient and reversible upon drug removal, we hypothesize that the introduction of drug holidays could significantly improve treatment response and recurrence times. Intermittent application of vemurafenib has proved to be successful in melanoma xenograft models (20), and ongoing clinical trials are testing intermittent versus continuous dosing of a combination of BRAF and MEK inhibitors (NCT02196181). However, we believe that a mechanistic and quantitative approach to treatment scheduling can improve the success of the otherwise empirical approach that these studies offer. We therefore systematically explored the space of holiday versus treatment days of an intermittent schedule treatment with BRAFi, combined with continuous FAKi.

Specifically, the targeted inhibitor is administered during the time windows |$[\scriptstyle {t_i^{TT,k},t_f^{TT,k}} ]{\rm{, }}$| for |$\scriptstyle k \in {ℕ} ^0}$| with |$\scriptstyle t_f^{TT,k} - t_i^{TT,k} = {T_{\rm{T}}},\ {\rm{ }}t_i^{TT,k + 1} - t_{f}^{TT,k} = H{\rm{, }}\ \forall {{k}}\, \ge\, {\rm{0}}{\rm{. }}$| That is, we consider treatments of fixed duration TT, with the time between the end of one treatment and the start of the next treatment being fixed at H. Figure 4 shows the treatment outcome in the holiday versus treatment space (H,TT), where the outcome of each treatment strategy over the time frame of [0, 70] days is quantified with Π, defined in (Eq. C). This reveals that the region corresponding to tumor burden minimization (Π maximization) is concentrated around the line H = 2TT. Intuitively, this means that the length of holiday needed to renormalize the system is proportional to the pulse of treatment. In addition, it indicates that longer treatment holidays are more effective at controlling tumor burden, while the total number of treatment days is reduced.

Figure 4.

Exploration of treat/holiday space for intermittent BRAFi combined with continuous FAKi to maximize control of tumor burden. Surface plot of Π (see Eq. C) in the treat/holiday space. Model parameterized on mouse IX of cell line 5555.

\begin{array}{@{}l@{}}{\rho _S} = 0.66325\ 1/{\rm day},}{\rm{\, }}{\rho _R} = 0.49543\ 1/{\rm day}, {{\rm{\, }}K = 4818.62\,\,\end{array}
\begin{array}{@{}l@{}}m{m^3},{\rm{ }}\tilde{\eta } = 26.876\ 1/{\rm day},{\rm{ \,}} \tilde{\alpha } = 14.4\ 1/{\rm day},{\rm{\, }}\theta = 0.03\ 1/{\rm day},{\rm{\, }}\xi = 0.01\ 1/{\rm day},{\rm{\, }}\varphi = 1\ 1/{\rm day}, \cr {S_0} = 48\,m{m^3},{\rm{ }}{R_0} = 12\,m{m^3},{\rm{ }}{F_0} = 60\,m{m^3},{\rm{ }}{A_0} = 0\,m{m^3}. \end{array}
The star indicates the treatment schedule simulated in Fig. 5.

Figure 4.

Exploration of treat/holiday space for intermittent BRAFi combined with continuous FAKi to maximize control of tumor burden. Surface plot of Π (see Eq. C) in the treat/holiday space. Model parameterized on mouse IX of cell line 5555.

\begin{array}{@{}l@{}}{\rho _S} = 0.66325\ 1/{\rm day},}{\rm{\, }}{\rho _R} = 0.49543\ 1/{\rm day}, {{\rm{\, }}K = 4818.62\,\,\end{array}
\begin{array}{@{}l@{}}m{m^3},{\rm{ }}\tilde{\eta } = 26.876\ 1/{\rm day},{\rm{ \,}} \tilde{\alpha } = 14.4\ 1/{\rm day},{\rm{\, }}\theta = 0.03\ 1/{\rm day},{\rm{\, }}\xi = 0.01\ 1/{\rm day},{\rm{\, }}\varphi = 1\ 1/{\rm day}, \cr {S_0} = 48\,m{m^3},{\rm{ }}{R_0} = 12\,m{m^3},{\rm{ }}{F_0} = 60\,m{m^3},{\rm{ }}{A_0} = 0\,m{m^3}. \end{array}
The star indicates the treatment schedule simulated in Fig. 5.

Close modal

Figure 5 shows the temporal dynamics for one of the best combination treatment schedules predicted by our model. FAKi is continuously administered, and helps control the tumor burden when EMDR sets in, whereas BRAFi is given periodically for 1 day, then off for 2 days. This treatment induces only minimal stromal activation and delays progression by approximately 10 days when compared with the untreated tumor. When compared with the continuous treatment, this intermittent treatment delays progression by approximately 20 days, while using a third of the amount of BRAFi. Although this study does not explicitly account for drug toxicity, total dose reduction is a desirable outcome, especially in the case of combination therapy, where resulting toxicity might be a significant issue.

Figure 5.

Example of combination therapy schedule (BRAFi + FAKi) to exploit tumor–stroma interactions. The model is parameterized as reported in Fig. 4. Bands above the graph indicate the BRAFi and FAKi administration windows. The BRAFi is intermittently administered for 1 day with 2-day holiday. The FAKi is continuously administered. This treatment schedule delays the disease progression by approximately 10 days (compare solid and dashed S(t) + R(t) curves) while using a third of the BRAFi dose.

Figure 5.

Example of combination therapy schedule (BRAFi + FAKi) to exploit tumor–stroma interactions. The model is parameterized as reported in Fig. 4. Bands above the graph indicate the BRAFi and FAKi administration windows. The BRAFi is intermittently administered for 1 day with 2-day holiday. The FAKi is continuously administered. This treatment schedule delays the disease progression by approximately 10 days (compare solid and dashed S(t) + R(t) curves) while using a third of the BRAFi dose.

Close modal

Molecularly targeted therapies for cancers with known driver mutations are extremely effective for 6 to 8 months (e.g., vemurafenib for BRAF V600E melanoma; ref. 39) and are accompanied by lower toxicity when compared with cytotoxic chemotherapeutic agents (3). However, with continuous and prolonged treatment, the emergence of drug resistance seems to be inevitable. Upon removal of the targeted drug, due to relapse, a typical disease flare is observed (e.g., EGFR-mutated lung cancer treated with a combination of tyrosine kinase inhibitors; ref. 40), suggesting that the treatment has somehow selected for a more aggressive clonal population in the tumor. However, subsequent treatment with the same inhibitor often leads to an additional response (41, 42), suggesting that selection of resistant clones alone cannot explain this disease etiology. The environment is now considered an important source of nonintrinsic drug resistance mechanisms (43), collectively referred to as EMDR. As the changes accompanying EMDR are considered transient and therefore reversible, the possibility of regulating EMDR dynamics with smarter treatment scheduling is promising. However, a necessary first step toward the design of such treatment strategies is a more quantitative understanding of the interactions and dynamics occurring between the tumor and the stroma.

In vitro model systems can accurately quantify temporal tumor growth and treatment response in controlled environments, whereas in vivo models more readily capture the native environment that is directly relevant to patients. However, both of these are models of human disease and only capture specific aspects of reality over very specific spatial and temporal scales. The ODE model we develop here bridges between these experimental scales, to integrate relevant information from each of them.

Starting from analysis of BRAF-mutated melanoma cell lines, we quantified the baseline dynamics of cancer cells in a uniform nutrient-rich environment. By comparing the growth rates of cells untreated and treated with the BRAF inhibitor, we were able to quantify the overall reduction of growth under drug application. Our model facilitates this analysis by classifying the cancer into two separate populations, growing with or without drug (R and S). Then, using approximate Bayesian computation, we calculate plausible regions of parameter values. This type of estimation can be particularly useful when assessing the error in fitting. Our parameter estimation method does not make assumptions on the initial conditions, and R0 and S0 are included in the parameters to be estimated. Consequently, the model is agnostic to the mechanisms producing the initial resistant population, R0. These mechanisms could be EMDR related as well as epigenetic or nonautonomous. However, at this stage, no data are available to distinguish between these instances of resistance, and we group all cells that grow under drug treatment in the R compartment, irrespective of the underlying mechanisms of resistance.

Subsequent analysis of data from mice xenografts implanted with the same cell lines allowed us to identify the relative contribution of the environment to drug resistance. This analysis revealed heterogeneity in both the local tissue-carrying capacity and in the stromal promotion of tumor growth. This heterogeneity may be one possible explanation for the spectrum of response observed across patients. In the context of metastatic disease, with tumors seeded across a variety of tissues, heterogeneity in stromal composition could be an important discriminating factor in the success of a systemic treatment. Therefore, quantifying this variation in individual patients could have a significant impact on the design of future treatment strategies that target both the tumor and stroma.

Within the current experimental and modeling framework, assessing the strength of stromal protection is nontrivial. This quantity is dependent on the abundance of the interacting stroma (we could only estimate the overall promotion rate |$\tilde{\eta } = \eta {F_0}$|⁠) as well as the speed of drug-induced stromal activation (we found high sensitivity to parameter θ). At the same time, with the available data, we can explain the variability of responses across mice by variation in carrying capacity and/or stromal promotion (Supplementary Fig. S4). Further investigation of the heterogeneity that our study revealed would require additional experimental quantification of these stromal-related processes. This would, in turn, allow us to address the main shortcoming of the current model, namely the high sensitivity of the estimate of stromal protection to the parameter θ (Supplementary Fig. S1).

Analysis of our ODE model revealed that the degree of stromal protection η, and cancer proliferation ρR under drug treatment, are key in discriminating between responses to the combined action of inhibitors targeting tumor and stromal processes (BRAFi and FAKi, respectively). We found that for slower growing tumors, it is possible to keep growth in check provided treatment with BRAFi is applied for a sufficient period of time. Conversely, for fast growing tumors or elevated stromal protection, the tumor burden increases, despite the administration of the inhibitors. However, for these tumors, we can exploit the transient nature of the EMDR-associated mechanisms and delay progression. Specifically, scheduling treatment holidays for the mouse- (patient-) specific calibrated model would allow for renormalization of the system directly translating into better disease burden control. We used our parameterized ODE model to explore the space of intermittent treatment strategies, with the hope of improving response in cancers falling into the treatment refractory category. Neglecting toxicity of targeted drugs, we searched the space of holiday versus treatment length for intermittent BRAFi application, combined with continuous FAKi. We found that most effective tumor control is achieved with short BRAFi treatment pulses and longer holidays, requiring significantly less inhibitor, when compared with the continuous treatment.

It is worth noting that in optimizing the treatment schedule for these inhibitors, we are only modulating the dynamics by reducing the emergence of EMDR. This allows us to delay recurrence by approximately 10 days. If we were to combine this strategy with a cytotoxic treatment, such as chemotherapy, which provides additional reduction of the tumor burden, then recurrence could be further delayed (44). However, to consider additional treatments for combination therapies, it is necessary to account for toxicity of the single agents, as well as toxicity resulting from their combination. The latter would impose an additional constraint in the optimization problem. Here, we made no assumption regarding the toxicity of both inhibitors and therefore allowed any length of continuous targeted drug administration. Nevertheless, it is worth noting that the intermittent drug treatment we propose not only delays progression but also uses only a third of the drug, when compared with continuous treatment.

The heterogeneity our study revealed from the in vivo experiments highlights the importance in accounting for mouse- (human-) specific microenvironmental parameters to accurately capture response dynamics. This heterogeneity is often ignored in preclinical models, as they aim at establishing general relationships of causality between biological mechanisms. However, as our study suggests, heterogeneity can be key in explaining the variation observed across replicates of an experimental system. Furthermore, models that exploit the transient nature of EMDR must rely on individually calibrated dynamics to propose effective and improved treatment strategies.

Although this study has been focused on melanoma, our model is also applicable to the treatment of other molecularly targeted tumors, such as non–small cell lung cancer. Within the practical constraints of frequency in monitoring a patient's systemic tumor burden and tissue characteristics, our simple model could be used to drive patient- (and tumor-) specific treatment strategies that target both the tumor and stroma. In addition, our approach is ideally suited to directly inform the design of adaptive therapies (45).

E. Sahai has received speakers bureau honoraria from Novartis. No potential conflicts of interest were disclosed by the other authors.

Conception and design: N. Picco, A.R.A. Anderson

Development of methodology: N. Picco, P.K. Maini, A.R.A. Anderson

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

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): N. Picco, A.R.A. Anderson

Writing, review, and/or revision of the manuscript: N. Picco, E. Sahai, P.K. Maini, A.R.A. Anderson

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

Study supervision: P.K. Maini, A.R.A. Anderson

The authors gratefully acknowledge useful discussions with Dr. Gary Mirams, Ross Johnston, and Dr. Robert Jenkins. We also thank Dr. Mark Robertson-Tessi for critical reading and feedback on earlier versions of this manuscript, and Dr. Eishu Hirata for useful insight into the data and experimental techniques.

N. Picco and A.R.A. Anderson were supported by NCI grant U01CA151924. N. Picco was supported by UK Engineering and Physical Sciences Research Council (EPSRC grant number EP/G037280/1).

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

1.
Sawyers
C
. 
Targeted cancer therapy
.
Nature
2004
;
432
:
294
7
.
2.
Sosman
JA
,
Kim
KB
,
Schuchter
L
,
Gonzalez
R
,
Pavlick
AC
,
Weber
JS
, et al
Survival in BRAF V600-mutant advanced melanoma treated with vemurafenib
.
N Engl J Med
2012
;
366
:
707
14
.
3.
Chapman
PB
,
Hauschild
A
,
Robert
C
,
Haanen
JB
,
Ascierto
P
,
Larkin
J
, et al
Improved survival with vemurafenib in melanoma with BRAF V600E mutation
.
N Engl J Med
2011
;
364
:
2507
16
.
4.
Zhou
C
,
Wu
YL
,
Chen
G
,
Feng
J
,
Liu
XQ
,
Wang
C
, et al
Erlotinib versus chemotherapy as first-line treatment for patients with advanced EGFR mutation-positive non-small-cell lung cancer (OPTIMAL, CTONG-0802): a multicentre, open-label, randomised, phase 3 study
.
Lancet Oncol
2011
;
12
:
735
42
.
5.
Larkin
J
,
Ascierto
PA
,
Dreno
B
,
Atkinson
V
,
Liszkay
G
,
Maio
M
, et al
Combined vemurafenib and cobimetinib in BRAF-mutated melanoma
.
N Engl J Med
2014
;
373
:
23
4
.
6.
Larkin
J
,
Chiarion-Sileni
V
,
Gonzalez
R
,
Grob
JJ
,
Cowey
CL
,
Lao
CD
, et al
Combined nivolumab and ipilimumab or monotherapy in untreated melanoma
.
N Engl J Med
2015
;
371
:
1867
76
.
7.
Menzies
AM
,
Long
GV
. 
Systemic treatment for BRAF-mutant melanoma: where do we go next?
Lancet Oncol
2014
;
15
:
e371
381
.
8.
Straussman
R
,
Morikawa
T
,
Shee
K
,
Barzily-Rokni
M
,
Qian
ZR
,
Du
J
, et al
Tumor microenvironment induces innate RAF-inhibitor resistance through HGF secretion
.
Nature
2012
;
487
:
500
4
.
9.
Obenauf
AC
,
Zou
Y
,
Ji
AL
,
Vanharanta
S
,
Shu
W
,
Shi
H
, et al
Therapy-induced tumor secretomes promote resistance and tumor progression
.
Nature
2015
;
520
:
368
72
.
10.
Wagle
N
,
Emery
C
,
Berger
MF
,
Davis
MJ
,
Sawyer
A
,
Pochanard
P
, et al
Dissecting therapeutic resistance to RAF inhibition in melanoma by tumor genomic profiling
.
J Clin Oncol
2011
;
29
:
3085
96
.
11.
Van Allen
EM
,
Wagle
N
,
Sucker
A
,
Treacy
D
,
Johannessen
CM
,
Goetz
EM
, et al
The genetic landscape of clinical resistance to RAF inhibition in metastatic melanoma
.
Cancer Discov
2014
;
4
:
94
109
.
12.
Meads
MB
,
Gatenby
RA
,
Dalton
WS
. 
Environmental-mediated drug resistance: a major contributor to minimal residual disease
.
Nat Rev Cancer
2009
;
9
:
665
74
.
13.
Hirata
E
,
Girotti
MR
,
Viros
A
,
Hooper
S
,
Spencer-Dene
B
,
Matsuda
M
, et al
Intravital imaging reveals how BRAF inhibition generates drug-tolerant microenvironments with high integrin β1/FAK signaling
.
Cancer Cell
2015
;
27
:
574
88
.
14.
Marusyk
A
,
Tabassum
DP
,
Janiszewska
M
,
Place
AE
,
Trinh
A
,
Rozhok
AI
, et al
Spatial proximity to fibroblasts impacts molecular features and therapeutic sensitivity of breast cancer cells influencing clinical outcomes
.
Cancer Res
2016
;
76
:
6495
506
15.
Shain
KH
,
Landowski
TH
,
Dalton
WS
. 
Adhesion-mediated intracellular redistribution of c-Fas-associated death domain-like IL-1-converting enzyme-like inhibitory protein-long confers resistance to CD95-induced apoptosis in hematopoietic cancer cell lines
.
J Immunol
2002
;
168
:
2544
53
.
16.
Hazlehurst
LA
,
Argilagos
RF
,
Dalton
WS
. 
β1 integrin mediated adhesion increases Bim protein degradation and contributes to drug resistance in leukemia cells
.
Br J Haematol
2007
;
136
:
269
75
.
17.
White
DE
,
Rayment
JH
,
Muller
WJ
. 
Addressing the role of cell adhesion in tumor cell dormancy
.
Cell Cycle
2006
;
5
:
1756
9
.
18.
Bhowmick
NA
,
Neilson
EG
,
Moses
HL
. 
Stromal fibroblasts in cancer initiation and progression
.
Nature
2004
;
432
:
332
7
.
19.
Fedorenko
IV
,
Abel
EV
,
Koomen
JM
,
Fang
B
,
Wood
ER
,
Chen
YA
, et al
Fibronectin induction abrogates the BRAF inhibitor response of BRAF V600E/PTEN-null melanomacells
.
Oncogene
2016
;
35
:
1225
35
.
20.
Das Thakur
M
,
Salangsang
F
,
Landman
AS
,
Sellers
WR
,
Pryer
NK
,
Levesque
MP
, et al
Modelling vemurafenib resistance in melanoma reveals a strategy to forestall drug resistance
.
Nature
2013
;
494
:
251
5
.
21.
Mumenthaler
SM
,
Foo
J
,
Choi
NC
,
Heise
N
,
Leder
K
,
Agus
DB
, et al
The impact of microenvironmental heterogeneity on the evolution of drug resistance in cancer cells
.
Cancer Inform
2015
;
14
:
19
31
.
22.
Sun
X
,
Bao
J
,
Shao
Y
. 
Mathematical modeling of therapy-induced cancer drug resistance: connecting cancer mechanisms to population survival rates
.
Sci Rep
2016
;
6
:
22498
.
23.
Silva
AS
,
Gatenby
RA
,
Gillies
RJ
,
Yunes
JA
. 
A quantitative theoretical model for the development of malignancy in ductal carcinoma in situ
.
J Theor Biol
2010
;
262
:
601
13
.
24.
Robertson-Tessi
M
,
Gillies
RJ
,
Gatenby
RA
,
Anderson
ARA
. 
Impact of metabolic heterogeneity on tumor growth, invasion and treatment outcomes
.
Cancer Res
2015
;
75
:
1567
79
.
25.
Lavi
O
,
Gottesman
MM
,
Levy
D
. 
The dynamics of drug resistance: a mathematical perspective
.
Drug Resist Updat
2012
;
15
:
90
7
.
26.
Greene
JM
,
Levy
D
,
Fung
KL
,
de Souza
PS
,
Gottesman
MM
,
Lavi
O
. 
Modeling intrinsic heterogeneity and growth of cancer cells
.
J Theor Biol
2015
;
367
:
262
77
.
27.
Pennisi
M
. 
A mathematical model of immune-system-melanoma competition
.
Comput Math Methods Med
2012
;
2012
:
850754
.
28.
DePillis
LG
,
Radunskaya
AE
,
Wiseman
CL
. 
A validated mathematical model of cell-mediated immune response to tumor growth
.
Cancer Res
2005
;
65
:
7950
8
.
29.
Isaeva
OG
,
Osipov
VA
. 
Different strategies for cancer treatment: mathematical modeling
.
Comput Math Methods Med
2009
;
10
:
253
72
.
30.
Eikenberry
S
,
Thalhauser
C
,
Kuang
Y
. 
Tumor-immune interactions, surgical treatment, and cancer recurrence in a mathematical model of melanoma
.
PLoS Comput Biol
2009
;
5
:
e1000362
.
31.
Flach
EH
,
Rebecca
VW
,
Herlyn
M
,
Smalley
KSM
,
Anderson
ARA
. 
Fibroblasts contribute to melanoma tumor growth and drug resistance
.
Mol Pharm
2011
;
8
:
2039
49
.
32.
Kim
E
,
Rebecca
V
,
Fedorenko
IV
,
Messina
JL
,
Mathew
R
,
Maria-Engler
SS
, et al
Senescent fibroblasts in melanoma initiation and progression: an integrated theoretical, experimental, and clinical approach
.
Cancer Res
2013
;
73
:
6874
85
.
33.
Basanta
D
,
Strand
DW
,
Lukner
RF
,
Franco
OE
,
Cliffel
DE
,
Ayala
GE
, et al
The role of TGF-β mediated tumor-stroma interactions in prostate cancer progression: an integrative approach
.
Cancer Res
2009
;
69
:
7111
20
.
34.
Kim
Y
,
Wallace
J
,
Li
F
,
Ostrowski
M
,
Friedman
A
. 
Transformed epithelial cells and fibroblasts/myofibroblasts interaction in breast tumor: a mathematical model and experiments
.
J Math Biol
2010
;
61
:
401
21
.
35.
Bollag
G
,
Hirth
P
,
Tsai
J
,
Zhang
J
,
Ibrahim
PN
,
Cho
H
, et al
Clinical efficacy of a RAF inhibitor needs broad target blockade in BRAF-mutant melanoma
.
Nature
2010
;
467
:
596
9
.
36.
Veening
JW
,
Smits
WK
,
Kuipers
OP
. 
Bistability, epigenetics and bet-hedging in bacteria
.
Annu Rev Microbiol
2008
;
62
:
193
210
37.
Sharma
SV
,
Lee
DY
,
Li
B
,
Quinlan
MP
,
Takahashi
F
,
Maheswaran
S
, et al
A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations
.
Cell
2010
;
141
:
69
80
.
38.
Ramirez
M
,
Rajaram
S
,
Steininger
RJ
,
Osipchuk
D
,
Roth
MA
,
Morinishi
LS
, et al
Diverse drug-resistance mechanism can emerge from drug-tolerant cancer persister cells
.
Nat Commun
2016
;
7
:
10690
.
39.
McArthur
GA
,
Chapman
PB
,
Robert
C
,
Larkin
J
,
Haanen
JB
,
Dummer
R
, et al
Safety and efficacy of vemurafenib in BRAFV600E and BRAFV600K mutation-positive melanoma (BRIM-3): extended follow-up of a phase 3, randomized, open-label study
.
Lancet Oncol
2014
;
15
:
323
32
.
40.
Chaft
JE
,
Oxnard
GR
,
Sima
CS
,
Kris
MG
,
Miller
VA
,
Riely
GJ
. 
Disease flare after tyrosine kinase inhibitor discrontinuation in patients with EGFR-mutant lung cancer and acquired resistance to erlotinib or gefitinib: implications for clinical trial design
.
Clin Cancer Res
2011
;
17
:
6298
303
.
41.
Lee
JC
,
Jang
SH
,
Lee
KY
,
Kim
Y
. 
Treatement of non-small cell lung carcinoma after failure of epidermal growth factor receptor tyrosine kinase inhibitor
.
Cancer Res Treat
2013
;
45
:
79
85
.
42.
Mackiewicz-Wysocka
M
,
Krokowicz
L
,
Kocur
J
,
Mackiewicz
J
. 
Resistance to vemurafenib can be reversible after treatment interruption
.
Medicine
2014
;
93
:
e157
.
43.
Mueller
MM
,
Fusenig
NE
. 
Friends or foes – bipolar effects of the tumor stroma in cancer
.
Nat Rev Cancer
2004
;
4
:
839
49
.
44.
Kim
E
,
Rebecca
VW
,
Smalley
KSM
,
Anderson
ARA
. 
Phase I trials in melanoma: a framework to translate preclinical findings to the clinic
.
Eur J Cancer
2016
;
67
:
213
22
.
45.
Enriquez-Navas
PM
,
Kam
Y
,
Das
T
,
Hassan
S
,
Silva
A
,
Foroutan
P
, et al
Exploiting evolutionary principles to prolong tumor control in preclinical models of breast cancer
.
Sci Transl Med
2016
;
8
:
327ra24
.