Mathematical Modeling of Therapy-induced Cancer Drug Resistance: Connecting Cancer 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 dose-dependent 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 proposed model might facilitate the quantitative evaluation and optimization of combination therapeutics and cancer clinical trial design.

Text S2. Sensitivity analysis methods.    Table S1. Values of parameters in the model.

Text S1. Parameter calibration
Biological descriptions and values of the model parameters are listed in Table S1. Most of the parameters involved in the cellular dynamics and pharmacokinetics modules were collected from previous studies [1][2][3], while other parameters were calibrated according to recent experimental [4,5] and clinical data [6]. We provide details regarding the calibration of the parameter values below.
Parameters involved in DIRF secretion. Recent studies [4] (extended data, Figure 4d in Ref. [4]) 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 equation (7) in the main text.
The DIRF basal level was assumed to be 0. 1 Parameters involved in drug-resistant cell growth. The experimental data in Ref. [4] 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 equation (9)  Parameters involved in metastasis. The drug-resistant cancer cell dissemination rate coefficient that was regulated by DIRFs is described by equation (10) in the main text. We set the values of  and 3_DIRF K based on experimental data [4] ( Figure 2D in Ref. [4]). 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 3_DIRF K was set to 0.01.
Parameters fitted to the clinical data. The clinical data of progression-free survival percentages ( Figure   1A in Ref. [6]) were used to estimate parameters including d S (maximal inhibition rate of drug on drugsensitive cells), K Tm (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 (see Supplementary Fig. S1). It should be noted that another independent set of clinical data ( Figure 1A in Ref. [7]) was used for validating our model predictions at patient population level, as presented in Fig. 3 in the main text.

Text S2. 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 [8] for total tumor cell number with respect to parameter j p was calculated as follows: T , with L =360 and T =360 days in the simulation.  The sensitivity coefficient of the survival percentage with respect to parameter j p was also calculated according to the following formula:

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 equations (S1-S2). 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 two-     Remarks: "_" in the above Table S1 represents dimensionless unit.