 Research Article
 Open Access
 Published:
An in silico exploration of combining Interleukin12 with Oxaliplatin to treat livermetastatic colorectal cancer
BMC Cancer volume 20, Article number: 26 (2020)
Abstract
Background
Combining anticancer 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 Interleukin12 (IL12) can eliminate preexisting 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 Mifepristoneinduced IL12 and chemotherapy agent OXP.
Methods
A multiscale 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 IL12 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, partialresponders, and nonresponders.
Results
The multiscale model captures tumor growth patterns of the three phenotypic responses observed in mice in response to the combination therapy against a tumor rechallenge and was used to explore the impacts of changing the dose and timing of the mixed immunechemotherapy on tumor growth subjected to a tumor rechallenge 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 IL12 therapy worked more efficiently in responders by increased priming of T cells, enhanced CD8 + T cellmediated killing, and functional inhibition of regulatory T cells. In a virtual cohort that mimics nonresponders and partialresponders, simulations show that an increased dose of OXP alone would improve the response. In addition, enhanced IL12 expression alone or an increased number of treatment cycles of the mixed immunechemotherapy can barely improve tumor control for nonresponders 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 cancerrelated 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].
Interleukin12 (IL12) 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 IL12 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 readministration is impaired by the rapid emergence of neutralizing antibodies [3]. To allow a good control of the strength and duration of IL12 expression, highcapacity adenoviral vectors containing a liverspecific, Mifepristoneinducible system for the expression of murine IL12 (HCAd/RUmIL12) were recently designed to control primary or metastatic liver cancer [6]. Since standalone chemo or radiotherapeutic regimens are insufficient (with a few notable exceptions) to eradicate neoplastic lesion [7], HCAd/RUmIL12 was combined with chemotherapy agent Oxaliplatin (OXP) to treat liverimplanted colon cancer cells [6]. As a consequence of the combination therapy, preexisting liver metastases of colorectal cancer were eradicated, and enhanced establishment of a protective immune response against tumor rechallenge 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 antiPDL1 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 ODEbased 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 immunochemotherapy 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 liverspecific expression of IL12 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 tumorassociated antigenspecific T cell response.
The structure of this paper is as follows. First, we present a multiscale mechanistic model of antitumor immunity and tumor growth in response to a combined immunochemotherapy 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 tumorfree and hightumor 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 nonresponders and partial responders.
Methods
Our method was to develop a multiscale 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 RungeKutta 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 nonresponder) 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 longterm 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 nonresponder mice after the combination treatment, and local stability analysis are discussed in the following subsections.
Results
A multiscale model of tumor growth subject to IL12 and OXP therapy
Our mathematical model is based on the experimental data presented by Manuela GonzalezAparicio and collaborators in [6] using the MC38Luc1 cell lines for murine metastatic colorectal cancer. Using this mouse model, OXP and IL12 combination therapy eradicated preexisting liver metastases, established a protective immune response against tumor rechallenge, and increased overall survival of animals. To better understand the dynamics of the primary response to adenovirusmediated induction of an antitumor immune response, we developed a threecompartment mathematical model to quantify the cytotoxic CD8 ^{+} T cell response to IL12 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.
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^{3}). We assume that naïve CD8 ^{+} T cells are produced at a constant rate c_{1} 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_{1} (antigenpresenting 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^{3}). 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^{3} mm^{3} is the volume of the blood compartment [35] and Vol_{ln}=0.25 mm^{3} 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_{12}·T_{E1} is the efflux rate. We assume that T_{E1} cells are killed by chemotherapy agent OXP_{1} (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_{1}, units: cells per mm^{3}). We assume that APCs in the lymph node have a natural death rate of k_{d3}·APC_{1} 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_{3} is the concentration of APCs in tumor, \({Vol}_{t}= \frac {\epsilon +C(t)}{1V_{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^{3}. 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_{1}, units: mg per kg). We assume that OXP decays naturally at a rate k_{d4}·OXP_{1} 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_{2} 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^{3}). 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_{21}·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_{23}·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_{2}, units: cells per mm^{3}). 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_{2} is the growth rate constant and K is the carrying capacity. We assume a b_{23}·APC_{2} 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_{2}, units: mg per kg). We assume that OXP decays naturally at a rate k_{d4}·OXP_{2} and the efflux rates of OXP from blood to lymph node and from blood to tumor are c_{21}·OXP_{2} and c_{23}·OXP_{2}, 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^{3}). 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 IL12 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 IL12 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_{3}) 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^{3}). We assume that IFN_{γ} is secreted solely by effector CD8 ^{+} T cells within the tumor with stimulation from IL12 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_{3}, units: cells per mm^{3}). We assume that APCs in the tumor microenvironment have a natural death rate of k_{d3}·APC_{3}, 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_{31}·APC_{3}.
$$ {}\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)Interleukin12 (IL, units: ng per ml). Interleukin12 (IL12) is produced by APCs at a rate of c_{5}·APC_{3} and decays naturally at a rate of k_{d8}·IL. The extra IL12 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_{3}, units: mg per kg). We assume that OXP decays naturally at a rate k_{d4}·OXP_{3} 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^{3}). Regulatory T cells are produced at a constant rate c_{6} 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^{3}). 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 cellmediated 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^{+}}\\ &\quadk_{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}}\\ &\quadk_{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^{3}). 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 rechallenge 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 IL12 and OXP caused by the therapies, respectively. The sudden change in the volume of MHC class I negative tumor cells due to tumor rechallenge 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_{2}(t)=OXP_{2}(t^{+})−OXP_{2}(t^{−}) reflect the abrupt changes of IL12 and oxaliplatin at administration times t_{k1} and t_{k2}, while IL_{k1} and {OXP_{2}}_{k2} are the dosages of IL12 and oxaliplatin at the administration times t_{k1} and t_{k2} with k1=1,2,⋯,n_{1} and k2=1,2,⋯,n_{2}, respectively; ΔC_{MHCI}(t)=C_{MHCI}(t^{+})−C_{MHCI}(t^{−}) represents the sudden changes in tumor size due to tumor rechallenge with implantation size C_{k3} mm^{3} at time t_{k3} with k3=1,2,⋯,n_{3}.
A schematic diagram summarizing this threecompartment model is shown in Fig. 1. Model parameters and their meanings are listed in Table 1.
Nonnegativity of solutions to the model
For any mathematical model that has biological implications, it is important to make sure that solutions with nonnegative initial conditions remain nonnegative all the time. For the model comprised of Eqs. 1)  18), we can see that, with positive impulsive inputs IL_{k1}, {OXP_{2}}_{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 nonnegative initial conditions will remain nonnegative 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 IL12 requires efficient methods to control the plasma concentration of this potent immunostimulatory cytokine in order to avoid toxicity [6]. It was determined in an MC38 syngeneic tumor model that a blood concentration of IL12 less than 20 ng/ml has no antitumor effect, while concentrations greater than 700 ng/ml are associated with toxicity [43]. For comparison, the normal range of IL12 in humans is around 7.5 pg/ml [44]. GonzalezAparicio and colleagues designed a new induction protocol to keep IL12 within this therapeutic range [6]. Once the liver of a group of C57BL/6 mice was transduced with the vector (typically 2.5* 10^{8} IU), a suboptimal amount of Mif (125 μg/kg) is administered for the first 2 days in order to prevent toxicity. The concentration of IL12 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 IL12 expression in mice treated with the HCAd/RUmIL12 vector (Fig. 2). Before we start the calibration of model parameters, we first quantified the IL12 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 IL12 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 Mifinduced IL12 treatment effects. Since the authors [6] verified the Mifinduction 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 Mifinduced treatment cycles in our model for simplicity.
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 (IL12 + OXP group).
Concentration of Interferon gamma was obtained from Fig. 4b (IL12 + OXP group).
The ratio between CD8 ^{+} T lymphocytes and T regulatory cells was obtained from Fig. 5b (experimental results for IL12 + OXP group in tumor).
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^{5} MC38Luc1) and the immunological protection against cancer cells treated with the combination therapy after a tumor rechallenge (10^{6} 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}=0 is the day that 5×10^{5} cells /mouse MC38Luc1 tumor cells were inoculated in the liver of C57BL /6 mice [45, 46], T_{N}(0)=0.0714 cells per mm^{3} (\(=\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^{3} mm^{3}), \({APC}_{2}(0)=214.2857\;(=\frac {3*10^{5}}{1.4*10^{3}})\) cells per mm^{3} 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_{2}(s)=5 mg/ kg with s=10, 34, 100; ΔIL(t) follows \(f(t)=\frac {13.6127*(t12)^{2}+0.8606*(t12)+1}{0.313*(t12)^{2}0.6216*(t12)+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 longterm 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 10day IL12 induction starting day 103, which follows the adjusted protocol as described in Fig. 2 in [6], were administered after a tumor rechallenge on day 75. This treatment occurred about two weeks after the mice survived two cycles of combined treatments with 10day Mifinduced IL12 (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 nonresponders (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_{2}(9)=5 mg/kg; ΔIL(t) follows \(IL(t)=\frac {13.6127*(ta)^{2}+0.8606*(ta)+1}{0.313*(ta)^{2}0.6216*(ta)+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: nonresponders versus responders While both the responders and nonresponders survived two cycles of combination therapy treatment before tumor rechallenge and then underwent the third cycle following the tumor rechallenge on day 75, the simulations in Fig. 3b  d suggest that the nonresponders 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^{3} cells per mm^{3} 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 + antiCTLA4 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 highdimensional 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 IL12 plus OXP. Experimental results for longterm management of colorectal cancer by observing cooperation of IL12 and OXP for the control of experimental relapses in distant locations are compared against model predictions in Fig. 3b, c, and d for responders, partialresponders and nonresponders, 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 tumorfree equilibrium \(\vec {X}_{0}\), a second tumor free equilibrium \(\vec {X}_{1}\) when APC_{2} growth rate is larger than the rate constant for APC_{2} flowing from blood to tumor (i.e., r_{2}>b_{23}), 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
where \(T_{E_{1}}\) satisfies the following polynomial equation
(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
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
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_{23}+r_{2}, −k_{d4}, −k_{d4}, −k_{d4}−c_{21}−c_{23}, −k_{d6}, −k_{d7}, −b_{31}−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 λ^{2}+(k_{d5}+a_{12}+a_{21}+a_{23})λ+a_{12}(k_{d5}+a_{23})=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_{2}<b_{23} (i.e., APC_{2} growth rate is smaller than the rate constant for APC_{2} 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_{2}>b_{23} (i.e., APC_{2} growth rate is larger than the rate constant for APC_{2} 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 lefthalf 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_{23}−r_{2}, and the other two satisfies the following quadratic equation
It is found that both eigenvalues are in the lefthalf of the complex plane if k_{p4}<k_{d10}. Hence this tumorfree equilibrium point is locally stable if r_{2}>b_{23} 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_{23}+r_{2},−k_{d4}−c_{21}−c_{23},−k_{d7},−b_{31}−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(λ)=λ^{3}+α_{1}λ^{2}+α_{2}λ+α_{3}, where
Based on the list of the first twelve eigenvalues, the hightumor equilibrium is unstable if either r_{2}>b_{23} or k_{d10}>k_{p4} is satisfied. Suppose r_{2}<b_{23} and k_{d10}<k_{p4}, it is easy to see that α_{i}>0, i=1,2,3. According to the RouthHurwitz 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 hightumor equilibrium is stable if r_{2}<b_{23} (i.e., APC_{2} growth rate is smaller than the rate constant for APC_{2} 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_{2}<b_{23} 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 IL12 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_{2}>b_{23}), MHCI negative tumor cells (\(\mathrm {C}_{\text {MHCI}^{}}\)) will either be completely removed in which case solutions of this dynamic approach the second tumorfree 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, partialresponders, and nonresponders, 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 rechallenge) 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^{3} using 0.1×α_{i} and other parameters in the ith calibrated parameter set and \(\bar {y_{i}}\) is the predicted tumor size in mm^{3} 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 nonresponders. In addition, there are 10 (out of 48) parameters whose value changes greatly affect tumor size of responders but not the size of nonresponders and partialresponders. These potentially OXP and IL12 treatment important parameters include c_{23} (OXP flow rate from blood to tumor), K (APC carrying capacity), c_{4}  c_{6} (IFN_{γ}, IL12, and T_{R} production rate constants, respectively), g_{10}  g_{14} (IL12, IFN_{γ}, OXP_{3}, T_{R}, and \({\mathrm {C}}_{MHCI}^{+}\) killing by OXP_{3} saturation rate constants, respectively). In addition, changes in the value of the following 7 parameters cause from zero for nonresponders to increasing changes in normalized tumor size from partialresponders to responders: T cell flow rates from blood to lymph node and from tumor to blood, a_{21} and a_{32}, respectively; APC flow rates from tumor to lymph node b_{31}; production rate constant of naive T cells c_{1}; transfer rate constant of naive T cell to T effector cell in lymph node c_{2}; IL12 natural death rate constant k_{d8}, and APC growth rate constant r_{2}. We also note that no change in tumor size for all three mice group (nonresponders, partialresponders, and responders) results from the value changes of following 5 parameters: \({\mathrm {C}}_{MHCI}^{}\) killing by OXP_{3} saturation rate constant g_{15}, 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 partialresponders and nonresponders, we simulated the following alternative treatment scenarios: changing the dose and frequency of chemotherapy drug OXP administration, changing the strength of Mifinduced IL12 expression, and changing the number of combined treatment cycles.
Partialresponders
We note from Fig. 5 that increased number of treatment cycles in the IL12 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 IL12 expression alone would slightly reduce tumor size more rapidly after the tumor reaches its maximum size.
Nonresponders
From Fig. 6, we found that neither increased strength of IL12 expression nor moderately increased OXP dose alone in the IL12 and oxaliplatin (OXP) combination therapy seems to improve tumor control for the median, 90th percentile lower and 90th percentile upper responses for the 30 nonresponder 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 IL12 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.
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 threecompartment 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 IL12 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 IL12 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 nonresponders. 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 nonresponders, partialresponders, 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 nonresponders, numerical simulation of multiple OXP and IL12 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 antigenspecific immunochemotherapies 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
 IL12:

Interleukin12
 MHC:

Major histocompatibility complex
 NK:

Natural killer cells
 ODE:

Ordinary differential equations
 OXP:

Oxaliplatin
 RF:

Random forest
 TIL:

Tumor infiltrating lymphocyte
References
 1
Cancer.org: Colon and Rectum Cancer. http://www.cancer.org/cancer/colonandrectumcancer/detailedguide/. Accessed 1 July 2019.
 2
Martinet O, Ermekova V, Qiao JQ, Sauter B, Mandeli J, Chen L, Chen SH. Immunomodulatory gene therapy with interleukin 12 and 41bb ligand: longterm remission of liver metastases in a mouse model. J Natl Cancer Inst. 2000; 92(11):931–6.
 3
HernandezAlcoceba R, Berraondo P. Immunochemotherapy against colon cancer by gene transfer of interleukin12 in combination with oxaliplatin. Oncoimmunology. 2012; 1(1):97–99.
 4
Colombo MP, Trinchieri G. Interleukin 12 in antitumor immunity and immunotherapy. Cytokine Growth Factor Rev. 2002; 13(2):155–68.
 5
Klinke DJ. A multiscale systems perspective on cancer, immunotherapy, and interleukin 12. Mol Cancer. 2010; 9(242):18.
 6
GonzalezAparicio M, Alzuguren P, Mauleon I, MedinaEcheverz J, HervasStubbs S, Mancheno U, Berraondo P, Crettaz J, GonzalezAseguinolaza G, Prieto J, HernandezAlcoceba 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
Vacchelli E, Prada N, Kepp O, Galluzzi L. Current trends of anticancer immunochemotherapy. Oncoimmunology. 2013; 2(6):3.
 8
Kosinsky Y, Dovedi SJ, Peskov V. KndVoronova, Chu L, Tomkinson H, AlHuniti 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
Wodarz D, Komarova NL. Computational Biology of Cancer: Lecture Notes and Mathematical Modeling. New Jersey: World Sci; 2005.
 10
Klinke DJ, Wang Q. Understanding immunology via engineering design: the role of mathematical prototyping. Comput Math Methods Med. 2012; 9:9.
 11
Klinke DJ. In silico modelbased inference: A contemporary approach for hypothesis testing in network biology. Biotechnol Prog. 2014; 30:1247–61.
 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
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
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
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
Zheng Y, Bresson D, Fradkin M, Chan JR, von Herrath M, C W. Biosimulations predict optimal oral insulin ∖ anticd3 and oral insulin/exendin4 combination treatment regimens for the reversal of diabetes in the nonobese diabetic (nod) mouse. Clin Immunol. 2008; 127:101–2.
 17
Pappalardo F, MartinezForero I, Pennisi M, Palazon A, Melero I, Motta S. Simb16: modeling induced immune system response against b16melanoma. PLoS ONE. 2011; 6(10):12.
 18
Dovedi SJ, Adlard AL, LipowskaBhalla G, McKenna C, et al. Acquired resistance to fractionated radiotherapy can be overcome by concurrent pdl1 blockade. Cancer Res. 2014; 74(19):5458–68.
 19
Dovedi SJ, Cheadle EJ, Popple A, Poon E, Morrow M, Stewart R, et al. Fractionated radiation therapy stimulates antitumor immunity mediated by both resident and infiltrating polyclonal tcell populations when combined with pd1 blockade. Clin Cancer Res. 2017; 23(18):5514–26.
 20
Kirschner D, Panetta JC. Modeling immunotherapy of the tumor immune interaction. J Math Biol. 1998; 37(3):235–52.
 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
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
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
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
Eladdadi A, P K, Mallet D. Mathematical Models of Tumor Immune System Dynamics. Switzerland: Springer; 2014.
 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
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
de Pillis LG, Eladdadi A, Radunskaya AE. Modeling cancer immune responses to therapy. J Pharmacokinet Pharmacodyn. 2014; 41(5):461–78.
 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
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
Jenkins MK, Chu HH, McLachlan JB, Moon JJ. On the composition of the preimmune repertoire of t cells specific for peptidemajor histocompatibility complex ligands. Ann Rev Immunol. 2010; 28:275–94.
 32
Miller MJ, Wei SH, Parker I, Cahalan MD. Twophoton imaging of lymphocyte motility and antigen response in intact lymph node. Science. 2002; 296(5574):1869–73.
 33
Krummel M, Bartumeus F, A G. Tcell migration, search strategies and mechanisms. Nat Rev Immunol. 2016; 16(3):193–201.
 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
C, 57BL ∖ 6 Inbred Mouse. https://www.jax.org/newsandinsights/2005/october/howmuchbloodcanitakefromamousewithoutendangeringitshealth. Accessed 30 Dec 2019.
 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
Pufnock JS, Cigal M, Rolczynski LS, AndersenNissen 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
Abbas AK, Lichtman AH. Cellular and Molecular Immunology, 5th edn. Philadelphia: Saunders; 2003.
 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 crosspresentation, crosspriming rather than crosstolerizing host tumorspecific cd8 t cells. J Immunol. 2003; 170(10):4905–13.
 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
Cao X, Leonard K, Collins LI, Cai SF, Mayer JC, Payton JE, Walter MJ, PiwnicaWorms 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
McGray AJ, Hallett R, Bernard D, Swift SL, Zhu Z, Teoderascu F, Vanseggelen H, Hassell JA, Hurwitz AA, Wan Y, Bramson JL. Immunotherapyinduced CD8+ T cells instigate immune suppression in the tumor. Mol Ther. 2014; 22(1):206–18.
 43
Wang L, HernandezAlcoceba 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
Lichtenauer M, Franz M, Fritzenwanger M, Figulla HR, Gerdes N, Jung C. Elevated plasma levels of interleukin12p40 and interleukin16 in overweight adolescents. Biomed Res Int. 2015; 2015:940910.
 45
Zabala M, Alzuguren P, Benavides C, Crettaz J, GonzalezAseguinolaza G, Ortiz de Solorzano C, GonzalezAparicio M, Kramer MG, Prieto J, HernandezAlcoceba 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
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
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
Ishwaran H, Kogalur UB, Gorodeski EZ, Minn AJ, Lauer MS. Highdimensional variable selection for survival data. J Am Stat Assoc. 2010; 105:205–17.
 49
Ishwaran H, Kogalur UB, Chen X, Minn AJ. Random survival forests for highdimensional data. Stat Anal Data Min. 2011; 4:115–32.
 50
TwymanSaint Victor C, Rech AJ, Maity A, Rengan R, Pauken KE, Stelekati E, Benci JL, et al. Radiation and dual checkpoint blockade activate nonredundant immune mechanisms in cancer. Nature. 2015; 520(67547):373–7.
Acknowledgements
The authors would like to thank a variety of undergraduate research assistants for assisting with calibrating parameter values of functions describing relationships between IL12 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
Affiliations
Contributions
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.
Corresponding author
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.
About this article
Cite this article
Wang, Q., Wang, Z., Wu, Y. et al. An in silico exploration of combining Interleukin12 with Oxaliplatin to treat livermetastatic colorectal cancer. BMC Cancer 20, 26 (2020). https://doi.org/10.1186/s1288501965009
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s1288501965009
Keywords
 Adenoviral vector
 Combination therapy
 Mathematical modeling
 Impulsive ordinary differential equation
 Stability analysis