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

BACKGROUND: Combining anti-cancer therapies with orthogonal modes of action, such as direct cytotoxicity and immunostimulatory, hold promise for expanding clinical benefit to patients with metastatic disease. For instance, a chemotherapy agent Oxaliplatin (OXP) in combination with Interleukin-12 (IL-12) can eliminate pre-existing liver metastatic colorectal cancer and protect from relapse in a murine model. However, the underlying dynamics associated with the targeted biology and the combinatorial space consisting of possible dosage and timing of each therapy present challenges for optimizing treatment regimens. To address some of these challenges, we developed a predictive simulation platform for optimizing dose and timing of the combination therapy involving Mifepristone-induced IL-12 and chemotherapy agent OXP.
METHODS: A multi-scale mathematical model comprised of impulsive ordinary differential equations was developed to describe the interaction between the immune system and tumor cells in response to the combined IL-12 and OXP therapy. An ensemble of model parameters were calibrated to published experimental data using a genetic algorithm and used to represent three different phenotypes: responders, partial-responders, and non-responders.
RESULTS: The multi-scale model captures tumor growth patterns of the three phenotypic responses observed in mice in response to the combination therapy against a tumor re-challenge and was used to explore the impacts of changing the dose and timing of the mixed immune-chemotherapy on tumor growth subjected to a tumor re-challenge in mice. An increased ratio of CD8 + T effectors to regulatory T cells during and after treatment was key to improve tumor control in the responder cohort. Sensitivity analysis indicates that combined OXP and IL-12 therapy worked more efficiently in responders by increased priming of T cells, enhanced CD8 + T cell-mediated killing, and functional inhibition of regulatory T cells. In a virtual cohort that mimics non-responders and partial-responders, simulations show that an increased dose of OXP alone would improve the response. In addition, enhanced IL-12 expression alone or an increased number of treatment cycles of the mixed immune-chemotherapy can barely improve tumor control for non-responders and partial responders.
CONCLUSIONS: Overall, this study illustrates how mechanistic models can be used for in silico screening of the optimal therapeutic dose and timing in combined cancer treatment strategies.

Background

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

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

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

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

Methods

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

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

Results

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

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

Fig. 1

Schematic diagram illustrating the interactions among species present in the three compartments. State variables and transport relations are shown in black. Parameters are in red while influence relationships are in blue. Naïve CD8 + T cells (T) are activated and become CD8 + T effectors (T) 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 (T) and enter tumor microenvironment (T) 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 ( and ) through increased presentation of tumor-associated antigens by Major Histocompatibility Complex protein class I (MHCI). During this process, IL-12 (IL) helps promote T cell proliferation and suppresses regulatory T (T) cells’ proliferation and immunosuppressive action on effector CD8 + T cells. In addition, the chemotherapy drug Oxaliplatin in the lymph node and tumor (OXP where i=1,3) will kill fast-proliferating cells such as T effectors and tumor cells

Table 1

List of parameters in the model

Parameter	Units	Description
kd1	day−1	Naïve CD8 + T cell natural death rate constant
kd2	day−1	T effector (in lymph node or tumor) or tumor cell death rate constant due to OXP
kd3	day−1	APC (in lymph node or tumor) natural death rate constant
kd4	day−1	OXP natural decay rate constant
kd5	day−1	T effector (in blood) natural death rate constant
kd6	day−1	T effector (in tumor) natural death rate constant
kd7	day−1	Interferon γ natural death rate constant
kd8	day−1	IL-12 natural death rate constant
kd9	day−1	Regulatory T cell natural death rate constant
kd10	day−1	Tumor cell natural death rate constant
kd11	mm3∗day−1	MHC class I positive tumor cell death rate constant due to T effector (in tumor) lysis
kp1	day−1	T effector (in lymph node) proliferation rate constant due to tumor antigens presented by APC in lymph node
kp2	day−1	T effector (in tumor) proliferation rate constant
kp3	ng∗mm−3·day−1	Regulatory T cell proliferation rate constant due to tumor growth and proliferation of T effector in tumor
kp4	day−1	Tumor cell proliferation rate constant
a12	day−1	Rate constant for T cell flow from lymph node to blood
a21	day−1	Rate constant for T cell flow from blood to lymph node
a23	day−1	Rate constant for T cell flow from blood to tumor
a32	day−1	Rate constant for T cell flow from tumor to blood
b23	day−1	Rate constant for APC flow from blood to tumor
b31	day−1	Rate constant for APC flow from tumor to lymph node
c1	cell·mm−3·day−1	Naïve T cell natural production rate constant
c2	day−1	Naïve T cell to T effector (in lymph node) transfer rate constant
c21	day−1	Rate constant for OXP flow from blood to lymph node
c23	day−1	Rate constant for OXP flow from blood to tumor
c4	pg·mm−3·day−1	Interferon γ secretion constant
c5	pg·cell−1·day−1	IL-12 production rate constant by APC in tumor
c6	cells·mm−3·day−1	Regulatory T cell production rate constant
c7	day−1	MHC class I negative to positive tumor cells transfer rate constant
α	(cell·mm−3)2	T effector (in lymph node) saturation constant
K	cell·mm−3	Carrying capacity of APC (in blood)
g1	cell·mm−3	APC (in lymph node) saturation constant
g2	cell·mm−3	APC (in lymph node) saturation constant
g3	mg·kg−1	OXP (in lymph node) saturation constant
g4	ng∗ml	IL-12 saturation rate constant
g5	cell·mm−3	Regulatory T cell saturation constant
g6	mg·kg−1	OXP (in tumor) saturation constant
g7	ng·ml	IL-12 saturation rate constant
g8	cell·mm−3	Regulatory T cell saturation rate constant
g9	cell·mm−3	T effector (in tumor) saturation rate constant
g10	ng·ml	IL-12 saturation constant
g11	pg·mm−3	Cellular Interferon γ saturation constant
g12	mg·kg−1	OXP (in tumor) saturation rate constant
g13	cell·mm−3	Regulatory T cell saturation rate constant
g14	mg·kg−1	OXP (in tumor) killing MHC class I positive tumor cells saturation rate constant
g15	mg·kg−1	OXP (in tumor) killing MHC class I negative tumor cells saturation rate constant
r1	cell−1·day−1	Constant in tumor logistic growth
r2	day−1	Growth rate constant of APC (in blood)

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 (, units: cells per mm). We assume that naïve CD8 + T cells are produced at a constant rate c1 from thymus and die naturally at a rate k·T [31]. Naïve T cells are recruited and activated by tumor antigens presented by APC1 (antigen-presenting cells in lymph node) at a rate [32-34].

Effector CD8 T cells in lymph node (, units: cells per mm). The increase in the rate of concentration of effector CD8 + T cells in the lymph node due to activation of naïve CD8 + T cells from the blood stream is given by , where Vol=1.4∗103 mm3 is the volume of the blood compartment [35] and Vol=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 T, a saturable term that represents antigen presenting cells (APC) and defined by , and an immune checkpoint term defined by , where α is the square root of the saturation constant of T [37]. We also assume that influx rate of effector T cells from blood to lymph node is , where T is the concentration of T effectors in blood, and a12·T is the efflux rate. We assume that T cells are killed by chemotherapy agent OXP1 (Oxaliplatin in lymph node) at the rate .

Antigen Presenting Cells in lymph node (APC, units: cells per mm). We assume that APCs in the lymph node have a natural death rate of k·APC1 and the influx rate of APCs from tumor to lymph node is , where APC3 is the concentration of APCs in tumor, (since Vol=ε+C(t)+V·T·Vol, where T 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 () and MHC class I negative tumor cells (). The average size of a T effector cell (V) is equal to [38].

Chemotherapy agent Oxaliplatin in lymph node (OXP, units: mg per kg). We assume that OXP decays naturally at a rate k·OXP1 and the influx rate of Oxaliplatin (OXP) from blood to lymph node is , where OXP2 is the concentration of OXP in blood.

Effector CD8 T cells in blood (, units: cells per mm). We assume the effector CD8 + T cells die naturally in blood at a rate k·T. The influx rate of effector CD8 + T cells from lymph node to blood is equal to and the efflux rate of effector CD8 + T cells from blood to lymph node is equal to a21·T. The influx rate of CD8 + T effectors from the tumor to blood is , where T is the concentration of T effectors in tumor and the efflux rate of CD8 + T effectors from blood to tumor is a23·T.

Antigen Presenting Cells in blood (APC, units: cells per mm). A logistic growth pattern 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.

Chemotherapy agent Oxaliplatin in blood (OXP, units: mg per kg). We assume that OXP decays naturally at a rate k·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).

Effector CD8 T cells in tumor microenvironment (, units: cells per mm). We assume that effector CD8 + T cells can proliferate locally upon recognition of the corresponding tumor antigen presented by MHC class I positive tumor cells upon IL-12 stimulation and subject to suppression from T regulatory cells at a saturable rate equal to , where IL is the concentration of IL-12 and T 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 . The efflux rate of effector CD8 + T cells from the tumor to blood is . Effector T cells have a finite lifespan and die within the tumor microenvironment at a rate equal to k·T. T effector cells are assumed to be killed by chemotherapy agent OXP (OXP3) at the rate in the tumor microenvironment.

Interferon gamma (IFN, units: pg per mm). We assume that IFN is secreted solely by effector CD8 + T cells within the tumor with stimulation from IL-12 and inhibition from regulatory T cellsat a rate of [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.

Antigen Presenting Cells in tumor (APC, units: cells per mm). We assume that APCs in the tumor microenvironment have a natural death rate of k·APC3, the influx rate of APCs from blood to tumor is 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.

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 k·IL. The extra IL-12 expression obtained through the combined therapy is approximated using the discrete Eq. (16).

Chemotherapy agent Oxaliplatin in tumor (OXP, units: mg per kg). We assume that OXP decays naturally at a rate k·OXP3 and the influx rate of OXP from blood to tumor is .

Regulatory T cells (, units: cells per mm). Regulatory T cells are produced at a constant rate c6 from thymus and die naturally at a rate k·T. We assume that regulatory T cells proliferate at a rate of [3, 4, 6].

MHC class I positive tumor cells (, units: mm). MHC class I positive tumor cells are converted from MHC class I negative tumor cells () with the assistance of IFNγ at a rate and the rate of effector CD8 + T cell-mediated killing of MHC class I positive tumor cells is [6, 31]. We assume that the dilution rate of MHC class I positive tumor cells due to proliferation is . The natural death rate of MHC class I positive tumor cells is assumed to be and MHC class I positive tumor cells are killed by chemotherapy agent OXP in tumor at a rate .

MHC class I negative tumor cells (, units: mm). MHC class I negative tumor cells are converted to MHC class I positive tumor cells with the assistance of IFN at a rate of . We assume that the proliferation rate of MHC class I positive tumor cells is equal to . 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 and these cells are killed by chemotherapy agent OXP in tumor at a rate . The difference in tumor volume caused by tumor re-challenge is described by the discrete Eq. (18).

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).

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 t and t, while IL and {OXP2} are the dosages of IL-12 and oxaliplatin at the administration times t and t 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 C mm3 at time t 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 IL, {OXP2}, and C at impulsive moments t, t and t in Eqs. (16) - (18), all solutions of the system of ODEs (1) - (15) with non-negative initial conditions will remain non-negative because for x=0 and x≥0, where i,j=1,2,⋯,15 and i≠j.

Model calibration

Therapeutic use of IL-12 requires efficient methods to control the plasma concentration of this potent immuno-stimulatory cytokine in order to avoid toxicity [6]. It was determined in an MC38 syngeneic tumor model that a blood concentration of IL-12 less than 20 ng/ml has no anti-tumor effect, while concentrations greater than 700 ng/ml are associated with toxicity [43]. For comparison, the normal range of IL-12 in humans is around 7.5 pg/ml [44]. Gonzalez-Aparicio and colleagues designed a new induction protocol to keep IL-12 within this therapeutic range [6]. Once the liver of a group of C57BL/6 mice was transduced with the vector (typically 2.5* 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

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

Comparison of model predictions with experimental measures of therapeutic response upon tumor re-challenge. a. Comparison of model predictions with experimental measures of tumor volume, IFN and T/ T 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 T/ T (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

Parameter	Figure 3a.	Figure 3b. Par.-Res.	Figure 3c. Res.	Figure 3d. Non-Res.
kd1	8.060∗10−2	4.613∗10−5	7.892∗10−6	7.629∗10−5
kd2	3.954∗10−2	9.923	3.286	8.370
kd3	5.738∗10−1	3.863	5.982	5.727∗10−1
kd4	1.760	2.658	1.924∗10−3	2.358
kd5	6.293∗10−3	3.312∗10−4	5.385∗10−1	6.170∗10−5
kd6	9.849∗10−3	5.026∗10−6	2.258∗10−2	5.988∗10
kd7	8.159∗10−2	3.390∗10−2	7.620∗10−2	5.603
kd8	6.097∗10−1	5.085∗10−6	4.692∗10−4	2.513∗10−7
kd9	4.382∗10−6	5.310∗10−6	2.206∗10−2	9.706∗10−7
kd10	2.116∗10−5	2.140∗10−6	7.795∗10−6	9.434∗10−6
kd11	7.956∗10−3	6.953∗10−6	7.047∗10−5	6.539∗10−5
kp1	2.188∗10	4.770∗104	9.808∗10	5.515∗102
kp2	6.445∗10−6	8.367∗10−7	1.687∗10−12	9.255∗10−7
kp3	5.015∗10−7	7.807∗10−9	4.276∗10−7	8.696∗10−7
kp4	5.800∗10−2	3.952∗10−1	3.297∗10−1	2.186∗10−1
a12	5.497∗10−1	9.636	3.031∗10−2	8.575∗10−2
a21	1.133∗10−4	9.984∗10−1	9.011∗10−4	1.614∗10−1
a23	7.254∗10−7	9.567∗10−6	2.431∗10−3	4.776∗10−11
a32	4.575∗10−3	9.246∗10−1	2.573∗10−3	7.638∗10−3
b23	9.581∗105	6.334∗10−9	1.406∗10−2	5.710∗10−13
b31	6.872∗10−7	8.489∗10−8	5.700∗10−10	3.980∗10−9
c1	8.563∗10−4	9.611	4.366∗10−4	4.108∗10−2
c2	7.360∗10−2	8.751∗10−2	3.373∗10	6.011∗10−2
c21	7.211∗10−8	5.684∗10−4	5.475∗10−2	9.907∗10−1
c23	3.648∗10−6	5.873	3.355	9.905∗10−1
c4	6.878∗10−2	8.671∗104	7.547∗102	2.327∗102
c5	5.263∗105	9.844∗10−9	7.279∗10−10	9.722∗10−11
c6	1.490∗10−2	3.260∗10−2	6.253∗10−4	7.635∗102
c7	8.842∗102	7.546	3.412∗102	9.289∗103
α	5.530∗109	3.647∗106	1.556∗107	9.834∗103
K	7.539∗104	9.292∗1012	3.116∗1011	4.964∗105
g1	7.971∗109	4.578	9.304∗102	3.445∗104
g2	2.441∗10−2	3.568∗10−10	6.108∗10−13	5.108∗10−12
g3	3.561∗10−7	1.081∗10−11	9.294∗10−12	5.325∗10−5
g4	6.740∗102	9.535∗104	1.518∗10	3.713∗10
g5	9.702∗107	2.440∗105	9.131∗10−3	8.235∗1011
g6	2.318∗105	8.034∗103	3.843∗106	6.932∗104
g7	9.387∗107	9.307∗108	9.040∗108	7.528∗108
g8	8.412∗10−5	7.079∗10−9	7.897∗10−7	5.385∗10−8
g9	3.757∗10−5	5.234∗10−3	5.185∗10−8	9.378∗10
g10	6.242∗10−5	9.338∗10−6	5.537∗10−6	5.110∗10−8
g11	2.183∗108	5.064∗1010	7.624∗106	4.565∗109
g12	3.229∗10−8	1.542∗10−11	1.930∗10−10	7.116∗10−9
g13	8.733∗106	2.273∗107	4.689∗107	9.747∗105
g14	5.981∗103	4.809∗109	4.799∗1010	3.855∗104
g15	3.068∗109	7.927∗105	8.642∗103	5.392∗106
r1	2.140∗10−10	2.448∗10−2	6.718∗10−2	7.260∗10−5
r2	8.065	2.780∗10−4	3.394∗10−9	2.793∗10−6

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

Fig. 4

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

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

Fig. 5

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

The model was calibrated against these data to reflect the effects of the combination therapy with one treatment cycle after tumor cell liver implantation (5×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 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], T(0)=0.0714 cells per mm3 ( 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∗103 mm3), cells per mm3 according to [47]. Other initial values are zero: for i=1,2,3. , ΔOXP2(s)=5 mg/ kg with s=10, 34, 100; ΔIL(t) follows 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 , ΔOXP2(9)=5 mg/kg; ΔIL(t) follows for a≤t≤a+9 with a=13, 37, 103. A sample set of parameter values for each of the response groups of mice used in the simulations are listed in Table 2.

Difference in treatment efficacy: non-responders versus responders While both the responders and non-responders survived two cycles of combination therapy treatment before tumor rechallenge and then underwent the third cycle following the tumor re-challenge on day 75, the simulations in Fig. 3b - d suggest that the non-responders show near zero concentration of IFN and near zero ratio of T effectors to regulatory T cells in the tumor all the time comparing to a stable concentration of IFN and ratio of T effectors to regulatory T cells in the tumor (at least 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 [48-50]. In [48] and [49], Random forest (RF) machine learning analysis of tumor infiltrating lymphocytes (TILs) demonstrated that the top predictor of resistance, as measured by variable importance scores and selection, was the CD8+CD44+ to Treg (CD8/Treg) ratio. The authors in [50] indicated that the CD8/Treg ratio in resistant tumors failed to increase after RT + anti-CTLA4 as it did in sensitive tumors because CD8 + CD44 + T cells did not significantly expand despite reduction in Tregs.

Genetic algorithms

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

Goodness of fit

Descriptions of model parameters are shown in Table 1. A couple of sample sets of estimated values of parameters obtained through fitting predictions of the impulsive ODE model 1) - 18) to data from a group of experiments in [6] are listed in Table 2. Figure 3a illustrates the comparison between model solutions and experimental measurements on tumor control and immunological protection against cancer cells in animals treated with IL-12 plus OXP. Experimental results for long-term management of colorectal cancer by observing cooperation of IL-12 and OXP for the control of experimental relapses in distant locations are compared against model predictions in Fig. 3b, c, and d for responders, partial-responders and non-responders, respectively. Trajectories of tumor growth, IFN and ratio of T to T 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 , a second tumor free equilibrium 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 when proliferation rate of tumor cells is higher than natural death rate of tumor cells (i.e., k>k).

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

(Based on the Descartes’ rule, there is only one positive solution from the equation), , T satisfies the following quadratic equation which has only one positive solution , , , , and

By simple calculation, thirteen of the fifteen eigenvalues of the Jacobian matrix of linearized system at the equilibrium are given by {−k, −k, −b23+r2, −k, −k, −k−c21−c23, −k, −k, −b31−k, −k, −k, −k−k, and −k+k} with the rest of the eigenvalues satisfying the quadratic equation λ2+(k+a12+a21+a23)λ+a12(k+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 is stable if both r2APC2 growth rate is smaller than the rate constant for APC2 flowing from blood to tumor) and ktumor proliferation rate is less than tumor natural death rate), otherwise it is unstable.

The second tumor free equilibrium 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 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 b23−r2, and the other two satisfies the following quadratic equation

It is found that both eigenvalues are in the left-half of the complex plane if ktumor-free equilibrium point is locally stable if r2>b23 and k

The high tumor equilibrium exists when k>k. Twelve of the fifteen eigenvalues of the Jacobian matrix are found as {−k,−k,−k,−k,−b23+r2,−k−c21−c23,−k,−b31−k,−k,−k,−k−k, and k−k}.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 high-tumor equilibrium is unstable if either r2>b23 or k>k is satisfied. Suppose r20, 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 r2APC2 growth rate is smaller than the rate constant for APC2 flowing from blood to tumor) and ktumor proliferation rate is larger than tumor natural death rate), otherwise it is unstable.

The stability conditions for tumor free equilibrium indicate that small tumors may not grow into a threatening size when tumor proliferation rate is less than tumor natural death rate (kOXP, IL and T, i=1,2,3 in be to zero). Under the combined IL-12 and OXP treatment, MHC class I positive tumor cells will ultimately be eliminated. Depending on effects of the combined treatment (reflected by the remaining level of APC and T, i=1,2,3 and whether r2>b23), MHCI negative tumor cells () will either be completely removed in which case solutions of this dynamic approach the second tumor-free equilibrium 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 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 is the predicted tumor size in mm3 using 0.1×α and other parameters in the ith calibrated parameter set and is the predicted tumor size in mm3 using α and other parameters in the i calibrated parameter set and α is the calibrated value for parameter α in the i 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 T production rate constants, respectively), g10 - g14 (IL-12, IFN, OXP3, T, and 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 k, 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: killing by OXP3 saturation rate constant g15, natural death rate constant of naive T cells k, killing rate constant of T effectors or tumor cells by OXP k, APCs natural death rate constant k, and natural decay rate constant of OXP k (see Fig. 4).

Model simulations

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

Partial-responders

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

Non-responders

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

Fig. 6

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

Discussion

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

Conclusion

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