An in silico exploration of combining Interleukin-12 with Oxaliplatin to treat liver-metastatic colorectal cancer

Abstract

Background

Combining anti-cancer therapies with orthogonal modes of action, such as direct cytotoxicity and immunostimulatory, hold promise for expanding clinical benefit to patients with metastatic disease. For instance, a chemotherapy agent Oxaliplatin (OXP) in combination with Interleukin-12 (IL-12) can eliminate pre-existing liver metastatic colorectal cancer and protect from relapse in a murine model. However, the underlying dynamics associated with the targeted biology and the combinatorial space consisting of possible dosage and timing of each therapy present challenges for optimizing treatment regimens. To address some of these challenges, we developed a predictive simulation platform for optimizing dose and timing of the combination therapy involving Mifepristone-induced IL-12 and chemotherapy agent OXP.

Methods

A multi-scale mathematical model comprised of impulsive ordinary differential equations was developed to describe the interaction between the immune system and tumor cells in response to the combined IL-12 and OXP therapy. An ensemble of model parameters were calibrated to published experimental data using a genetic algorithm and used to represent three different phenotypes: responders, partial-responders, and non-responders.

Results

The multi-scale model captures tumor growth patterns of the three phenotypic responses observed in mice in response to the combination therapy against a tumor re-challenge and was used to explore the impacts of changing the dose and timing of the mixed immune-chemotherapy on tumor growth subjected to a tumor re-challenge in mice. An increased ratio of CD8 + T effectors to regulatory T cells during and after treatment was key to improve tumor control in the responder cohort. Sensitivity analysis indicates that combined OXP and IL-12 therapy worked more efficiently in responders by increased priming of T cells, enhanced CD8 + T cell-mediated killing, and functional inhibition of regulatory T cells. In a virtual cohort that mimics non-responders and partial-responders, simulations show that an increased dose of OXP alone would improve the response. In addition, enhanced IL-12 expression alone or an increased number of treatment cycles of the mixed immune-chemotherapy can barely improve tumor control for non-responders and partial responders.

Conclusions

Overall, this study illustrates how mechanistic models can be used for in silico screening of the optimal therapeutic dose and timing in combined cancer treatment strategies.

Background

Carcinomas of the colon or rectum, termed colorectal cancer, are the third most common cancer diagnosed in both men and women in the United States. The American Cancer Society estimates the number of new cases of colorectal cancer in the United States for 2019 at 145,600 [1]. With 60,000 fatalities per year, colorectal cancer is second only to lung cancer as a cause of cancer-related deaths in the United States. Upon diagnosis, 10%-20% of patients have already developed liver metastases while 70% of patients with colorectal cancer ultimately develop liver metastases. Unfortunately, the prognosis for patients with liver metastatic colorectal cancer is poor because hepatectomy, palliative chemotherapy and symptomatic treatments are the only available options [2].

Interleukin-12 (IL-12) is a potent immunostimulatory cytokine that activates the proliferation and function of key cellular effectors of innate and adaptive immunity such as T lymphocytes and natural killer (NK) cells [3–5]. While toxicity is a serious obstacle for use of IL-12 as a systemic therapy in humans, an attractive alternative is to use adenoviral vectors to induce expression in specific tissues. However, transgene expression tends to be transient and the efficacy of re-administration is impaired by the rapid emergence of neutralizing antibodies [3]. To allow a good control of the strength and duration of IL-12 expression, high-capacity adenoviral vectors containing a liver-specific, Mifepristone-inducible system for the expression of murine IL-12 (HC-Ad/RUmIL-12) were recently designed to control primary or metastatic liver cancer [6]. Since stand-alone chemo- or radiotherapeutic regimens are insufficient (with a few notable exceptions) to eradicate neoplastic lesion [7], HC-Ad/RUmIL-12 was combined with chemotherapy agent Oxaliplatin (OXP) to treat liver-implanted colon cancer cells [6]. As a consequence of the combination therapy, pre-existing liver metastases of colorectal cancer were eradicated, and enhanced establishment of a protective immune response against tumor re-challenge and increased overall survival of animals were observed. In addition, a dramatic increase in the ratio of cytotoxic CD8+ T lymphocytes to immunosuppressive cell populations was detected in the tumor microenvironment [3].

Mathematical modeling using systems of ordinary differential equations (ODEs) can improve the design and administration of cancer treatments, especially when experimental data are incorporated [8–13]. In silico screening of parameter regions that seem most promising for optimal timing and dosage of therapy can be suggested using calibrated mathematical models and clinical trials can focus on those regions [13–17]. For instance, a quantitative systems pharmacology model [8] was developed to reproduce experimental data of CT26 tumor dynamics upon administration of radiation therapy and an anti-PD-L1 agent in [18] and [19]. The calibrated model was further used as an in silico tool to predict the best treatment combination schedules and sequences. Over the past decades, a variety of ODE-based mathematical models have been developed to better understand cancer progression and response to immunotherapy (see details in [20–25]). In exploring immunotherapy in combination with other treatment modalities, de Pillis et. al developed an ODE model governing cancer growth on a cell population level with a combination of immuno-chemotherapy treatments [26–29]. In addition, Kim and colleagues formulated a mathematical model of therapy with oncolytic viruses that simultaneously express immunostimulatory cytokines and costimulatory molecules [12]. Inspired by these studies, we developed an impulsive ODE model to represent the interaction between tumor and immune system in response to the chemotherapy drug OXP combined with liver-specific expression of IL-12 therapy to explore therapeutic options in the context of liver metastatic colorectal cancer. The current model extends the impulsive ODE model in [13] that only considered an immunotherapy initiated by an adenovirus vaccination to stimulate a tumor-associated antigen-specific T cell response.

The structure of this paper is as follows. First, we present a multi-scale mechanistic model of anti-tumor immunity and tumor growth in response to a combined immuno-chemotherapy using a set of impulsive ODEs. Second, we describe how we calibrated the model parameters against published experimental data using a genetic algorithm. Next we investigate the stability of tumor-free and high-tumor equilibria based on the linearized system. Then we study how alter parameter values may change the tumor growth dynamics. Finally, we used the simulation platform to explore potential ways to improve treatment regimes for non-responders and partial responders.

Methods

Our method was to develop a multi-scale impulsive ODE model based on our understanding of the corresponding biology, which is described in the following paragraphs. Numerical solutions of the model were obtained using simulators generated by C Sharp. The resulting mechanistic mathematical model was calibrated against existing experimental data.

A genetic algorithm was used to find parameter sets that closely match the experimental data in [6]. Each parameter set was modeled using an individual chromosome in order to apply the algorithm to search in the parameter space. For each generation, the impulsive ODE set was solved using the Runge-Kutta method of order four for each individual or parameter set [30]. The fitness function value, or variance, was calculated using the sum of error squared between experimental data and corresponding model predictions. To reduce the dependence of our model predictions about optimal treatment strategies for the combined therapy on any individual calibrated set of parameter values, we generated an ensemble of 30 parameter sets for each phenotypic cohort (i.e., responder, partial responder, and non-responder) that generated similar good fits against the experimental data. The simulation results using these ensemble of parameter sets were characterized in terms of the median and upper and lower bounds that enclose 90% of the predicted responses. Simulations start on day 0, which corresponds to the time of tumor implantation. At the initial time point, we assume that there is no activated tumor specific effector T cells present in the blood and at the site of the tumor. The calibrated mechanistic model was then used to investigate the long-term behavior through stability analysis. Details of model development, parameter calibration, goodness of fit, difference between major variables of immune response of responder mice, partial responder mice, and non-responder mice after the combination treatment, and local stability analysis are discussed in the following subsections.

Results

A multi-scale model of tumor growth subject to IL-12 and OXP therapy

Our mathematical model is based on the experimental data presented by Manuela Gonzalez-Aparicio and collaborators in [6] using the MC38Luc1 cell lines for murine metastatic colorectal cancer. Using this mouse model, OXP and IL-12 combination therapy eradicated pre-existing liver metastases, established a protective immune response against tumor re-challenge, and increased overall survival of animals. To better understand the dynamics of the primary response to adenovirus-mediated induction of an anti-tumor immune response, we developed a three-compartment mathematical model to quantify the cytotoxic CD8 ⁺ T cell response to IL-12 and OXP combined therapy and subsequent inhibition of tumor cell growth, as shown schematically in Fig. 1. Model parameters and their meanings are listed in Table 1.

Fig. 1

Schematic diagram illustrating the interactions among species present in the three compartments. State variables and transport relations are shown in black. Parameters are in red while influence relationships are in blue. Naïve CD8 ⁺ T cells (T_N) are activated and become CD8 ⁺ T effectors (T_E1) when they encounter tumor antigen presented by the antigen presenting cells (APC₁) in the lymph node. Once activated, effector CD8 ⁺ T cells circulate within the blood (T_E2) and enter tumor microenvironment (T_E3) where they are retained upon recognition of the corresponding tumor-associated antigen. Effector CD8 ⁺ T cells secrete Interferon gamma (IFN_γ) which assist with the CD8 ⁺ T cell-mediated killing of tumor cells ($C_{\text {MHCI}^{+}}$ and $C_{\text {MHCI}^{-}}$) through increased presentation of tumor-associated antigens by Major Histocompatibility Complex protein class I (MHCI). During this process, IL-12 (IL) helps promote T cell proliferation and suppresses regulatory T (T_R) cells’ proliferation and immunosuppressive action on effector CD8 ⁺ T cells. In addition, the chemotherapy drug Oxaliplatin in the lymph node and tumor (OXP_i where i=1,3) will kill fast-proliferating cells such as T effectors and tumor cells

Table 1 List of parameters in the model

Among these three compartments, we consider the dynamics of fifteen state variables that are regulated by the following governing biological processes and assumptions:

• Naïve CD8⁺ T cells (T_N, units: cells per mm³). We assume that naïve CD8 ⁺ T cells are produced at a constant rate c₁ from thymus and die naturally at a rate k_d1·T_N [31]. Naïve T cells are recruited and activated by tumor antigens presented by APC₁ (antigen-presenting cells in lymph node) at a rate $c_{2}\cdot T_{N}\cdot \frac {\text {APC}_{1}}{\text {APC}_{1}+g_{1}}$ [32–34].

  $$ {}\frac{dT_{N}}{dt}=c_{1}-k_{d1}\cdot T_{N}-c_{2}\cdot T_{N}\cdot\frac{\text{APC}_{1}}{\text{APC}_{1}+g_{1}} $$

  (1)

• Effector CD8⁺ T cells in lymph node (T_E1, units: cells per mm³). The increase in the rate of concentration of effector CD8 ⁺ T cells in the lymph node due to activation of naïve CD8 ⁺ T cells from the blood stream is given by $c_{2} \cdot \frac {T_{N}\text {Vol}_{b}}{\text {Vol}_{ln}}\cdot \frac {\text {APC}_{1}}{\text {APC}_{1}+g_{1}}$, where Vol_b=1.4∗10³ mm³ is the volume of the blood compartment [35] and Vol_ln=0.25 mm³ is the volume of the lymph node compartment [36]. We assume that the natural death of effector T cells in the lymph node is negligible. Effector CD8 ⁺ T cells in the lymph node proliferate at a rate proportional to T_E1, a saturable term that represents antigen presenting cells (APC) and defined by $\frac {\text {APC}_{1}}{\text {APC}_{1}+g_{2}}$, and an immune checkpoint term defined by $\frac {\alpha }{\alpha +T^{2}_{E1}}$, where α is the square root of the saturation constant of T_E1 [37]. We also assume that influx rate of effector T cells from blood to lymph node is $a_{21}\cdot \frac {T_{E2}\text {Vol}_{b}}{\text {Vol}_{ln}}$, where T_E2 is the concentration of T effectors in blood, and a₁₂·T_E1 is the efflux rate. We assume that T_E1 cells are killed by chemotherapy agent OXP₁ (Oxaliplatin in lymph node) at the rate $k_{d2}\cdot \frac { T_{E1}\cdot \text {OXP}_{1}}{\text {OXP}_{1}+g_{3}}$.

  $$ {}\begin{aligned} \frac{dT_{E1}}{dt}&\!=c_{2} \cdot \frac{T_{N}\text{Vol}_{b}}{\text{Vol}_{ln}}\cdot\frac{\text{APC}_{1}}{\text{APC}_{1}+g_{1}} +k_{p1}\cdot T_{E1}\cdot \frac{\text{APC}_{1}}{\text{APC}_{1}+g_{2}}\\ &\quad \cdot\frac{\alpha}{\alpha+T^{2}_{E1}}+a_{21}\cdot\frac{T_{E2}\text{Vol}_{b}}{\text{Vol}_{ln}}-a_{12}\cdot T_{E1}-k_{d2}\\ &\quad\cdot \frac{ T_{E1}\cdot\text{OXP}_{1}}{\text{OXP}_{1}+g_{3}} \end{aligned} $$

  (2)

• Antigen Presenting Cells in lymph node (APC₁, units: cells per mm³). We assume that APCs in the lymph node have a natural death rate of k_d3·APC₁ and the influx rate of APCs from tumor to lymph node is $b_{31}\cdot \text {APC}_{3}\cdot \frac {\text {Vol}_{t}}{\text {Vol}_{ln}}$, where APC₃ is the concentration of APCs in tumor, ${Vol}_{t}= \frac {\epsilon +C(t)}{1-V_{i}\cdot T_{E3}}$ (since Vol_t=ε+C(t)+V_i·T_E3·Vol_t, where T_E3 is the concentration of T effectors in tumor) is the volume of the tumor compartment, ε is a small positive constant representing a small volume of tissue that excludes tumor and effector CD8 ⁺ T cells in the tumor compartment, where C(t) is the volume of tumor cells in mm³. The total volume of tumor cells is comprised of the volumes of major histocompatibility complex (MHC) class I positive tumor cells ($\phantom {\dot {i}\!}C_{MHCI^{+}}$) and MHC class I negative tumor cells ($\phantom {\dot {i}\!}C_{MHCI^{-}}$). The average size of a T effector cell (V_i) is equal to $\phantom {\dot {i}\!}10^{-7}\,mm^{3}$ [38].

  $$\begin{array}{@{}rcl@{}} {}\frac{d \text{APC}_{1}}{dt}&=&-k_{d3}\cdot\text{APC}_{1}+b_{31}\cdot \text{APC}_{3}\cdot \frac{\text{Vol}_{t}}{\text{Vol}_{ln}} \end{array} $$

  (3)

• Chemotherapy agent Oxaliplatin in lymph node (OXP₁, units: mg per kg). We assume that OXP decays naturally at a rate k_d4·OXP₁ and the influx rate of Oxaliplatin (OXP) from blood to lymph node is $c_{21}\cdot \frac {\text {OXP}_{2}\text {Vol}_{b}}{\text {Vol}_{ln}}$, where OXP₂ is the concentration of OXP in blood.

  $$\begin{array}{@{}rcl@{}} {}\frac{d \text{OXP}_{1}}{dt}&=&-k_{d4}\cdot \text{OXP}_{1}+c_{21}\cdot \frac{\text{OXP}_{2}\text{Vol}_{b}}{\text{Vol}_{ln}} \end{array} $$

  (4)

• Effector CD8⁺ T cells in blood (T_E2, units: cells per mm³). We assume the effector CD8 ⁺ T cells die naturally in blood at a rate k_d5·T_E2. The influx rate of effector CD8 ⁺ T cells from lymph node to blood is equal to $a_{12}\cdot \frac {T_{E1}\text {Vol}_{ln}}{\text {Vol}_{b}}$ and the efflux rate of effector CD8 ⁺ T cells from blood to lymph node is equal to a₂₁·T_E2. The influx rate of CD8 ⁺ T effectors from the tumor to blood is $a_{32}\cdot \frac {C_{MHCI^{-}}}{\epsilon +C(t)}\cdot \frac {T_{E3}\text {Vol}_{t}}{\text {Vol}_{b}}$, where T_E3 is the concentration of T effectors in tumor and the efflux rate of CD8 ⁺ T effectors from blood to tumor is a₂₃·T_E2.

  $$\begin{array}{@{}rcl@{}} {}\frac{dT_{E2}}{dt}&\,=\,&-k_{d5} \cdot\! T_{E2} +a_{12}\cdot \frac{T_{E1}\text{Vol}_{ln}}{\text{Vol}_{b}} \,-\, a_{21}\cdot T_{E2}\,-\,a_{23}\cdot T_{E2} \\&&+a_{32}\cdot \frac{C_{MHCI^{-}}}{\epsilon+C(t)}\cdot \frac{T_{E3}\text{Vol}_{t}}{\text{Vol}_{b}} \end{array} $$

  (5)

• Antigen Presenting Cells in blood (APC₂, units: cells per mm³). A logistic growth pattern $r_{2}\cdot \text {APC}_{2}\cdot \bigg (1-\frac {\text {APC}_{2}}{K}\bigg)$ is used for APCs in blood where r₂ is the growth rate constant and K is the carrying capacity. We assume a b₂₃·APC₂ efflux rate of APCs from blood to tumor.

  $$\begin{array}{@{}rcl@{}} {}\frac{d \text{APC}_{2}}{dt}&=&r_{2}\cdot \text{APC}_{2}\cdot \bigg(1-\frac{\text{APC}_{2}}{K}\bigg)-b_{23}\cdot \text{APC}_{2} \end{array} $$

  (6)

• Chemotherapy agent Oxaliplatin in blood (OXP₂, units: mg per kg). We assume that OXP decays naturally at a rate k_d4·OXP₂ and the efflux rates of OXP from blood to lymph node and from blood to tumor are c₂₁·OXP₂ and c₂₃·OXP₂, respectively. The source of OXP is provided by each administration whose dose and time are reflected by the discrete Eq. (17).

  $$\begin{array}{@{}rcl@{}} {}\frac{d \text{OXP}_{2}}{dt}&=&-k_{d4}\cdot \text{OXP}_{2}-c_{21}\cdot \text{OXP}_{2}-c_{23}\cdot \text{OXP}_{2} \end{array} $$

  (7)

• Effector CD8⁺ T cells in tumor microenvironment (T_E3, units: cells per mm³). We assume that effector CD8 ⁺ T cells can proliferate locally upon recognition of the corresponding tumor antigen presented by MHC class I positive tumor cells upon IL-12 stimulation and subject to suppression from T regulatory cells at a saturable rate equal to $k_{p2}\cdot \frac {C_{MHCI^{+}}}{\epsilon +C(t)}\cdot \frac {\text {IL}}{\text {IL}+g_{4}}\cdot \frac {T_{E3}}{T_{R}+g_{5}}$, where IL is the concentration of IL-12 and T_R is the concentration of regulatory T cells [39, 40]. The influx rate of effector CD8 ⁺ T cells from the blood to tumor is defined by $a_{23}\cdot \frac {T_{E2}\text {Vol}_{b}}{\text {Vol}_{t}}$. The efflux rate of effector CD8 ⁺ T cells from the tumor to blood is $a_{32}\cdot T_{E3}\cdot \frac {C_{MHCI^{-}}}{\epsilon +C(t)}$. Effector T cells have a finite lifespan and die within the tumor microenvironment at a rate equal to k_d6·T_E3. T effector cells are assumed to be killed by chemotherapy agent OXP (OXP₃) at the rate $k_{d2}\cdot T_{E3}\cdot \frac {\text {OXP}_{3}}{\text {OXP}_{3}+g_{6}}$ in the tumor microenvironment.

  $$ {}\begin{aligned} \frac{dT_{E3}}{dt}&=a_{23}\cdot \frac{T_{E2}\text{Vol}_{b}}{\text{Vol}_{t}}-a_{32}\cdot T_{E3}\cdot \frac{C_{MHCI^{-}}}{\epsilon+C(t)} + k_{p2}\\ & \quad \cdot \frac{C_{MHCI^{+}}}{\epsilon+C(t)}\cdot\frac{\text{IL}}{\text{IL}+g_{4}}\cdot\frac{T_{E3}}{T_{R}+g_{5}}-k_{d6}\cdot T_{E3}-k_{d2}\\ & \quad\cdot T_{E3}\cdot\frac{ \text{OXP}_{3}}{\text{OXP}_{3}+g_{6}} \end{aligned} $$

  (8)

• Interferon gamma (IFNγ, units: pg per mm³). We assume that IFN_γ is secreted solely by effector CD8 ⁺ T cells within the tumor with stimulation from IL-12 and inhibition from regulatory T cellsat a rate of $c_{4}\cdot \frac {\text {IL}}{\text {IL}+g_{7}}\cdot \frac {T_{E3}}{T_{R}+g_{8}}$ [41]. While this assumption may not hold in all model systems, the presence of IFN_γ in the tumor was dependent on CD8 ⁺ T cell activation [42]. IFN_γ decays at a rate proportional to its concentration with a rate constant k_d7.

  $$\begin{array}{@{}rcl@{}} {}\frac{d \text{IFN}_{\gamma}}{dt}&=&-k_{d7}\cdot \text{IFN}_{\gamma}+c_{4}\cdot \frac{\text{IL}}{\text{IL}+g_{7}}\cdot\frac{T_{E3}}{T_{R}+g_{8}} \end{array} $$

  (9)

• Antigen Presenting Cells in tumor (APC₃, units: cells per mm³). We assume that APCs in the tumor microenvironment have a natural death rate of k_d3·APC₃, the influx rate of APCs from blood to tumor is $b_{23}\cdot \text {APC}_{2}\cdot \frac {\text {Vol}_{b}}{\text {Vol}_{t}}$ and APCs take tumor antigen in tumor microenvironment and migrate to the lymph node to present tumor antigens to T cells at the rate of b₃₁·APC₃.

  $$ {}\frac{d \text{APC}_{3}}{dt}=b_{23}\cdot \text{APC}_{2} \cdot \frac{\text{Vol}_{b}}{\text{Vol}_{t}}-b_{31}\cdot \text{APC}_{3}-k_{d3}\cdot\text{APC}_{3} $$

  (10)

• Interleukin-12 (IL, units: ng per ml). Interleukin-12 (IL-12) is produced by APCs at a rate of c₅·APC₃ and decays naturally at a rate of k_d8·IL. The extra IL-12 expression obtained through the combined therapy is approximated using the discrete Eq. (16).

  $$\begin{array}{@{}rcl@{}} {}\frac{d\text{IL}}{dt}&=&c_{5}\cdot \text{APC}_{3}-k_{d8}\cdot\text{IL} \end{array} $$

  (11)

• Chemotherapy agent Oxaliplatin in tumor (OXP₃, units: mg per kg). We assume that OXP decays naturally at a rate k_d4·OXP₃ and the influx rate of OXP from blood to tumor is $c_{23}\cdot \frac {\text {OXP}_{2}\cdot \text {Vol}_{b}}{\text {Vol}_{t}}$.

  $$\begin{array}{@{}rcl@{}} {}\frac{d \text{OXP}_{3}}{dt}&=&-k_{d4} \cdot \text{OXP}_{3}+c_{23}\cdot \frac{\text{OXP}_{2} \cdot \text{Vol}_{b}}{\text{Vol}_{t}} \end{array} $$

  (12)

• Regulatory T cells (T_R, units: cells per mm³). Regulatory T cells are produced at a constant rate c₆ from thymus and die naturally at a rate k_d9·T_R. We assume that regulatory T cells proliferate at a rate of $k_{p3}\cdot \frac {C_{MHCI^{-}}}{\epsilon +C(t)}\cdot \frac {T_{E3}}{T_{E3}+g_{9}}\cdot \frac {T_{R}}{\text {IL}+g_{10}}$ [3, 4, 6].

  $$ {}\frac{dT_{R}}{dt}=c_{6}-k_{d9}\cdot T_{R}+k_{p3} \cdot\frac{C_{MHCI-}}{\epsilon+C(t)}\cdot\frac{T_{E3}}{T_{E3}+g_{9}}\cdot\frac{T_{R}}{\text{IL}+g_{10}} $$

  (13)

• MHC class I positive tumor cells ($C_{MHCI^{+}}\phantom {\dot {i}\!}$, units: mm³). MHC class I positive tumor cells are converted from MHC class I negative tumor cells ($C_{MHCI^{-}}\phantom {\dot {i}\!}$) with the assistance of IFNγ at a rate $c_{7}\cdot \frac {\text {IFN}{\gamma }}{\text {IFN}{\gamma }+g_{11}}\cdot C_{MHCI^{-}}$ and the rate of effector CD8 + T cell-mediated killing of MHC class I positive tumor cells is $k_{d11} \cdot \bigg (1+\frac {\text {OXP}_{3}}{\text {OXP}_{3}+g_{12}}\bigg)\cdot \frac {C_{MHCI^{+}}}{\epsilon +C(t)}\cdot \frac {T_{E3}}{T_{R}+g_{13}}$ [6, 31]. We assume that the dilution rate of MHC class I positive tumor cells due to proliferation is $k_{p4} \cdot C_{MHCI^{+}}\phantom {\dot {i}\!}$. The natural death rate of MHC class I positive tumor cells is assumed to be $k_{d10}\cdot C_{MHCI^{+}}\phantom {\dot {i}\!}$ and MHC class I positive tumor cells are killed by chemotherapy agent OXP in tumor at a rate $k_{d2}\cdot C_{MHCI^{+}}\cdot \frac { \text {OXP}_{3}}{\text {OXP}_{3}+g_{14}}\phantom {\dot {i}\!}$.

  $$ {}\begin{aligned} \frac{dC_{MHCI^{+}}}{dt}&=c_{7}\cdot \frac{\text{IFN}_{\gamma}}{\text{IFN}_{\gamma}+g_{11}}\cdot C_{MHCI^{-}}-k_{d10} \cdot C_{MHCI^{+}}\\ &\quad-k_{d11}\cdot\frac{ C_{MHCI^{+}}}{\epsilon+C(t)}\cdot\bigg(1+\frac{\text{OXP}_{3}}{\text{OXP}_{3}+g_{12}}\bigg)\cdot\frac{T_{E3}}{T_{R}+g_{13}}\\ &\quad-k_{p4}\cdot C_{MHCI^{+}}- k_{d2}\cdot C_{MHCI^{+}}\cdot \frac{\text{OXP}_{3}}{\text{OXP}_{3}+g_{14}} \end{aligned} $$

  (14)

• MHC class I negative tumor cells ($C_{MHCI^{-}}\phantom {\dot {i}\!}$, units: mm³). MHC class I negative tumor cells are converted to MHC class I positive tumor cells with the assistance of IFN_γ at a rate of $c_{7}\cdot \frac {\text {IFN}_{\gamma }}{\text {IFN}_{\gamma }+g_{11}}\cdot C_{MHCI^{-}}\phantom {\dot {i}\!}$. We assume that the proliferation rate of MHC class I positive tumor cells is equal to $\phantom {\dot {i}\!}2\cdot k_{p4} \cdot C_{MHCI^{+}}$. As MHC class I positive tumor cells proliferate, they lose MHC class I expression and become MHC class I negative cells. A logistic growth pattern is assumed for the number of MHC class I negative tumor cells in the absence of treatment. We assume that the natural death rate of MHC class I negative tumor cells is $\phantom {\dot {i}\!}k_{d10}\cdot C_{MHCI^{-}}$ and these cells are killed by chemotherapy agent OXP in tumor at a rate $\phantom {\dot {i}\!}k_{d2}\cdot C_{MHCI^{-}}\cdot \frac { \text {OXP}_{3}}{\text {OXP}_{3}+g_{15}}$. The difference in tumor volume caused by tumor re-challenge is described by the discrete Eq. (18).

  $$ {}\begin{aligned} \frac{d C_{MHCI^{-}}}{dt}&=-c_{7}\cdot \frac{\text{IFN}_{\gamma}}{\text{IFN}_{\gamma}+g_{11}}\cdot C_{MHCI^{-}}-k_{d10} \cdot C_{MHCI^{-}}\\ &+k_{p4} \cdot C_{MHCI^{-}} -r_{1} \cdot C^{2}_{MHCI^{-}}+2\cdot k_{p4} \cdot C_{MHCI^{+}}\\ &-k_{d2}\cdot C_{MHCI^{-}}\cdot \frac{ \text{OXP}_{3}}{\text{OXP}_{3}+g_{15}} \end{aligned} $$

  (15)

• We use the difference Eqs. (16) and (17) to reflect the abrupt changes in the concentrations of IL-12 and OXP caused by the therapies, respectively. The sudden change in the volume of MHC class I negative tumor cells due to tumor re-challenge is described by the difference Eq. (18).

  $$\begin{array}{*{20}l} {}\Delta \text{IL}(t)&=\text{IL}_{k1}, \;\;\;t=t_{k1},\;\;k1=1,2, \cdots, n_{1} \end{array} $$

  (16)
  $$\begin{array}{*{20}l} {}\Delta \text{OXP}_{2}(t)&=\{{\text{OXP}}_{2}\}_{k2}, \ t=t_{k2},\;\;k2=1,2,\cdots, n_{2} \end{array} $$

  (17)
  $$\begin{array}{*{20}l} {}\Delta C_{\text{MHCI-}}(t)&=C_{k3}, \;\;\;t=t_{k3},\;\;k3=1,2,\cdots, n_{3} \end{array} $$

  (18)

  where ΔIL(t)=IL(t⁺)−IL(t⁻) and ΔOXP₂(t)=OXP₂(t⁺)−OXP₂(t⁻) reflect the abrupt changes of IL-12 and oxaliplatin at administration times t_k1 and t_k2, while IL_k1 and {OXP₂}_k2 are the dosages of IL-12 and oxaliplatin at the administration times t_k1 and t_k2 with k1=1,2,⋯,n₁ and k2=1,2,⋯,n₂, respectively; ΔC_MHCI-(t)=C_MHCI-(t⁺)−C_MHCI-(t⁻) represents the sudden changes in tumor size due to tumor re-challenge with implantation size C_k3 mm³ at time t_k3 with k3=1,2,⋯,n₃.

A schematic diagram summarizing this three-compartment model is shown in Fig. 1. Model parameters and their meanings are listed in Table 1.

Non-negativity of solutions to the model

For any mathematical model that has biological implications, it is important to make sure that solutions with non-negative initial conditions remain non-negative all the time. For the model comprised of Eqs. 1) - 18), we can see that, with positive impulsive inputs IL_k1, {OXP₂}_k2, and C_k3 at impulsive moments t_k1, t_k2 and t_k3 in Eqs. (16) - (18), all solutions of the system of ODEs (1) - (15) with non-negative initial conditions will remain non-negative because $\frac {dx_{i}}{dt}\geq 0$ for x_i=0 and x_j≥0, where i,j=1,2,⋯,15 and i≠j.

Model calibration

Therapeutic use of IL-12 requires efficient methods to control the plasma concentration of this potent immuno-stimulatory cytokine in order to avoid toxicity [6]. It was determined in an MC38 syngeneic tumor model that a blood concentration of IL-12 less than 20 ng/ml has no anti-tumor effect, while concentrations greater than 700 ng/ml are associated with toxicity [43]. For comparison, the normal range of IL-12 in humans is around 7.5 pg/ml [44]. Gonzalez-Aparicio and colleagues designed a new induction protocol to keep IL-12 within this therapeutic range [6]. Once the liver of a group of C57BL/6 mice was transduced with the vector (typically 2.5* 10⁸ IU), a suboptimal amount of Mif (125 μg/kg) is administered for the first 2 days in order to prevent toxicity. The concentration of IL-12 is measured in serum 10 h after the first induction, and based on this information, a stepwise increase in Mif is scheduled to allow several cycles of sustained IL-12 expression in mice treated with the HC-Ad/RUmIL-12 vector (Fig. 2). Before we start the calibration of model parameters, we first quantified the IL-12 concentration as a function of time in days (Fig. 2a) and Mifepristone in μg/kg (Fig. 2b), respectively. Empirical functions were used to represent the IL-12 as a function of time and as a function of Mifepristone dose. These calibrated functions are shown in Fig. 2, where they are compared against the experimental data reported in [6]. Overall, the curves show a good match between experimental data and model predictions used to describe Mif-induced IL-12 treatment effects. Since the authors [6] verified the Mif-induction system is functional for more than 5 months with a slow decrease in the intensity of expression in each cycle, we used the same relationships for repeated Mif-induced treatment cycles in our model for simplicity.

Fig. 2

Quantified Mif-induced IL-12 expression. Simulated IL-12 expression as a function of time and Mif (the green curves) was calibrated to (mean + s.d.) experimental data reported in Fig. 1a and b in [6] in (a) and (b), respectively.The HC-Ad/RUmIL-12 vector was administered at 2.5* 10⁸ IU/mouse in C57BL/6 mice by intrahepatic injection. A set of 8 mice received an adjusted protocol (red circles, n =8) that consisted of 125 μg/kg Mifepristone days 1-2; 250 μg/kg days 3-5; 500 μg/kg days 5-7 and 1000 μg/kg days 9-11 in (a) and a single dose of Mifepristone (125; 250; 1000; 2000 or 4000 μg/kg) was administered intraperitoneally to different groups of animals (n =5) after 2 weeks in (b). The concentration of IL-12 in serum was determined 10 h after induction at the indicated days. Experimental data in error bars represent mean + s.d

We then calibrated model parameters in system described by Eqs. 1) - 15) using two sets of experimental data from [6]. The first sets of data are listed below:

• Total volume of MC38Luc1 tumor cells was calibrated against data shown in Figs. 2b, c, and 3 (IL-12 + OXP group).

  Fig. 3

  

  Comparison of model predictions with experimental measures of therapeutic response upon tumor re-challenge. a. Comparison of model predictions with experimental measures of tumor volume, IFN _γ and T_E3/ T_R of mice subjected to tumor re-challenge after one cycle of IL-12 and OXP treatment at day 57. The experimental data were acquired for a group of C57BL/6 mice with 5* 10⁵ MC38Luc1 cells inoculated in the liver on day 0 and subjected to one cycle of OXP (on day 9) and Mif-induced IL-12 (started on day 12 and continued 10 days) treatment. To check the immunological protection against cancer cells in treated animals, the cured mice had a tumor re-challenge of 10⁶ MC38Luc1 cells about one month after completion of previous treatment. Experimental measures of tumor volume, IFN _γ, and T_E3/ T_R (crosses, represent average of n = 16) from Figs. 2 - 5 in [6] were compared to the model predictions (blue curve) generated using a genetic algorithm. b - d. The experimental data were acquired for a group of C57BL/6 mice bearing hepatic tumors treated with the HC-Ad/RUmIL-12 vector and received two cycles of Mifepristone (Mif) induction preceded by OXP (5 mg/kg, intraperitoneally). Animals cured from their hepatic tumors were subjected to a subcutaneous challenge with the same tumor cells (MC38Luc1), and received a third cycle of IL-12 and OXP treatment starting on day 103. Experimental measures of tumor volume for individual mice (squares, triangles, and crosses) from Fig. 7 in [6] were compared to the model predictions (blue curve) generated using a genetic algorithm. Model predictions calibrated to tumor volume for responder, partial-responder, and non-responder mice treated with one cycle of combined therapy after tumor re-challenge are shown in panels b, c, d, respectively. Each graph displays a collection of 30 good fits of model predictions against experimental data. The solid blue curve provides the median model prediction of the 30 good fits, and the dashed purple and green curves indicate the 90% upper and lower boundaries in the model predictions of 30 good fits, respectively. Example parameter values of good fits in each panel are included in Table 2

  Table 2 Examples of calibrated parameter values against experimental data

• Concentration of Interferon gamma was obtained from Fig. 4b (IL-12 + OXP group).

  Fig. 4

  

  Violin plots of normalized tumor size changes with 30 good fits of parameter sets for responders, partial-responders, and non-responders on day 120. The sample set of parameter values for each group used in the plots are listed in Table 2

• The ratio between CD8 ⁺ T lymphocytes and T regulatory cells was obtained from Fig. 5b (experimental results for IL-12 + OXP group in tumor).

  Fig. 5

  

  Treatment strategies for partial-responders. The distribution in responses of 30 patients were sketched for each treatment strategy using 30 sets of good fits of calibrated parameters for partial responders (solid line represents median while dotted lines enclose 90% of the predictions). See a sample set of parameter values used in the plots in Table 2. a. Effects of increased interleukin-12 (IL-12) dose from 1 time (1X, control) to 3, 5, 10 times (3X, 5X, 10X, respectively). b. Effects of moderately increased OXP dose from 1 time (1X, control) to 2, 4, 6 times (2X, 4X, 6X, respectively). c. Effects of aggressively increased OXP dose from 1 time (1X, control) to 10, 100, 200 times (10X, 100X, 200X, respectively). d. Effects of increased number of treatment cycles from 3 cycles to 4, 5, and 6 cycles

The model was calibrated against these data to reflect the effects of the combination therapy with one treatment cycle after tumor cell liver implantation (5×10⁵ MC38Luc1) and the immunological protection against cancer cells treated with the combination therapy after a tumor re-challenge (10⁶ MC38Luc1 cells 35 days after the completion of the previous treatment). Calibration results, including the median (solid blue curve), 90th percentile upper (dashed purple curve), and 90th percentile lower responses (dashed green curve) of 30 good fits, are included in Fig. 3a), where $C_{{MHCI}^{-}}(0)=1 \,mm^{3}$ since t₀=0 is the day that 5×10⁵ cells /mouse MC38Luc1 tumor cells were inoculated in the liver of C57BL /6 mice [45, 46], T_N(0)=0.0714 cells per mm³ ($=\frac {100}{1.4*10^{3}}$ as we assume that the number of naïve CD8 ⁺ T cells in a mouse is 100 and the volume of the blood system of a mature mouse is 1.4∗10³ mm³), ${APC}_{2}(0)=214.2857\;(=\frac {3*10^{5}}{1.4*10^{3}})$ cells per mm³ according to [47]. Other initial values are zero: $\phantom {\dot {i}\!}T_{Ei}(0)=C_{MHCI^{+}}(0)=\text {OXP}_{i}(0)=\text {IFN}_{\gamma }(0)=\text {APC}_{1}(0)=\text {APC}_{3}(0)=T_{R}(0)=IL(0)=0$ for i=1,2,3. $\phantom {\dot {i}\!}\Delta C_{MHCI^{-}}(57)=2 \,mm^{3}$, ΔOXP₂(s)=5 mg/ kg with s=10, 34, 100; ΔIL(t) follows $f(t)=\frac {13.6127*(t-12)^{2}+0.8606*(t-12)+1}{0.313*(t-12)^{2}-0.6216*(t-12)+1}$ for 12≤t≤21 (see details in Fig. 2a). The parameter values used in the simulations are listed in Table 2 with biological meanings of each parameter listed in Table 1.

To show the long-term management of colorectal cancer using the combined therapy, tumor growth of a group of mice subjected to one cycle of treatment was calibrated to data from Fig. 7(D) in [6]. In the combined therapy, 5 mg/kg OXP on day 100 and 10-day IL-12 induction starting day 103, which follows the adjusted protocol as described in Fig. 2 in [6], were administered after a tumor re-challenge on day 75. This treatment occurred about two weeks after the mice survived two cycles of combined treatments with 10-day Mif-induced IL-12 (induction started on days 13 and 37) and OXP treatments (5 mg/kg on days 10 and 34), which in turn started on day 14. The experimental results were split into responders (n= 2, Fig. 3b), partial responders (n = 2, Fig. 3c) and non-responders (n = 3, Fig. 3d) groups. The mathematical model was calibrated separately to these different response groups. Figure 3b - d with the median (solid blue curve), 90th percentile upper (dashed purple curve), and 90th percentile lower responses (dashed green curve) of 30 good fits illustrate the results of our simulations compared with the corresponding experimental data. Here, we have $\Delta C_{MHCI^{-}}(75)=2 \,mm^{3}\phantom {\dot {i}\!}$, ΔOXP₂(9)=5 mg/kg; ΔIL(t) follows $IL(t)=\frac {13.6127*(t-a)^{2}+0.8606*(t-a)+1}{0.313*(t-a)^{2}-0.6216*(t-a)+1}$ for a≤t≤a+9 with a=13, 37, 103. A sample set of parameter values for each of the response groups of mice used in the simulations are listed in Table 2.

Difference in treatment efficacy: non-responders versus responders While both the responders and non-responders survived two cycles of combination therapy treatment before tumor rechallenge and then underwent the third cycle following the tumor re-challenge on day 75, the simulations in Fig. 3b - d suggest that the non-responders show near zero concentration of IFN_γ and near zero ratio of T effectors to regulatory T cells in the tumor all the time comparing to a stable concentration of IFN_γ and ratio of T effectors to regulatory T cells in the tumor (at least 10³ cells per mm³ after the combination therapy treatment) in responders and partial responders. The simulations also indicate that whether the immune system can maintain a high ratio of T effectors to regulatory T cells in the tumor as well as generating a moderate but stable concentration of IFN_γ might be crucial to control tumor growth. This finding is consistent with results from experimental studies [48–50]. In [48] and [49], Random forest (RF) machine learning analysis of tumor infiltrating lymphocytes (TILs) demonstrated that the top predictor of resistance, as measured by variable importance scores and selection, was the CD8+CD44+ to Treg (CD8/Treg) ratio. The authors in [50] indicated that the CD8/Treg ratio in resistant tumors failed to increase after RT + anti-CTLA4 as it did in sensitive tumors because CD8 + CD44 + T cells did not significantly expand despite reduction in Tregs.

Genetic algorithms

In this research, genetic algorithms are applied to find good parameter sets in a high-dimensional parameter space. To apply genetic algorithms, chromosomes are encoded using parameter sets. For each generation, the ODE set is solved for each individual and the fitness function value, or variance, is calculated. For each generation, the chromosomes are sorted using fitness function values. The top half is passed to the next generation. The best individual in each generation is free from mutation. Pairing is done in the sorted sequence then crossover is conducted to generate the next generation. Adaptive versions of GAs were designed and implemented. Besides variable initial values and parameter ranges, step size and GA control parameters can also be specified. Adaptive GAs are able to change the number of crossover points, use random crossover points, use random mutation rates, and dynamically evolve based on the simulation progress, which is evaluated by the best fitness function value achieved at the moment. Therefore, the search process can move more efficiently in the large parameter space and have smaller chances to get trapped in local minimums.

Goodness of fit

Descriptions of model parameters are shown in Table 1. A couple of sample sets of estimated values of parameters obtained through fitting predictions of the impulsive ODE model 1) - 18) to data from a group of experiments in [6] are listed in Table 2. Figure 3a illustrates the comparison between model solutions and experimental measurements on tumor control and immunological protection against cancer cells in animals treated with IL-12 plus OXP. Experimental results for long-term management of colorectal cancer by observing cooperation of IL-12 and OXP for the control of experimental relapses in distant locations are compared against model predictions in Fig. 3b, c, and d for responders, partial-responders and non-responders, respectively. Trajectories of tumor growth, IFN_γ and ratio of T_E3 to T_R arising from the model are extremely close to corresponding data from the experiments. For each calibration, an excess of data points (57, 65, 58, 84 for Fig. 3a - d, respectively) relative to the number of parameters (48) suggests that the model is identifiable in theory.

Model stability analysis

In this section, we discuss local stability of equilibria of the model using linearized system evaluated at these points. We found that system comprised of equations (1) - (15) has a tumor-free equilibrium $\vec {X}_{0}$, a second tumor free equilibrium $\vec {X}_{1}$ when APC₂ growth rate is larger than the rate constant for APC₂ flowing from blood to tumor (i.e., r₂>b₂₃), and a high tumor equilibrium $\vec {X}_{2}$ when proliferation rate of tumor cells is higher than natural death rate of tumor cells (i.e., k_p4>k_d10).

By setting the right hand sides of the equation system (1) - (15) to zero and solving the equations simultaneously, we obtain

$${}\begin{aligned} \vec{X}_{0}&=\left(\;T_{N},{\vphantom{C_{MHCI}^{-}}}\;T_{E1},\;\text{APC}_{1},\,\text{OXP}_{1},\,T_{E2},\;\text{APC}_{2}, \,\text{OXP}_{2},T_{E3},\right.\\& \quad\left.\,\text{IFN}_{\gamma},\text{APC}_{3},\,\text{IL},\,\text{OXP}_{3}, \,T_{R},\ C_{MHCI}^{+}, \;C_{MHCI}^{-}\,\right)^{T}\\ &=\left(\;\frac{c_{1}}{k_{d1}},\,0,\,0,\,0,\,0,\,0,\,0,\,0,\,0,\,0,\,0,\,0,\ \frac{c_{6} }{k_{d9} },\,0,\,0\;\right)^{T}, \end{aligned} $$

$${}\begin{aligned} \vec{X}_{1}&=\left({\vphantom{C_{MHCI}^{-}}}T_{N},\, T_{E1},\;\text{APC}_{1},\,\text{OXP}_{1},\,T_{E2},\;\text{APC}_{2}, \,\text{OXP}_{2},\,T_{E3},\right.\\ &\quad\left. \text{IFN}_{\gamma},\, \text{APC}_{3},\,\text{IL},\,\; \text{OXP}_{3}, \,T_{R}, \, C_{MHCI}^{+}, \;C_{MHCI}^{-}\right)^{T}\\ &=\left(\!T_{N},\,T_{E_{1}},\, \frac{b_{23} \text{Vol}_{b} K(r_{2} -b_{23})}{r_{2} k_{d3} \text{Vol}_{\ln} },\,0,\, T_{E2},\frac{K(r_{2} -b_{23})}{r_{2} },\right.\\ &\quad\left.0,\,T_{E3},\, \text{IFN}_{\gamma},\,\text{APC}_{3},\, \text{IL},\,0,\,\frac{c_{6} }{k_{d9} },\,0,\,0\right)^{T}, \end{aligned} $$

where $T_{E_{1}}$ satisfies the following polynomial equation

$${}\begin{aligned} \frac{a_{12} k_{d5} }{a_{21} +k_{d5}} T_{E1}^{3} &+\frac{\text{Vol}_{b}}{\text{Vol}_{\ln} } (k_{d1} T_{N} -c_{1})T_{E1}^{2}\\ &+\alpha \left(\frac{a_{12} k_{d5} }{a_{21} +k_{d5}} -k_{p1} \frac{APC_{1} }{APC_{1} +g_{2}} \right)T_{E1}\\ & +\alpha \frac{\text{Vol}_{b} }{\text{Vol}_{\ln} } (k_{d1} T_{N} -c_{1})=0, \end{aligned} $$

(Based on the Descartes’ rule, there is only one positive solution from the equation), $T_{E2}=\frac {a_{12} \text {Vol}_{\ln } T_{E1}}{(k_{d5} +a_{21})\text {Vol}_{b} }$, T_E3 satisfies the following quadratic equation

$$k_{d6} V_{i} T_{E3}^{2} +k_{d6} \varepsilon T_{E3} -a_{23} \text{Vol}_{b} T_{E2} =0$$

which has only one positive solution $T_{E3}=\frac {-\varepsilon k_{d6} +\sqrt {(\varepsilon k_{d6})^{2} +4a_{23} V_{i} k_{d6} \text {Vol}_{b} T_{E2}} }{2V_{i} k_{d6} }$, $\text {IFN}_{\gamma }=\frac {c_{4} T_{E3} IL}{k_{d7} (IL+g_{7})(T_{R} +g_{8})}$, $\text {APC}_{3}=\frac {b_{23} {APC}_{2} \text {Vol}_{b} }{(k_{d3} +b_{31})\text {Vol}_{t} }$, $\text {IL}=\frac {c_{5} {APC}_{3} }{k_{d8} }$, and

$${}\begin{aligned} \vec{X}_{2}&=\left({\vphantom{C_{MHCI}^{-}}}T_{N},\;T_{E1},\;\text{APC}_{1},\,\text{OXP}_{1},\,T_{E2},\;\text{APC}_{2}, \,\text{OXP}_{2},\,T_{E3},\right.\\ &\qquad\left.\text{IFN}_{\gamma}, \,\text{APC}_{3},\,\text{IL},\,\text{OXP}_{3}, \,T_{R},\, C_{MHCI}^{+}, \;C_{MHCI}^{-}\right)^{T}\\ &=\left(\frac{c_{1}}{k_{d1}},\,0,\,0,\,0,\,0,\,0,\,0,\,0,\,0,\,0,\,0,\,0,\,\frac{c_{6} }{k_{d9} },\,0,\right.\\ &\qquad\left.\frac{k_{p4}-k_{d10}}{r_{1}}\right)^{T}. \end{aligned} $$

By simple calculation, thirteen of the fifteen eigenvalues of the Jacobian matrix of linearized system at the equilibrium $\vec {X}_{0}$ are given by {−k_d1, −k_d3, −b₂₃+r₂, −k_d4, −k_d4, −k_d4−c₂₁−c₂₃, −k_d6, −k_d7, −b₃₁−k_d3, −k_d8, −k_d9, −k_d10−k_p4, and −k_d10+k_p4} with the rest of the eigenvalues satisfying the quadratic equation λ²+(k_d5+a₁₂+a₂₁+a₂₃)λ+a₁₂(k_d5+a₂₃)=0. It is easy to see that both eigenvalues are in the left half of the complex plane for a wide range of the parameters. Thus the first tumor free equilibrium $\vec {X}_{0}$ is stable if both r₂<b₂₃ (i.e., APC₂ growth rate is smaller than the rate constant for APC₂ flowing from blood to tumor) and k_p4<k_d10 (i.e., tumor proliferation rate is less than tumor natural death rate), otherwise it is unstable.

The second tumor free equilibrium $\vec {X}_{1}$ exists when r₂>b₂₃ (i.e., APC₂ growth rate is larger than the rate constant for APC₂ flowing from blood to tumor). Similar to the previous case, most of the eigenvalues of the Jacobian matrix at $\vec {X}_{1}$ are in the left-half of the complex plane. It is easy to see that twelve of the fifteen eigenvalues are negative. With respect to the remaining three, one is b₂₃−r₂, and the other two satisfies the following quadratic equation

$${}\begin{aligned} \lambda^{2} &+\!\left(\frac{c_{7} {IFN}_{\gamma} }{IFN_{\gamma} +g_{11}} +\frac{k_{d11} T_{E3} }{\varepsilon (T_{R} +g_{13})} +2k_{d10} \right)\lambda\\ & +\!\left(\!-k_{d10} -k_{p4} -\frac{k_{d11} T_{E3} }{\varepsilon (T_{R} +g_{13})} \right)\\ &\times\!\left(\!-k_{d10} +k_{p4} -\frac{c_{7} {IFN}_{\gamma} }{IFN_{\gamma} +g_{11}}\! \right)\!-2k_{p4} \frac{c_{7} {IFN}_{\gamma} }{IFN_{\gamma} +g_{11}} \!=0. \end{aligned} $$

It is found that both eigenvalues are in the left-half of the complex plane if k_p4<k_d10. Hence this tumor-free equilibrium point is locally stable if r₂>b₂₃ and k_p4<k_d10; otherwise it is unstable.

The high tumor equilibrium $\vec {X}_{2}$ exists when k_p4>k_d10. Twelve of the fifteen eigenvalues of the Jacobian matrix are found as {−k_d1,−k_d3,−k_d4,−k_d4,−b₂₃+r₂,−k_d4−c₂₁−c₂₃,−k_d7,−b₃₁−k_d3,−k_d8,−k_d9,−k_d10−k_p4, and k_d10−k_p4}.The other three eigenvalues are the roots of the polynomial p(λ)=λ³+α₁λ²+α₂λ+α₃, where

$${}\begin{aligned} \alpha_{1} &= a_{12} +a_{21} +a_{23} +k_{d5} +k_{d6} +\frac{a_{32} (k_{p4} -k_{d10})}{r_{1} \varepsilon +k_{p4} -k_{d10} }, \\ \alpha_{2} &= k_{d6} (k_{d5} +a_{12} +a_{21} +a_{23})+a_{12} (k_{d5} +a_{23})\\ & +(k_{d5} +a_{12} +a_{21})\frac{a_{32} (k_{p4} -k_{d10})}{r_{1} \varepsilon +k_{p4} -k_{d10} }, \\ \alpha_{3} &= a_{12} k_{d5} k_{d6} +a_{12} a_{23} k_{d6} +k_{d5} a_{12} \frac{a_{32} (k_{p4} -k_{d10})}{r_{1} \varepsilon +k_{p4} -k_{d10} }. \end{aligned} $$

Based on the list of the first twelve eigenvalues, the high-tumor equilibrium is unstable if either r₂>b₂₃ or k_d10>k_p4 is satisfied. Suppose r₂<b₂₃ and k_d10<k_p4, it is easy to see that α_i>0, i=1,2,3. According to the Routh-Hurwitz criterion, p(λ) is a Hurwitz polynomial if and only if α₁α₂−α₃>0, which is indeed the case after simplifying the expression. Therefore, the high-tumor equilibrium is stable if r₂<b₂₃ (i.e., APC₂ growth rate is smaller than the rate constant for APC₂ flowing from blood to tumor) and k_d10<k_p4 (i.e., tumor proliferation rate is larger than tumor natural death rate), otherwise it is unstable.

The stability conditions for tumor free equilibrium $\vec {X}_{0}$ indicate that small tumors may not grow into a threatening size when tumor proliferation rate is less than tumor natural death rate (k_p4<k_d10) without any treatment (r₂<b₂₃ with all the APC_i, OXP_i, IL and T_Ei, i=1,2,3 in $\vec {X}_{0}$ be to zero). Under the combined IL-12 and OXP treatment, MHC class I positive tumor cells will ultimately be eliminated. Depending on effects of the combined treatment (reflected by the remaining level of APC_i and T_Ei, i=1,2,3 and whether r₂>b₂₃), MHCI negative tumor cells ($\mathrm {C}_{\text {MHCI}^{-}}$) will either be completely removed in which case solutions of this dynamic approach the second tumor-free equilibrium $\vec {X}_{1}$ or the MHC class I negative tumor cells eventually approach the carrying capacity. This can occur when MHC class I positive tumor cells are all killed by tumor infiltrating lymphocytes, which results in the exhaustion of effector CD8 ⁺ T cells and cytokines while naïve T cells and MHC class I negative tumor cells remaining at constant levels.

Sensitivity of parameters

To test the impact of how the change of a certain parameter value (e.g. α) would affect tumor growth pattern for responders, partial-responders, and non-responders, normalized differences of tumor sizes $y_{i}=\frac {|\hat {y_{i}}-\bar {y_{i}}|}{\bar {y_{i}}}, \,(i=1,2,\cdots,30)$ on day 120 (8 days post the second treatment cycle of the combination therapy after tumor re-challenge) were used to draw the violin plots in Fig. 4 for each of the 48 parameters for all three patient groups, where $\hat {y_{i}}$ is the predicted tumor size in mm³ using 0.1×α_i and other parameters in the ith calibrated parameter set and $\bar {y_{i}}$ is the predicted tumor size in mm³ using α_i and other parameters in the i_th calibrated parameter set and α_i is the calibrated value for parameter α in the i_th calibrated parameter set.

In general, we found that changing the value of each of the 48 parameters barely affected the tumor growth for non-responders. In addition, there are 10 (out of 48) parameters whose value changes greatly affect tumor size of responders but not the size of non-responders and partial-responders. These potentially OXP and IL-12 treatment important parameters include c₂₃ (OXP flow rate from blood to tumor), K (APC carrying capacity), c₄ - c₆ (IFN_γ, IL-12, and T_R production rate constants, respectively), g₁₀ - g₁₄ (IL-12, IFN_γ, OXP₃, T_R, and ${\mathrm {C}}_{MHCI}^{+}$ killing by OXP₃ saturation rate constants, respectively). In addition, changes in the value of the following 7 parameters cause from zero for non-responders to increasing changes in normalized tumor size from partial-responders to responders: T cell flow rates from blood to lymph node and from tumor to blood, a₂₁ and a₃₂, respectively; APC flow rates from tumor to lymph node b₃₁; production rate constant of naive T cells c₁; transfer rate constant of naive T cell to T effector cell in lymph node c₂; IL-12 natural death rate constant k_d8, and APC growth rate constant r₂. We also note that no change in tumor size for all three mice group (non-responders, partial-responders, and responders) results from the value changes of following 5 parameters: ${\mathrm {C}}_{MHCI}^{-}$ killing by OXP₃ saturation rate constant g₁₅, natural death rate constant of naive T cells k_d1, killing rate constant of T effectors or tumor cells by OXP k_d2, APCs natural death rate constant k_d3, and natural decay rate constant of OXP k_d4 (see Fig. 4).

Model simulations

In order to investigate potential ways to improve treatment regimes for partial-responders and non-responders, we simulated the following alternative treatment scenarios: changing the dose and frequency of chemotherapy drug OXP administration, changing the strength of Mif-induced IL-12 expression, and changing the number of combined treatment cycles.

Partial-responders

We note from Fig. 5 that increased number of treatment cycles in the IL-12 and OXP combination therapy does not seem to improve tumor control in the first 8 months post treatment while increased dose of OXP alone would achieve better tumor control and enhanced strength of IL-12 expression alone would slightly reduce tumor size more rapidly after the tumor reaches its maximum size.

Non-responders

From Fig. 6, we found that neither increased strength of IL-12 expression nor moderately increased OXP dose alone in the IL-12 and oxaliplatin (OXP) combination therapy seems to improve tumor control for the median, 90th percentile lower and 90th percentile upper responses for the 30 non-responder patients. Meanwhile, aggressively increased OXP dose (for instance, 10 + times) in the combination therapy shows reduced tumor size and delayed time of tumor reaching its carrying capacity only for the 90th percentile lower responses for the 30 patients. The reduction of tumor size slows greatly when OXP dose is increased to more than 100 times. In addition, increased number of treatment cycles in the IL-12 and OXP combination therapy reduced tumor size and delayed the time of tumor reaching its carrying capacity only for the 90th percentile lower responses for the 30 patients.

Fig. 6

Treatment strategies for non-responders. The distribution in responses of 30 patients were sketched for each treatment strategy using 30 sets of good fits of calibrated parameters for non-responders (solid line represents median while dotted lines enclose 90% of the predictions). A sample set of parameter values used in the plots is listed in Table 2. a. Effects of increased interleukin-12 (IL-12) dose from 1 time (1X, control) to 3, 5, 10 times (3X, 5X, 10X, respectively). b. Effects of moderately increased OXP dose from 1 time (1X, control) to 2, 4, 6 times (2X, 4X, 6X, respectively). c. Effects of aggressively increased OXP dose from 1 time (1X, control) to 10, 100, 200 times (10X, 100X, 200X, respectively). d. Effects of increased number of treatment cycles from 3 cycles to 4, 5, and 6 cycles

Discussion

Developing mathematical models that represent known features of the biological system and that are calibrated to experimental studies can help improve understanding of the underlying biology targeted by drugs and enables exploring therapeutic scenarios that may be difficult or costly to test experimentally. In this paper, we developed a three-compartment mechanistic mathematical model to describe the clonal expansion of CD8 + T cells in a mouse model of metastatic colorectal cancer in response to a combined therapy of IL-12 plus the chemotherapy drug Oxaliplatin. Based on the collective knowledge of the underlying biology, the model represents the primary CD8 + T cell response under a boosting effect of IL-12 and OXP and the subsequent impact on the growth of a tumor based on the syngeneic MC38Luc1 mouse model for metastatic colorectal cancer, where the observed response was characterized by three phenotypes: responders, partial responders, and non-responders. Model parameters were calibrated against published experimental data that describes the primary response for these three phenotypes. The sensitivity analysis of parameters helped explain the differences in calibrated values of parameters between non-responders, partial-responders, and responders. To reduce the dependence of our model predictions on any single calibrated set of parameter values, we generated an ensemble of 30 parameter sets for each phenotype that provided a similar good fit against the experimental data and show the distribution in phenotypic responses for those virtual cohorts. Using the corresponding ensemble of model predictions for non-responders, numerical simulation of multiple OXP and IL-12 combination therapy suggest that aggressively increasing the dose (between 10 and 100 times of the control) of OXP will improve tumor control results while increasing the number of treatment cycles of the combined therapy can decrease the tumor size as well. We also found that only increasing the OXP dose in the combination therapy can dramatically decrease the tumor size for partial responders.

Conclusion

Overall, these results illustrate how mechanistic models can be used to predict tumor growth response to antigen-specific immuno-chemotherapies and screen in silico for optimal therapeutic dosage and timing in treating patients with metastatic colorectal cancer.