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][4][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][9][10][11][12][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][14][15][16][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][21][22][23][24][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][27][28][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 tumorfree 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.

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 preexisting 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.
Among these three compartments, we consider the dynamics of fifteen state variables that are regulated by the following governing biological processes and assumptions:

1). 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 (antigen-presenting cells in lymph node) at a rate c 2 · T N · APC 1 APC 1 +g 1 [32][33][34].
2). 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 Vol ln · APC 1 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 APC 1 APC 1 +g 2 , and an immune checkpoint term defined by α α+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 · T E2 Vol b 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 · T E1 ·OXP 1 OXP 1 +g 3 .
3). 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 · APC 3 · Vol t Vol ln , where APC 3 is the concentration of APCs in tumor, Vol t = +C(t)  During this process, IL-12 (IL) helps promote T cell proliferation and suppresses regulatory T (T R ) cells' proliferation and immunosuppressive action on effector CD8 + T cells. In addition, the chemotherapy drug Oxaliplatin in the lymph node and tumor (OXP i where i = 1, 3) will kill fast-proliferating cells such as T effectors and tumor cells 4). 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 · OXP 2 Vol b Vol ln , where OXP 2 is the concentration of OXP in blood.

5). 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 · T E1 Vol ln 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 · , 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 .

6). Antigen Presenting Cells in blood (APC 2 , units: cells per mm 3 ). A logistic growth pattern
Carrying capacity of APC (in blood) Regulatory T cell saturation rate constant 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.
. 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).

8). 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 IL-12 stimulation and subject to suppression from T regulatory cells at a saturable rate equal to k p2 · C MHCI + +C(t) · IL IL+g 4 · T E3 T R +g 5 , where IL is the concentration of IL-12 and T R is the concentration of regulatory T cells [39,40]. The influx rate of effector CD8 + T cells from the blood to tumor is defined by a 23 · T E2 Vol b Vol t . The efflux rate of effector CD8 + T cells from the tumor to blood is a 32 · T E3 · C MHCI − +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 · T E3 · OXP 3 OXP 3 +g 6 in the tumor microenvironment.

9). 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 IL-12 and inhibition from regulatory T cellsat a rate of c 4 · IL IL+g 7 · 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 .
10). 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 · APC 2 · Vol b 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 .
11). Interleukin-12 (IL, units: ng per ml). Interleukin-12 (IL-12) is produced by APCs at a rate of c 5 · APC 3 and decays naturally at a rate of k d8 · IL. The extra IL-12 expression obtained through the combined therapy is approximated using the discrete Eq. (16).
12). 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

13). 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 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 · C MHCI + . The natural death rate of MHC class I positive tumor cells is assumed to be k d10 · C MHCI + and MHC class I positive tumor cells are killed by chemotherapy agent OXP in tumor at a rate k d2 · C MHCI + · OXP 3 OXP 3 +g 14 .
We assume that the proliferation rate of MHC class I positive tumor cells is equal to 2 · k p4 · 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 k d10 · C MHCI − and these cells are killed by chemotherapy agent OXP in tumor at a rate k d2 · C MHCI − · OXP 3 OXP 3 +g 15 . The difference in tumor volume caused by tumor re-challenge is described by the discrete Eq. (18).
16). 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).
A schematic diagram summarizing this threecompartment model is shown in Fig. 1. Model parameters and their meanings are listed in Table 1.

Non-negativity of solutions to the model
For any mathematical model that has biological implications, it is important to make sure that solutions with non-negative initial conditions remain non-negative all the time. For the model comprised of Eqs. 1) -18), we can see that, with positive impulsive inputs IL k1 , {OXP 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 non-negative initial conditions will remain nonnegative because dx i dt ≥ 0 for x i = 0 and x j ≥ 0, where i, j = 1, 2, · · · , 15 and i = j.

Model calibration
Therapeutic use of IL-12 requires efficient methods to control the plasma concentration of this potent immunostimulatory cytokine in order to avoid toxicity [6]. It was determined in an MC38 syngeneic tumor model that a blood concentration of IL-12 less than 20 ng/ml has no anti-tumor effect, while concentrations greater than 700 ng/ml are associated with toxicity [43]. For comparison, the normal range of IL-12 in humans is around 7.5 pg/ml [44]. Gonzalez-Aparicio and colleagues designed a new induction protocol to keep IL-12 within this therapeutic range [6]. Once the liver of a group of C57BL/6 mice was transduced with the vector (typically 2.5*10 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 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 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). • Concentration of Interferon gamma was obtained from Fig. 4b (IL-12 + OXP group). • The ratio between CD8 + T lymphocytes and T regulatory cells was obtained from Fig. 5b (experimental results for IL-12 + 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 re-challenge (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] 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  [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 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 Table 2.
Difference in treatment efficacy: non-responders versus responders While both the responders and non-responders survived two cycles of combination therapy treatment before tumor rechallenge and then underwent the third cycle following the tumor re-challenge on day 75, the simulations in Fig. 3b -d suggest that the non-responders show near zero concentration of IFN γ and near zero ratio of T effectors to regulatory T cells in the tumor all the time comparing to a stable concentration of IFN γ and ratio of T effectors to regulatory T cells in the tumor (at least 10 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][49][50]. In [48] and [49], Random forest (RF) machine learning analysis of tumor infiltrating lymphocytes (TILs) demonstrated that the top predictor of resistance, as measured by variable importance scores and selection, was the CD8+CD44+ to Treg (CD8/Treg) ratio. The authors in [50] indicated that the CD8/Treg ratio in resistant tumors failed to increase after RT + anti-CTLA4 as it did in sensitive tumors because CD8+CD44+ T cells did not significantly expand despite reduction in Tregs.

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

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

Model stability analysis
In this section, we discuss local stability of equilibria of the model using linearized system evaluated at these points. We found that system comprised of equations (1) - (15) has a tumor-free equilibrium X 0 , a second tumor free equilibrium 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 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 , T E3 satisfies the following quadratic equation , and , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, c 6 k d9 , 0, By simple calculation, thirteen of the fifteen eigenvalues of the Jacobian matrix of linearized system at the equilibrium X 0 are given by , 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 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 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 X 1 are in the left-half of the complex plane. It is easy to see that twelve of the fifteen eigenvalues are negative. With respect to the remaining three, one is b 23 − r 2 , 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 k p4 < k d10 . Hence this tumor-free equilibrium point is locally stable if r 2 > b 23 and k p4 < k d10 ; otherwise it is unstable.
The high tumor equilibrium X 2 exists when k p4 > k d10 . Twelve of the fifteen eigenvalues of the Jacobian matrix are found as 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 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 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 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 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 X 0 be to zero). Under the combined IL-12 and OXP treatment, MHC class I positive tumor cells will ultimately be eliminated. Depending on effects of the combined treatment (reflected by the remaining level of APC i and T Ei , i = 1, 2, 3 and whether r 2 > b 23 ), MHCI negative tumor cells (C MHCI − ) will either be completely removed in which case solutions of this dynamic approach the second tumor-free equilibrium X 1 or the MHC class I negative tumor cells eventually approach the carrying capacity. This can occur when MHC class I positive tumor cells are all killed by tumor infiltrating lymphocytes, which results in the exhaustion of effector CD8 + T cells and cytokines while naïve T cells and MHC class I negative tumor cells remaining at constant levels.

Sensitivity of parameters
To test the impact of how the change of a certain parameter value (e.g. α) would affect tumor growth pattern for responders, partial-responders, and non-responders, normalized differences of tumor sizes y i = |ŷ i −ȳ i | y i , (i = 1, 2, · · · , 30) on day 120 (8 days post the second treatment cycle of the combination therapy after tumor re-challenge) were used to draw the violin plots in Fig. 4 for each of the 48 parameters for all three patient groups, whereŷ i is the predicted tumor size in mm 3 using 0.1 × α i and other parameters in the ith calibrated parameter set and 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 non-responders. In addition, there are 10 (out of 48) parameters whose value changes greatly affect tumor size of responders but not the size of non-responders and partial-responders. These potentially OXP and IL-12 treatment important parameters include c 23 (OXP flow rate from blood to tumor), K (APC carrying capacity), c 4c 6 (IFN γ , IL-12, and T R production rate constants, respectively), g 10g 14 (IL-12, IFN γ , OXP 3 , T R , and 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 non-responders to increasing changes in normalized tumor size from partial-responders to responders: T cell flow rates from blood to lymph node and from tumor to blood, a 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 ; IL-12 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 (non-responders, partial-responders, and responders) results from the value changes of following 5 parameters: 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 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.

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  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 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 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 antigenspecific immuno-chemotherapies and screen in silico for optimal therapeutic dosage and timing in treating patients with metastatic colorectal cancer.