Skip to main content

Advertisement

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

  • 136 Accesses

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 [35]. 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 [813]. 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 [1317]. 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 [2025]). 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 [2629]. 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
figure1

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 (TN) are activated and become CD8 + T effectors (TE1) when they encounter tumor antigen presented by the antigen presenting cells (APC1) in the lymph node. Once activated, effector CD8 + T cells circulate within the blood (TE2) and enter tumor microenvironment (TE3) 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 (TR) cells’ proliferation and immunosuppressive action on effector CD8 + T cells. In addition, the chemotherapy drug Oxaliplatin in the lymph node and tumor (OXPi 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 (TN, units: cells per mm3). We assume that naïve CD8 + T cells are produced at a constant rate c1 from thymus and die naturally at a rate kd1·TN [31]. Naïve T cells are recruited and activated by tumor antigens presented by APC1 (antigen-presenting cells in lymph node) at a rate \(c_{2}\cdot T_{N}\cdot \frac {\text {APC}_{1}}{\text {APC}_{1}+g_{1}}\) [3234].

    $$ {}\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 (TE1, units: cells per mm3). 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 Volb=1.4103 mm3 is the volume of the blood compartment [35] and Volln=0.25 mm3 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 TE1, 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 TE1 [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 TE2 is the concentration of T effectors in blood, and a12·TE1 is the efflux rate. We assume that TE1 cells are killed by chemotherapy agent OXP1 (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 (APC1, units: cells per mm3). We assume that APCs in the lymph node have a natural death rate of kd3·APC1 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 APC3 is the concentration of APCs in tumor, \({Vol}_{t}= \frac {\epsilon +C(t)}{1-V_{i}\cdot T_{E3}}\) (since Volt=ε+C(t)+Vi·TE3·Volt, where TE3 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 mm3. 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 (Vi) 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 (OXP1, units: mg per kg). We assume that OXP decays naturally at a rate kd4·OXP1 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 OXP2 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 (TE2, units: cells per mm3). We assume the effector CD8 + T cells die naturally in blood at a rate kd5·TE2. 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 a21·TE2. 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 TE3 is the concentration of T effectors in tumor and the efflux rate of CD8 + T effectors from blood to tumor is a23·TE2.

    $$\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 (APC2, units: cells per mm3). 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 r2 is the growth rate constant and K is the carrying capacity. We assume a b23·APC2 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 (OXP2, units: mg per kg). We assume that OXP decays naturally at a rate kd4·OXP2 and the efflux rates of OXP from blood to lymph node and from blood to tumor are c21·OXP2 and c23·OXP2, 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 (TE3, units: cells per mm3). 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 TR 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 kd6·TE3. T effector cells are assumed to be killed by chemotherapy agent OXP (OXP3) 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 mm3). 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 kd7.

    $$\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 (APC3, units: cells per mm3). We assume that APCs in the tumor microenvironment have a natural death rate of kd3·APC3, 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 b31·APC3.

    $$ {}\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 c5·APC3 and decays naturally at a rate of kd8·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 (OXP3, units: mg per kg). We assume that OXP decays naturally at a rate kd4·OXP3 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 (TR, units: cells per mm3). Regulatory T cells are produced at a constant rate c6 from thymus and die naturally at a rate kd9·TR. 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: mm3). 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: mm3). 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 ΔOXP2(t)=OXP2(t+)−OXP2(t) reflect the abrupt changes of IL-12 and oxaliplatin at administration times tk1 and tk2, while ILk1 and {OXP2}k2 are the dosages of IL-12 and oxaliplatin at the administration times tk1 and tk2 with k1=1,2,,n1 and k2=1,2,,n2, respectively; ΔCMHCI-(t)=CMHCI-(t+)−CMHCI-(t) represents the sudden changes in tumor size due to tumor re-challenge with implantation size Ck3 mm3 at time tk3 with k3=1,2,,n3.

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 ILk1, {OXP2}k2, and Ck3 at impulsive moments tk1, tk2 and tk3 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 xi=0 and xj≥0, where i,j=1,2,,15 and ij.

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* 108 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
figure2

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* 108 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
    figure3

    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 TE3/ TR 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* 105 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 106 MC38Luc1 cells about one month after completion of previous treatment. Experimental measures of tumor volume, IFN γ, and TE3/ TR (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
    figure4

    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
    figure5

    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×105 MC38Luc1) and the immunological protection against cancer cells treated with the combination therapy after a tumor re-challenge (106 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 t0=0 is the day that 5×105 cells /mouse MC38Luc1 tumor cells were inoculated in the liver of C57BL /6 mice [45, 46], TN(0)=0.0714 cells per mm3 (\(=\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.4103 mm3), \({APC}_{2}(0)=214.2857\;(=\frac {3*10^{5}}{1.4*10^{3}})\) cells per mm3 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}\), ΔOXP2(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}\!}\), ΔOXP2(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 ata+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 103 cells per mm3 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 [4850]. 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 TE3 to TR 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 APC2 growth rate is larger than the rate constant for APC2 flowing from blood to tumor (i.e., r2>b23), 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., kp4>kd10).

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} }\), TE3 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 {−kd1, −kd3, −b23+r2, −kd4, −kd4, −kd4c21c23, −kd6, −kd7, −b31kd3, −kd8, −kd9, −kd10kp4, and −kd10+kp4} with the rest of the eigenvalues satisfying the quadratic equation λ2+(kd5+a12+a21+a23)λ+a12(kd5+a23)=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 r2<b23 (i.e., APC2 growth rate is smaller than the rate constant for APC2 flowing from blood to tumor) and kp4<kd10 (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 r2>b23 (i.e., APC2 growth rate is larger than the rate constant for APC2 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 b23r2, 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 kp4<kd10. Hence this tumor-free equilibrium point is locally stable if r2>b23 and kp4<kd10; otherwise it is unstable.

The high tumor equilibrium \(\vec {X}_{2}\) exists when kp4>kd10. Twelve of the fifteen eigenvalues of the Jacobian matrix are found as {−kd1,−kd3,−kd4,−kd4,−b23+r2,−kd4c21c23,−kd7,−b31kd3,−kd8,−kd9,−kd10kp4, and kd10kp4}.The other three eigenvalues are the roots of the polynomial p(λ)=λ3+α1λ2+α2λ+α3, 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 r2>b23 or kd10>kp4 is satisfied. Suppose r2<b23 and kd10<kp4, 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 α1α2α3>0, which is indeed the case after simplifying the expression. Therefore, the high-tumor equilibrium is stable if r2<b23 (i.e., APC2 growth rate is smaller than the rate constant for APC2 flowing from blood to tumor) and kd10<kp4 (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 (kp4<kd10) without any treatment (r2<b23 with all the APCi, OXPi, IL and TEi, 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 APCi and TEi, i=1,2,3 and whether r2>b23), 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 mm3 using 0.1×αi and other parameters in the ith calibrated parameter set and \(\bar {y_{i}}\) is the predicted tumor size in mm3 using αi and other parameters in the ith calibrated parameter set and αi is the calibrated value for parameter α in the ith 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 c23 (OXP flow rate from blood to tumor), K (APC carrying capacity), c4 - c6 (IFNγ, IL-12, and TR production rate constants, respectively), g10 - g14 (IL-12, IFNγ, OXP3, TR, and \({\mathrm {C}}_{MHCI}^{+}\) killing by OXP3 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, a21 and a32, respectively; APC flow rates from tumor to lymph node b31; production rate constant of naive T cells c1; transfer rate constant of naive T cell to T effector cell in lymph node c2; IL-12 natural death rate constant kd8, and APC growth rate constant r2. 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 OXP3 saturation rate constant g15, natural death rate constant of naive T cells kd1, killing rate constant of T effectors or tumor cells by OXP kd2, APCs natural death rate constant kd3, and natural decay rate constant of OXP kd4 (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
figure6

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.

Availability of data and materials

All data was obtained from the originally cited publications. No new data was generated as part of the research for this manuscript.

Abbreviations

APC:

Antigen presenting cells

GA:

Genetic algorithms

IL-12:

Interleukin-12

MHC:

Major histocompatibility complex

NK:

Natural killer cells

ODE:

Ordinary differential equations

OXP:

Oxaliplatin

RF:

Random forest

TIL:

Tumor infiltrating lymphocyte

References

  1. 1

    Cancer.org: Colon and Rectum Cancer. http://www.cancer.org/cancer/colonandrectumcancer/detailedguide/. Accessed 1 July 2019.

  2. 2

    Martinet O, Ermekova V, Qiao JQ, Sauter B, Mandeli J, Chen L, Chen SH. Immunomodulatory gene therapy with interleukin 12 and 4-1bb ligand: long-term remission of liver metastases in a mouse model. J Natl Cancer Inst. 2000; 92(11):931–6.

  3. 3

    Hernandez-Alcoceba R, Berraondo P. Immunochemotherapy against colon cancer by gene transfer of interleukin-12 in combination with oxaliplatin. Oncoimmunology. 2012; 1(1):97–99.

  4. 4

    Colombo MP, Trinchieri G. Interleukin 12 in antitumor immunity and immunotherapy. Cytokine Growth Factor Rev. 2002; 13(2):155–68.

  5. 5

    Klinke DJ. A multiscale systems perspective on cancer, immunotherapy, and interleukin 12. Mol Cancer. 2010; 9(242):18.

  6. 6

    Gonzalez-Aparicio M, Alzuguren P, Mauleon I, Medina-Echeverz J, Hervas-Stubbs S, Mancheno U, Berraondo P, Crettaz J, Gonzalez-Aseguinolaza G, Prieto J, Hernandez-Alcoceba R. Oxaliplatin in combination with liver specific expression of interleukin 12 reduces the immunosuppressive microenvironment of tumours and eradicates metastatic colorectal cancer in mice. Gut. 2011; 60(3):341–9.

  7. 7

    Vacchelli E, Prada N, Kepp O, Galluzzi L. Current trends of anticancer immunochemotherapy. Oncoimmunology. 2013; 2(6):3.

  8. 8

    Kosinsky Y, Dovedi SJ, Peskov V. KndVoronova, Chu L, Tomkinson H, Al-Huniti N, Stanski DR, G H. Radiation and pd-(l)1 treatment combinations: immune response and dose optimization via a predictive systems model. J Immunother Cancer. 2018; 6(1):17.

  9. 9

    Wodarz D, Komarova NL. Computational Biology of Cancer: Lecture Notes and Mathematical Modeling. New Jersey: World Sci; 2005.

  10. 10

    Klinke DJ, Wang Q. Understanding immunology via engineering design: the role of mathematical prototyping. Comput Math Methods Med. 2012; 9:9.

  11. 11

    Klinke DJ. In silico model-based inference: A contemporary approach for hypothesis testing in network biology. Biotechnol Prog. 2014; 30:1247–61.

  12. 12

    Kim PS, Crivelli JJ, Choi IK, Yun CO, Wares JR. Quantitative impact of immunomodulation versus oncolysis with cytokine expressing virus therapeutics. Math Biosci Eng. 2015; 12(4):841–58.

  13. 13

    Wang Q, Klinke DJ, Wang Z. Cd8 t cell response to adenovirus vaccination and subsequent suppression of tumor growth: modeling, simulation and analysis. BMC Syst Biol. 2015; 9(27):19.

  14. 14

    Eftimie R, Dushoff J, Bridle BW, Bramson JL, Earn DJ. Multistability and multiinstability phenomena in a mathematical model of tumor immune virus interactions. Bull Math Biol. 2011; 73(12):2932–61.

  15. 15

    Gadkar KG, Shoda LK, Kreuwel HT, Ramanujan S, Zheng Y, Whiting CC, Young DL. Dosing and timing effects of anticd40l therapy: predictions from a mathematical model of type 1 diabetes. Ann New York Acad Sci. 2007; 1103:63–68.

  16. 16

    Zheng Y, Bresson D, Fradkin M, Chan JR, von Herrath M, C W. Biosimulations predict optimal oral insulin anti-cd3 and oral insulin/exendin-4 combination treatment regimens for the reversal of diabetes in the non-obese diabetic (nod) mouse. Clin Immunol. 2008; 127:101–2.

  17. 17

    Pappalardo F, Martinez-Forero I, Pennisi M, Palazon A, Melero I, Motta S. Simb16: modeling induced immune system response against b16-melanoma. PLoS ONE. 2011; 6(10):12.

  18. 18

    Dovedi SJ, Adlard AL, Lipowska-Bhalla G, McKenna C, et al. Acquired resistance to fractionated radiotherapy can be overcome by concurrent pd-l1 blockade. Cancer Res. 2014; 74(19):5458–68.

  19. 19

    Dovedi SJ, Cheadle EJ, Popple A, Poon E, Morrow M, Stewart R, et al. Fractionated radiation therapy stimulates anti-tumor immunity mediated by both resident and infiltrating polyclonal t-cell populations when combined with pd1 blockade. Clin Cancer Res. 2017; 23(18):5514–26.

  20. 20

    Kirschner D, Panetta JC. Modeling immunotherapy of the tumor immune interaction. J Math Biol. 1998; 37(3):235–52.

  21. 21

    Kim PS, Lee PP. Modeling protective antitumor immunity via preventative cancer vaccines using a hybrid agent based and delay differential equation approach. PLoS Comput Biol. 2012; 8(10):16.

  22. 22

    Eftimie R, Bramson JL, Earn DJ. Interactions between the immune system and cancer: a brief review of nonspatial mathematical models. Bull Math Biol. 2011; 73(1):2–32.

  23. 23

    Pappalardo F, Palladini A, Pennisi M, F C, Motta S. Mathematical and computational models in tumor immunology. Math Model Nat Phenom. 2012; 7(1):25.

  24. 24

    Tsygvintsev A, Marino S, Kirschner DE. Mathematical Methods and Models in Biomedicine Lecture Notes on Mathematical Modelling in the Life Sciences. New York: Springer; 2013.

  25. 25

    Eladdadi A, P K, Mallet D. Mathematical Models of Tumor Immune System Dynamics. Switzerland: Springer; 2014.

  26. 26

    de Pillis LG, Gu W, Radunskaya AE. Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations. J Theor Biol. 2006; 238(4):841–62.

  27. 27

    de Pillis LG, Fister KR, Gu W, Collins C, M D, D G, Moore J, Preskill B. Seeking bangbang solutions of mixed immunochemotherapy of tumor. Electonic J Differ Equ. 2007; 2007(171):1–24.

  28. 28

    de Pillis LG, Eladdadi A, Radunskaya AE. Modeling cancer immune responses to therapy. J Pharmacokinet Pharmacodyn. 2014; 41(5):461–78.

  29. 29

    Mahasa KJ, Ouifki R, Eladdadi A, de Pillis LG. Mathematical model of tumor immune surveillance. J Theor Biol. 2016; 404(171):312–30.

  30. 30

    Wang Z, Wang Q. Numerical simulation of a tumor growth dynamics model using particle swarm optimization. J Comput Sci Syst Biol. 2015; 9(1):1–5.

  31. 31

    Jenkins MK, Chu HH, McLachlan JB, Moon JJ. On the composition of the preimmune repertoire of t cells specific for peptide-major histocompatibility complex ligands. Ann Rev Immunol. 2010; 28:275–94.

  32. 32

    Miller MJ, Wei SH, Parker I, Cahalan MD. Two-photon imaging of lymphocyte motility and antigen response in intact lymph node. Science. 2002; 296(5574):1869–73.

  33. 33

    Krummel M, Bartumeus F, A G. T-cell migration, search strategies and mechanisms. Nat Rev Immunol. 2016; 16(3):193–201.

  34. 34

    Tang M, Diao J, Cattral MS. Molecular mechanisms involved in dendritic cell dysfunction in cancer. Cell Mol Life Sci. 2017; 74(5):761–76.

  35. 35

    C, 57BL 6 Inbred Mouse. https://www.jax.org/news-and-insights/2005/october/how-much-blood-can-i-take-from-a-mouse-without-endangering-its-health. Accessed 30 Dec 2019.

  36. 36

    Linderman JJ, Riggs T, Pande M, Miller M, Marino S, Kirschner DE. Characterizing the dynamics of cd4+ t cell priming within a lymph node. J Immunol. 2010; 184(6):2873–85.

  37. 37

    Pufnock JS, Cigal M, Rolczynski LS, Andersen-Nissen E, Wolfl M, McElrath MJ, Greenberg PD. Priming cd8+ t cells with dendritic cells matured using tlr4 and tlr7/8 ligands together enhances generation of cd8+ t cells retaining cd28. Blood. 2011; 117(24):6542–51.

  38. 38

    Abbas AK, Lichtman AH. Cellular and Molecular Immunology, 5th edn. Philadelphia: Saunders; 2003.

  39. 39

    Nowak AK, Lake RA, Marzo AL, Scott B, Heath WR, Collins EJ, Frelinger JA, Robinson BW. Induction of tumor cell apoptosis in vivo increases tumor antigen cross-presentation, cross-priming rather than cross-tolerizing host tumor-specific cd8 t cells. J Immunol. 2003; 170(10):4905–13.

  40. 40

    Wei LZ, Xu Y, Nelles EM, Furlonger C, Wang JC, Di Grappa MA, Khokha R, Medin JA, Paige CJ. Localized interleukin 12 delivery for immunotherapy of solid tumours. J Cell Mol Med. 2013; 17(11):1465–74.

  41. 41

    Cao X, Leonard K, Collins LI, Cai SF, Mayer JC, Payton JE, Walter MJ, Piwnica-Worms D, Schreiber RD, Ley TJ. Interleukin 12 stimulates ifn gamma mediated inhibition of tumor induced regulatory t cell proliferation and enhances tumor clearance. Cancer Res. 2009; 69(22):8700–9.

  42. 42

    McGray AJ, Hallett R, Bernard D, Swift SL, Zhu Z, Teoderascu F, Vanseggelen H, Hassell JA, Hurwitz AA, Wan Y, Bramson JL. Immunotherapy-induced CD8+ T cells instigate immune suppression in the tumor. Mol Ther. 2014; 22(1):206–18.

  43. 43

    Wang L, Hernandez-Alcoceba R, Shankar V, Zabala M, Kochanek S, Sangro B, Kramer MG, Prieto J, C Q. Prolonged and inducible transgene expression in the liver using gutless adenovirus: a potential therapy for liver cancer. Gastroenterology. 2004; 126(1):278–89.

  44. 44

    Lichtenauer M, Franz M, Fritzenwanger M, Figulla HR, Gerdes N, Jung C. Elevated plasma levels of interleukin-12p40 and interleukin-16 in overweight adolescents. Biomed Res Int. 2015; 2015:940910.

  45. 45

    Zabala M, Alzuguren P, Benavides C, Crettaz J, Gonzalez-Aseguinolaza G, Ortiz de Solorzano C, Gonzalez-Aparicio M, Kramer MG, Prieto J, Hernandez-Alcoceba R. Evaluation of bioluminescent imaging for noninvasive monitoring of colorectal cancer progression in the liver and its response to immunogene therapy. Mol Cancer. 2009; 8(2):13.

  46. 46

    Choy G, O’Connor S, Diehn FE, Costouros N, Alexander HR, Choyke P, Libutti SK. Comparison of noninvasive fluorescent and bioluminescent small animal optical imaging. Biotechniques. 2003; 35(5):1022–6.

  47. 47

    Duriancik DM, Hoag KA. The identification and enumeration of dendritic cell populations from individual mouse spleen and peyers patches using flow cytometric analysis. Cytometry A. 2009; 75(11):951–9.

  48. 48

    Ishwaran H, Kogalur UB, Gorodeski EZ, Minn AJ, Lauer MS. High-dimensional variable selection for survival data. J Am Stat Assoc. 2010; 105:205–17.

  49. 49

    Ishwaran H, Kogalur UB, Chen X, Minn AJ. Random survival forests for high-dimensional data. Stat Anal Data Min. 2011; 4:115–32.

  50. 50

    Twyman-Saint Victor C, Rech AJ, Maity A, Rengan R, Pauken KE, Stelekati E, Benci JL, et al. Radiation and dual checkpoint blockade activate non-redundant immune mechanisms in cancer. Nature. 2015; 520(67547):373–7.

Download references

Acknowledgements

The authors would like to thank a variety of undergraduate research assistants for assisting with calibrating parameter values of functions describing relationships between IL-12 concentrations and Mif or time in days to experimental data and the sensitivity analysis.

Funding

This work was supported by grants from the National Science Foundation (CAREER 1053490 to DJK), the National Cancer Institute (R01CA193473 to DJK), and the National Institute of General Medical Sciences of the National Institutes of Health grant as part of the West Virginia IDeA Network of Biomedical Research Excellence (P20GM103434 to QW). The content is solely the responsibility of the authors and does not necessarily represent the official views of the NSF, the NCI, or the National Institutes of Health.

Author information

DJK and QW conceived the study. DJK and QW developed the model and drafted the manuscript, ZW wrote the computer simulation code and designed the online simulators, YW performed stability analysis, QW did the numerical simulations and analyzed simulation data. All authors have read and approved the manuscript.

Correspondence to David J. Klinke II.

Ethics declarations

Ethics approval and consent to participate

The authors declare that no animal or human experiments have been performed as part of the research for this manuscript.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Wang, Q., Wang, Z., Wu, Y. et al. An in silico exploration of combining Interleukin-12 with Oxaliplatin to treat liver-metastatic colorectal cancer. BMC Cancer 20, 26 (2020) doi:10.1186/s12885-019-6500-9

Download citation

Keywords

  • Adenoviral vector
  • Combination therapy
  • Mathematical modeling
  • Impulsive ordinary differential equation
  • Stability analysis