Multiscale Stochastic Modeling Connects Cancer Drug Resistance Mechanisms to Population Survival Rates

Drug resistance significantly limits the long-term effectiveness of targeted therapeutics for cancer patients. Recent experimental studies have demonstrated that cancer cell heterogeneity and microenvironment adaptations to targeted therapy play important roles in promoting the rapid acquisition of drug resistance and in increasing cancer metastasis. The systematic development of effective therapeutics to overcome drug resistance mechanisms poses a major challenge. In this study, we used a modeling approach to connect cellular mechanisms underlying cancer drug resistance to population-level patient survival. To predict progression-free survival in cancer patients with metastatic melanoma, we developed a set of stochastic differential equations to describe the dynamics of heterogeneous cell populations while taking into account micro-environment adaptations. Clinical data on survival and circulating tumor cell DNA (ctDNA) concentrations were used to confirm the effectiveness of our model. Moreover, our model predicted distinct patterns of dosedependent synergy when evaluating a combination of BRAF and MEK inhibitors versus a combination of BRAF and PI3K inhibitors. These predictions were consistent with the findings in previously reported studies. The impact of the drug metabolism rate on patient survival was also discussed. The proposedmodel might facilitate the quantitative evaluation and optimization of combination therapeutics and cancer clinical trial design.


Introduction
Drug resistance places an often inevitable limit on the long-term effectiveness of targeted therapeutics for cancer patients [1,2]. Considerable efforts have been made to combat drug resistance and improve patient survival. Although the underlying molecular and cellular mechanisms are complex, some paradigms of drug resistance mechanisms have been established [3][4][5][6][7][8].
It is widely acknowledged that the inherent heterogeneity [9,10] of cancer cell populations, which is assumed containing both drug-sensitive and drug-resistant cells, contributes to drug resistance and metastasis [11][12][13][14]. A recent study [15] revealed a novel drug resistance mechanism in which drug-sensitive cancer cells secrete various soluble factors (e.g., IGF and HGF) into the tumor microenvironment in response to targeted therapy. These secreted factors can promote the growth, dissemination and metastasis of drug-resistant cancer cells and support the survival of drug-sensitive cells. Therefore, microenvironment adaptation [16] plays an important role in the rapid emergence of acquired drug resistance.
Evaluating cancer therapeutics in the context of tumor heterogeneity and microenvironment adaptation is very complex. In traditional in vitro and in vivo experiments, multiple cell types and multiple drug dosages must be considered, in addition to other experimental conditions and challenges in human population studies. As such, these studies are expensive and time consuming. Therefore the systematic development of effective therapeutics to overcome drugresistance mechanisms has posed a major challenge. Mathematical modeling may potentially serve to bridge molecular/cellular mechanisms of drug resistance and population-level patient survival, and facilitates the quantitative evaluation and optimization of combination therapeutics and cancer clinical trial design.
Many mathematical and computational models have been developed to simulate tumor growth and drug response. For example, the cellular automata model [17,18] or agentbased model [19][20][21], continuum partial differential equations model [22,23] and hybrid discrete-continuum model [24,25] have all been applied to evaluate tumor growth at the molecular, cellular and/or tissue level. These models have substantially advanced our understanding of tumor initiation and progression. However, due to their complexity and/or intensive computing burden, these models have rarely been applied to predict populationscale patient survival. Haeno et al. [26] developed a mathematical framework to describe pancreatic metastasis using a branching process to help understand cancer growth dynamics during metastasis and identify optimal therapeutic interventions. However, this framework focused on genetic mutation-induced drug resistance and did not address the role of targeted therapy-induced microenvironment adaptations in drug resistance. The use of combination therapy has been suggested in cases of drug resistance, such as in advanced melanoma patients with BRAF mutations [15,16]. Therefore, the development of mathematical models capable of quantitatively evaluating synergism in combination drug therapy is desirable.
In this study, we created a multiscale model comprising a set of stochastic differential equations driven by both the Wiener process and Poisson process to describe pharmacokinetics, cellular dynamics, and progression-free survival at the patient level while accounting for microenvironment adaptations (Figure 1). Our model was subsequently verified using population-and cellular-scale clinical data. Then, we evaluated the efficiency and synergy of different combination therapies (combinations of BRAF, MEK and PI3K inhibitors). Our modeling revealed that different patterns of synergy existed for these combinations. Finally, sensitivity analysis revealed several key parameters that may combine with each other to affect the cancer cellular dynamics and patient survival, and facilitates the quantitative evaluation and design of combination therapeutics. In addition, we examined the impacts of different treatment schedules and drug metabolism rates on patient survival.

Cellular dynamics
Tumors are heterogeneous (e.g., in their mutation profiles), resulting in some tumor cells possessing sensitivity to drug therapy and others in the same tumor exhibiting resistance to it. To model growth, transition and dissemination dynamics in drug-sensitive and drugresistant cancer cells in patients with metastatic disease, we employed the following set of stochastic differential equations (SDEs): Targeted therapy (e.g., treatment of melanoma using BRAF kinase inhibitors) is effective on drug-sensitive cells; however, a small number of pre-existing drugresistant cancer cells are unaffected by treatment. In response to drug treatment, drug-sensitive cancer cells secrete various compounds (e.g., IGF and HGF) into the tumor microenvironment. These compounds were termed drug-induced resistance factors (DIRFs) in this study. The secreted DIRFs enhance the growth, dissemination and metastasis of cancer cells [15]. In our mathematical model, the metastatic cells refer to the new metastatic cells after the initiation of drug treatment. (B) A flowchart of our work. We constructed a stochastic model comprised of a set of stochastic differential equations driven by both Wiener and Poisson processes to model pharmacokinetics, DIRF secretion and cellular dynamics based on the recent experiments and clinical data. This enabled us to calculate progression-free survival in silico. Our model was then verified by comparing its predictions against clinical patient survival data. We used the model to quantitatively evaluate the efficiency and synergy of two combination therapies (BRAF inhibitor plus MEK inhibitor and BRAF inhibitor plus PI3K inhibitor). Furthermore, sensitivity analysis revealed several important parameters in the model that may provide implications for the design of combination therapies. In addition, we also examined cellular-and patient-level responses to different drug treatment schedules and investigated the impact of heterogeneity in drug metabolism rates on patient survival.
where C S and C R represent relative numbers (assumed in the unit of 10 8 [26][27][28]) of drugsensitive cancer cells and drug-resistant cancer cells, respectively. The first terms on the righthand side of Eqs. (1) and (2) describe the growth of sensitive cells and resistant cells, respectively.r S andr R are growth rate coefficients associated with these two cell types. The growth of drug-sensitive and drug-resistant tumor cells was assumed to follow a logistic growth law [29,30]. T max represents maximal carrying capacity. The second terms in equation (1 and 2) describe the transition from sensitive cells to resistant cells, e.g., due to genetic/epigenetic mutations. u represents the mutation rate in drug-sensitive cells as they convert to drugresistant cells (i.e., mutation-driven drug resistance). The third term in Eq. (1) describes the drug-induced death of drug-sensitive cells.d S is the death rate of drug-sensitive tumor cells following treatment (e.g., BRAF inhibitors for V600 mutated melanoma) and depends on drug where d S represents maximal death rate, and K Drug is a Michaelis constant representing the drug concentration associated with reaching the half-maximal inhibition effect. The fourth term (also called a diffusion term) in Eq. (1) simulates the stochastic fluctuation of cell numbers and is modeled by the standard Wiener process W that is described as where N(0,1) is a unit normal distribution. The third term in Eq. (2) is similar. σ i (i = 1, 2) represents the diffusion rate. The last terms in equation (1-2) describe the dissemination of existing cancer cells.
Independently of W i , N t represents a Poisson process with intensity λ and describes the count of metastasis within a cancer cell population [31,32]. Specifically, the Poisson process N t (t ≥ 0) is characterized by where λ is the expectation of disseminating cell number within per unit time (Day). In addition, N t has independent increments, and N 0 ¼ 0. In the above equations (1 and 2), both drugsensitive and drug-resistant cancer cells were assumed to have the potential to further metastasize. q M andq M represent the dissemination rates of drug-sensitive and drug-resistant cells, respectively.q M is regulated by drug-induced resistance factors as described below. It should be noted that the metastasized cells in patients before therapy were considered to be included in these sensitive or resistant cells, and a new variable was introduced to account for new metastasis after the initiation of targeted therapy as follows.
Therapy-induced drug resistance can intensify tumor metastasis [15,16]. The growth of new metastatic tumor cells following the drug treatment was modeled using a SDE driven by a jump process as follows: Metastasis from resistant cells Second metastasis (5) where C M represents the number of new metastatic cells after the initiation of new therapy. The first term in Eq. (3) describes the growth of the metastatic cells, and r M is a metastatic cell growth rate coefficient. M max is the maximal carrying capacity of metastatic cell growth. The second term (diffusion term) simulates fluctuation of metastatic cell population as mentioned above. Metastasis from existing cancer and metastatic emissions by the metastases themselves (i.e., secondary metastasis) [33] were taken into account, which were modeled in the last three terms of Eq. (3). q M andq M respectively represent dissemination rates of drug-sensitive and drug-resistant cancer cells as described above. Metastatic rates were assumed to depend on existing tumor size (i.e. drug-sensitive/resistant cancer cell numbers C S and C R ) [34] and angiogenic cell number (C K ) [29]. The drug effect on drug-sensitive metastatic cells was incorporated by using 1 Àd S C S in the third term. By assuming that newly developed metastasis sites are more supportive of the growth of invasive cancer cells, a positive net increase rate of new metastatic cells due to the secondary metastasis was introduced in the last term of the above equation.
Angiogenic growth in tumors is induced by the secretion of angiogenic growth factors (e.g., VEGF). We modeled angiogenesis based on previously established work [35] with the following equation: where C K represents the number of angiogenic cells. The first term describes the growth of angiogenic cells induced by tumor cells.r K is a growth rate coefficient associated with angiogenic cells, and K max is the maximal carrying capacity for blood vessel growth. The second term describes the growth inhibition of angiogenic cells by tumor cells with a coefficient d K .
Newly grown blood vessels can provide tumor cells with nutrients (such as oxygen and glucose) and thus influence the maximal carrying capacity of tumor cells [29,30] as follows:

Pharmacokinetics
Pharmacokinetics describe the dynamics of drug absorption, metabolism and elimination by the body [36]. These processes are often modeled as follows [37,38]: where D represents drug concentration in the body. The first term in the above equation models the first-order elimination rate of drugs, and d drug is a metabolic rate coefficient of patients. U t ð Þ in the second term is the rate of drug delivery. Brownian motion was also assumed in the above equation to accommodate stochasticity [37,38].
The initial conditions of the Eqs. (1-2, 5, 6) were set to C S = 0.2, C R = 0.001 and C K = 0.1, simulating the relative cell number (in the unit of 10 8 [26][27][28]) in patients (e.g., patients with metastatic melanoma and BRAF V600 mutations) before initiation of the new therapy. Starting from the initiation of the drug treatment, the number of new metastatic cells was counted, with an initial value C M = 0. The initial concentration of drug was set to 0. The uniqueness of the solution to the above SDEs (Eqs. (1-7)) were easily obtained, since their coefficients satisfies the appropriate growth conditions and are locally Lipschitz continuous [39]. We employed a timeadapted Euler scheme [40] to provide numerical solutions to the SDEs driven by both diffusion (Brownian motion) and jump (Poisson process).

Microenvironment adaptations to drug treatment
As demonstrated by recent experimental and preclinical studies [15], when BRAF inhibitors (BRAF-I), such as Vemurafenib and Dabrafenib, are administered to cancer patients with BRAF mutations, they can induce drug-sensitive cancer cells to secrete resistance factors (e.g., IGF and HGF) into the tumor microenvironment. We modeled the secretion of drug-induced resistance factors (DIRFs) by drug-sensitive tumor cells according to Michaelis-Menten kinetics [41] as follows: where V DS is the maximal secretion rate of DIRFs from drug-sensitive cells, and K DS is a Michaelis constant representing the drug concentration at which a half-maximal secretion rate is achieved. d DIRF represents DIRF degradation rate.
It was assumed that DIRF secretion and/or degradation in the microenvironment are much faster processes than cellular phenotype switching and shifts in cellular population dynamics. Therefore, using a quasi-steady state assumption we can express the secreted DIRF concentration as follows: Secreted DIRFs can promote outgrowth, dissemination and metastasis in drug-resistant cells and enhance survival in drug-sensitive cells [15,16]. The effects of DIRFs on cellular dynamics were modeled using the following functions, which correlate DIRF concentration to the growth, dissemination and metastasis rates of three types of cancer cells: Eq. (8)  In this way, we linked the short-term timescale (minutes) associated with intercellular signaling to the long-term timescale (days) necessary for cellular dynamics [42][43][44].

Progression-free survival analysis of a patient population
Cancer progression is often clinically evaluated using radiographic imaging. In the below analysis, if a patient's total tumor cell number or tumor volume exceeded a pre-set threshold, C Th (assumed to be 1.6 in this work), then we considered the patient's cancer to be progressing. Therefore, progression-free survival (PFS) time (T PFS ) was defined as the length of time between initiation of therapy (t = 0, the starting time in our model) and initiation of cancer progression or death as follows: where T PFS is a random variable due to the stochastic nature of cancer progression. In our simulations, N, which represents number of patients, was set to 100. We calculated the progression-free survival time for each patient in the simulation and then computed overall survival percentages and survival frequencies for the entire patient population under different treatment schedules. In the following text, patient survival refers to progression-free survival unless stated otherwise.

Incorporation of the effects of MEK and PI3K inhibitors into the model
Currently, in addition to BRAF inhibitors (e.g., Vemurafenib and Dabrafenib), several other targeted inhibitors, such as MEK inhibitors (e.g., Trametinib and Cobimetinib) and PI3K inhibitors (e.g., BEZ235) are in clinical trials for melanoma cancer patients. In this study, we investigated the synergy between BRAF inhibitors in conjunction with each of these other two inhibitor types by incorporating their effects into our model based on their different signaling mechanisms, respectively.
MEK is a downstream effector protein of RAF signaling [45]; therefore, we assumed that MEK inhibitors would produce similar effects to BRAF inhibitors when inducing the death of drugsensitive cancer cells. It has been shown that both BRAF inhibitors and MEK inhibitors can increase the death rate of drug-sensitive cells. MEK inhibitors can also promote the secretion of drug-induced resistance factors (DIRFs) from drug-sensitive cells [15] because RAF and MEK share the same downstream effector, transcription factor FRA1, which has been identified as a major regulatory factor of DIRF secretion. Therefore, the effects of MEK inhibitor use were incorporated into the model using Hill functions as in Refs. [46,47]: where d S represents basal death rate. K BRAFi and K MEKi are Michaelis constants representing the BRAF inhibitor and MEK inhibitor concentrations at which half-maximal inhibition effects are reached. K DS1 and K DS2 are also BRAF inhibitor and MEK inhibitor Michaelis constants for DIRF secretion. C S is the number of drug-sensitive cells.
In drug-resistant cells, the PI3-Kinase (PI3K)/AKT pathway is over-activated by DIRFs [15]; therefore, a PI3K inhibitor may repress DIRF-stimulated PI3K/AKT pathway activation and thus reduce DIRF effects. We used an inhibition Hill function [47] to include the effects of PI3K inhibitor-mediated DIRF signal modification into our model: where DIRF ½ Eff represents effective action of DIRF inhibited by PI3K inhibitor, and DIRF ½ BRAFi is the DIRF steady-state concentration following stimulation with a BRAF inhibitor. K PI3Ki is the Michaelis constant for the PI3K inhibitor's half-maximal inhibition concentration. In this simplified way, we incorporated the effects of PI3K inhibitor into the model without considering complex intracellular signaling networks. This strategy enabled us to reduce the complexity of the model.
It should be noted that we used dimensionless concentrations of BRAF inhibitor, MEK inhibitor and PI3K inhibitor in the simulation by respectively normalizing them to the Michaelis constants K BRAFi , K MEKi and K PI3Ki , as in our previous study [47]. As such, we did not introduce any additional parameters into the model to further reduce the number of unknown parameters.

Parameter estimation
Biological descriptions and values of the model parameters are listed in Table 1. Most of the parameters involved in the cellular dynamics and pharmacokinetics modules were collected from previous studies [26,35,37], while other parameters were calibrated according to recent experimental [15,16] and clinical data [48]. We provide details regarding the calibration of the parameter values below.

Parameters involved in DIRF secretion
Recent studies [15] (extended data, Figure 4d in Ref. [15]) have reported that the levels of some murine stroma-derived cytokines, such as HGF, IGF, MFGE8 and TREM1, were upregulated by 1.6-to 10-fold in response to treatment with Vemurafenib (a BRAF inhibitor) compared to vehicle treatment. The newly secreted DIRFs are described in Eq. (7) in the main text. The DIRF basal level was assumed to be 0.1. To ensure that the change in total DIRF (basal plus newly secreted DIRFs) remained within a range of 1 to 10 according to the experiments mentioned above,Ṽ DS and K DS1 (i = 1, 2) were set to 100 and 10, respectively, as the concentration of BRAF inhibitor changed from 0.1 to 1 in the simulation. Note that inṼ DS ¼ V DS =d DIRF , the degradation rate (d DIRF ) is assumed to be 0.01 per day. Thus, the maximal DIRF secretion rate from drug-sensitive cells (V DS ) was accordingly set to 1.

Parameters involved in drug-resistant cell growth
The experimental data in Ref. [15] also indicated that the relative number of drug-resistant cancer cells in drug-sensitive tumors increased by approximately 2-to 4-fold after 3 days of drug treatment (Vemurafenib, Crizotinib and Erlotinib). The growth rate coefficient corresponding to drug-resistant cells regulated by DIRFs is described in Eq. (9) in the main text. Because the minimal and maximal DIRF concentrations were set to 0.1 and 0.9, as stated in the above section, α and K 2_DIRF were set to 0.1 and 0.3, respectively, to ensure an approximate 2-to 4-fold change in the relative number of drug-resistant cancer cells induced by the drug treatment on day 3, in agreement with the above-mentioned experimental data. The value of K 1_DIRF , the DIRFassociated Michaelis constant, was set in a similar way for drug-sensitive cancer cells.

Parameters involved in metastasis
The drug-resistant cancer cell dissemination rate coefficient that was regulated by DIRFs is described by Eq. (10) in the main text. We set the values of α and K 3_DIRF based on experimental data [15] ( Figure 2D in Ref. [15]). The relative migration of drug-resistant cancer cells increased by 2-to 3-fold following BRAF inhibitor treatment compared to vehicle treatment. Therefore, to reproduce this change in our simulation, α was set to 0.3, and K 3_DIRF was set to 0.01. λ (intensity of Poisson process) in Eq. (11) was scaled and estimated from Ref. [26].

Parameters fitted to the clinical data
The clinical data of progression-free survival percentages ( Figure 1A in Ref. [48]) were used to estimate parameters including d S (maximal inhibition rate of drug on drug-sensitive cells), K Tm Each parameter is listed with its symbol, value, biological description and reference. Remarks: "_" in the above Table 1 represents dimensionless unit.

(Michaelis constant for maximal carrying capacity of drug-sensitive/resistant cells), M max (maximal carrying capacity of new metastatic cells) and σ (diffusion coefficient in Wiener process).
The clinical data include progression-free survival percentages of a patient cohort treated with daily dosing of BRAF inhibitor vemurafenib and MEK inhibitor cobimetinib administered for 21 days, followed by 7 days off. By iteratively calibration, we chose the values of the above parameters which produced good agreement with the clinical data ( Figure 2).

Sensitivity analysis
Sensitivity analysis was used to quantitatively explore the critical parameters in the model that affected cellular dynamics and survival percentage. A sensitivity coefficient [49] for total tumor cell number with respect to parameter p j was calculated as follows: where t i ; i ¼ 1; ⋯L f gis an equal partition of 0; T ½ , with L = 360 and T = 360 days in the simulation. C t; p j and C t;p j represent the median total tumor cell numbers with parameter p j and varied parameterp j at time t, respectively. ΔC T t i ð Þ ¼ C T t i ;p j À C T t i ; p j is the change of median total tumor cell number with the parameter perturbation Δp j ¼p j À p j at time t i (day). Then, we obtained the relative changes for the area under curve of total tumor cell number with respect to the examined parameters.
The sensitivity coefficient of the survival percentage with respect to parameter p j was also calculated according to the following formula: where SP t; p j and SP t;p j represent survival percentages with parameter p j and varied parameterp j at time t respectively. ΔSP t i ð Þ ¼ SP t i ;p j À SP t i ; p j is the change of survival percentage with the parameter perturbation Δp j ¼p j À p j at time t i (day). Then, we obtained the relative changes for the area under curve of survival percentage with respect to the examined parameters.

Single-parameter sensitivity analysis
Each parameter was increased by 50% from its estimated value and we then obtained the relative changes in area under curve of total tumor cell number and area under survival curve, as defined by Eqs. (16) and (17). The computations were repeated 20 times, and the mean value and standard deviation of the sensitivity coefficients were calculated.

Two-parameter sensitivity analysis
To assess the combinatorial effects of parameters involved in the model, we performed a twoparameter sensitivity analysis. The values of each pair of two different parameters were increased by 50% from their original values simultaneously. All other parameters remained at their base values. Taking into account the stochasticity of the model, the computations were repeated 20 times, and the mean value of the sensitivity coefficients was calculated.

In silico prediction of cellular kinetics and patient survival following drug treatment
We investigated cellular response kinetics following drug treatment in silico. In Figure 3, a typical simulation of BRAF-I treatment is shown. Figure 3A details BRAF-I kinetics for 100 cancer patients. The inhibitor was administered daily for 3 weeks followed by 1 week of no treatment, concordant with the drug schedule of a previous study [48]. Figure 3B and 2E show the time courses of all 100 samples with respect to numbers of drug-sensitive cancer cells, drug-resistant cancer cells, metastatic cells and total tumor cells. Drug-sensitive cancer cell growth was repressed following drug administration, but it periodically rebounded during no treatment weeks. Interestingly, metastatic cell growth ( Figure 3D) showed a similar pattern among the patients: an initial slow growth period (the length of which varied in different patients) followed by a rapid increase within~1 month. The metastatic cell populations in different patients exhibited different transition times, resulting in heterogeneous sizes of cancer cells among patient population ( Figure 3E). Interestingly, this "all-or-no" metastasis causes a bimodal distribution for the number of total tumor cells after 12 months ( Figure 3F), indicating that cancer in some patients progressed but not yet in others.
For comparison, cellular dynamics without drug treatment are shown in Figure 4. In this case, the distribution ( Figure 4D) of the number of total tumor cells was changed to an asymmetric bimodal distribution with decreased frequencies of cell numbers at high level ( Figure 4E) due to the lack of the drug-induced metastasis ( Figure 4C). We also applied our model to investigate the effects of targeted therapy on patient survival at the population level. Figure 4F shows the survival percentage of 100 cancer patients undergoing BRAF-I treatment from 0 to 360 days compared with no treatment controls. Our simulation demonstrated that treatment with BRAF-I significantly prolonged progression-free survival in melanoma patients harboring BRAF mutations. This result is consistent with clinical studies of melanoma patients harboring the BRAF V600E mutation [48,50].

Validation of cellular kinetics and patient survival using clinical data
We next validated our model using clinical data. In Figure 5, a comparison of our model predictions with clinical population-scale survival data is shown [51]. The clinical data included distributions of progression-free survival times for 54 patients in each group treated either with Dabrafenib monotherapy (a BRAF-I) ( Figure 5A) or with a combination of BRAF-I and Trametinib (a MEK inhibitor, or MEK-I) at doses of 150 and 1 mg (1X; Figure 5B) or 150 and 2 mg (2X; Figure 5C). To simulate the conditions used to produce the clinical data, our model (red lines in Figure 5A-C) examined the following three treatment strategies: BRAF-I alone, BRAF-I combined with 1 mg MEK-I, and BRAF-I combined with 2 mg MEK-I. It should be noted that the drug doses in the simulation have been normalized (refer to the Model section). We computed progression-free survival times using our model and compared them against the clinical data; the predicted and experimental results were in good agreement. Furthermore, as shown in Figure 5D, our model predicted that the combination therapies enhanced progression-free survival more than the monotherapy, consistent with the clinical data. Figure 6 shows a validation of the sudden increase observed in metastatic cell number in the model using clinical values of circulating tumor DNA (ctDNA). ctDNA has been proposed as a promising biomarker for monitoring metastatic cancers [52]. We used clinical data consisting of plasma ctDNA concentrations from nine patients [53] treated with BRAF-I and MEK-I in combination to verify the predicted metastatic cell growth. Of the nine evaluated patients, four showed elevated ctDNA levels while undergoing combination therapy. In Figure 6A, a comparison of the simulated growth pattern of metastatic cells with the clinical ctDNA data is shown. Both the clinical data and the model prediction showed a pattern of explosive metastatic cell growth in several patients. In addition, the model predicted that new metastatic cell numbers would either remain at an undetectable, low level (48%) or significantly increase (52%) ( Figure 6B). A threshold of metastatic cell number was set to 1, separating these two distinct levels. This prediction agrees with the clinical plasma ctDNA concentrations, in which ctDNA levels are either undetectable (5/9) or elevated (4/9).

Evaluation of drug combination synergy
Currently, in addition to BRAF inhibitors (e.g., Vemurafenib and Dabrafenib), several other targeted inhibitors, including MEK inhibitors (e.g., Trametinib and Cobimetinib) and PI3K inhibitors (e.g., BEZ235), are being evaluated in clinical trials for treatment of melanoma patients. In the following study, we investigated whether the co-administration of BRAF-I with either MEK-I or PI3K inhibitor (PI3K-I) produces synergistic effects. To accomplish this, we The predicted progression-free survival percentages showed improved survival over time following the administration of combination therapeutics compared to BRAF-I monotherapy. The predicted survival curve shown here has a highly similar pattern to that in Figure 1A of Ref. [51].
incorporated the effects of these inhibitors into our model (see details in Model section) based on their specific signaling mechanisms.
A Bliss combination index [54,55] was used to quantitatively evaluate the synergy produced by BRAF-I and MEK-I co-treatment and BRAF-I and PI3K-I co-treatment as follows: where R i x i ð Þ ¼ C T À C T x i ð Þ À Á =C T is the reduction of the median total tumor cell number (C T ) induced by the single BRAF inhibitor (i ¼ 1), or either MEK inhibitor or PI3K inhibitor (i ¼ 2) alone at dose Importantly, the predicted differences in the synergies of these two combinations are consistent with experimental studies [56,57], which have reported stronger synergy for BRAF-I & PI3K-I co-treatment than BRAF-I & MEK-I co-treatment. Our model also predicted that the maximal CI value for the BRAF-I & PI3K-I combination (up to 0.3, Figure 7A) is greater than that of the BRAF-I & MEK-I combination (less than 0.15, Figure 7B). This demonstrates good agreement between our model's predictions and the experimental data. In addition, the dose-dependent synergy predicted by our model might also help explain contradictory experimental observations that different PI3K/AKT inhibitors (e.g., MK2206 and Perifosine) produce opposing effects when combined with BRAF-I (PLX4032) [58].

Sensitivity analysis
We conducted parameter sensitivity analysis to examine whether the model is robust to parameter variations and to identify parameters with critical effects on both cellular dynamics and patient survival. Figure 8A, B show a single-parameter sensitivity analysis. Here, a parameter was regarded to be sensitive/critical if its sensitivity coefficient is greater than 0.2. The results showed that the model was rather robust with respect to the variations of most parameters including those calibrated. Moreover, the following parameters were critical to model outputs: growth rate and death rate in angiogenic cells (r K , d K ) and metastatic rate (λ). Since angiogenesis positively supports the sustained growth of the tumor cells and provides an avenue for dissemination and translocation of metastatic cancer cells (especially drugresistant metastatic cells), therefore the growth/death rates in angiogenic cells as well as metastatic rate of cancer cells were shown to significantly affect cellular dynamics and patient survival.
We further conducted a two-parameter sensitivity analysis to investigate how the parameters combine with each other to affect cellular dynamics and patient survival. The values of each pair of two different parameters were increased by 50% from their original values simultaneously. The computations were repeated 20 times, and the mean value of the sensitivity coefficients was calculated. Figure 8C, D show the relative changes of the areas under the curves of the total tumor cell number ( Figure 8C) and patient survival percentage ( Figure 8D) with respect to the combinatorial variations in parameter values. This global sensitivity analysis result revealed some interesting phenomena. For example, the combination of r S (growth rate of sensitive cells) and r K , the combination of d S (drug-induced death rate of sensitive cells) and d K , the combination of r R (growth rate of resistant cells) and r K , and the combination of r K and λ show significant impacts on the tumor growth and patient survival. This also suggests that combining anti-angiogenic therapy with targeted therapy to combat drug resistance and cancer progression may improve cancer patient survival.

Discussion
In this study, to examine therapy-induced drug resistance and cancer metastasis, we designed a novel stochastic model that connects the biological mechanisms underlying cancer drug resistance to population-level patient survival. A set of stochastic differential equations (SDEs) was used to model the dynamics of heterogeneous cellular populations containing drug-sensitive, drug-resistant, and metastatic cancer cells. An associated random variable characterizing progression-free survival time was subsequently defined. Our approach revealed several interesting features associated with cancer metastasis and progression kinetics. For example, dosedependent synergy was evident in both BRAF-I and MEK-I co-treatment and BRAF-I and PI3K-I co-treatment. This result suggests that combination therapy with optimized dosages of inhibitors might reduce drug resistance.
Our model demonstrated that cancer metastasis and progression occur in bursts ( Figure 3D, F and Figure 6A). As such, metastasis may occur earlier than can be detected using common radiographic imaging approaches [59]. Furthermore, metastatic cancers may quickly progress after being detected. Based on these phenomena, new prognostic methods that offer much earlier detection of metastasis and progression are needed. The identification of biomarkers that can be easily and regularly measured in cancer patients to detect cancer metastasis would be promising. Recently, ctDNA [60] has been acknowledged as a promising biomarker for several types of cancer [52,[61][62][63][64]. Using new PCR technologies, ctDNA can be quantitatively measured in cancer patients [65,66]. Dynamic changes in ctDNA levels might serve as a biomarker of early relapse in cases of surgically resected metastatic melanoma.
We also investigated the impact of heterogeneity in patient drug-metabolism rates on survival percentage. To include the effects of population heterogeneity in our mechanistic models and examine the impact of patient drug metabolism heterogeneity on survival percentage, we considered three patient subclasses with different metabolic rates (low, medium, and high). As shown in Figure 9, metabolic rate differences significantly affected the probability of patient survival. The patient subclass with a high rate of drug metabolism showed a lower survival probability compared to the patient subclass with a low rate of drug metabolism. This was true under both the 3 weeks on/1 week off treatment schedule ( Figure 9A) and the daily treatment schedule ( Figure 9C). In Figure 9B, D, the effects of different drug combination dosages on survival percentage in high-metabolism patients are shown. An appropriately elevated dosage of BRAF-I combined with MEK-I (2-fold of normal dosages) improved progression-free survival in cancer patients with high metabolic rates. As a strategy for personalized therapy, optimizing dosing based on patient metabolic rate might be beneficial. Furthermore, our model indicated that daily treatment of cancer patients ( Figure 9C, D) resulted in higher survival rates versus the 3 weeks on/1 week off treatment schedule ( Figure 9A, B). These results suggest that daily drug treatment might produce better results than therapies with drug discontinuance.
Our model predicted that combined BRAF and MEK inhibitions produced a different pattern of synergy from that created by combined BRAF and PI3K inhibition. As demonstrated in [15], the use of a BRAF-I can induce drug-sensitive cells to secrete DIRFs that promote PI3K/AKT/ mTOR pathway activation in drug-resistant cells. Thus, BRAF-I and PI3K-I co-treatment represses the actions of DIRFs on drug-resistant cells and reduces growth, dissemination and metastasis in drug-resistant cells. Therefore, at appropriate dosages, this combination therapy produced a synergistic effect, reducing the number of tumor cells ( Figure 5C, D). Conversely, because MEK is a downstream target of BRAF signaling, MEK inhibition does not necessarily strengthen BRAF inhibition. As such, BRAF-I and MEK-I co-treatment produced only weak synergistic effects at relatively low doses. This lack of cooperation led to a smaller reduction in tumor cell number ( Figure 5A, B).
Drug resistance is an obstacle often encountered during oncoprotein-targeted therapy and develops by very complex mechanisms. Tumors frequently display a substantial amount of spatial heterogeneity in both cell population composition and microenvironment factors, such as drug, oxygen and glucose concentrations [2,67]. Recent studies have suggested that the emergence of drug resistance is driven by the tumor microenvironment [2]. In our study, based on recent clinical and preclinical data, we modeled feedback in drug-sensitive cancer cells in response to drug treatment and interactions between heterogeneous cell populations and their microenvironments to understand drug resistance and metastasis dynamics in tumors. In future work, we aim to investigate spatial heterogeneity in tumor cells by developing an agent-based model [21,43] using partial differential equations to assess spatial-temporal changes in the tumor microenvironment.
In summary, our study utilized a set of SDEs to model the dynamics of targeted therapyinduced drug resistance and metastasis. The model predicted that different patterns of synergy exist for the combination of BRAF-I and MEK-I and for the combination of BRAF-I and PI3K-I. Our study provides insight into the microenvironment-mediated mechanisms underlying drug resistance and delineates the implications associated with optimal dosing and scheduling of combination therapy for melanoma patients with BRAF mutations. It is our hope that this predictive model will facilitate cancer clinical trial design and accelerate the design of effective and robust tumor therapeutics.