Using immune checkpoint modulators in the clinic to increase the number and activity of cytotoxic T lymphocytes that recognize tumor antigens can prolong survival for metastatic melanoma. Yet, only a fraction of the patient population receives clinical benefit. In short, these clinical trials demonstrate proof-of-principle but optimizing the specific therapeutic strategies remains a challenge. In many fields, CAD (computer-aided design) is a tool used to optimize integrated system behavior using a mechanistic model that is based upon knowledge of constitutive elements. The objective of this study was to develop a predictive simulation platform for optimizing anti-tumor immunity using different treatment strategies.

Methods

To better understand the therapeutic role that cytotoxic CD8 ^{+} T cells can play in controlling tumor growth, we developed a multi-scale mechanistic model of the biology using impulsive differential equations and calibrated it to a self-consistent data set.

Results

The multi-scale model captures the activation and differentiation of naïve CD8 ^{+} T cells into effector cytotoxic T cells in the lymph node following adenovirus-mediated vaccination against a tumor antigen, the trafficking of the resulting cytotoxic T cells into blood and tumor microenvironment, the production of cytokines within the tumor microenvironment, and the interactions between tumor cells, T cells and cytokines that control tumor growth. The calibrated model captures the modest suppression of tumor cell growth observed in the B16F10 model, a transplantable mouse model for metastatic melanoma, and was used to explore the impact of multiple vaccinations on controlling tumor growth.

Conclusions

Using the calibrated mechanistic model, we found that the cytotoxic CD8 ^{+} T cell response was prolonged by multiple adenovirus vaccinations. However, the strength of the immune response cannot be improved enough by multiple adenovirus vaccinations to reduce tumor burden if the cytotoxic activity or local proliferation of cytotoxic T cells in response to tumor antigens is not greatly enhanced. Overall, this study illustrates how mechanistic models can be used for in silico screening of the optimal therapeutic dosage and timing in cancer treatment.

Background

Cytotoxic CD8 ^{+} T cells are important effectors in the adaptive immune response against intracellular pathogens and play an important role in immunosurveillance against malignancy [1, 2]. Modulating an immune checkpoint to increase cytotoxic T lymphocytes (CTLs) that target malignant cells can cure patients of metastatic melanoma [3, 4]. While this clinical success demonstrates proof-of-principle, the clinical response is limited to a subset of patients. Yet, these results encourage alternative approaches to direct host immunity against tumors, including adoptive transfer of autologous T cells extracted from a patient’s own tumor (e.g., [5]), engineering of T cell receptors to recognize tumor antigens (e.g., [6, 7]), or vaccination against tumor antigens (e.g., [8-11]). Cancer vaccines based on patient-specific material is attractive as it would enable personalized treatments that enhance CTL response to the specific antigens expressed by a patient’s tumor [12]. One approach is to use adenoviruses that were initially developed as vehicles for gene therapy. Attempts to replace missing or faulty genes by adenoviral gene transfer were largely unsuccessful in experimental animals and human volunteers alike due to innate and adaptive immune responses induced by the adenoviral antigens ([13]). Replication-deficient adenovirus vectors have been pursued as vaccine carriers in the clinic as they showed high efficiency in some rodent and simian preclinical models [13, 14]. The profile of the immune response elicited by adenovirus vaccines against tumor antigens in murine models was investigated by some research groups (see [15-17]). While the approach seems promising, the results are suboptimal as similarly observed for the immune checkpoint modulators. In exploring one treatment variation, sequential treatments involving adenovirus and oncolytic viruses may lead to improved antitumor response [18]. However a more systematic approach to explore treatment variants may be helpful to improve overall response.

As illustrated by [19-22], a variety of mathematical models based on ordinary differential equations (ODEs) have been developed to better understand cancer progression and response to immunotherapy in the last couple of decades. Early work employed Lotka-Volterra equations to describe the interactions between tumor and the immune system where effector cells acted as predators and tumor cells as prey ([21, 23]). The immune surveillance phenomena was described qualitatively in [23] where low doses of tumor cells can escape immune defenses and grow into a larger tumor whereas larger doses of tumor cells are eliminated. The simple predator-prey model was generalized by Kirschner [21, 24], de Pillis et al. [25], Eftimie et al. [26], Wilson and Levy [27], and other researchers where different components of the immune system, such as particular cytokines or natural killer cells, were introduced into the model depending on different cancer treatment strategies. The effect of time delay in the immune response was considered in [20, 28, 29] where authors found that impact of time delay on tumor growth is almost negligible.

Mathematical modeling and computer simulations can be powerful tools in optimizing therapeutic strategies. Mathematical modeling and simulations can be used to screen in silico parameter regions that seem most promising for optimal timing and dosage of therapy and clinical trials can be focused on those regions [30-32]. In [33], the authors explore how the timing of oral insulin delivery and immunomodulatory drugs can be optimized for maximum effect. Moreover, an in silico approach can suggest targeted experiments and then minimize the number of needed experiments [34]. It can also be applied to combine in a virtual way different modes of actions that are well characterized in isolation, such as immunotherapy and chemotherapy, and see how they may be combined to maximum benefit. For instance, Eftimie et al. explored how vaccination using two different viruses that carry the same tumor antigen achieves a greater therapeutic response than if one virus is used alone [26]. In this paper, we use simulations to investigate the impact of multiple adenovirus vaccinations on CD8 ^{+} T cell proliferation and recruitment to the tumor microenvironment and to identify important parameter ranges that control tumor growth through vaccination-induced anti-tumor immunity.

The structure of this paper is as follows. First, we present a multi-scale mechanistic model of anti-tumor immunity and tumor growth based on a set of coupled impulsive ODEs. Second, we describe how we calibrated the parameters of the model 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. Finally, we used the simulation platform to explore the impact of multiple adenovirus vaccinations on T cell proliferation and recruitment to the tumor microenvironment to control tumor growth.

Methods

Here, we developed a multi-scale impulsive ODE model based on our mechanistic understanding of underlying biology and calibrated the model using existing experimental data. This multi-scale mathematical model represents the cytotoxic T cell response to adenovirus vaccination against a tumor antigen and subsequent control of the growth of B16F10 tumors. For the reported experiments, the B16F10 cell line was purchased from American Tissue Culture Collection (ATCC, Bethesda, MD). Numerical solutions of the model were obtained using simulators generated by C Sharp. Simulations start on day 0, the time of tumor implantation and conclude on day 49. 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. A genetic algorithm was used to find parameter sets that closely match the experimental data [15, 16]. 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 parameter set. The fitness function value, or variance, was calculated using a linear combination of sum of error squared and sum of differences between slopes of lines of experimental data and corresponding model predictions. The calibrated mechanistic model was then used to investigate the long-term behavior through stability analysis. Finally, we used the calibrated model to explore the impact of multiple vaccinations on tumor growth to improve anti-tumor immunity, a scenario that is difficult to test experimentally using pre-clinical mouse models but could be potentially used in the clinic to treat patients. Details of model development, parameter calibration, stability analysis, and numerical simulations of multiple vaccinations are described in the following sections.

Results

A multi-scale model of CD8 ^{+} T cell control of tumor growth

Our mathematical model is based on the experimental data presented by Bramson and coworkers [15, 16] using the B16F10 model for metastatic melanoma. The B16F10 model is one of a number of transplantable models of cancer that have been used as pre-clinical models to test anti-tumor immunotherapies [35]. In these models, a malignant cell line derived from a spontaneous mouse cancer is cultured in vitro and then injected back into syngeneic immunocompetent mice. The B16 model was used to help demonstrate the efficacy of anti-CTLA4 therapy [36], a drug that has revolutionized the treatment of metastatic melanoma in humans [4, 37]. As an alternative approach to enhance anti-tumor immunity, Bramson and coworkers examined how a cytotoxic CD8 ^{+} T cell response directed against a tumor antigen using a recombinant adenovirus vector can help control tumor growth. This adenovirus vector induces both the transient expression of a defined tumor antigen and triggers innate immunity to initiate a primary adaptive immune response against this tumor antigen. A primary adaptive immune response is organized spatially: presentation of antigen and initial activation of naïve CD8 ^{+} T cells occurs in the secondary lymphoid organs, activated effector T cells circulate in the blood and peripheral tissues in search of tumor antigens, and effector T cells remain in tissues that express tumor antigens and selectively kill cells that express the cognate tumor antigen [38]. In order 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 mode to quantify the cytotoxic CD8 ^{+} T cell response to adenovirus vaccination and subsequent inhibition of tumor cell growth, as shown schematically in Fig. 1. Among these three compartments, we consider the dynamics of nine 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}). As the immunization protocol induces the clonal expansion of small subset of CD8 ^{+} T cell clones rather than globally changing T cell numbers, we assumed that naïve CD8 ^{+} T cells expressing the T cell receptor that recognizes the epitope derived from the immunized tumor antigen are produced at a constant rate c_{1} from thymus and die naturally at a rate k_{
d1}T_{
N
}. Naïve CD8 ^{+} T cells are maintained at a constant level in the absence of adenovirus, i.e., c_{1}=k_{
d1}T_{
N
}(0). Naïve CD8 ^{+} T cells are recruited to the lymph node and activated by adenovirus vaccination and become effector CD8 ^{+} T cells (T_{
E1}) when they encounter adenovirus-induced antigen expression (LV) at a rate proportional to T_{
N
} and a saturable adenovirus-induced antigen (LV) term defined by \(\frac {\text {LV}}{\text {LV}+\gamma }\).

Effector CD8^{+} T cells in lymph node (T_{
E
1
}, 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} \frac {T_{N}Vol_{b}}{Vol_{\textit {ln}}}\frac {\text {LV}}{\text {LV}+\gamma }\), where Vol_{
b
}=1.4∗10^{3}mm^{3} is the volume of the blood compartment ([39]) and Vol_{
ln
}=0.25 mm^{3} is the volume of the lymph node compartment ([40]). 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 adenovirus-induced antigen term defined by \(\frac {\text {LV}}{\text {LV}+\gamma }\), 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}. We also assume that influx rate of effector CD8 ^{+} T cells from blood to lymph node is \(a_{21}\frac {T_{E2}Vol_{b}}{Vol_{\textit {ln}}}\) and a_{12}·T_{
E1} is the efflux rate.

Adenovirus in lymph node (LV, units: Relative Light Units (RLU) per mm^{3}). Since the adenovirus used in the calibration experiments are replicate-defective and include a GFP expression plasmid, we assume an exponential decay model for LV. We also used a difference equation \(\Delta \text {LV}(t)=\text {LV}(t^{+})-\text {LV}(t^{-})=\text {LV}_{k}\) to reflect the abrupt change of the concentration of adenovirus during vaccination at time t_{
k
}, where LV_{
k
} represents the dosage of vaccination at t_{
k
} with k=1,2,…,n.

Effector CD8^{+} T cells in blood (T_{
E
2
}, units: cells per mm^{3}). We assume the effector CD8 ^{+} T cells die naturally at a rate k_{
d3}T_{
E2} in blood. The influx rate of effector CD8 ^{+} T cells from lymph node to blood is equal to \(a_{12}\frac {T_{E1}Vol_{\textit {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}\frac {C_{MHCI^{-}}}{\epsilon +C(t)}\frac {T_{E3}Vol_{t}}{Vol_{b}}\) and the efflux rate of CD8 ^{+} T effectors from blood to tumor is a_{23}T_{
E2}, where \(C(t)=C_{MHCI^{-}}+C_{MHCI^{+}}\phantom {\dot {i}\!}\) is the number of tumor cells, \(C_{MHCI^{+}}\phantom {\dot {i}\!}\) is the number of major histocompatibility complex (MHC) class I positive tumor cells, \(\phantom {\dot {i}\!}C_{MHCI^{-}}\) is the number of MHC class I negative tumor cells, and ε is a small positive constant representing a small volume of tissue that excludes tumor and effector CD8 ^{+} T cells in the tumor compartment.

MHC class I positive tumor cells (\(\phantom {\dot {i}\!}C_{MHCI^{+}}\), units: cell number). MHC class I positive tumor cells are converted from MHC class I negative tumor cells (\(C_{MHCI^{-}}\phantom {\dot {i}\!}\)) with the assistance of Interferon γ (IFNγ) at a rate \(c_{3}\frac {\text {IFN}{\gamma }}{k_{1} +\text {IFN}{\gamma }}C_{MHCI^{-}}\) and the effector CD8 ^{+} T cell-mediated MHC class I positive tumor cells death rate is \(c_{4} T_{E3}\frac {C_{MHCI^{+}}}{\epsilon +C(t)}\phantom {\dot {i}\!}\). We assume that the dilution rate of MHC class I positive tumor cells due to proliferation is \(k_{p2} C_{MHCI^{+}}\phantom {\dot {i}\!}\). The natural death rate of MHC class I positive tumor cells is assumed to be \(k_{d4}C_{MHCI^{+}}\phantom {\dot {i}\!}\).

MHC class I negative tumor cells (\(\phantom {\dot {i}\!}C_{MHCI^{-}}\), units: cell number). MHC class I negative tumor cells are converted to MHC class I positive tumor cells with the assistance of Interferon gamma (IFN_{
γ
}) at a rate \(c_{3}\frac {\text {IFN}{\gamma }}{k_{1} +\text {IFN}{\gamma }}C_{MHCI^{-}}\). We assume that the proliferation rate of MHC class I positive tumor cells is equal to \(2k_{p2} C_{MHCI^{+}}\phantom {\dot {i}\!}\). 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 vaccination treatment.

Effector CD8^{+} T cells in tumor microenvironment (T_{
E
3
}, units: cells per mm^{3}). We assume that effector CD8 ^{+} T cells can proliferate locally upon recognition of the corresponding tumor antigen presented by MHCI positive tumor cells at a saturable rate equal to \(k_{p3}\frac {C_{MHCI^{+}}}{\epsilon +C(t)}T_{E3}\phantom {\dot {i}\!}\). Effector CD8 ^{+} T cells have a finite lifespan and die within the tumor microenvironment as a rate equal to k_{
d5}·T_{
E3}. The influx rate of effector CD8 ^{+}T cells from the blood to tumor is defined by \(a_{23}\frac {T_{E2}Vol_{b}}{Vol_{t}}\), where Vol_{
t
}=ε+s_{
t
}C(t)+V_{
i
}T_{
E3}mm^{3} 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, s_{
t
}=6∗10^{−7}mm^{3} is the average size of a B16F10 tumor cell ([41]), and V_{
i
}=10^{−7}mm^{3} is the average size of a T effector cell ([42]). The efflux rate of effector CD8 ^{+} T cells from the tumor to blood is \(a_{32}T_{E3}\frac {C_{MHCI^{-}}}{\epsilon +C(t)}\).

Interferon gamma (IFNγ, units: moles per mm^{3}). We assume that Interferon γ is secreted solely by effector CD8 ^{+} T cells within the tumor at a rate proportional to the concentration of effector CD8 ^{+} T cells within the tumor microenvironment and decays at a rate proportional to its concentration. While this assumption may not hold in all model systems, the presence of IFN γ in the tumor was dependent on CD8 ^{+} T cell activation [43].

Tumor Necrosis Factorα (TNFα, units: moles per mm^{3}). We assume that Tumor Necrosis Factor α decays naturally at a rate proportional to its concentration and is secreted solely by effector CD8 ^{+} T cells in the tumor at a rate that includes both autocrine and constitutive production terms: \(\left (k_{c2}\frac {\text {TNF}{\alpha }}{k_{2}+\text {TNF}{\alpha }}+k_{3}\right)T_{E3}\). While this assumption may not hold in all model systems, the presence of TNF α in the tumor was also dependent on CD8 ^{+} T cell activation [43].

Based on the governing biological processes and assumptions described above, the dynamics of these cytotoxic T cell, tumor cell, and cytokine state variables are represented by a mass-action formalism and encoded by the following impulsive ordinary differential equations:

$$ \begin{aligned} \Delta \text{LV}(t)=\text{LV}_{k}, \;\;\;t=t_{k},\;\;k=1,2, 3,\cdots, n, \end{aligned} $$

((10))

where ΔLV(t)=LV(t^{+})−LV(t^{−}) reflects the abrupt change of adenovirus at vaccination time t and LV_{
k
} is the dosage of the adenovirus vaccination at the administration time t_{
k
} with k=1,2,3,⋯,n.

Model calibration

Next, we calibrated model parameters against self-consistent experimental data. These data were acquired from two papers. The first paper described the general dynamics of a CD8 ^{+} T cell response to vaccination with a recombinant human adenovirus serotype 5 (rHuAd5) vector that can be used as a general delivery vehicle to express human tumor antigens [16]. The second paper describes using this adenovirus vector to induce a CD8 ^{+} T cell response to the human dopachrome tautomerase antigen (hDCT; vector: rHuAd5-hDCT) [15]. In contrast, the same adenovirus vector engineered to vaccinate against the glycoprotein gp100 (rHuAd5-hgp100) was unable to control the growth of B16F10 in prophylactic and neo-adjuvant settings. The B16F10 cell line exhibits a defect in the processing and presentation of peptides derived from gp100 through the Major Histocompatibility Complex class I pathway [44]. Together these results suggest that the control of tumor growth induced by rHuAd5-hDCT is through tumor-specific CD8 ^{+} T cells.

The experimental data used in calibrating the mathematical model are listed as follows:

CD8 ^{+} T cells in the secondary lymph nodes (T_{
E1}) and effector CD8 ^{+} T cells in the blood (T_{
E2}) are obtained from Figure 1(A) of Yang’s paper ([16]).

Antigen expression derived from adenovirus vaccination (LV) corresponds to data presented in Figure 3(B) of Yang’s paper ([16]).

Total volume of B16F10-derived tumors was calibrated against data shown in Figure 1(B) of McGray’s paper ([15]).

The concentration of effector CD8 ^{+} T cells present within the tumor (T_{
E3}) are found in Figure 4(A) of McGray’s paper ([15]).

Expression of Interferon gamma (\(\overline {\text {IFN}} \gamma \)) and Tumor Necrosis Factor alpha (\(\overline {\text {TNF}}\alpha \)) genes within the tumor are obtained from Figure 1(E) of McGray’s paper ([15]).

As there are more data points (93) than parameters (27) parameters, the mechanistic model is identifiable in theory.

Simulation results for the modeled variables along with their corresponding experimental measurements are shown in Fig. 2, where t_{0}=0 is the day of tumor inoculation, t=5 is the day of adenovirus immunization. The starting concentration of naïve CD8 ^{+} T cells (\(T_{N}(0)=0.0714 \,\text {cells}\,\text {per}\,\text {mm}^{3}\)) was estimated by assuming 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}). Initially, 2×10^{6} tumor cells were injected into mice and assumed to not express tumor antigens (\(C_{MHCI^{-}}\phantom {\dot {i}\!}\)). The remaining state variables were initially zero. Vaccination was simulated by abruptly changing the concentration of adenovirus (LV) at the administration time (t_{1}=5, the 5th day after tumor implantation) using an impulse dose equal to LV_{1} (\(\text {LV}_{1}=1.100\times 10^{6} \;\text {RLU}\,\text {per}\, \text {mm}^{3}\)). The calibrated parameter values obtained using the genetic algorithm are listed in Table 1.

Stability analysis

The dynamics of the nonlinear ODE model comprised of equations (1) - (9) are complicated. To gain insight into the behavior of the system, we explored the steady state solutions of the system using stability analysis. By setting the right hand sides of the ODE system (1) - (9) to 0 and solving the equations simultaneously, we notice that the system of ODEs (1) - (9) only has a tumor-free equilibrium \(\overrightarrow {X}_{0}\):

We note that the tumor-free equilibrium has only one none-zero element: the naïve T cells T_{
N
}, this occurs when there are no tumor cells present and no adenovirus immunization treatment is administered and also corresponds to tumor-specific effector CD8 ^{+} T cells and cytokines being equal to zero. The high tumor equilibrium \(\overrightarrow {X}_{1}\) has two non-zero elements: the naïve T cells T_{
N
} and the MHC class I negative tumor cells \(\phantom {\dot {i}\!}C_{\textit {MHCI}}^{-}\), which reflects the status of the steady state when the one-time adenovirus vaccination treatment failed to completely eradicate the tumor cells. This situation occurs when adenovirus LV decays to zero and the MHC class I positive tumor cells are all killed by tumor infiltrating lymphocytes, which causes exhaustion of effector CD8 ^{+} T cells in three compartments and cytokines decay to zero. The rest of the MHC class I negative tumor cells then approach the carrying capacity and the naïve T cells return to their original constant level.

By simple calculation, we obtain the Jacobian matrix \(J(\overrightarrow {X})\) of the ODE system (1)- (9):

where J_{44}=−k_{
d3}−a_{21}−a_{23} and J_{55}=−k_{
p2}−k_{
d4}.

It can be shown that \(J(\overrightarrow {X}_{0})\) has the following eigenvalues: λ_{1}=−k_{
d1}, λ_{2}=−k_{
d2}, λ_{3}=−k_{
d5}, λ_{4}=−k_{
d6}, λ_{5}=−k_{
d7}, λ_{6}=−k_{
p2}−k_{
d4}, λ_{7}=k_{
p2}−k_{
d4}, while λ_{8} and λ_{9} are determined by the quadratic equation

In both cases, λ_{8} and λ_{9} have negative real parts.

Thus when k_{
p2}>k_{
d4}, the tumor-free equilibrium \(\overrightarrow {X}_{0}\) is unstable and when k_{
p2}<k_{
d4}, the tumor-free equilibrium \(\overrightarrow {X}_{0}\) is stable since all eigenvalues of the Jacobian matrix have negative real parts.

The Jacobian matrix evaluated at the high tumor equilibrium \(J(\overrightarrow {X}_{1})\) is given by

It can be shown that \(J(\overrightarrow {X}_{1})\) has the following eigenvalues: λ_{1}=−k_{
d1}, λ_{2}=−k_{
d2}, λ_{3}=−k_{
d6}, λ_{4}=−k_{
d7}, λ_{5}=−k_{
p2}−k_{
d4}, λ_{6}=−(k_{
p2}−k_{
d4}), and λ_{7}, λ_{8} and λ_{9} are determined by

where a_{2}=a_{12}+a_{21}+a_{23}+a_{32}+k_{
d3}+k_{
d5}, a_{1}=a_{12}(k_{
d5}+a_{32}+k_{
d3}+a_{23})+k_{
d3}(k_{
d5}+a_{32})+a_{21}(k_{
d5}+a_{32})+k_{
d5}a_{23}, and a_{0}=a_{12}(k_{
d3}k_{
d5}+a_{32}k_{
d3}+a_{23}k_{
d5}). It is easy to see that a_{
i
}>0 for i=0,1,2 since all parameters are positive. By simple calculation, we can obtain a_{2}a_{1}−a_{0}>0 which implies that roots of equation (13) (λ_{7}, λ_{8}, λ_{9}) all have negative real parts by the Routh-Hurwitz criterion. Thus when the proliferation rate of tumor cells is greater than the natural death rate of tumor cells (i.e., k_{
p2}>k_{
d4}), the high tumor equilibrium \(\overrightarrow {X}_{1}\) is stable and when the proliferation rate of tumor cells is less than the natural death rate of tumor cells (i.e., k_{
p2}<k_{
d4}), the high-tumor equilibrium \(\overrightarrow {X}_{1}\) is unstable.

Therefore, using the parameter values obtained from model calibration (k_{
p2}>k_{
d4}), the tumor-free equilibrium \(\overrightarrow {X}_{0}\) is unstable and the high tumor equilibrium \(\overrightarrow {X}_{1}\) is stable. This implies that under the current status of the mouse immune system, a small tumor will keep growing to its carrying capacity because of the fast proliferation of tumor cells without adenovirus vaccination treatment. On the other hand, the one-time adenovirus immunization as applied in the experiment was not very successful in completely eliminating tumor cells due to limited effects on enhancing CTL immune response. It seems that the CTL response falls to zero before all MHC negative tumor cells are converted to MHC positive tumor cells and killed by the cytotoxic CD8 ^{+} T cells. Then, MHC positive tumor cells approach zero while MHC negative tumor cells approach the carrying capacity and all T cell effectors and cytokines drop back to zero soon after the one time vaccination treatment (see Fig. 2).

Effect of multiple vaccinations

Next, we investigated in silico the impact of multiple adenovirus vaccinations on T cell proliferation and recruitment, cytokine secretion, and tumor growth using the calibrated model in conjunction with an impulsive control mechanism where control laws are discrete in time, as represented by Eq. (10).

To explore the impact of dose-dependence, adenovirus vaccinations were simulated with an increasing dosage of 1×10^{6} RLU per mm^{3}, 1×10^{8} RLU per mm^{3}, and 1×10^{10} RLU per mm^{3} and delivered every 5 days for a total of 10 times following the initial vaccination (see Fig. 3). The simulation results indicate that with the same frequency and number of times of immunizations, the increase in immunization dose increases the length and the maximum magnitude of the adenovirus vaccinations and CTL response. As we can see in Fig. 3, a 100-fold increase in immunization dose results in a 15-day (or 10 %) increase in the length of the adenovirus vaccinations as well as CTL response. However, the impact of increase in vaccination dose on the strength of CTL response is very limited: a 100-fold increase in immunization dose generates less than a 1 % increase in the collective magnitude of T cell concentration in all three compartments. The increased immunization dosage shows a significant impact on the length but not the collective magnitude of Interferon gamma (IFNG) and Tumor Necrosis Factor alpha (TNF α) in the tumor microenvironment. Surprisingly, the increase in dose with multiple immunizations does not appear to affect the peak time and maximum magnitude of tumor growth. In addition, the time period of tumor staying at its carrying capacity barely changes with increased dose of vaccinations.

Increasing the time period between successive vaccinations was the next variable that we explored (see Fig. 4). While the dose and number of immunization were kept the same, the biological responses were simulated for adenovirus vaccinations that were administered every 5 days, 10 days, 15 days, and 20 days with a dose of 1×10^{6} RLU per mm^{3} for 10 times in total. Simulation results suggest that shorter time periods between successive immunizations generate higher maximum immunizations magnitudes. In addition, immunizations last longer with longer time periods between successive immunizations. Length of T cell responses in all compartments would increase with the increase in time periods between successive vaccinations and then drop to zero in about 60 to 100 days after the last vaccination for all four cases with a smaller duration of response for longer time between successive immunizations. In general, the increased time periods between successive immunizations result in a 4 % to 18 % decrease in the maximum magnitude for T cell concentration in blood and almost no impact on the maximum magnitudes of T cell concentrations in lymph node and tumor microenvironment. The expression of IFNG and TNF α do not seem to be affected although shorter time periods between successive immunizations correspond to shorter durations of both cytokine responses. In summary, simulation results suggest that with the same dosage and number of times of vaccinations, the shorter the time between successive vaccinations, the higher T cell concentration in blood and the shorter the T cell responses in all three compartments. However, the increase of time periods between successive vaccinations does not seem to affect tumor growth.

We also explored the impact of increasing the number of times of immunizations with the same dosage and time period between successive vaccinations (see Fig. 5). Model predictions were obtained with adenovirus vaccinations with dose of 1×10^{10} RLU per mm^{3} every 30 days administered once, 5 times, and 10 times following the initial immunization. We found that, when dosage and time period between successive vaccinations were kept the same, increasing the number of times of immunizations increased the duration of the T cell. In addition, the increase in number of times of vaccinations only increased the maximum T cell concentration in blood by 8 % to 34 % but did not impact the maximum T cell concentrations in lymph node and tumor microenvironment. Figure 5 also shows an interesting result that the time period of the tumor staying in its maximum value and then drop to its carrying capacity is increased with the increase of number of times of multiple immunizations while increasing the length of the T cell immune response (i.e., \(C_{\textit {MHCI}}^{+}\doteq 0\) and \(C_{\textit {MHCI}}^{-}= \frac {k_{p2}-k_{d4}}{r_{2}}\doteq 1.497\times 10^{9}\) in number of tumor cells or equivalently, \(s_{t}*C_{\textit {MHCI}}^{-}\doteq 898.2 \,mm^{3}\) in tumor size).

Impact of the single adenovirus vaccination with enhanced T cell cytotoxicity or proliferation

Finally, we explored the impact of a single adenovirus vaccination with enhanced T cell cytolytic ability or T cell proliferation on tumor growth by changing the corresponding parameters. Model predictions with enhanced T cell cytolytic ability (i.e., c_{4}=2.49×10^{5}) and improved T cell proliferation rate in the tumor compartment (i.e., k_{
p3}=5.73) were obtained and compared to the prediction using calibrated parameters reported in Table 1 (see Fig. 6). The cell concentrations of naïve CD8 ^{+} T cells and adenovirus are not included in the figure as the model predictions were exactly the same for all three cases. With an enhanced T cell cytolytic ability, the concentrations of effector CD8 ^{+} T cells in the lymph node and blood as well as the total number of T cells in the tumor were almost the same relative to simulations using the calibrated parameter values. In contrast, the concentration of effector CD8 ^{+} T cells in the tumor, IFNG and TNF α were increased by 100 to 1000 times compared to the prediction obtained using the c_{4} from model calibration. The increase in CD8 ^{+} T cell response within the tumor decreased the tumor size to effectively zero (i.e., less than the size of a single tumor cell) within 100 days. To check whether tumor would relapse during the average life span (1 to 2 years) of a mouse, we have run the simulation for more than 800 days and found that tumor cells continued to decrease to near zero after 100 days as shown in Fig. 6. Collectively, the simulation results suggest that tumor cells can be completely eliminated by a single adenovirus vaccination with greatly strengthened T cell cytolytic ability (c_{4}). In the model prediction with enhanced T cell proliferation rate k_{
p3}, we see that T cell concentrations in all three compartments were greatly enhanced and tumor cells were reduced to near zero within 100 days. To summarize, the simulation results suggest that the tumor may be completely eradicated by a single adenovirus vaccination by enhancing either the CD8 ^{+} T cell cytolytic ability or the CD8 ^{+} T cell proliferation rate.

Discussion

Mathematical modeling and simulation are increasingly being used in the pharmaceutical industry to better understand the underlying biology targeted by a drug and to explore therapeutic scenarios that may be difficult to test experimentally [45]. Here, we developed a three-compartment mechanistic mathematical model to describe the clonal expansion of CD8 ^{+} T cells in a mouse model of metastatic melanoma in response to adenovirus vaccination against a defined tumor antigen. Based on the collective knowledge of this pre-clinical mouse model, the model represents the primary CD8 ^{+} T cell response to adenovirus immunization and the subsequent impact on the growth of a tumor derived from the B16F10 cell line. Using the mechanistic model as a framework to integrate different experimental studies, model parameters were calibrated against published experimental data that describes the primary response. As shown in Fig. 2, our model predictions of adenovirus concentration, tumor size, concentrations of CD8 ^{+} T effectors in blood and tumor, gene expression of IFNG and TNF α matched the experimental data. The proposed model structure reflects a trade-off between biological realism, parameter identifiability, and a fitness-for-purpose. As mentioned above, an excess of data points (93) relative to the number of parameters (27) suggests that the model is identifiable in theory. Efforts are underway to identify the appropriate topology of the network, given the available data [46].

In terms of the trade-off between biological realism and fitness-for-purpose, we settled on the proposed model structure to facilitate stability analysis. Stability analysis of tumor-free and high tumor equilibria was conducted based on the linearized system. Impulsive stabilization using the Lyapunov method will be considered in the future to provide conditions on parameters such that the high-tumor equilibrium may be stabilized using impulsive control through manipulation of strength and frequency of the multiple vaccinations. However, the proposed model structure imposes some limitations in how the model represents the system and interpreting the model predictions. In particular, we note that the model predicted an earlier peak time for the concentration of effector CD8 ^{+} T cells in the lymph node compared to experimental data, which may suggest a more complicated model structure for the lymph node compartment than proposed here. While additional model structure may help in capturing the dynamics of T cells within the lymph node, the current structure is sufficient to capture the dynamics of CD8 ^{+} T cells within the blood, which is the pool that gets recruited to the tumor compartment. Additional lymph node structure would then have limited impact on our conclusions. We also note that the effective concentration of effector CD8 ^{+} T cells within the tumor compartment (e.g., number of effector CD8 ^{+} T cells per weight of tumor) peaked at 10 days despite the blood population of CD8 ^{+} T cells peaking at day 20. As CD8 ^{+} T cell recruitment from the blood into the tumor compartment was assumed to be independent of tumor size and the parameter values suggest that proliferation of CD8 ^{+} T cells within the tumor was negligible, this decline in CD8 ^{+} T cell concentration was due to dilution of recruited effector CD8 ^{+} T cells into an exponentially growing tumor mass. While a direct measure of tumor infiltrating CD8 ^{+} lymphocytes was reported at a single time point, expression of IFNG and TNF α are implicit surrogate markers for CD8 ^{+} T cell infiltration as these two cytokines are directly proportional to the concentration of CD8 ^{+} T cells within the tumor compartment. Measuring the number of tumor infiltrating lymphocytes in this mouse model at additional time points would help confirm these assumptions. This would be interesting, as tumors, like the B16 model, are known to develop immunosuppressive mechanisms that are proportional to tumor size that could alter the relationship between the presence of tumor infiltrating lymphocytes and cytokine production [47].

Given the rapid growth of the B16F10 model and ethical limitations of animal studies, studies using this pre-clinical mouse model is limited typically to a single round of therapy. Yet, the treatment of human cancers typically involves multiple rounds of therapy to control tumor growth. Here we used a calibrated mechanistic model coupled with computer simulation to explore clinically relevant treatment options in silico. In exploring the impact of multiple vaccinations, our model indicates that increasing the dose of adenovirus vaccination, the time period between successive adenovirus vaccinations, or the number of adenovirus vaccinations results in a prolonged lifespan of effector CD8 + T cells in all three compartments and extended length of secretion of the cytokines IFNG and TNF α within the tumor microenvironment. However these changes in multiple vaccinations have little impact on the magnitude of the clonal CD8 ^{+} T cell immune response and therefore have very little impact on reducing tumor growth. If technically and ethically feasible, additional animal experiments using multiple vaccinations may be helpful to confirm these predictions. As the adenovirus vector promotes clonal expansion of CD8 ^{+} T cells that recognize a small number of epitopes derived from tumor antigens, the number and diversity of effector CD8 ^{+} T cells might not be sufficient to eliminate tumor cells completely. The results are consistent with recent findings in literature, where Budhu et al. reported that CD8 ^{+} T cell concentration determines their efficiency in killing melanoma cells [48]. As reported in [49, 50], immunotherapy of patients with cancer requires the in vivo generation of large numbers of highly reactive anti-tumor lymphocytes that are not restrained by normal tolerance mechanisms and are capable of sustaining immunity against solid tumors. Immunization of melanoma patients with a broader array of cancer antigens can increase the number of circulating effector CD8 ^{+} T cells (eCTLs), but to date this has not correlated with clinical tumor regression, suggesting a defect in function of the eCTLs. In contrast, a clinical benefit has been observed in patients with metastatic melanoma using antibodies against CTLA4, which globally increase the number of circulating CD8 ^{+} T cells irrespective of antigen specificity [4]. The one-time adenovirus immunization experimental data [15, 16] and the computational simulations exploring the efficacy of multiple adenovirus vaccinations using our calibrated model all indicate the limited impact of a single antigen-specific therapy to eliminate tumors. Our simulation results also suggest that increasing the cytotoxic activity of effector CD8 ^{+} T cell or the local proliferation of CD8 ^{+} T cells within the tumor microenvironment, as observed following anti-PD1 therapy [51], may completely eliminate tumor cells.

Conclusions

In summary, we present a multi-scale mechanistic model of CD8 ^{+}-mediated control of tumor growth in response to adenovirus-vaccination mediated T cell stimulation using a system of impulsive ordinary differential equations. The model parameters were calibrated against experimental data employing a genetic algorithm whose fitness function is given by a linear combination of sum of normalized error squared and sum of normalized difference of slopes squared. With the calibrated parameter values, our model predictions match experimental data very well. Stability analysis via linearization implies that, in the case of no vaccination treatment, a small tumor will grow to its carrying capacity as a result of a stable tumor-free equilibrium and a unstable high tumor equilibrium. Using the calibrated model, numerical simulation of multiple adenovirus vaccinations suggest that this treatment strategy will significantly prolong T cell immune response but not necessarily enhance a cytotoxic CD8 ^{+} T cell response to a tumor antigen that noticeably reduces tumor size. A reduction in tumor size can be obtained if the cytotoxic activity or proliferation of effector CD8 ^{+} T cells present within the tumor microenvironment are enhanced. Along those lines, simulation results also show that a tumor may be completely eliminated by a single adenovirus vaccination that creates a highly enhanced cytotoxic T cell efficacy or with enhanced local proliferation of cytotoxic T cells within the tumor. Overall, the results illustrate how mechanistic models can be used to predict tumor growth response to antigen-specific immunotherapies and screen in silico for optimal therapeutic dosage and timing in treating patients with cancer.

Ethics

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

References

Parish IA, Kaech SM. Diversity in cd8(+) t cell differentiation. Curr Opin Immunol. 2009; 21(3):291–7.

Topalian SL, Hodi FS, Brahmer JR, Gettinger SN, Smith DC, McDermott DF, et al. Safety, activity, and immune correlates of anti-PD-1 antibody in cancer. N Engl J Med. 2012; 366(26):2443–454.

Rosenberg SA, Yang JC, Sherry RM, Kammula US, Hughes MS, Phan GQ, et al. Durable complete responses in heavily pretreated patients with metastatic melanoma using T-cell transfer immunotherapy. Clin Cancer Res. 2011; 17(13):4550–557.

Robbins PF, Morgan RA, Feldman SA, Yang JC, Sherry RM, Dudley ME, et al. Tumor regression in patients with metastatic synovial cell sarcoma and melanoma using genetically engineered lymphocytes reactive with NY-ESO-1. J Clin Oncol. 2011; 29(7):917–24.

Porter DL, Levine BL, Kalos M, Bagg A, June CH. Chimeric antigen receptor-modified T cells in chronic lymphoid leukemia. N Engl J Med. 2011; 365(8):725–33.

Kenter GG, Welters MJ, Valentijn AR, Lowik MJ, Berends-van der Meer DM, Vloon AP, et al. Vaccination against HPV-16 oncoproteins for vulvar intraepithelial neoplasia. N Engl J Med. 2009; 361(19):1838–47.

McGray AJ, Bernard D, Hallett R, Kelly R, Jha M, Gregory C, et al. Combined vaccination and immunostimulatory antibodies provides durable cure of murine melanoma and induces transcriptional changes associated with positive outcome in human melanoma patients. Oncoimmunology. 2012; 7(4):419–31.

Yang TC, Millar J, Groves T, Grinshtein N, Parsons R, Takenaka S, et al. The CD8 ^{+} t cell population elicited by recombinant adenovirus displays a novel partially exhausted phenotype associated with prolonged antigen presentation that nonetheless provides long-term immunity. J Immunol. 2006; 176(1):200–10.

Yang TC, Millar J, Groves T, Zhou W, Grinshtein N, Parsons R, et al. On the role of cd4 ^{+} t cells in the cd8 ^{+} t-cell response elicited by recombinant adenovirus vaccines. Mol Ther. 2007; 15(5):997–1006.

Bridle BW, Stephenson KB, Boudreau JE, Koshy S, Kazdhan N, Pullenayegum E, et al. Potentiating cancer immunotherapy using an oncolytic virus. Mol Ther. 2010; 18(8):1430–9.

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

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

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.

Kuznetsov V, Makalkyn M, Taylor M, Perelson A, Chaplain MA. Nonlinear dynamics of immunogenic tumors: parameter estimation and global bifurcation analysis. Bull Math Biol. 1994; 56(2):295–321.

Eftimie R, Dushoff J, Bridle B, Bramson JL, Earn DJD. Multi-stability and multi-instability phenomena in a mathematical model of tumor-immune-virus interactions. Bull Math Biol. 2011; 73(12):2932–961.

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

Bagheri N, Shiina M, Lauffenburger DA, Korn WM. A dynamical systems model for combinatorial cancer therapy enhances oncolytic adenovirus efficacy by mek-inhibition. PLoS Comput Biol. 2011; 17(2):10.

Gadkar KG, Shoda LK, Kreuwel HT, Ramanujan S, Zheng Y, Whiting CC, et al. Dosing and timing effects of anti-cd40l therapy: predictions from a mathematical model of type 1 diabetes. Ann N Y Acad Sci. 2007; 1103:63–8.

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

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

van Elsas A, Hurwitz AA, Allison JP. Combination immunotherapy of B16 melanoma using anti-cytotoxic T lymphocyte-associated antigen 4 (CTLA-4) and granulocyte/macrophage colony-stimulating factor (GM-CSF)-producing vaccines induces rejection of subcutaneous and metastatic tumors accompanied by autoimmune depigmentation. J Exp Med. 1999; 190(3):355–66.

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–885.

Moy BC, Petzold GL, Badiner GJ, Kelly RC, Aristoff PA, Adams EG, et al. Characterization of B16 melanoma cells resistant to the CC-1065 analogue U-71,184. Cancer Res. 1990; 50(8):2485–492.

McGray AJ, Hallett R, Bernard D, Swift SL, Zhu Z, Teoderascu F, et al. Immunotherapy-induced CD8+ T cells instigate immune suppression in the tumor. Mol Ther. 2014; 22(1):206–18.

Leitch J, Fraser K, Lane C, Putzu K, Adema GJ, Zhang QJ, et al. CTL-dependent and -independent antitumor immunity is determined by the tumor not the vaccine. J Immunol. 2004; 172(9):5200–205.

Wen FT, Thisted RA, Rowley DA, Schreiber H. A systematic analysis of experimental immunotherapies on tumors differing in size and duration of growth. Oncoimmunology. 2012; 1(2):172–8.

Budhu S, Loike JD, Pandolfi A, Han S, Catalano G, Constantinescu A, et al. Cd 8^{+} t cell concentration determines their efficiency in killing cognate antigen-expressing syngeneic mammalian cells in vitro and in mouse tissues. J Exp Med. 2010; 207(81):223–35.

Dudley ME, Wunderlich JR, Robbins PF, Yang JC, Hwu P, Schwartzentruber DJ, et al. Cancer regression and autoimmunity in patients after clonal repopulation with antitumor lymphocytes. Science. 2002; 298:850–4.

Rosenberg SA, Yang JC, Schwartzentruber DJ, Hwu P, Marincola FM, Topalian SL, et al. Immunologic and therapeutic evaluation of a synthetic peptide vaccine for the treatment of patients with metastatic melanoma. Nat Med. 1998; 4:321–7.

Tumeh P, Harview C, Yearley J, Shintaku I, Taylor E, Robert L, et al. PD-1 blockade induces responses by inhibiting adaptive immune resistance. Nature. 2014; 515:568–71.

This work was supported by grants from the National Science Foundation (CAREER 1053490 to DJK), the National Cancer Institute (R15CA132124 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. The authors thank Logan Lyda and Darryl Johnson for assistance with data analysis of the simulation results.

Author information

Author notes

Authors and Affiliations

Department of Computer Sciences, Mathematics, and Engineering, Shepherd University, Shepherdstown, 25443, WV, USA

Qing Wang & Zhijun Wang

Department of Chemical Engineering and Mary Babb Randolph Cancer Center, West Virginia University, Morgantown, 25606, WV, USA

David J Klinke II

Department of Microbiology, Immunology, & Cell Biology, West Virginia University, Morgantown, 25606, WV, USA

The authors declare that they have no competing interests.

Authors’ 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, QW performed stability analysis and analyzed simulation data. All authors approved the final version of the manuscript.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.

The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

The Creative Commons Public Domain Dedication waiver (https://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Wang, Q., J Klinke, D. & Wang, Z. CD8 ^{+} T cell response to adenovirus vaccination and subsequent suppression of tumor growth: modeling, simulation and analysis.
BMC Syst Biol9, 27 (2015). https://doi.org/10.1186/s12918-015-0168-9