Adaptive Therapy for Metastatic Melanoma: Predictions from Patient Calibrated Mathematical Models

Simple Summary

Tumors are composed of different cancer cells with varying degrees of treatment resistance, which compete for a shared resource. Adaptive therapy is an evolution-based treatment approach that exploits this competition between heterogeneous cancer cells. The approach permits a significant number of drug-sensitive cells to survive, with less dose or with treatment breaks, so that they suppress the proliferation of drug-resistant cells via competition. How can one decide when to stop or resume treatment for each patient? This study presents two mathematical models that guide therapy on and off cycles in a patient-specific manner. The models were applied to melanoma patients and predicted patient-specific adaptive therapy schedules that significantly delayed disease progression with far less therapy (in terms of time on treatment) than the current standard of care. The benefits of adaptive therapy varied between patients. Model-based predictive factors were identified to predict the clinical time gain of individual patients.

Abstract

Adaptive therapy is an evolution-based treatment approach that aims to maintain tumor volume by employing minimum effective drug doses or timed drug holidays. For successful adaptive therapy outcomes, it is critical to find the optimal timing of treatment switch points in a patient-specific manner. Here we develop a combination of mathematical models that examine interactions between drug-sensitive and resistant cells to facilitate melanoma adaptive therapy dosing and switch time points. The first model assumes genetically fixed drug-sensitive and -resistant popul tions that compete for limited resources. The second model considers phenotypic switching between drug-sensitive and -resistant cells. We calibrated each model to fit melanoma patient biomarker changes over time and predicted patient-specific adaptive therapy schedules. Overall, the models predict that adaptive therapy would have delayed time to progression by 6–25 months compared to continuous therapy with dose rates of 6–74% relative to continuous therapy. We identified predictive factors driving the clinical time gained by adaptive therapy, such as the number of initial sensitive cells, competitive effect, switching rate from resistant to sensitive cells, and sensitive cell growth rate. This study highlights that there is a range of potential patient-specific benefits of adaptive therapy and identifies parameters that modulate this benefit.

1. Introduction

Current targeted therapies in melanoma are based on continuous treatment. Patients with advanced BRAFV600E mutant melanoma are eligible for targeted therapy with BRAF and MEK-inhibitors, with objective tumor response rates in up to 75% of the patients [1]. Despite impressive initial responses, a majority of patients with metastatic melanoma experience disease progression. Median progression-free survival ranges from 11–15 months [2,3]. In melanoma, a major driver of this resistance is intratumor heterogeneity [4]. Adaptive therapy is an evolutionarily inspired treatment strategy that exploits this heterogeneity, specifically harnessing competition between different cancer cell types for limited resources [5,6]. Standard treatments using maximum tolerated dose (MTD) rapidly remove drug-sensitive populations. If drug-resistant cells are present before treatment [7,8,9], then this aggressive MTD treatment alters the competition between drug-sensitive cells with drug-resistant populations. MTD treatments release resistant cells from the competition with their sensitive counterparts via a competitive release. As a result, the resistant cells rapidly come to dominate the tumor [10,11]. Several studies have shown that phenotypic switching between sensitive and resistance types provides another avenue for acquired resistance during therapy [8,12,13]. Once resistant cells dominate, treatment fails since it no longer prevents tumor growth. The goal of adaptive therapy is to control each patient’s tumor volume by allowing a significant number of sensitive cells to survive through treatment breaks or dose-modulation. The central mechanism of adaptive therapy is competition, in that any remaining sensitive cells can compete with the resistant population, effectively suppressing their growth and significantly extending progression-free survival [14,15,16].

Adaptive therapy has successfully controlled cancers in preclinical xenograft model systems [17,18]. An ongoing metastatic castrate-resistant prostate cancer clinical trial at Moffitt Cancer Center (NCT02415621) has shown that adaptive therapy can delay disease progression for 27 months using just a 53% cumulative drug rate compared to standard of care (MTD) [19]. For melanoma, previous work has shown that intermittent dosing of a BRAF inhibitor, vemurafenib, given to patient-derived mouse model systems, controlled tumor volumes significantly better than continuous therapy [20]. However, this fixed intermittent strategy is less likely to be effective in patients due to inter-patient heterogeneity. In fact, a trial that explored fixed 5 weeks on/3 weeks off dosing of a BRAF/MEK inhibitor regimen, dabrafenib + trametinib, demonstrated a lower median progression-free survival than standard, continuous dosing in patients with metastatic melanoma [21]. Our group has previously shown that a mathematically driven adaptive approach can work in vivo using xenograft mouse models of melanoma [22]. We posit that this personalized treatment approach, whereby mathematical models facilitate optimizing a patient’s treatment regime based on their current tumor state and historic response, will be far more effective.

Several challenges exist in designing adaptive therapies. First, optimizing the timing of treatment withdrawal and re-challenge for each patient is difficult to implement in clinical practice. Second, it is critical to identify predictive factors for selecting patients who will likely benefit the most from this adaptive therapy. The effectiveness of adaptive therapy will vary among patients, as observed in the prostate cancer trial [19]. Several mathematical models have already been developed to address these challenges. Combined with in vitro experiments, an agent-based model predicted that spatially constrained competition drove an effective control of resistant cell growth [14]. Another agent-based model developed by Gallaher et al. assumed a cost of resistance, with resistant cells having a slower growth rate in the absence of treatment [15]. They modeled resistance as a continuous phenotypic trait and revealed three key factors for successful adaptive therapy: the proportion of resistant cells before the start of treatment, the rate of cancer cell migration, and the speed of evolution towards more resistant phenotypes. An ordinary differential equation model developed by Hansen et al. identified thresholds for the level of initial resistance for successful adaptive therapy [16]. West et al. developed an evolutionary game theory model to determine the optimal timing of multi-drug adaptive therapies [23]. Brady-Nicholls et al. developed patient-specific mathematical models that explain prostate cancer inter-tumor dynamics that can guide intermittent androgen deprivation therapy [24]. Mathematical modeling has shown the consequences of having spontaneous versus induced resistance responses to therapy [25]. Strobl et al. showed how cancer cells’ turnover rate within the tumor influences benefits derived from adaptive therapy [26]. More general and more in-depth mathematical analyses by Viossat and Noble predicted the clinical benefits of tumor containment therapy [27].

Here we present two different ordinary differential equation models that assume different modes of competition between sensitive and resistant populations. We then calibrate both models to a cohort of eight patients with metastatic melanoma. By simulating adaptive therapy schedules with these models, we identify predictive factors that correlate with the largest clinical gains compared to continuous MTD therapy. The first model is a Lotka-Volterra competition model (LV) [28], where genetically fixed drug-sensitive and -resistant populations compete for limited resources. The model assumes resistant cell growth inhibition by a drug-sensitive population. The second model considers phenotypic switching (SW). In addition to competition within and between cell phenotypes, this model allows for switching from sensitive to resistant states or vice-versa when treatment is on or off, respectively. We estimated model parameters by minimizing the difference between model prediction and patient data. For each patient, we used longitudinal data from a serologic marker that is used to monitor advanced melanoma.

The cohort of patients had advanced/metastatic melanoma. All were treated with continuous therapy at MTD. Their therapy consisted of BRAF/MEK inhibitors (either vemurafenib + cobimetinib, or dabrafenib + trametinib). Several had disease progression within 6 months of treatment (Figure 1). While melanoma does not have an ideal biomarker of burden, LDH (Lactate dehydrogenase) is clinically used in melanoma treatment decision making as a correlate of tumor burden and cancer dynamics. LDH is the only serologic marker used for monitoring advanced melanoma in the US [29]. Elevated serum LDH is associated with worse outcomes in patients treated with BRAF/MEK inhibitors, based on the results of a pooled analysis of three trials involving dabrafenib/trametinib with over 600 patients [30]. In the cohort of our current study, all patients had an elevated LDH at the start of treatment, and serial LDH levels were measured in blood at baseline and during routine blood draws approximately every 2–4 weeks (Figure S1, Table S1). LDH was used as a biomarker to correlate with melanoma tumor burden, which could only be measured directly in patients with computed tomography (CT) imaging every 2 to 3 months. It is worth noting that an elevated level of LDH appeared to correlate with disease response and progression determined by CT imaging (Figure S1, Table S1). Temporal LDH profiles for each of the eight patients highlight heterogeneous responses, with some rapidly developing therapy resistance (P6–8), while others took longer to progress (P1–5). Please note that the initial tumor burden level is also different across the patients (LDH range: ~300 to 1500).

Model calibration (see methods) to the melanoma patient LDH data provided a suite of parameter sets that fit the patient data equally well, defining a virtual cohort of patients [31]. Using this virtual cohort, we used the models to predict what might have been the patient’s responses to different adaptive therapy schedules. Our results show that adaptive therapies can delay the time to progression up to several months with less cumulative drug dose rate compared to the continuous MTD regime. Furthermore, we identified key model parameters that determine the benefit of adaptive therapy that could be used to select patients suitable for this evolutionary therapeutic approach.

2. Materials and Methods

2.1. Mathematical Modeling

The first model uses Lotka-Volterra (LV) competition Equations (1) and (2) to describe the competition between two distinct cancer cell populations, drug-sensitive (S) and drug-resistant (R) cells (Figure 2A):

$$\frac{dS}{dt}=r_S\left(1-\frac{S+R}{K}\right)S-\delta S,$$

(1)

$$\frac{dR}{dt}=r_R\left(1-\frac{C\ast S+R}{K}\right)R,$$

(2)

where $r_{S, R}$ indicate the intrinsic growth rates of S and R, respectively. The term $\delta$ > 0 imposes a death rate on S due to therapy. In the absence of treatment, we set $\delta$ = 0. Furthermore, we assume that treatment stops any proliferation by sensitive cells, and so we set $r_S=0$ when treatment is on and $r_S>0$ when treatment is off. We made this assumption since the drugs (BRAF/MEK inhibitors) in this study have shown growth inhibition as well as apoptosis promotion [22]. The two populations, S and R, share the same carrying capacity $K$, the maximum size of the tumor due to nutrient and space constraints. The coefficient $C$ scales the degree to which sensitive cells inhibit the growth rate of resistant cells. If C > 1 (or C < 1), then sensitive cells have a greater (or smaller) competitive effect on resistant cells than resistant cells have on themselves.

As a second variant of the therapy model, we consider the case where phenotypic plasticity (3) and (4) allows for switching between drug-sensitive and resistant cell types (SW, Figure 2B):

$$\frac{dS}{dt}=r_S\left(1-\frac{S+R}{K}\right)S-\delta S-\alpha S+\beta R,$$

(3)

$$\frac{dR}{dt}=r_R\left(1-\frac{S+R}{K}\right)R+\alpha S-\beta R.$$

(4)

This model is a direct extension of the first model, and as such, all overlapping parameters are the same. There are two key differences. First, sensitive cells no longer have a scaled effect on resistant cells, and so C = 1. Second, sensitive cells can switch to resistant ones at rate $\alpha$, or resistant cells to sensitive ones at rate $\beta$, depending on whether treatment is on or off. Please note that $\alpha$ is non-zero when treatment is on and zero otherwise, whereas $\beta$ is non-zero when treatment is off and zero when treatment is on.

To more easily visualize and understand the dynamics of these two models, we non-dimensionalized them ($\tau =r_{R}t,\overline{S}=S/K,\overline{R}=R/K)$ and examined their dynamics in the phase plane of $\overline{S} vs. \overline{R}$. In the LV model, the competition coefficient $C$ and initial populations ($\overline{S}\left(0\right), \overline{R}\left(0\right)$) determine the number of intermittent therapy cycles before resistance dominates (Figure 2A). Intermittent therapy results in more on-off treatment cycles when both $C$ and $\overline{S}\left(0\right)$ are large (Figure 2A bottom panel). In the SW model (Figure 2B), intermittent therapy results in more on-off treatment cycles when both the initial sensitive population ($\overline{S}\left(0\right)$) and the switching rate from resistant to sensitive populations ($\beta$) are large (Figure 2B bottom panel). These analyses help focus our attention on the key parameters that might determine the efficacy of adaptive therapy.

2.2. Parameter Estimation

We identified model parameters that minimized the difference between model predictions and patient data (5) (Figure 1). The cost function for this optimization is to minimize $L_2$ norm of the difference,

$$\underset{\underset{\theta }{\to }}{\min }\Vert V\left(t;\underset{\theta }{\to }\right)-L\left(t\right){\Vert }_2^2 ,$$

(5)

where $\underset{\theta }{\to }$ is the entire model parameter set, and $V$ is the predicted total tumor burden ($V=S+R$) at time $t,$ and $L$ is the actual tumor burden of each patient at time $t$. Here, we assume that LDH is directly proportional to S + R, and the basal (normal) level of LDH is 100 units per liter, as LDH is also produced by normal tissue [32]. The LV model parameter set includes intrinsic growth rates, carrying capacity, death rate, and the competition coefficient. Note, we assume that sensitive cells do not divide when treatment is on ($r_S=0,$ with therapy on). All patient data is for continuous treatment (at the time of writing, we do not have intermittent therapy results for such patients). Thus, a parameter set for the LV model is $\theta =\left\{S_0, K, \delta ,r_R, C\right\}$, where $S_0$ is the initial level of LDH produced by sensitive cells and $R_0=LD{H}_0-S_0$, where $LD{H}_0$ is LDH level at a time $t_0,$ is the initial level of LDH produced by resistant cells. In the SW model, the transition rate from resistant to sensitive is assumed to be zero when treatment is on. The parameter set for the SW model is $\theta =\left\{S_0,K, \delta ,r_R, \alpha \right\}$. We employed a steepest descent optimization algorithm with implicit filtering [33] to identify best-fit parameters for both models (see Supplementary Excel file). Parameter estimations were conducted in MATLAB using the implicit filtering algorithm [33]. We generated patient-specific parameter estimates for each model. The fits to the individual patient data are shown in Figures S2 and S3, Figure 3A,D. The distributions of estimated parameters, variability assessment by interquartile ranges, and L₂ norm of the errors are presented in Figures S4−S7. In addition, we perturbed one of the estimated parameter sets 10% and made model predictions with the perturbed parameters (Figures S8 and S9). A few predictions still fit to the data well (close to the data, small error in Figures S8 and S9), while most predictions failed to explain the data (Figures S8 and S9). Across patients, under the LV model, time to progression is negatively correlated with the intrinsic growth rate of resistant cells ($r_R$) and the competition coefficient ($C$) (Figure 3B). Under the SW model, time to progression is negatively correlated with the intrinsic growth rate of resistant cells (Figure 3E).

3. Results

3.1. Adaptive Therapy Delays Time to Progression

Both the LV and SW models provide fits to the patient LDH data (Figure 3A,D, R² = 0.84 and R² = 0.81, respectively, when comparing actual to predicted across all patient data points). Model parameterization generated multiple similar fits to the data with a similar error. Since there is no single ideal fit, we used the top 50 fits for each model and each patient, resulting in a cohort of 400 virtual patients (8 patients; 50 virtual patients for each real patient for each model, LV and SW, respectively). Using this virtual cohort, we simulated both continuous and adaptive therapy for each of the virtual patients. For adaptive therapy simulations, we assumed treatment dose per treatment day is the same as the dose in continuous therapy. We stopped treatment when a patient’s LDH level dropped to 50% of their initial LDH. Treatment remained off until the LDH grew back to its initial level. We estimated time to disease progression as the moment when a virtual patient’s LDH reached 150% of the patient’s initial LDH. For each virtual patient, we determined the time gained from adaptive therapy relative to continuous therapy by subtracting the time to progression under continuous therapy from under adaptive therapy. The LV model predicts that adaptive therapy delayed the progression of patient #1 by 4.5 months on average (Figure 3C. continuous (cyan) vs. adaptive (green)) with an average cumulative dose rate of ~54% of continuous therapy. It is worth noting that the cumulative dose rate in this study refers to the percentage of time on therapy, not the actual dose, since we are simulating treatment on and off days with the same dose (effect of the drug) per unit time (day). Of note, the one free parameter, the intrinsic growth rate of the sensitive population, was set to be the same as that of the resistant population. The SW model predicts a significant improvement from adaptive therapy (~6.8 months) with an average cumulative dose rate of ~46% of continuous therapy (Figure 3F). We fixed two of the free parameters by setting the growth rates of sensitive cells to be the same as the estimated resistant cell growth rates ($r_S=r_R$, when the treatment is off), and the transition rate to be 45% of the resistant cell-intrinsic growth rate ($\beta =0.45r_R$, when the treatment is off). To assess the robustness of model predictions, we simulated adaptive therapy for patient 1 with new choices of $r_R$ and $\delta$. Specifically, we fixed all the other parameters and re-estimated the two parameters $r_R$ and $\delta$. We varied the free parameters $r_S$ ($r_S=$ 5–200% of the new $r_R$) for the LV model and ($r_S\in \left[0.05r_R{}^{new},2r_R{}^{new}\right],\beta \in \left[0.05r_R{}^{new},0.95r_R{}^{new}\right],$ when the therapy is off) for the SW model. The LV model predicted that adaptive therapy delayed progression 3–7.2 months. The SW model predicted up to 22 months of progression delay (Figure S10).

3.2. Key Parameters That Determine Clinical Gains

Both forms of the mathematical model predicted that adaptive therapy would be beneficial, but the patient-specific benefits were substantially different across the eight patients (Figure 4A). The LV model predicted that the average time gained for patient #1 was about 3.5 months, while that for patient #8 was only about one month (Figure 4A).

To further explore the significant differences across the eight patients in the LV model, we simulated patient-specific adaptive therapies using the cohort of 400 virtual patients over a broad range of intrinsic growth rates ($r_S=$ 5–200% of the estimated $r_R$). The estimated competition coefficient, C, and the estimated initial size of the sensitive population ($S_0/K$) are key for determining the benefit (time gained) of adaptive therapy relative to continuous therapy (Figure 4B). In Figure 4B, each colored dot is an average virtual patient of $r_S=$ 5–200% of the estimated $r_R$ per each estimated parameter set obtained from the real patients. If the initial sensitive cell percentage is smaller than 10%, time to tumor progression under adaptive therapy is not significantly different from continuous therapy (time gained $\le$ 1 month, patient #3, #7, #8). The larger the competition coefficient ($C$), the more time was gained (Figure 4B), provided that the initial sensitive cell population ranges from 10% to 40% of carrying capacity. Adaptive therapy is most successful for tumors with a high initial number of sensitive cells (>40% of K, see patients #1). For adaptive therapy, the percent time on treatment varied between patients (i.e., the cumulative dose rate of 12% to 100% of continuous therapy), with an overall average of 65%. In the simulations, we assume treatment dose per treatment day in adaptive therapy is the same as continuous therapy. The drug dose implies time on the drug. The cumulative dose rate of 12% means 12% of the time on the drug compared to continuous therapy.

To explore the significant differences of adaptive therapy response across the eight patients in the SW model, we simulated adaptive therapy using the cohort of 400 virtual patients generated from fitting the SW model to the patient data (Figure 1 and Figure 3D). Adaptive therapy is predicted to delay progression even more (Figure 4C vs. Figure 4A) with this model. Among the eight patients, adaptive therapy increased patient #1’s time to progression by up to 7 months while patient #8 still gained only one month on average when we fixed two of the free parameters $r_S\left(=r_R\right)$ and $\beta \left(=0.45r_R\right)$ when the therapy is off ($r_S=0, \beta =0$, when the therapy is on). We further ran a suite of simulations with the virtual cohort and compared the average time gained. It is worth noting again that most model parameters were estimated from each patient LDH (Figure 1 except for two parameters, the intrinsic growth rate of sensitive cells ($r_S$) and the switching rate from resistant to sensitive during treatment holidays ($\beta$). For these two free parameters ($r_S,\beta$), we considered a range of values ($r_S\in \left[0.05r_R,2r_R\right],\beta \in \left[0.05r_R,0.95r_R\right]$). A larger switching rate ($\beta )$ results in a larger time gained from adaptive therapy (Figure 4D). This effect occurred for all patients. Interestingly, changing the intrinsic growth rate of sensitive cells had little impact on the efficacy of adaptive therapy if the switching rate is low (time gained ≤ 2 months for all $r_S$ if $\beta <0.1r_R$). When the switching rate was larger than $0.1r_R$, a tumor composed of slower growing sensitive cells responded better to adaptive therapy (increasing time gained relative to continuous therapy). The cumulative dose rates of adaptive therapy varied substantially over the ranges of parameters. For large switching and slow-growing tumors, the cumulative dose rate was 20% of continuous therapy, while for low switching rates and fast-growing tumors, the rate was 74% of continuous therapy.

Both model formulations demonstrate that, in general, adaptive therapy can be effective in delaying tumor progression using significantly less drug than continuous therapy. For the LV model, adaptive therapy is most effective if tumors are composed of a sufficiently large number of sensitive cells initially (>40% of carrying capacity K). Adaptive therapy gains superiority as the sensitive cells have a higher competition coefficient (C) and exert more significant inhibition on resistant cells through competition. The SW model highlights that adaptive therapy could be beneficial if drug-sensitive and -resistant states are phenotypically plastic. The predicted time gained is about 20 months if the switching rate from drug-resistant to drug-sensitive cells is large, and the sensitive cell growth rate is slow during treatment holidays. Even if the switching rate is very low, adaptive therapy can delay time to progression from 2–5 months. Taken together, our simulations show the potential benefits of adaptive therapy and identify the key parameters and conditions for adaptive therapy to be superior in controlling tumor volume relative to continuous therapy.

3.3. A Different Treatment-Stopping Criterion

So far, for adaptive therapy, we used 50% of the initial tumor burden as the treatment-stopping criterion. This criterion was implemented in the first clinical trial of adaptive therapy for castrate-resistant prostate cancer [19]. Would a less aggressive treatment stop criterion be better? We simulated a suite of adaptive therapy simulations with a new stopping criterion of 20% reduction to address this question. Therapy is held off if the tumor burden decreases to below 20% of the initial burden. The tumor is re-challenged with therapy when the burden returns to its initial value.

The LV model predicts that a 20% stopping criterion is indeed a better threshold if tumors are already responding well to adaptive therapy. In Figure 4B, we showed that tumors respond well if the initial population size of sensitive cells is high, and the competitive effect of sensitive on resistant cells is large. Under these conditions, the 20% threshold is more effective in controlling tumor burden than the 50% threshold (Figure 5A, Patient #1 under the 20% threshold adaptive therapy gained 7–20 months over continuous therapy vs. Figure 4B, Patient #1 under the 50% threshold gained 3–7 months over continuous therapy). We observed a similar pattern in the SW model. A lower threshold is more effective if tumors are already responding well. If the switching rate from resistant to sensitive cells is large and the intrinsic growth rate of sensitive cells is small, adaptive therapy with either the 50% or 20% threshold substantially delays progression relative to continuous therapy. Under these conditions, an adaptive therapy with the 20% threshold delays progression up to 5 months more (up to 25 months gained Figure 5B vs. up to 20 months gained with 50% threshold Figure 4D). Taken together, a smaller tumor burden reduction criterion may be more effective in delaying tumor progression when the tumors underlying dynamics satisfy specific conditions, such as a large number of initially sensitive cells, strong competitive inhibition of resistant cells by sensitive ones, a high switching rate from R to S, and a slow-growing sensitive population.

3.4. Progression-Free Survival

Finally, we compared the probability of progression-free survival of all virtual patients under continuous, −50%, and −20% stopping adaptive therapies. Since we developed two mathematical models, LV and SW, we generated two sets of Kaplan-Meier curves (Figure 6, A: LV, B: SW). For these K-M curves, we used all virtual patients with either the LV or the SW model. This “trial” group was then subjected to the two adaptive trial procedures (−50% and −20% stopping threshold) and continuous therapy. Here, we assume that tumor progression occurred when LDH levels reached 150% of their initial LDH (Figure 6, A: LV, B: SW). The two models predicted a longer progression-free survival under both adaptive therapies (p < 0.001), with a 20% threshold being superior to the 50% stopping criterion.

4. Discussion

The adaptive abiraterone trial for metastatic castration-resistant prostate cancer (NCT02415621) extended time to progression over 27 months [19]. The success of this trial has inspired three more clinical trials testing for the feasibility of adaptive therapy in patients with castrate sensitive prostate cancer (NCT03511196), advanced melanoma (NCT03543969), and thyroid cancer (NCT03630120). The actual real-world benefit from adaptive therapy would likely vary among patients, as already observed in [19]. Understanding the underlying mechanism for this variability may be crucial for patient selection in future clinical trials. As previous mathematical models have demonstrated [15,16,19,22,23,34], the interaction between the drug-sensitive population and resistant population drives the outcome of adaptive therapy. These competitive aspects have been further investigated experimentally [35], but more work needs to be done to better understand how resistance emerges and is maintained under different treatment strategies. A sensitive cell population can inhibit the growth of a resistant population via competition. Sensitive cells can acquire resistance, and resistant cells can be re-sensitized to the drug. To examine these interactions in the context of metastatic melanoma, we developed two different mathematical models, a standard Lotka-Volterra competition model (LV) with a varying competition coefficient (effect of sensitive cells on resistant cell population) and an extension with phenotype switching (SW) but with all competition coefficients set equal to 1. The LV model describes the competition between distinct sensitive and resistant populations for a limited resource. In contrast, the SW model assumes individual cells can acquire resistance when therapy is on and become re-sensitized when therapy is off.

Our analyses and simulations of these two models demonstrate that adaptive therapy, with timed treatment holidays, prolongs survival compared to standard of care (MTD continuous therapy). We identified key tumor parameters where this benefit is most significant. The LV model shows that adaptive therapy can extend the time to progression significantly whenever the initial population of sensitive cells is large enough, and when sensitive cells have a large competitive effect on resistant cells. When the initial number of a resistant cell population is large, the clinical time gained is less significant (e.g., P4, time gained: ~3 months). An even higher resistant cell initial percentage (e.g., R (0) = 50%) required an increased competition coefficient or a higher transition rate from R to S for a meaningful clinical gain (Figure S4). Under these conditions, we predict 1~2 years of time gained with a relatively large fraction of the R population (~80% at the time of tumor progression, see Figure S11). The SW model illustrates that adaptive therapy is most beneficial whenever the re-sensitization rate of the resistant cell population to the sensitive population is high, and the growth rate of sensitive cells is low. Under these constraints, adaptive therapy can more than double the time to progression (max time of delay is 20 months in adaptive vs. 6 months with continuous therapy). Interestingly, tumors without these properties did not respond well to the standard of care either (<3 months of progression-free survival). In these cases, switching to alternative therapies as soon as possible might be more desirable. For example, in metastatic melanoma, adding or switching to immunotherapy may provide a clinical benefit, particularly for patients expected to respond poorly to either adaptive therapy or continuous therapy. Taken together, our study identified a patient group that may benefit the most from adaptive therapy, with predicted clinical time gains for such patients. This study raises the challenge and the opportunity of directing clinical research towards swiftly measuring the key parameters for a given patient to facilitate truly personalized medicine and deciding on the efficacy of adaptive therapy.

A key assumption of our analysis is that a patient’s initial tumor burden is tolerable and not immediately life-threatening. Maintaining a potentially large tumor burden, as opposed to focusing on shrinking tumors as fast as possible, can make both clinicians and patients uneasy. There are cancers where an initial tumor burden is not tolerable and may be life-threatening (e.g., aggressive glioblastoma or very painful tumors that may be bleeding, etc.). There are, however, situations where a patient’s initial tumor burden is tolerable, and the patient is asymptomatic. In fact, maintaining the sum of all tumor diameters (tumor burden) is one of the response evaluation criteria in the measurement of solid tumors in clinical trials. Stable disease is defined as a <20% increase in the sum of target lesions per the widely used RECIST 1.1 criterion (Response Evaluation Criteria in Solid Tumors) [36]. There is increasing evidence that chronic control of disease burden is more effective at improving patient survival, not only in cancers [19], but also in bacterial infections [37,38,39,40]. Our proof of concept analysis on patients with advanced melanoma treated with BRAF/MEK inhibitor targeted therapy illustrates that maintaining a tolerable tumor burden may delay progression significantly.

The mathematical models presented here are simplified representations of what may be happening in actual cancers under treatment. The models rest on the assumption that two key cancer cell populations compete and interact with each other in a well-mixed environment. In reality, a tumor is not a well-mixed population of cancer cells, but is a spatially stratified population [41,42,43,44] that drives different pharmacokinetics [45] and is impacted by the influence of a heterogeneous and dynamic microenvironment [46,47,48,49]. Such aspects could be incorporated in future studies; however, adding more complexity does not guarantee a better understanding. For now, the clinical data is non-spatial (blood levels of LDH). Therefore, while a spatial model may deliver a better representation of the tumor in a more relevant context, it would be more complex and burdened with additional assumptions that cannot be evaluated by blood biomarkers alone.

We chose our modeling approach as a starting point to better understand how different modes of drug-sensitive and -resistant interactions, distinct sensitive and resistant subpopulations vs. phenotypic switching, impact the outcomes of adaptive therapy. The results presented show the need to better understand the properties of the sensitive and resistant populations. If they are genetically distinct and pre-existing, then for our cohort of eight patients, adaptive therapy would provide less benefits to fewer of the patients than if drug sensitivity and resistance are reversible phenotypically plastic states. These properties of the sensitive and resistant populations must be determined empirically. We have previously examined the behavior of the metastatic melanoma cell line WM164 and found it spears to exhibit this more phenotypically plastic behavior, whereas another melanoma cell line 1205LU did not [22]. Therefore, there may be significant heterogeneity across patients in terms of which resistance mechanisms are at play and how they emerge and are maintained.

Equally as important as the phenotypic heterogeneity is the actual tumor burden—accurate and frequent measurement of tumor burden is key to calibrating mathematical models so that they can better reproduce tumor dynamics and make reliable predictions. In the ongoing adaptive therapy prostate cancer trials, prostate-specific antigen (PSA) is used as a marker for disease burden (NCT02415621 and NCT03511196). LDH is a standardized systemic biomarker for melanoma, although it can also increase in certain conditions not related to tumor burden, such as in patients with liver toxicity from drug treatment [50]. The modeling approach here assumed that this LDH level is linearly correlated with tumor burden. In the cohort of this study, all patients had an elevated LDH level initially, which decreased while responding to therapy, and then increased due to progression (Table S1 and Figure S1). Newer blood biomarkers, such as circulating tumor DNA levels, may correlate better with tumor response and progression in melanoma patients [51]. However, these biomarkers may only highlight the sensitive cell population, and all of them lack spatial information. Development of novel imaging technologies is urgently needed to allow for non-invasive serial assessment of tumor burden; not only would this provide greater temporal resolution but would offer the opportunity to understand the spatial dynamics better.

Our analysis of patients with advanced melanoma identifies tumor-specific conditions (parameters) and resulting dynamics where adaptive therapy may lead to significant clinical gains. Identifying such patients before treatment would accelerate clinical translation. Here, we identify such parameters by calibrating our models to individual patient tumor burden dynamics. This study highlights another potential benefit of using mathematical models in the clinic by supporting patient selection in clinical trials based on tumor parameter identification. Having chosen the appropriate patients, the mathematical models can then be tailored to each patient treatment response to predict and drive their next treatment decision. This mathematical model driven treatment decision paradigm is especially critical for adaptive therapy (and more generally, personalized medicine) since it directly facilitates both therapy dosing and switch points. We therefore advocate for greater integration of predictive mathematical models in the clinical decision process.