Modelling combined virotherapy and immunotherapy: strengthening the antitumour immune response mediated by IL-12 and GM-CSF expression

Combined virotherapy and immunotherapy has been emerging as a promising and effective cancer treatment for some time. Intratumoural injections of an oncolytic virus instigate an immune reaction in the host, resulting in an influx of immune cells to the tumour site. Through combining an oncolytic viral vector with immunostimulatory cytokines an additional antitumour immune response can be initiated, whereby immune cells induce apoptosis in both uninfected and virus infected tumour cells. We develop a mathematical model to reproduce the experimental results for tumour growth under treatment with an oncolytic adenovirus co-expressing the immunostimulatory cytokines interleukin 12 (IL-12) and granulocyte-monocyte colony stimulating factor (GM-CSF). By exploring heterogeneity in the immune cell stimulation by the treatment, we find a subset of the parameter space for the immune cell induced apoptosis rate, in which the treatment will be less effective in a short time period. Therefore, we believe the bivariate...

While the cytotoxic effects of viruses are most commonly viewed in terms of pathogenicity, it is possible to harness their activity for therapeutic purposes. Oncolytic viruses are genetically engineered viruses that selectively infect, self-replicate and lyse cancer cells (Parato, Senger, Forsyth, & Bell, 2005). They may be fortuitously tumour selective in wildtype; however, in most cases, attenuated forms are engineered to provide tumour selectivity. In this article, we will examine the usefulness of a gene-attenuated oncolytic adenovirus genetically engineered by Choi, Zhang, Choi, Kim, and Yun (2012a) to selectively infect and replicate within tumour cells and release the immunostimulatory cytokines interleukin 12 (IL-12) and granulocyte monocyte stimulating factor (GM-CSF). In this way, Choi et al. manufactured a virus that not only eradicates tumour cells, but also activates immune cells to do the same.
Human immune systems have a very powerful and effective way of eliminating an invading pathogen such as a virus. The front line of defence includes three key innate immune cells: macrophages, natural killer cells (NK) and dendritic cells (DCs) (Janeway, Travers, Walport, & Shlomchik, 2005). Macrophages present at the tumour site are activated by virus-infected cells as well as cytokines released by cell lysis (Janeway, Travers, Walport, & Shlomchik (2005)). Macrophages are known primarily to be antigen presenting cells (APCs) and are a critical part of the immune system activation as they serve as stimulates for the naïve immune cells. Dendritic cells (DCs) are one of the most important APCs as they have the ability to activate both helper T cells (Th cells) and killer T cells (CTLs) (Janeway et al., 2005). Th cells are known to activate CTLs and produce the cytokine interleukin 2 (IL-2) which is required for CTLs to proliferate. The primary job of a CTL is to induce apoptosis in cells that have been infected by the virus. CTLs present at the tumour site can be stimulated to induce apoptosis in uninfected tumour cells by the release of immunostimulatory cytokines or tumour antigen in the micro-environment. This is known as the antitumour immune response and is the immune pathway that combined oncolytic virotherapy and immunotherapy looks to exploit.
Cancer cells are known to suppress the immune system by the induction of anergy or tolerance in the host, enabling the cancer cells to proliferate uncontrollably (Janeway et al., 2005). The focus of immunotherapy is to overcome this suppression by stimulating an antitumour immune response. Oncolytic viruses can be used, not only to target and lyse cancer cells, but also as cytokine delivery and generating vectors. IL-12 and GM-CSF are cytokines commonly used as immunotherapeutic agents in cancer gene therapy due to their immunostimulatory ability (Choi et al., 2012a). IL-12 is known to have potent anti-tumour effects through promotion of the immunity of Th1 cells and activation of CTLs (Choi et al., 2012a). IL-12 stimulates Th cells to produce Th1 cytokines: TNF, IFN-γ and IL-2 which in turn stimulate the proliferation of CTLs and Th cells (Janeway et al., 2005). Choi et al. found that intratumoural doses of adenovirus expressing IL-12 strongly induced the activation and recruitment of T cells, including helper T cells and killer T cells (Choi et al., 2012a).
While IL-12 is a known stimulator of T cells the cytokine GM-CSF is known primarily as an attractant for APCs. Choi et al. (2012a) found that GM-CSF expressed in the tumour tissue strongly recruited APCs to the tumour site. GM-CSF is also known to enhance the processing and presentation of antigen on antigen presenting cells (Heystek, Mudde, Ohler, & Kalthoff, 2000). In combination, IL-12 and GM-CSF have been shown to overcome cancer immune suppression and stimulate a highly active antitumour immune response (Choi et al., 2012a).
The key to improving combined virotherapy and immunotherapy lies in quantifying the dependence of treatment outcome on immune stimulation. Here, we develop a system of ordinary differential equations (ODEs) to model the interaction between the immune system, oncolytic adenovirus expressing immunostimulatory cytokines and a population of tumour cells. Parameters are obtained by optimizing the model to the murine tumour time-series measurements under treatment with an oncolytic adenovirus co-expressing IL-12 and GMCSF (Choi et al., 2012a). The model focuses on the specific groups of immune cells that would be most affected by the use of these immunostimulatory cytokines and proposes innovative insights into the key dynamics of this interaction. Using our model, we demonstrate how the up-regulation of the antitumour immune response, stimulated by an oncolytic virus co-expressing IL-12 and GM-CSF, creates more effective cancer treatment. Additionally, a range of immune cell induced apoptosis rates for which the treatment will ultimately be more effective are quantified and the bivariate nature of the treatment outcome is explained.

Experimental data
To determine whether modification of an oncolytic adenovirus with either IL-12 or GM-CSF could improve treatment efficacy, Choi et al. (2012a) investigated the antitumour effect of an oncolytic adenovirus (Ad) co-expressing IL-12 and GM-CSF (Ad/IL12/GMCSF) compared to an oncolytic adenovirus expressing IL-12 (Ad/IL12) or GM-CSF (Ad/ GMCSF). From their experiments Choi et al. determined that oncolytic viruses were most effective when co-expressing the cytokines IL-12 and GM-CSF. In this work, we will help to further their understanding of why this cytokine combination was the most effective, what interactions between the immune system, virus particles and the tumour are most affected by this cytokine combination and how this treatment could be improved.
Oncolytic adenovirus expressing IL-12 and GM-CSF were generated by Choi et al. through inserting specific murine genes into the adenovirus E1 and E3 gene region. Tumours in 6-8 mice were injected with either phosphate-buffered saline (PBS), Ad, Ad/IL12, Ad/GMCSF or Ad/IL12/GMCSF every other day for three total injections, beginning when the average size of the tumour was 80-100 mm 3 . The length and width of each tumour was measured using a calliper and the tumour volume was estimated as .523×length×width 2 . Here we assume the tumour volume is proportional to the number of tumour cells, and the population to be 10 6 cells/mm 3 . The tumour volume was monitored for an experiment-specific number of days.

Model development
To create a deeper understanding of the antitumour effect of an oncolytic adenovirus expressing IL-12 and GM-CSF we develop a mathematical model, drawing on current modelling work in the field. We extend the model developed by Jenner, Yun, Kim, and Coster (in press) that investigated the interaction between an oncolytic virus and tumour cell population in the absence of the immune response. Following this, to model the immune response to oncolytic virotherapy, we draw on the modelling work of Kim, Crivelli, Choi, Yun, and Wares (2015) which looked at the effect of an oncolytic virus expressing 4-1BBL and IL-12 on the immune response. Whilst this model is similar to ours, it does not explicitly account for the helper cell population or limited growth of a tumour (due to space limitations), both factors we incorporate into our model.
We consider a deterministic compartmental model of a tumour cell population interacting with an oncolytic adenovirus co-expressing the immunostimulatory cytokines IL-12 and GM-CSF, Figure 1. The model is used to measure the antitumour effect by anticipating changes in the tumour size with respect to intratumoural administration of oncolytic adenoviruses with differing expression of immunostimulatory cytokines.
We have used mass action in our model as a mean-field approximation of the geometric and spatial effects of the virus-tumour-immune interaction. Frequency-dependent rates have been included to model cell-cell and cell-virus interactions as we assume these occur at a rate proportional to the total number of cells at the tumour site. As is shown below, the model can replicate the observed experimental results under these assumptions.
There are six state variables considered, modelled using the system of ODEs and initial conditions: where t is time, U is the population of uninfected tumour cells, I is the population of infected tumour cells and V is the population of virus particles. As the model was developed for an adenovirus expressing IL-12 and GM-CSF the populations of immune cells considered here are those most affected by these cytokines: antigen presenting cells To model tumour growth we have assumed uninfected cells to be the only population of tumour cells undergoing proliferation. We have modelled tumour cell proliferation using a Gompertz growth function, r log L U U, where r is the growth constant and L is the carrying capacity (Laird, 1964). We have chosen to use Gompertzian growth as opposed to exponential growth as we believe it is crucial to account for the carrying capacity of the tumour due to space and nutrient limitations. While the logistic growth function does account for a tumour carrying capacity, the Gompertz function has the inbuilt characteristic of more rapid increase away from the start point compared to the approach of the carrying capacity. We believe this is what makes the Gompertz function a good approximation for tumour growth, as tumours grow exponentially fast at the start and taper off slowly towards the end.
Virus particles infect uninfected tumour cells at a frequency-dependent rate with constant β. For the purposes of this study, we do not include multiple infections of tumour cells. Virus-infected cells undergo lysis at a rate d I , generating α new virus particles. Once a virus is injected into the body, there are specific immune agents that work to eliminate virus particles, such as those from the complement immune system, encoded in the virus decay term d V appropriately. Following the experiments by Choi et al., an amount V 0 of virus is injected intratumourally on days 0, 2 and 4, this is modelled using a delta function in Equation (7). Antigen presenting cells (APCs) include both dendritic cells and macrophages. These cells are stimulated by infected cells at a rate s A and decay at a rate d A . Helper T cells are then stimulated by APCs at a rate s H and decay at a rate d H . Both APCs and helper T cells then activate CTLs K at a rate s KA and s KH respectively. CTLs induce apoptosis in uninfected and infected tumour cells at a frequency-dependent rate with constant k. CTLs decay at a rate d K . We have assumed that initially there are no stimulated immune cells as these are only generated through the presence of virus infected tumour cells I.

Model fitting
To reduce the degrees of freedom in our model, certain parameters were fixed to those found in the literature. The average time taken for an infected tumour cell to undergo lysis is one day, so we fixed d 1 = 1 (Ganly, Mautner, & Balman, 2000). The number of viral particles created through lysis was estimated as α = 3500 (Chen et al., 2001). We then based our estimation for the rate that the virus leaves the tumour site d V on laboratory observations of Kim, Lee, and Levy (2011), Li et al. (2008) and Wang, Wang, Li, and Yuan (2006), assuming 90% of the virus population decays in one day. Using the half-life decay formula N(t) = N 0 ( 1 2 ) t/t 1/2 we can calculate d V = log (.1) = 2.3/day. Helper T cells are known to have a half-life of 3 days , which gives d H = − log (2)/3 = .23/days. For the immune cell death rates, it was assumed that APCs and helper T cells die or exit the system at a similar rate; therefore, d A = d H = .23/day . For the CTL population, these cells have a shorter half life, estimated at 48 hours, giving d K = .35/day . These parameter estimates are summarized in Table 1 The remaining parameters in the model were obtained by fitting parameters for submodels of Equations (1)-(7) to the data described in Section 2.1, and fixing their values for higher level models in accordance with gradual modifications of the base adenovirus, see Figure 2. The first experiment was the control (PBS). Since there were no viral particles present, the model was reduced to the uninfected tumour population U, by fixing I = A = H = K = 0 in Equations (1)-(7). The remaining parameter values, describing the tumour replication constant r and carrying capacity L, were optimized to the data and fixed for all subsequent optimizations. See Table 2 for the summary of the PBS model specifics.
The first virus-based experiment was the oncolytic adenovirus (Ad) with no immunostimulatory cytokines. We assumed the effects on the populations of immune cells were negligible, considering the absence of IL-12 and GM-CSF, and fixed A = H = K = 0 in Equations (1)-(7). This resulted in the model of Jenner et al. (in press). The remaining parameters of the model describing the infection rate of the virus β and initial tumour size U 0 were optimized and their values were fixed for all subsequent optimizations. See Table  2 for the summary of the Ad experiment model specifics.
The last three viruses tested were modifications of the adenovirus with the different cytokines: Ad/GMCSF, Ad/IL12 and Ad/IL12/GMCSF. Choi et al. found that intratumoural doses of adenovirus expressing IL-12 strongly induced the activation and recruitment of T cells, including helper T cells and killer T cells (Choi et al., 2012a). Hence, to optimize the tumour time-series measurements under treatment with Ad/IL12, we considered the population of APCs negligible and fixed A = 0 in Equations (1)-(7). Similarly for the adenovirus expressing GM-CSF, it was assumed the effect on the helper cell population was negligible as GM-CSF primarily stimulates the antigen presenting cells, (Choi et al., 2012a). Therefore, for this experiment the model was adjusted to exclude the helper T cells, H = 0, and the remaining model was optimized to the data. For the Ad/IL12/GMCSF virus the full model Equations (1)-(7) was used to optimize the model as both cytokines were present.
Due to the overlap in the cytokines expressed by the three viruses, we assumed the stimulation rate of the APCs (s A ) and helper T cells (s H ) could be determined specifically from optimization to the Ad/GMCSF and Ad/IL12 data, respectively. Once the values for s A and s H were obtained, they were fixed in the optimization to the Ad/IL12/GMCSF experiment. The remaining parameters, s KA , s KH and k were then allowed to vary between the three experiments and were used to quantify the major differences in the outcome of treatment from the cytokine expression of the three viruses. A full summary of the experiment-specific hierarchical optimization for the five data sets can be found in Table  2.
The model was optimized to the data using a least-squares non-linear fitting algorithm 'lsqnonlin' in MATLAB2017a. The termination tolerance, which is the minimum change in the objective function, was 1×10 −6 . The maximum number of function evaluations was fixed as 100× number of parameters. The maximum number of iterations for each fit was 400. The solver ode45 was used to solve the model in Equations (1)-(7) for each iteration of the fitting algorithm and to simulate the model. The model was fit to the mean of the data with normalization using the standard error. When solving Equations (1)- (7) numerically, N was replaced by N + for = .001, to avoid the singularity occurring as T → 0. As a second and third injection of treatment was given on days 2 and 4, the model was solved piecewise to account for the addition of V 0 virus particles into the interaction.

Optimized parameter values for the PBS (control) experiment
To assess the antitumour effectiveness of the immunostimulatory adenovirus, Choi et al. (2012a) first conducted a control (PBS) experiment that monitored tumour growth in the absence of treatment, previously discussed in Section 2.1. The model parameters r, the replication rate of tumour cells, and L, the carrying capacity of the tumour, were optimized using the tumour time-series measurements.  (1) Equation (1) Equation (1) Equation (1) Equation (1)  equations Equation (2) Equation (2) Equation (2) Equation (2) Equation (3) Equation (3) Equation (3) Equation (3) Equation (4) Equation (4) Equation (5) Equation (5) Equation (6) Equation (6) Equation (6) Variables Figure 2. Output of the optimized tumour growth model, Table 2, for the PBS (control) case. The individual mouse data is plotted as grey circles, with the mean and standard error bar at each time point in blue. The model output is plotted as a solid black line.
The trajectory of tumour growth arising from the optimized model is close to the tumour growth data from the experiment, Figure 2. The estimates obtained for the parameters are presented in Table 3 with the corresponding goodness of fit estimates in Table 4. These values were then used when optimizing the model parameters using the other, virus-based, experiments.

Optimized parameter values for the adenovirus experiment with no immunostimulatory cytokines (Ad)
To create a baseline for the effectiveness of oncolytic adenoviruses without IL-12 or GM-CSF, Choi et al. (2012a) monitored the growth of pre-established tumours in eight mice after treatment with an adenovirus, previously discussed in Section 2.1. The tumour timeseries measurements exhibited high variability and illustrate the heterogeneity in response to treatment. Multiple mice did not survive the entire experimental duration. Within this data, however, there were three clear subgroups of treatment responses: those that died early, Figure 3(a); low responders, those whose tumours grew slowly until about day 10, after which point the tumours grew exponentially, Figure 3(b); and high responders, those with small tumours over the whole duration of the experiment, Figure 3(c).
To determine whether the model could adequately represent the observed behaviour, we optimized the model and parameter values using each subgroup of data, Figure 3. The optimized values for the infection rate β and initial tumour size U 0 differed for the different subgroups. For the subset that died early we obtained β = 1.3 and U 0 = 220. For the low responder subgroup we obtained β = .92 and U 0 = 27 and for the high responder subset we obtained β = 1.1 and U 0 = 18.
The dynamics of the model optimized to each subgroup was qualitatively similar: each of the solutions rises to a maximum and then decays. Perturbations in β and U 0 alter the location and value of the turning point, not the existence. The large range of initial tumour sizes, U 0 , obtained is an accurate reflection of the initial tumour sizes observed in the experiment. The difference in the infection rates, β, between the three treatment response subgroups was less variable, and the model output was less sensitive to changes in β than U 0 .
Optimizing the parameters to all data simultaneously, Figure 3(d), resulted in β = 1.2 and U 0 = 85, inside the range obtained for the 3 subgroups. Due to the different trajectories, the mean trend of the data and the individual points diverge around day 11, and do not represent any given mouse in the observations. For mice undergoing different treatment protocols we cannot predict whether they would have been high or low responders if treated with adenovirus with no immunostimulatory cytokines. The estimates β = 1.2 and U 0 = 85 obtained using all the data simultaneously lie within the range of the other subgroup estimates, and all the model outputs exhibit the same features. Thus we use these values when optimizing the other model parameters using the data from the more highly modified treatments

Optimized parameter values for the Ad/IL12, Ad/GMCSF and Ad/GMCSF/IL12 experiments
The model parameters were optimized to each immunostimulatory adenovirus based experiment of Choi et al. (2012a) (i.e. Ad/IL12, Ad/GMCSF, Ad/IL12/ GMCSF) as detailed in Table 2. Figure 4 shows the tumour cell population as a function of time for each experiment overlaid with the optimized model. The parameter values obtained are presented in Table  3 and the goodness of fit measures in Table 4.
It can be seen that the model is a good representation of the features of the experimental tumour growth trajectories. As with the Ad experiments, some of the experiments show different response levels to the treatments. In these cases, the model presented reflects the mean behaviour of the data rather than that of any particular subgroup (for instance in Figure 4(c) the mean value straddles two subgroups of responders). Perturbing the parameter values about these mean estimates allows the model to better represent one or other of the subgroups observed.

Parameter values and goodness of fit statistics
In Table 3 we have summarized all the parameter values obtained in the optimization of our model to Choi et al. (2012a). To measure the accuracy of our model in describing the variability of the data, we have calculated the goodness of fit statistics for all five data-sets, see Table 4. A low residual norm and a coefficient of determination (R 2 ) and Pearson's r value close to unity representing a good fit. R 2 is a measurement for how well the fit approximates the data and Pearson's r is a quantity that gives the quality of a least squares fitting to data.

Simulating heterogeneity in immune efficacy
In evaluating treatment efficacy, it is imperative to consider how the outcome of treatment depends on the heterogeneity in immune characteristics. Using model parameters from our optimization to the Ad/IL12/GMCSF data (Table 3), we investigated the effects of perturbations in the rates of immune stimulation and apoptosis induction. We considered whether increasing the immunostimulatory capability of infected cells on APCs, or APCs      Table 3 column Ad/IL12/GMCSF, a detailed view for short times is shown inset.
alter treatment outcome, the killing rate of CTLs, k, was also separately perturbed, Figure  5(c).
In Figure 5(a), we see, as expected, that the higher the stimulation rate of APCs by infected tumour cells, the larger the number of tumour cells. Decreasing the stimulation rate of APCs, results in a much smaller tumour burden, smaller than even the initial tumour size. These findings suggest that increasing APC stimulation has a negative effect on the ability of the treatment to reduce tumour size, and this rate should actually be decreased for an optimal treatment to be obtained. Comparing this to the perturbation in the immunostimulatory rate of helper T cells, Figure 5(b), we see the opposite occurring, with the larger stimulation rates resulting in the smallest tumour size.
As the experiments of Choi et al. (2012a) showed significant tumour growth over the space of 33 days, we take this as the therapeutic window over which the treatment needs to be effective. Tumours were quite large by the end of 33 days for the control case in the experiments of Choi et al. (2012a), and we believe it is important to control the growth within this therapeutic window. We therefore simulate the model in this therapeutic window and look at the effect of perturbations in the immune cell-induced apoptosis rates of tumour cells on treatment outcome. In Figure 5(c), we see that larger k values result in effective early containment of the tumour growth. We also find that for very large values of k, close to k = 1.5, the tumour is completely eradicated (for this model we consider complete tumour eradication to occur if the total tumour population drops below 10 −3 ) in this window of time. However, for mid-range values of the apoptosis rate, e.g. k = 1.25, the treatment results in a large growth increase of the tumour around day 25. These two treatment responses (complete eradication or unbounded growth) mimic the results seen in Figure 4(c). Interestingly, when k is much smaller, e.g. k = 1, we see a lower maximum tumour count achieved within this time frame. These findings suggest the existence of a mid-range interval of k values for which the treatment is significantly less effective in the time frame of 33 days (the therapeutic window discussed earlier), than may have been anticipated outside of this interval. We also see that for large values of k complete tumour eradication can be obtained. This indicates that tumour cell apoptosis is a critical feature in the efficacy of treatment.

Discussion
The mathematical model presented in this article was used to identify the primary processes in the interaction between a population of tumour cells and an oncolytic adenovirus coexpressing IL-12 and GM-CSF. Using our model, we successfully replicated and embodied the experimental results from Choi et al. (2012a), see Figure 4. It is evident through visual inspection of Figure 4 that our model, along with the hierarchical fitting algorithm presented in Table 2, provides a reliable representation of the data. Goodness of fit measurements in Table 4 confirm that the models closely approximated the true system, with R 2 values greater than .98 and Pearson's correlation coefficient greater than .87.
Heterogeneity within individual mice tumour responses under treatment with an oncolytic adenovirus is clearly seen in the tumour time-series measurements of Choi et al. (2012a), Figure 3. There are three noticeable subgroups of treatment responses: those that died early; low responders, those whose tumours grew slowly until about day 10, after which point the tumours grew exponentially; and high responders, those with small tumours over the whole duration of the experiment. To examine the differences behind the heterogeneity between the subgroups, and to see whether the model was sufficiently flexible to embody all the observed behaviour, we optimized the model to each subgroup. All were well explained by the model. The subgroups had slightly different rates of infectivity of the treatment, but more importantly started with different initial tumour sizes. Examining the models for the subgroups and that for all the data, Figure 3, we can see that the long-term dynamics of these underlying subgroups are qualitatively similar.
The hierarchical modelling fixed the rate of infectivity and initial tumour size to the values determined by all of the data for Ad mouse treatment response. The other model parameters were then optimized using the data for the other treatment responses. When the rate of infectivity and initial tumour size were varied to reflect the different Ad subgroups and the model then re-optimized to the subsequent treatment data very similar fits resulted (not shown). Whilst the optimized model does not represent the treatment response of every mouse, perturbation of the model parameters allows for each variation to be observed, so the structure of the model of the model is flexible enough to encompass all the observed responses.
Quantifying the effects of cytokine combinations on treatment efficacy is possible through optimizing parameters in our model to the experiments of Choi et al. (2012a). The nature of the hierarchical experiments allows us concentrate on the primary differences between the immunostimulatory oncolytic adenoviruses: Ad/GMCSF Ad/IL12 and Ad/GMCSF/IL12. Comparing CTL induced apoptosis rate, k, for the Ad/IL12 and Ad/GMCSF viruses, Table 3, we clearly see that expression of cytokine IL-12 results in a higher immune cell killing rate. Therefore, the addition of IL-12 has the greater effect on improving immune cell killing rate k, and consequently tumour cell death. This is also evident when comparing the tumour time-series measurements obtained in Figure 4(a) and (d) where it is clear that Ad/IL12 has a greater antitumour potency then Ad/GMCSF. We found that the largest immune cell killing rate was obtained for co-expression of both cytokines -Ad/IL12/GMCSF. This suggests that it is only with both cytokines that the treatment reaches its maximal effectiveness in stimulating the immune system to attack the tumour cells.
Furthermore, the optimization results suggest a competition in IL-12 and GM-CSF's effectiveness on immune cell stimulation. Examining Figure 4(c) we see a much more advantageous result achieved under treatment by both cytokines, compared to all other oncolytic adenoviruses tested. However, when we examine the parameter values for the Ad/IL12, Ad/GMCSF and Ad/IL12/GMCSF experiments, we see a decrease in the rate of Th activation, s H , and an increase in the rate of APC activation, s A , when both cytokines are being expressed. Mechanistically, the model suggests that combining both cytokines reduces the number of T cells produced and increases the presence of APCs at the tumour site, suggesting the existence of negative feedback. Currently, however, there is no experimental evidence to support this hypothesis so we pose it as a motivation for future research. Future investigations in the field of immunotherapy could look into the effect of combining cytokines IL-12 and GM-CSF on immune stimulation to determine whether it does result in an increased immune cell killing response and increased stimulation of APCs and a decrease in Th cell activation as our model proposes, and whether this dynamic shift in immune behaviour is what results in a more effective treatment.
While the results of the Ad/IL12/GMCSF experiments reduce the tumour population most significantly out of the five experiments, the finding that Th cell activation is decreased requires further investigation. How exactly this might be hindering the immune interaction and the obtaining of optimal treatment efficacy is the subject of future work. In Figure 4, it can be seen that only one mouse in the Ad/IL12/GMCSF experiment had tumour growth after day 20. From our optimization results, we propose that tumour cells in this case may have escaped immune removal by down-regulation of Th cell activation.
To investigate further the immune heterogeneity, individual responses to changes in immune efficacy were simulated. The analysis in Section 3.5 suggests there is a counter intuitive relationship between treatment efficacy and immune stimulation rates. Using the model optimized for Ad/IL12/GMCSF as a platform for prediction, we see that the dependence of treatment efficacy on APC stimulation, s A , and Th cell stimulation, s H , differs significantly, Figure 5(a) and (b). Simulations show that increasing the stimulation rate of APCs has a negative effect on the treatment efficacy, allowing for the tumour cell population to escape the control of treatment and grow unbounded. However, increasing Th stimulation rates has a positive effect on treatment efficacy, allowing for the tumour cell population to be controlled for longer and, for certain parameters, to be completely eliminated. These results suggest that there is a sensitive threshold of APC stimulation, above which a negative effect on the immune response occurs. Biologically, this could signify an over-stimulation of immune cells that results in an incapacity of the treatment to contain and eliminate the tumour cell population. On the other hand, increased stimulation of Th cells consistently promoted tumour cell death. The results we present here are purely hypothetical and suggest that further investigations of this cancer treatment could examine how increasing the expression of IL-12 cytokine and decreasing GM-CSF expression has a downstream effect of the probable increase in Th cell stimulation and decrease in APC stimulation.
Our model suggests that reducing APC stimulation and increasing helper T cell stimulation could possibly improve treatment. Researchers have suggested the possibility that chemical inhibition of the MAPK ERK pathway in DCs reduces the maturation of these APCs, and therefore, the stimulation rate, (Liechtenstein, Dufait, Lanna, Breckpot, & Escors, 2012;Puig-Kröger et al., 2001). This is one possible avenue of investigation that could be undertaken to test the results in Figure 5(a). To increase the helper T cell activation, both cytokines IL-1 and IL-12 are known to heavily stimulate the differentiation of naive T cells (Liechtenstein et al., 2012;Macatonia et al., 1995). So to test the results seen in Figure  5(b), we could consider an additional intravenous injection of IL-12 or IL-1 as a possible way of increasing the stimulation rate of helper T cells.
Heterogeneity in immune-cell-induced apoptosis is a key determinant of treatment outcome. Perturbations in the rate of CTL induced apoptosis, k, for the Ad/IL12/GMCSF model, see Figure 5(c), demonstrate a very interesting phenomenon: the existence of a parameter window for which the treatment is relatively ineffective compared to parameter values outside this interval. We discovered an extremely sensitive non-linear relationship between treatment outcome and CTL induced apopotosis rate k. In Figure 5(c) we see that for lower values of k the tumour population is initially controlled with slow growth over time, reaching a turning point, after which tumour volume decreases. It may be that the immune system is able to control the tumour growth even with this smaller killing rate. What is interesting is if we increase k we notice that the tumour volume at the turning point increases and, for a range of k values, the tumour population is able to grow unbounded. Increasing k further, we find that the tumour population can be completely eradicated (in this model we consider complete eradication to be obtained when the total tumour population drops below 10 −3 ). From this result we suggest that the different responses in the mice in Figure 4(c) (i.e. tumour eradication or unbounded growth) could possibly be explained by a difference in the immune cell killing rate of tumour cells. These results also suggest that there may be a window of CTL induced apoptosis rates, for which the treatment is ineffective, but outside of which we can see either controlled tumour growth within the time period of 33 days or complete eradication.
Controlling the CTL induced apoptosis rates k could be achieved through the introduction of an experimental cancer treatment known as CTLA-4 blockades (Henson, Macaulay, Kiani-Alikhan, & Akbar, 2008;Parry et al., 2005). CTLA-4 is a well-known T cell inhibitory B7-receptor that is expressed by activated T cells. Using the CTLA-4 blockade, researchers have shown that this treatment can enhance T cell cytotoxic responses and induce the differentiation of cytotoxic CD4 T cells Leach et al. (1996). Therefore, one possible way to investigate whether perturbations in the CTL induced apoptosis rates k would result in the outcomes presented in Figure 3(e), would be to look at combining CTLA-4 blockades with oncolytic virus expressing IL-12 and GM-CSF.
In Figure 5(c), it is also evident that we have an exchange in the dominant processes acting as a function of the cell induced apoptosis rate, k . For high values of k, it is clear that tumour cells are predominantly removed by the immune system, which is why the tumour is eventually completely eradicated. However, reducing the value of k results in the initial decrease in tumour cell numbers due to viral interactions rather than the immune system. This result reinforces the importance of stimulating the correct mechanisms at the right stage of tumour growth when investigating improvements for combined oncolytic virotherapy and immunotherapy.

Conclusion
We as humans are incredibly heterogeneous and it is important as modellers to analyse and investigate this heterogeneity. Our model shows that further investigation into the immune stimulation rate and response rate is incredibly important if we wish to further the efficacy of this treatment. We propose the existence of certain CTL-induced apoptosis rates for which this treatment is ineffective and that stimulation rates of APCs and Th cells require further experimental optimization if we wish to improve the outcome of treatment.
Immunotherapy has become a promising new frontier in cancer treatment. Combined with oncolytic virotherapy, these two fields could hold the key to a curative treatment; however, there is still a long way to go before their complex nature is completely understood. The mathematical model presented in this article presents a simple yet elegant dissection of the virus-tumour-immune interaction for the interaction of an oncolytic adenovirus co-expressing IL-12 and GM-CSF and a population of tumour cells. While specifically derived for the experiments of Choi et al., we suggest that our model could be used for future combined oncolytic virotherapy and immunotherapy experiments. We have illustrated the effectiveness of a hierarchical optimization algorithm in obtaining parameter estimates and replicating and embodying experimental data in the model.
We believe that in devising future strategies for tumour eradication, we must first investigate the rates of immune stimulation and CTL induced apoptosis. Only for certain ranges of these parameters do we see complete eradication of the tumour cells. We believe this work opens further questions about these interactions that will lead us to develop more effective cancer treatments in this field.