Predicting response to combination evofosfamide and immunotherapy under hypoxic conditions in murine models of colon cancer

1.   Introduction

The aberrant growth of tumor cells can lead to an oxygen demand that is higher than the existing vasculature can provide ^[1]. This lack of sufficient oxygenation and nutrients in tumor cells may result in regions of hypoxia (defined as $\leq 10$ mm Hg ^[2,3]. When compared to normoxic tumors (i.e., $> 10$ mm Hg), hypoxic tumors have been linked with a poorer outcome for several anticancer therapies, including radiotherapy, chemotherapy, and immunotherapy ^[4,5,6,7]. However, it has been shown that resistance to immunotherapy can be alleviated by combining immunotherapy with hypoxia-activated prodrugs and altering the oxygenation within a tumor ^[8,9].

Hypoxia-activated prodrugs are activated only under low oxygen tension, resulting in a hypoxia-specific release of the active drug in areas of low oxygen, such as solid tumors, while sparing normal tissue that has adequate oxygen levels ^[10]. For example, under hypoxic conditions, evofosfamide (TH-302) is converted into the alkylating agent bromo-isophosphoramide mustard ^[10]. When combined with radiotherapy and immunotherapy, evofosfamide can improve treatment outcome ^[8,9,10,11]. Several clinical trials and in vivo preclinical experiments have investigated the effects of evofosfamide in combination with a range of treatment types in cancer with favorable results ^{[12,13,14,15]}. For example, Laubach et al. designed a phase Ⅰ/Ⅱ study for multiple myeloma, to investigate the effects of combination of evofosfamide and dexamethasone where disease was observed to be stabilized in over 80% of the patients ^[13]. Hegde et al. evaluated the tumor response to a combination of evofosfamide and immune checkpoint inhibitor, ipilimunab ^[9]. In this study, of 18 patients (with a variety of primary tumors), three achieved partial response and 12 achieved stable disease ^[9]. Compared to nonresponders, the responders exhibited improved peripheral T-cell proliferation and increased intratumoral T-cell infiltration into hypoxia regions ^[9]. In preclinical models, evofosfamide was shown to reduce or eliminate hypoxia in prostate cancer and the combination of evofosfamide and immunotherapy cured 80% of mice. Following treatment, mice treated with evofosfamide (alone or in combination with immunotherapy) demonstrated that less than 7% of the tumor area was hypoxic compared to 33% in control mice ^[15]. These studies indicate that it may be useful to characterize the oxygenation status of a tumor to determine if evofosfamide should be administered in combination with anti-cancer therapy.

There are several methods for noninvasively quantitatively interrogating tumor oxygen status via positron emission tomography (PET; ^[16]), such as [¹⁸F]FMISO-PET and FAZA-PET. For example, Reeves et al. quantified tumor hypoxia using [¹⁸F]-fluoromisonidazole ([¹⁸F]FMISO)-PET before and after immunotherapy ^[8]. The data was collected in experiments with mice injected with either MC38 (murine colorectal cancer cells) or E0771 (murine breast cancer cells). The authors demonstrated that within five days of starting immunotherapy, nonresponding colorectal and breast tumors exhibited a significant increase in hypoxia. After identifying the nonresponding tumors, the investigators treated them with a combination of immunotherapy and evofosfamide, resulting in significant improvement in tumor volume control and overall survival. This combination of evofosfamide and immunotherapy led to a significant reduction in hypoxia within hypoxic tumors. Further, Kumar et al. ^[17] have demonstrated that evofosfamide has a synergistic effect when combined with anti-VEGF inhibitors, indicating it's potential to enhance tumor response to immunotherapy. Given these promising results on identifying hypoxia as a potential imaging biomarker for immunotherapy response and utilizing that information to enhance the efficacy of subsequent treatment with immunotherapy, it is natural to work towards optimizing this process.

Mathematical models have been developed to study the combination of evofosfamide and combretastin in glioblastoma ^[18], evofosfamide and erlotinib in non-small cell lung cancer ^[19], and evofosfamide and radiotherapy in chondrosarcoma ^[20]. For example, Meaney et al. developed a system of reaction diffusion equations to model a glioblastoma treated with evofosfamide and combretastin ^[18]. They simulated patients receiving just one of the drugs or their combination. The authors also explored scenario where the drugs were released from the vessels or from nanocells. The parameters of the model were taken from previous mathematical models that were calibrated to experimental data, and their results indicate that the nanocells have the potential to increase the delivery and activation of both drugs. Lindsay et al. developed a stochastic mathematical model, parametrized using clinical and experimental data (e.g., experimental growth rates, cell response to evofosfamide, tumor oxygenation measurements, and pharmacokinetic data), and optimized the treatment protocol of evofosfamide and erlotinib ^[19]. They found that the optimal treatment protocol to reduce the tumor burden is to sequentially alternate a single dose of each drug while minimizing the time between the evofosfamide and erlotinib doses. Hamis et al. developed a hybrid, multiscale cellular automata model to study solid tumors subjected to evofosfamide and radiotherapy ^[20]. The model was parameterized by published data from an in vitro experiment with multicellular tumor spheroids of human chondrosarcoma. Numerical simulations demonstrated that evofosfamide increased the radiotherapy efficacy in hypoxic tumors. Despite these advances on modeling the combination of evofosfamide with other therapies, the combination of evofosfamide and immunotherapy is currently unexplored by mathematical modeling.

In this contribution, motivated by the results obtained in ^[8], we develop a mathematical model to reproduce the growth of normoxic and hypoxic tumors, as defined by [¹⁸F]FMISO-PET imaging data, under different treatment regiments of immunotherapy and evofosfamide ^[8]. The experimental data collected by Reeves et al. is used to calibrate the model. We employed a staggered Bayesian calibration method to isolate the model parameters within each calibration scenario. This approach is designed to systematically explore the main mechanisms of our model that allow it to fit every experimental scenario. The calibration is completed excluding the data from normoxic tumor cells treated with immunotherapy and evofosfamide, which will be used as a validation scenario. Thus, we develop a framework capable of identifying the effects of combining immunotherapy and evofosfamide, and their roles in reproducing the experimental data.

2.   Methods

2.1.   Murine model

As details are described elsewhere ^[8], here we present only the salient details on how we obtained tumor size changes as function of time and treatment in a murine model of colorectal cancer. MC38 murine colorectal cancer cells were obtained from Kerafast (September 2019). These cells are an immunogenic, grade Ⅲ adenocarcinoma of colorectal cancer and are characterized by microsatellite instability. They were cultured in Dulbecco's Modified Eagle Medium (DMEM) with 10% fetal bovine serum (FBS), 2 mmol/L L-glutamine, and 1 mmol/L sodium pyruvate in a humidified incubator with 5% CO2 at 37C. On day 0, $5\times 10^5$ cells were diluted in 40% Matrigel and 60% serum free DMEM and subcutaneously injected into the upper right shoulder of 6- to 12-week-old C57BL/6 mice (Charles River Laboratories). Tumors were grown until they reached a volume of a 100 mm$^3$ (7-10 days post inoculation).

In Figure 1, we present the experimental protocol and indicate the days that: 1) [¹⁸F]FMISO-PET data were collected, 2) mice received 200 $\mu$g anti-PD-1 (clone RPM1-14, Bio X Cell) + 100 $\mu$g anti-CTLA-4 (clone 9H10, Bio X Cell) combination therapy, 3) mice were treated with 50 mg/kg evofosfamide (TH-302 Selleckchem), and 4) the tumor volume was measured. Immunotherapy and evofosfamide doses were chosen based on concentrations that reflected human clinical studies. To separate mice into normoxic or hypoxic groups, the standard uptake value (SUV) of [¹⁸F]FMISO in the tumor and muscle were quantified and their ratio was used to stratify tumors as normoxic or hypoxic (SUV$>$1.86, sensitivity 100%) as previously defined ^[8]. Thus, mice were divided into six experimental groups (separated based on imaging metrics): Ⅰ) control, Ⅱ) hypoxic mice that received immunotherapy, Ⅲ) normoxic mice that received immunotherapy, Ⅳ) evofosfamide, Ⅴ) hypoxic mice that received immunotherapy and evofosfamide, and Ⅵ) normoxic mice that received immunotherapy and evofosfamide.

Figure 1.  C57BL/6 mice were subcutaneously implanted with MC38 cancer cells. Starting at day 6, the tumor volume was measured with calipers (days marked in magenta). [¹⁸F]FMISO-PET data was collected at days 12 and 17 (green arrows) to determine percent of the tumor that is hypoxic during treatment. In the days indicated by the yellow arrows, mice received immunotherapy. Mice that received evofosfamide were treated between days 13 and 17 (red arrows). The mice were divided in the six scenarios according to the treatments received and fraction of hypoxic voxels.

DownLoad: Full-Size Img PowerPoint

2.2.   Mathematical models

We developed a general mathematical model that characterizes the temporal dynamics of the tumor volume by directly accounting for the effects of anti-hypoxia prodrugs, evofosfamide, and checkpoint blockade immunotherapy, anti-PD1 and anti-CTLA4, in different scenarios. We propose a system of ordinary differential equations that captures the temporal evolution of tumor volume, $T(t)$, and the well-vascularized tumor fraction, $V(t)$. The well-vascularized fraction is defined as $V(t) = 1 - H(t)$, where $H(t)$ is the hypoxic tumor fraction. We assume that regions with high hypoxia are not vascularized. While hypoxia imaging furnishes valuable insights into a tumor's oxygenation status ^[8], our assumption may not fully capture the complexity of vascularization and other factors influencing hypoxia. Nevertheless, given these considerations, we have opted to utilize hypoxia imaging as a proxy for vascularization and have incorporated it into our model. We also assume that the tumor volume increases exponentially at the rate, $k_T$. This proliferation function is selected as the tumor is able to grow without any apparent limiting factor during the time frame from our experiments (e.g., without appearing to reach a maximum tumor volume).

The tumor volume decreases as the immune system eliminates the tumor cells, which we assume occurs more frequently in well-vascularized regions. The tumor has a constant death rate, $\mu_T$, but the tumor death rate can increase by $\gamma_i$ each time, $t_i$, the immunotherapy is delivered. The well-vascularized tumor fraction increases at the rate, $k_V$. This vascularization rate can then increase by a constant, $\gamma_e$, each time, $t_e$, the mouse receives evofosfamide. As the tumor increases, the hypoxic fraction also increases, and thus, the well-vascularized region decreases at a rate $\mu_V$. The pharmacodynamics and pharmacokinetics of the combination of evofosfamide and immunotherapy are defined through $\gamma_i$, $\gamma_e$, $r_i$, and $r_e$. The tumor death rate increase, $\gamma_i$, and vascularization rate increase, $\gamma_e$, describe the pharmacodynamics of the treatment, as they capture the impact of immunotherapy and evofosfamide on tumor cell death and tumor vascularization, respectively. Conversely, the decay rates of the effects of immunotherapy, $r_i$, and evofosfamide, $r_e$, are vital for describing both the pharmacokinetics of the drugs (and their wash out) in conjunction with the pharmacodynamics with respect to drug efficiency, as they govern the time-dependent decline in the efficacy of both treatments. With these assumptions, the model takes the following form:

where $r_i$ and $r_e$ are the decay rates due to checkpoint blockade anti-PD1+anti-CTLA4 immunotherapy and evofosfamide, respectively, and $N_i$ and $N_e$ are the number of days that the mice received immunotherapy and evofosfamide, respectively. In Eqs 2.1 and 2.2, $\mathcal{H}$ is the Heaviside function. (See Table 1 for a description of all model parameters.)

Table 1.  Parameter's definitions and uniform priors. $\mathcal{U}(a, b)$ represents a uniform distribution with bounds between $a$ and $b$ for the prior values of the corresponding parameters.

Parameter	Meaning	Prior
$k_T$	tumor growth rate	$\mathcal{U}(0, 0.4) h^{-1}$
$k_V$	vascularization rate	$\mathcal{U}(0, 0.5) h^{-1}$
$\mu_V$	well-vascularized tumor fraction decay rate	$\mathcal{U}(0, 0.5) h^{-1}$
$\mu_T$	tumor death rate by immune system	$\mathcal{U}(0, 0.5) h^{-1}$
$\gamma_i$	immunotherapy tumor death rate increase	$\mathcal{U}(0, 3.5) h^{-1}$
$r_i$	immunotherapy effect decay rate	$\mathcal{U}(0, 0.2) h^{-1}$
$\gamma_e$	evofosfamide vascularization rate increase	$\mathcal{U}(0, 3.5) h^{-1}$
$r_e$	evofosfamide effect decay rate	$\mathcal{U}(0, 0.2) h^{-1}$

---

 | Show Table

DownLoad: CSV

2.3.   Bayesian model calibration

The parameters within Eqs 2.1 and 2.2 are calibrated using the experimental data ^[8]. To account for the data uncertainties and model inadequacy, the models are calibrated via a Bayesian framework defined as:

where $\boldsymbol{D}$ is the experimental data, $\boldsymbol{\theta}$ is the vector of model parameters to be calibrated, $\pi(\boldsymbol{\theta})$ is the prior knowledge about the model parameters (see Table 1 for the prior used for each parameter), $\pi(\boldsymbol{D}|\boldsymbol{\theta})$ is the likelihood that the data is observed for a given set of parameters, $\pi(\boldsymbol{D})$ is a normalizing factor, and $\pi(\boldsymbol{\theta}|\boldsymbol{D})$ is the posterior distribution of the parameters. Assuming that the data is normally distributed, the log-likelihood is:

where $N_t$ and $N_v$ are the number of measurements from the experimental data of the tumor and the well-vascularized region, respectively, $\boldsymbol{Y}$ is the output of the model (i.e., $\boldsymbol{Y}^T$ is the tumor volume over time, and $\boldsymbol{Y}^V$ is the well-vascularized fraction over time), and $\sigma_T$ and $\sigma_V$ are the variance of the total error (i.e., the sum of the variance of the model uncertainties and the model inadequacy) of the tumor and the well-vascularized region, respectively. In our model, we have two quantities of interest, the tumor itself and the well-vascularized region within the tumor. As our goal is to calibrate the model to accurately reproduce the whole tumor volume data, and since we only have two time points to calibrate the well-vascularized region, we assign a higher weight to the tumor ($\alpha$) than to the well-vascularized region ($\beta$) in the calibration process. Thus, we take $\alpha = 2$ and $\beta = 1$ for the weights in Eq 2.5.

In Figure 2, we present the model calibration framework. We start by defining the complete tumor growth model (i.e., Eqs 2.1 and 2.2), which characterizes the effects of evofosfamide and the immune system. The parameters $k_T$, $k_V$, $\mu_V$, $\mu_T$, $\gamma_i$, $r_i$, and the initial well-vascularized fractions are calibrated to scenarios Ⅰ, Ⅱ, and Ⅲ (note that each scenario has its own initial well-vascularized fraction). Thus, we can calibrate just the effects of the immunotherapy and the tumor growth without the interaction with evofosfamide. Next, we calibrate the parameters $\gamma_e$ and $r_e$, and the initial well-vascularized fractions, to scenarios Ⅳ and Ⅴ while using the values obtained from the previous scenarios to assign the other parameters. After this step, all model parameters are calibrated and we can test the ability of the model to predict the tumor volume over time in scenario Ⅵ. If the prediction error (as quantified by, for example, the mean percent error, the Pearson correlation coefficient (PCC), or the concordance correlation coefficient (CCC)) is below a defined threshold, we can conclude that our model is able to accurately predict scenario Ⅵ with the parameter values obtained from scenarios Ⅰ - Ⅴ. This would mean that only knowledge of the evofosfamide effects on hypoxic cells is sufficient to accurately inform the temporal changes in normoxic tumor cells treated with immunotherapy and evofosfamide combination. However, if we find that the model fails to accurately predict the tumor growth (i.e, the prediction error is above the threshold), we find the smallest number of parameters to re-calibrate to capture the tumor growth observed in scenario Ⅵ. If we are able to find this set of parameters, it means that there are further differences between normoxic and hypoxic tumors that are not captured by the model when using the same parameter values obtained from scenario Ⅰ - Ⅴ. If we are not able to identify a parameter set to capture the tumor growth in scenario Ⅵ, the underlying mathematical model described by Eqs 2.1 and 2.2 need to be changed.

Figure 2.  Our model calibration framework begins with first defining an initial model. Next, we calibrate the model parameters to the data from scenarios Ⅰ, Ⅱ, and Ⅲ. The parameters regarding the effects of evofosfamide are excluded from this first set of calibrations, but are determined during calibration to the data from scenarios Ⅳ and Ⅴ. After all the parameters are calibrated, we predict the tumor dynamics in scenario Ⅵ and compare the prediction to the data. If the error is below a threshold, we stop the process because we have arrived at a suitable parameter set to describe all scenarios. If the error is above the threshold, we proceed to find which parameter(s) needs to be re-calibrated. If we can reproduce scenario Ⅵ with the new parameter value, we then stop. However, if we are not able to capture the tumor dynamics for all scenarios, the mathematical model of Eqs 2.1 and 2.2 would need to be adjusted.

DownLoad: Full-Size Img PowerPoint

2.4.   Sensitivity analysis

In this work, we utilize the Sobol method, which is a variance-based global sensitivity analysis method ^[21,22,23]. One notable aspect of this method is its ability to assess the contributions of individual parameters, as well as capture higher-order effects arising from parameter interactions. We now summarize the application of the Sobol method to our problem.

Let $\vec{M}(\vec{\theta})$ be a model parameterized by parameters $\vec{\theta}$ which belong to a parameter space $\vec{\Theta}\subset\mathbb{R}^k$. Among all possible methods for computing the total sensitivity index (which quantifies all effects of the parameter on the model output), we have selected the techniques presented in ^[24] due to their efficient convergence with smaller sample sizes compared to other approaches. Initially, we construct two $N\times K$ matrices, ${\bf A}$ and ${\bf B}$, where $K$ is the number of parameters, and $N$ is the sample size. In these matrices, $N$ random samples are drawn from a uniform distribution corresponding to the range of each parameter's uncertainty. Each row of these matrices represents a sampled value for $\vec{\theta}$. Additionally, $K$ matrices ${\bf A}_{{\bf B}}^{(k)}, \; k = 1, 2, \cdots, K$, are defined where all columns are from ${\bf A}$ except the $k^{th}$ column, which comes from ${\bf B}$. The model is then evaluated for each row of the matrices ${\bf A}$ and ${\bf A}_{{\bf B}}^{(k)}$ and the outputs are stored in the vectors $\boldsymbol{Y}_{\boldsymbol{A}}$ and $\boldsymbol{Y}_{\boldsymbol{A}\boldsymbol{B}}^{(k)}$, respectively. Finally, post-processing these results consists of computing the total sensitivity index, $S_{T_k}$, for each parameter $k$ ^[25]. The total effect, $S_{T_k}$, measures the contribution of the input $\theta_k$ to the model output variation. The total-effect index for each parameter, $\{S_{T_k}\}_{k = 1}^K$, can be approximated using the following estimator ^[24],

By employing the algorithm outlined in the previous paragraph, we can estimate multi-dimensional integrals with just $N(K+1)$ model evaluations ^[24]. For time-dependent processes, we repeat these steps for each time step to capture the importance of each model parameter over time. In the results section, we present the time evolution of $S_{T_k}$ for multiple treatment protocols, including control, immunotherapy, evofosfamide, and immunotherapy plus evofosfamide. The order of parameter importance can be used to define the order in which to change parameters, as shown in the "find the parameters to change" step presented in Figure 2. Specifically, we re-calibrate the parameter that has the highest ranking and keep the other parameters fixed. This process helps to identify the minimum number of parameters that need to be re-calibrated to accurately capture scenario Ⅵ.

2.5.   Numerical implementation

The model given by Eqs 2.1 and 2.2 and the calibration framework (presented in Figure 2) are implemented in C++. The model is solved via a fourth order Runge-Kutta method ^[26]. We employ a parallel, adaptive, multilevel Markov Chain Monte Carlo (MCMC) sampling method to compute the posterior density $\pi(\boldsymbol{\theta}|\boldsymbol{D})$ (see, e.g., ^[27,28,29]). The MCMC method is available in the C++ library QUESO (Quantification of Uncertainty for Estimation, Simulation, and Optimization) ^[30]. The code itself, as a well as a description of how to use it, is provided at https://github.com/Ernesto-Lima/CEIRPrediction.

3.   Results

Using the priors defined in Table 1 and following the framework presented in Figure 2, we calibrate the parameters $k_T$, $k_V$, $\mu_V$, $\mu_T$, $\gamma_i$, $r_i$, and the initial well-vascularized fractions to the data from scenarios Ⅰ (control), Ⅱ (hypoxic tumors that received immunotherapy), and Ⅲ (normoxic tumors that received immunotherapy). In Figure 3, the experimental data and the calibrated model are presented for these scenarios. The mean absolute percent error for the tumor volume is below 14% and the CCC is above 0.8 in all three calibrated scenarios. Although the calibration is weighted towards the tumor volume, the mean absolute percent error for the well-vascularized fraction is below 23% in all three scenarios. It is worth noting that, while the initial well-vascularized fraction is similar in scenarios Ⅰ and Ⅱ (panels B and D, respectively), the initial vascularization is higher in scenario Ⅲ (panel F), which is the normoxic scenario.

Figure 3.  Temporal dynamics of the experimental data (black) and the model (blue) in scenarios Ⅰ (control, panels A and B), Ⅱ (hypoxic tumors that received immunotherapy, panels C and D), and Ⅲ (normoxic tumors that received immunotherapy, panels E and F). Left column represents the tumor volume, and right column the well-vascularized fraction. The dashed lines in blue indicate the bounds of the 95% confidence interval of the model output. The mean and 95% confidence interval from the data is also indicated in black. The mean absolute percent error for the tumor volume is below 14% and the CCC is above 0.8 in the three scenarios calibrated. Even though the calibration is weighted towards the tumor volume, the mean absolute percent error for the well-vascularized fraction is below 23% in the three scenarios.

DownLoad: Full-Size Img PowerPoint

In the next step of our framework, we calibrate the parameters $\gamma_e$ and $r_e$, and the initial well-vascularized fractions, to the data from scenarios Ⅳ (evofosfamide alone) and Ⅴ (hypoxic tumors that received combination immunotherapy and evofosfamide). The other six parameters are fixed to the values found during the calibration of scenarios Ⅰ, Ⅱ, and Ⅲ. At this step of our framework, we are able to calibrate every parameter from our model (i.e., Eqs 2.1 and 2.2). In Table 2, we present the mean and the 16th and 84th percentiles (i.e., the standard deviation of a normal distribution) of the calibrated parameters. In Figure 4, the experimental data and the calibrated model are presented for scenarios Ⅳ (evofosfamide) and Ⅴ (hypoxic mice that received immunotherapy and evofosfamide). The mean absolute percent error for the tumor volume is below 28% and the CCC is above 0.8 in the two scenarios calibrated. These results indicated that we are able to capture the tumor dynamics of these two scenarios using the parameter values obtained in scenarios Ⅰ, Ⅱ, and Ⅲ.

Table 2.  Mean, and 16th and 84th percentiles (i.e., the standard deviation for a normal distribution) of the calibrated parameters.

Scenario	Parameter
	$\gamma_e$ (h$^{-1}$)	$r_e$ (h$^{-1}$)	$\gamma_i$ (h$^{-1}$)	$k_T$ (h$^{-1}$)	$k_V$ (h$^{-1}$)	$\mu_V$ (h$^{-1}$)	$\mu_T$ (h$^{-1}$)	$r_i$ (h$^{-1}$)
Ⅰ	n/a	n/a	0.022 [0.016, 0.029]	0.33 [0.31, 0.34]	0.094 [0.087, 0.10]	0.082 [0.071, 0.094]	0.36 [0.34, 0.39]	0.17 [0.14, 0.19]
Ⅱ
Ⅲ
Ⅳ	0.0091 [0.0052, 0.014]	0.14 [0.072, 0.18]
Ⅴ
Ⅵ			0.41 [0.37, 0.45]

---

 | Show Table

DownLoad: CSV

Figure 4.  Temporal dynamics of the experimental data (black) and the model (blue) in scenarios Ⅳ (evofosfamide alone, panels A and B), and Ⅴ (hypoxic tumors that received combination immunotherapy and evofosfamide, panels C and D). Left column represents the tumor volume, and right column the well-vascularized fraction. The dashed lines in blue indicate the bounds of the 95% confidence interval of the model output. The mean and 95% confidence interval from the data is also indicated in black. The mean absolute percent error for the tumor volume is below 28% and the CCC is above 0.8 in the two scenarios calibrated.

DownLoad: Full-Size Img PowerPoint

The last step from our framework is to capture the tumor temporal evolution of the data from scenario Ⅵ (normoxic tumors that received combination immunotherapy and evofosfamide). Using the parameter values calibrated to the data from scenarios Ⅰ-Ⅴ, our model was unable to accurately predict the tumor dynamics in scenario Ⅵ, achieving a CCC of only -0.08). To increase the prediction accuracy (i.e., increase the CCC) of the change in tumor size over time in scenario Ⅵ, we fixed seven parameters to the values obtained in calibration steps Ⅰ - Ⅴ and calibrated the remaining one to the data from scenario Ⅵ. We test each parameter from the highest total-effect index to the lowest. In Figure 5, we present the sensitivity analysis of the eight model parameters and the initial condition of the well-vascularized fraction ($V_0$) for four scenarios (i.e., immunotherapy plus evofosfamide (panel A), evofosfamide (panel B), immunotherapy (panel C), and control (panel D)). We computed the sensitivity index of each parameter for all time points. Thus, our results cover the same time range as the duration of the experimental data. The sensitivity analysis shows that, without the immunotherapy (i.e., panels B and D), the tumor growth rate is the most important parameter, followed by the death rate induced by the immune system. However, when immunotherapy is delivered (i.e., panels A and C), the increase in tumor death rate due to immunotherapy becomes the most important parameter, followed by the tumor growth rate, and the death rate due to the immune system. It is worth noting that the period during which the increase in tumor death rate due to immunotherapy, $\gamma_i$, is of greatest importance is prolonged when evofosfamide is administered. In panel A (i.e., immunotherapy plus evofosfamide), the interval of highest importance is between days 7.6 and 26.6 (19.0 days) in the simulation interval of ^[6,39], while in panel C (i.e., immunotherapy) the interval is between days 7.6 and 19.2 (only 11.6 days, nearly 40% less). Based on these results, we found that we only need to change the value of $\gamma_i$ to capture the observed changes in tumors size in scenario Ⅵ. (Table 2 presented all calibrated values from all scenarios). Note that, to capture the observed tumor volume changes in scenario Ⅵ, we need to increase the value of $\gamma_i$ by a factor of $18.6 \pm 9.3$ compared to the value used in the other scenarios. This indicates that, in normoxic tumors, the delivery of evofosfamide dramatically increases the efficacy of the immunotherapy. Interestingly, this increase in $\gamma_i$ does not translate into a significant change in the total tumor death rate. The total tumor death rate is the tumor death rate by the immune system, $\mu_T$, plus the increase in death rate by the immunotherapy, $\gamma_i$. Thus, the total tumor death rate increases by a factor of $2.52 \pm 0.53$ (see supplemental Figure S1) when adding evofosfamide into the normoxic tumors treated with immunotherapy. In Figure 6, the experimental data and the calibrated model are presented for scenario Ⅵ. The mean absolute percent error for the tumor volume is 30.94% and the CCC is 0.98. The high percent error is due to tumor volumes that reach zero after day 23.

Figure 5.  Sensitivity analysis of the model in the following scenarios: immunotherapy plus evofosfamide (panel A), evofosfamide (panel B), immunotherapy (panel C), and control (panel D). The sensitivity index is computed throughout the entire experimental timeline, offering insights into how different factors play a role at different stages of tumor growth and response to treatment. In panels B and D, the tumor growth rate has the highest total effect index ($k_T$, blue line), followed by the death rate caused by the immune system ($\mu_T$, red line). However, with the immunotherapy, in panels A and C, the increase in tumor death rate due to immunotherapy ($\gamma_i$, purple line) has the highest total effect index. This is due to the increase of immune cells when delivering immunotherapy in panels A and C.

DownLoad: Full-Size Img PowerPoint

Figure 6.  Temporal dynamics of the experimental data (black) and the model (blue) in scenario Ⅵ (normoxic tumors that received combination immunotherapy and evofosfamide). Left column represents the tumor volume, and right column the well-vascularized fraction. The dashed lines in blue indicate the bounds of the 95% confidence interval of the model output. The mean and 95% confidence interval from the data is also indicated in black. The mean absolute percent error for the tumor volume is below 31% and the CCC and PCC are both above 0.98.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

We have developed a system of ordinary differential equations capable of reproducing the tumor and its vascularization fraction growth and decline under the combined effects of immunotherapy and evofosfamide treatment. To the best of our knowledge, this represents the first effort to mathematically describe the interactions of these two drugs. Our model was able to recapitulate the experimental data showing that evofosfamide can rescue the immune response and sensitize hypoxic tumors to immunotherapy ^[8]. Following the framework presented in Figure 2, we were able to isolate the effects of immunotherapy and evofosfamide. We calibrated the parameters that model the tumor growth and immunotherapy using the data from scenarios Ⅰ, Ⅱ, and Ⅲ (i.e., scenarios in which evofosfamide was not delivered). In scenario Ⅱ, the hypoxic tumors did not respond to immunotherapy. The next step was to calibrate the evofosfamide terms using the data from scenarios Ⅳ and Ⅴ, while keeping the previously calibrated parameters fixed. In scenario Ⅴ, the hypoxic tumors responded to immunotherapy following the evofosfamide treatment. Thus, the response of hypoxic tumors to immunotherapy, observed in scenario Ⅴ, is only due to the terms of our model effected by evofosfamide. More specifically, the increase in tumor vascularization leads to an increase in immunotherapy efficacy. Reeves et al. have shown that evofosfamide can increase tumor vascularization ^[8], and this may facilitate the delivery of immunotherapy to the tumor cells. However, when predicting the normoxic tumor response to immunotherapy and evofosfamide (i.e., scenario Ⅵ), our model was unable to obtain a good agreement with the experimental data. Our model was only able to reproduce scenario Ⅵ when the death rate of tumor cells due to immunotherapy increases directly as a result of evofosfamide, rather than solely through an indirect effect due to the increase in tumor vascularization. This finding conflicts with our experimental data, which shows that evofosfamide alone has no effect on tumor size (i.e., scenario Ⅳ). However, it should be noted that previous studies have demonstrated that evofosfamide can improve the effectiveness of immunotherapy, in addition to increasing tumor vascularization ^[9,31]. Phase Ⅰ experiments conducted by ^[9] and ^[31] have shown that evofosfamide can resensitize tumors to immunotherapy, supporting the idea that it may have a direct effect on the death rate of tumor cells due to immunotherapy. Therefore, our model's results are consistent with other studies demonstrating that evofosfamide can enhance the effectiveness of immunotherapy. In our results, even though the death rate by immunotherapy, $\gamma_i$, increased by a factor of $18.6 \pm 9.3$ from the hypoxic tumor to the normoxic tumor, the total death rate (i.e., the tumor constant death rate, $\mu_T$, plus death rate by immunotherapy) increased by a factor of $2.52 \pm 0.53$, which is consistent with experimental data.

Sensitivity analysis is a valuable technique to help understand the complex behavior of biological systems. To select which parameters would be kept fixed with the values from previous scenarios, and which ones would be calibrated to the data from scenario Ⅵ, we tested the parameters according to their sensitivity analysis rank (from most to least influential). By measuring the effects of model parameters on system outputs, sensitivity analysis helps determine the underlying biological mechanisms that govern system behavior ^[21,22,23]. In time dependent biological processes, sensitivity analysis allows for assessing the relative importance of each model parameter over time. For example, in tumor growth models, sensitivity analysis can identify the most critical model parameters that influence the growth and temporal dynamics of cancer ^[32,33,34]. In our study, the results of the sensitivity analysis identify $\gamma_i$, death rate due to immunotherapy, as the key parameter to change to accurately reproduce the tumor growth dynamics in scenario Ⅵ, revealing that in normoxic tumors, there is heightened response to immunotherapy.

There are a few opportunities for further investigation that arise from this work. From an experimental perspective, our analysis is based on data from mice implanted with colon adenocarcinoma cells and may not necessarily be generalizable to other types of tumors. A possible extension would be experiments with different cell lines, and the inclusion of anti-angiogenic therapies (as done in ^[17]), to verify the synergy between evofosfamide and VEGF inhibitors. Another possible extension is the inclusion of different treatment protocols, altering the order, timing, and sequencing of the therapies. The current results are based on data from a limited number of scenarios and may not accurately capture the full range of possible responses to treatment. Future studies should test how switching the order and doses of the two treatments affect the temporal change in tumor size and hypoxic status. With a broader dataset that included more diverse treatment protocols, the model could be calibrated to predict optimal drug dosing and timing for maximizing the synergistic effects of evofosfamide and immunotherapy (as done in ^[35]). This predictive capability could be valuable in guiding treatment decisions and designing more effective therapeutic strategies. Additionally, it would be desirable to have measurements on vascularity as afforded by (for example) dynamic contrast enhanced magnetic resonance imaging ^[36]. While the imaging allowed for tumor information on hypoxia, additional [¹⁸F]FMISO-PET time points would have provided the opportunity to better characterize the longitudinal changes in hypoxia. Both of these measures would allow for the separation of the vascular and hypoxia contributions to treatment efficacy. This additional imaging data would come with a significant increase in resources and exposure of ionizing radiation to the animals. However, by applying the methods of optimal experimental design (e.g., ^[37,38,39]) to our mathematical model, we can determine the optimal time points for acquiring additional [¹⁸F]FMISO-PET scans to minimize the uncertainty of our models and increase model's prediction accuracy. This approach could help us maximize the information obtained from the additional imaging data while minimizing the associated costs. The additional imaging data could help to validate whether evofosfamide can maintain a reduction of hypoxia that can allow for long-term sustained responses to immunotherapy.

From a mathematical and computational perspective, the model is based on a set of ordinary differential equations, which may not fully capture the complexity of the interactions between evofosfamide, immunotherapy, and tumor cells. Investigating spatial heterogeneity and interactions through partial differential equations or agent-based models could improve the model's accuracy, predictive power, and generate new hypotheses on the underlying biology. Mechanistic models like the quantitative systems pharmacology platform developed by ^[40] incorporate a more extensive set of parameters (282 parameters) to offer detailed insights into underlying biological mechanisms and interactions. Such models enhance comprehension of the immune-oncology system. However, our platform has an advantage of using a minimal parameter set. This enables model calibration using limited data from a single patient, leading to personalized calibration and reduced uncertainty in predictions. Given the available vascularity measurements, we developed a simplified vascularization model that could still capture the observed phenomena to the best extent possible. While more sophisticated models might be constructed with larger datasets, our model adequately captures the phenomena observed in these experiments. It assumes an increase in vascularization after each evofosfamide treatment, effectively depicting tumor growth across different experimental scenarios. Notably, our model excludes the direct cytotoxic effects of evofosfamide. This assumption aligns with the specific experimental data used to calibrate the model ^[8]. The experiments indicated minimal direct cytotoxic impact on the tumor volume, although it did influence combination therapy at the cellular level. Consequently, our focus lay in developing a simple model calibrated to the available experimental data, capable of capturing the key observed phenomena. We opted to exclude the direct cytotoxic effects and concentrate on evofosfamide's role as an oxygenating and sensitizing agent, as evident in the experimental data. However, even in the absence of a term representing the direct cytotoxic effects, evofosfamide indirectly impacts tumors. It increases the well-vascularized fraction, increasing the tumor death rate via the immune system's action through the term $\mu_T T V$. Finally, the current model is based on the assumption that evofosfamide increases tumor vascularization, facilitating the subsequent delivery of immunotherapy to tumor cells. However, as seen in scenario Ⅵ, this is not enough to capture the effects on normoxic cells, necessitating an increase in the death rate due to immunotherapy. Further study is needed to determine the mechanisms behind this effect and how it may be incorporated into the model. Our current model could be adapted and calibrated with experimental data for other therapies that affect oxygenation (or vasculature) to predict their potential synergistic effects with immunotherapy. While the specific effects may vary between different drugs, the efficacy of any hypoxia-activated prodrug in reducing hypoxic regions and enhancing the effects of immunotherapy would be critical for its success.

The results from our model, in combination with the experiments presented in ^[8], demonstrate that it would be feasible to use of [¹⁸F]FMISO-PET in conjunction with hypoxia-activated prodrug therapy in human clinical trials. F[¹⁸F]FMISO-PET is noninvasive, clinically-available, and can provide information on temporal and spatial variations in hypoxia within tumors ^[41]. It has potential to open new opportunities for personalized and image-guided adaptive immunotherapy approaches.

5.   Conclusions

Our mathematical model is able to capture a tumor's response to evofosfamide and checkpoint blockade immunotherapy. The model also captures the interaction between these two types of treatments and their effects on a tumor's vascularization. Notably, the model was able to capture the tumor's temporal evolution observed in the five different treatment conditions (i.e., control, hypoxic and normoxic tumors treated with immunotherapy, treated with evofosfamide, and hypoxic tumors treated with immunotherapy and evofosfamide), with an average CCC of $0.86\pm 0.05$. Interestingly, the model calibrated with parameter values from those five scenarios was not able to predict the response of normoxic tumors to combination immunotherapy and evofosfamide. However, increasing the tumor death rate due to immunotherapy ($\gamma_i$) by a factor of $18.6$ allowed us to predict the tumor growth with a CCC of $0.98$. These results strongly indicate that evofosfamide not only increases the tumor vascularization, as considered in our model, but also increases efficacy of the immunotherapy, suggesting a synergistic interaction between the two treatments. This potential synergy warrants further investigation, and may have significant implications for the optimization of the dosing and scheduling of these drugs for maximal clinical benefit. This work provides a framework for identifying testable hypotheses for the effects and actions of evofosfamide and immunotherapy in the treatment of tumors.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

We thank the American Cancer Society Grant RSG-18-006-01-CCE and the National Institute of Health for funding via R01CA240589, R01CA276540, and U24CA226110. We thank the Cancer Prevention and Research Institute of Texas for support through CPRIT RR160005. T.E.Y. is a CPRIT Scholar in Cancer Research. We also thank the Texas Advanced Computing Center for providing high performance computing resources. Research reported in this publication was supported by the National Institute of General Medical Sciences of the National Institutes of Health under Award Number T32GM135028. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Conflict of interest

The authors declare there is no conflict of interest.

Supplementary

Figure S1.  Scaling factor of the total tumor death rate (i.e., the tumor death rate by the immune system, constant in time, plus the increase in death rate by the immunotherapy, changing in time as the treatment is delivered) over time, computed as the ratio of the total tumor death rate in one treatment scenario to the total tumor death rate in another scenario. The scaling factor is shown for: A) immunotherapy compared to control, B) immunotherapy plus evofosfamide compared to control, and C), immunotherapy plus evofosfamide compared to immunotherapy. The average and standard deviation of the scaling factor in panels A, B, and C are $1.10 \pm 0.04$, $2.79 \pm 0.66$, and $2.52 \pm 0.53$, respectively. The scaling factor provides a useful way to visualize and compare the efficacy of different treatment scenarios in inducing tumor death.

DownLoad: Full-Size Img PowerPoint