Optimising Hydrogel Release Proﬁles for Viro-Immunotherapy Using Oncolytic Adenovirus Expressing IL-12 and GM-CSF with Immature Dendritic Cells

: Sustained-release delivery systems, such as hydrogels, signiﬁcantly improve cancer therapies by extending the treatment efﬁcacy and avoiding excess wash-out. Combined virotherapy and immunotherapy (viro-immunotherapy) is naturally improved by these sustained-release systems, as it relies on the continual stimulation of the antitumour immune response. In this article, we consider a previously developed viro-immunotherapy treatment where oncolytic viruses that are genetically engineered to infect and lyse cancer cells are loaded onto hydrogels with immature dendritic cells (DCs). The time-dependent release of virus and immune cells results in a prolonged cancer cell killing from both the virus and activated immune cells. Although effective, a major challenge is optimising the release proﬁle of the virus and immature DCs from the gel so as to obtain a minimum tumour size. Using a system of ordinary differential equations calibrated to experimental results, we undertake a novel numerical investigation of different gel-release proﬁles to determine the optimal release proﬁle for this viro-immunotherapy. Using a data-calibrated mathematical model, we show that if the virus is released rapidly within the ﬁrst few days and the DCs are released for two weeks, the tumour burden can be signiﬁcantly decreased. We then ﬁnd the true optimal gel-release kinetics using a genetic algorithm and suggest that complex proﬁles present unnecessary risk and that a simple linear-release model is optimal. In this work, insight is provided into a fundamental problem in the growing ﬁeld of sustained-delivery systems using mathematical modelling and analysis.


Introduction
Controlled localised delivery systems have been replacing systemic administration of cancer therapies for some time [1]. These sustained-release injectable formulations are usually designed as microparticles, implants or gels [1]. They extend the presence of therapy at the tumour site by gradually releasing treatment as they degrade, improving treatment efficacy and avoiding excessive treatment loss to neighbouring tissue or the lymphatic system. Sustained-delivery systems have improved the effectiveness of a range of cancer treatments, including chemotherapy [2,3], oncolytic virotherapy [4][5][6][7] and immunotherapy [8][9][10].
Oncolytic virotherapy and immunotherapy are two promising cancer therapies that have gained significant attention over recent years [11]. Oncolytic virotherapy focuses on genetically engineering viruses to infect, replicate within, and lyse tumour cells. Many oncolytic viruses can accommodate gene insertions and thus be "armed" with therapeutic transgenes to produce immunostimulatory signals that stimulate an antitumour immune response, making them useful immunotherapy agents [12,13]. Immunotherapy treatments typically involve administration of additional cytokines, molecules or immature immune cells.
Interleukin 12 (IL-12) and granulocyte-monocyte colony-stimulating factor (GM-CSF) are cytokines used in immunotherapy clinical trials [14,15]. IL-12 is known to have potent antitumour effects through promotion of the immunity of helper T cells and activation of killer T cells, whereas GM-CSF is known to enhance the processing and presentation of antigen on antigen-presenting cells (APCs) [13]. Immature immune cells, such as dendritic cells (DCs), are also tested in immunotherapy trials with the goal of overcoming the induction of anergy or tolerance by the cancer [16,17]. Intratumoural administration of DCs increases the probability of tumour antigen recognition and subsequent activation of helper T cells and killer T cells, leading to heightened tumour cell killing [18].
An increasing number of studies has shown that combined virotherapy and immunotherapy (viro-immunotherapy) leads to a prolonged antitumour response, where both viruses and immune cells work simultaneously to eradicate cancer cells [10,13,[19][20][21]. Conventionally, oncolytic virotherapy and immunotherapy are administered intratumourally or intravenously. This can result in rapid clearance of the virus and functional impairment of immune cells in the tumour microenvironment [10]. To overcome this, there have been many investigations into the potential of treatment loaded onto injectable gels or polymer matrix systems [5][6][7]10]. Oh et al. [10] showed that an oncolytic adenovirus expressing IL-12 and GM-CSF loaded onto a hydrogel with immature DCs had a higher retention in tumour tissue. Additionally, the use of the gel resulted in greater antitumour activity due to the sustained release of therapy. Immune checkpoint blockade therapy and agonists were also shown to be considerably more effective when loaded onto an degradable gel matrices [8,9]. A major challenge for sustained-release mechanisms is determining a release profile that achieves optimal therapeutic benefit.
Mathematical models describing the biological response to disease can be used as an important tool in therapy planning and optimisation [22][23][24][25]. Deterministic systems of ordinary differential equations (ODEs) are used to describe the key biophysical interactions and are then analysed to determine optimal therapeutic protocols [23]. In oncology, these techniques have been used to successfully optimise chemotherapy and immunotherapy treatment [25][26][27][28][29][30]. Zhu et al. [31] optimised the scheduling of one cycle of chemotherapy drug VP-16 and found a new optimal drug regime which improved the clinical efficacy. For combined immunotherapy and chemotherapy, de Pillis et al. [32] developed an ODE system that predicted which treatment protocols would result in tumour elimination. There are also a handful of examples where deterministic models have also been used to improve chemotherapy delivery from sustained-release systems [33,34], Mathematical modelling has also been used to determine optimal dosage protocols for oncolytic virotherapy and immunotherapy [35][36][37][38]. Using ODEs calibrated with experimental data, Barish et al. [38] demonstrated that lower oncolytic virus doses combined with higher dendritic cell doses resulted in an optimal robust therapy. Optimising the scheduling of systemic doses of DCs and oncolytic virus, Wares et al. [37] found that tumour eradication could be achieved only when oncolytic viruses were administered before DCs. In work by Cassidy and Craig [29], an optimised schedule for combined virotherapy and chemotherapy treatment was obtained using a genetic algorithm. So far, however, no one has optimised the release profile of viro-immunotherapy from a sustained-delivery system. The challenges faced in developing such a system lie in de-complexifying the immune response to viro-immunotherapy, while still accurately representing a relevant experimental treatment.
In this article, we develop a novel system of ODEs to optimise the effectiveness of an oncolytic adenovirus genetically modified to express IL-12 and GM-CSF loaded onto a hydrogel with a population of immature DCs [10]. The parameter values in the model are optimised to in vitro and in vivo measurements relating to the release of this therapy from a gelatin-based hydrogel. With this calibrated model, we conduct a numerical investigation of the hydrogel release profile and investigate optimised sustained-delivery protocols. Implementing a genetic algorithm, we then determine the true optimal release curves and discuss how complex release curves may impose a less effective treatment.

Mathematical Model
Human immune systems have a very powerful and effective way of eliminating viruses. Macrophages and dendritic cells (DCs) are activated by virus-infected cells through the presentation of virus antigen [18,39]. Additionally, when virus-infected cells undergo lysis, cytokines and antigens are released, and these can activate DCs and macrophages [18]. These activated (mature) cells then stimulate helper T cells and killer T cells. Helper T cells also secrete specific cytokines (specifically type Th1) that provide continual stimulation of killer T cells.
Previously, we developed a system of ODEs that assumed the primary driver of the immune response was virus-infected tumour cells [36]. The virus-infected tumour cells stimulated APCs which in turn activated killer T cells and helper T cells. This system was used to model an oncolytic adenovirus expressing IL-12 and GM-CSF (Ad/IL12/GMCSF) created by Choi et al. [13]. Intratumoural doses of Ad/IL12 were shown to induce the activation and recruitment of helper T cells and killer T cells, whereas administration of Ad/GMCSF strongly recruited APCs to the tumour site.
In this work, we consider the injection of Ad/IL12/GMCSF, combined with an injection of immature DCs. Due to the presence of additional immature DCs, we believe the stimulation of the immune system by uninfected tumour cells would be crucial, and we extend the system of ODEs in Jenner et al. [36] to consider this, dD L dt " f u DC ptq´s AU D L U´s AI D L I´d L D L , dA I dt " r AI I´s AU A I U´s AI A I I´d AI A I , where t is time (days), U is the uninfected tumour population, I is the infected tumour population and V is the number 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: activated (or mature) antigen-presenting cells (APCs), A M ; helper T cells, H; and killer T cells, K. Additionally, as immature DCs are injected into the system, we now consider populations of short-lived and long-lived immature DCs that have been injected, D S and D L , and immature APCs that are already present in the body, A I . The total cell population at the tumour site at any time t is given by T " U`I`D L`DS`AI`AM`H`K . Figure 1 summarises the interactions modelled in Equations (1)- (9). The biological mechanisms described by each ODE in the system are described below.

•
In Equation (1), uninfected tumour cells, U, are growing at a rate described by a Gompertz function with proliferation rate r and carrying capacity L. Uninfected cells are infected by virus particles, V, at a frequency-dependent rate with rate constant β. Killer T cells, K, induce apoptosis in uninfected cells at a frequency-dependent rate with rate constant κ.

•
In Equation (2), tumour cells become infected cells, I, at rate βUV{N. Infected cells lyse at rate d I I, and, like uninfected cells, are killed at a frequency-dependent rate.

•
In Equation (3), virus particles enter the system at rate u V ptq, either by a single injection or release from the hydrogel. Virus particles decay at a rate d V , and α new viruses are created through lysis.

•
In Equation (4), short-term immature DCs, D S , enter the system at a rate p1´f qu DC ptq, either by a single injection or release from a hydrogel. They are then stimulated to become mature or activated APCs, A M , at rate s AU or s AI due to the interaction with either uninfected or infected tumour cells. They decay at a fast rate d S , where d S ąą d L .

•
In Equation (5), long-term immature DCs, D L , enter the system at rate f u DC ptq, where u DC ptq is the rate of injection or release from the hydrogel and f is the fraction of the initially loaded DCs that are long-term DCs, 0 ă f ă 1. Similarly to short-term DCs, they become activated APCs at rate s AU and s AI . These DCs decay slowly at rate d L .

•
In Equation (6), immature APCs, A I , are recruited to the tumour site by infected cells at rate r AI I. Immature APCs are then stimulated at the same rate as immature short-term and long-term DCs. The immature APCs decay at rate d AI .

•
In Equation (7), mature APCs are generated through the maturation of introduced long-term and short-term DCs and immature DCs recruited by uninfected tumour cells at rate s AU and infected tumour cells at rate s AI . These cells decay at rate d A .

•
In Equation (8), helper T cels, H, are stimulated by mature APCs at a rate s H . These cells decay at a rate d H . We assume the rate APCs stimulate helper cells and helper cells stimulate killer cells is independent of antigen.

•
In Equation (9), killer T cells, K, are activated by helper T cells and mature APCs at rates of s KH and s KA , respectively. These cells either die or leave the tumour site at a rate d K .
As the virus was administered either intratumourally or from the gel, there is no need to model the virus in the organs and blood (discussed by Jenner et al. [40]). Additionally, the influence of interferon-mediated antiviral immunity is not considered crucial. Note that this is the same killer T cell population that was considered by Jenner et al. [36], but this time the activation mechanism is modelled in more detail. The above model is similar to the model used by Kim et al. [35].
We consider two populations of immature dendritic cells (D L and D S ) in the initial administration as the in vitro release curve depicted a biphasic exponential decay profile for the DCs, see Section 3.1. We separated this population from the immature DCs already present in the mice, as we did not want to impose the assumption that the same separation in the survival rates would exist in this population.
If a gel is used to release the DCs and/or virus the initial conditions for system are where U 0 is the initial tumour size, whereas, for a single injection of virus or DCs, the initial conditions change to Vp0q " V 0 , D L p0q " f D 0 and D S p0q " p1´f qD 0 for the virus and DCs, where V 0 and D 0 are the initial number of viruses or DCs, respectively.

Calibrating Parameters to In Vitro and In Vivo Time-Series Measurements
To calibrate the parameters in Equations (1)-(9), we used the in vitro DC viability counts and in vivo tumour volume measurements of Oh et al. [10]. The release rate of DCs from the gel, u DC ptq, and the decay rate for the long-term and short-term DCs, d L and d S , were determined from DC viability counts with and without a degradable hydrogel [10]. We then fixed certain parameters to their values obtained in our previous work [36] where we optimised tumour volume measurements for the Ad/IL12/GMCSF virus based on experiments by Choi et al. [13]. The remaining parameters were then optimised to the in vivo tumour time-series measurements from [10] under treatment with a single injection of either DC, Ad/IL12/GMCSF or DC+Ad/IL12/GMCSF and treatment with a degradable hydrogel containing DC+Ad/IL12/GMCSF.

DC Release-Profile and Decay Rates
Oh et al. [10] counted the number of viable DCs after a single injection released from a gel over 6 days. To fit these measurements, we fixed U " I " V " A I " A M " H " K in Equations (1)-(9) as no tumour cells, immune cells or virus were present, and then optimised the decay rates of the long-term and short-term DCs, d L and d S , to the DC count after a single injection, see Figure 2a and Table 1. We initially considered a single decay rate for the injected immature DCs and found that the long-term model prediction was a significant underestimation. The data seems to suggest a fraction of the injected DCs live longer and contribute to a slower population decline after the first few days. As the DCs are derived from bone marrow, the difference in decay rates may be a result of heterogeneity imposed from the progenitor cell lineages that generate the DCs [41,42]. We, therefore, assumed that there are both long-term and short-term surviving immature DCs in the injected population. This biphasic decay rate is evident in other studies of bone marrow-or monocyte-derived DCs [43]. The decay rate could be modelled more accurately with a age-structured partial differential equation system, but due to the limited data, we decided this was not necessary.
From the measurements of the DCs released from the hydrogel [10], we then obtained the release-rate function u DC ptq. Oh et al. [10] used a gelatin-hydroxyphenyl propionic acid (GHPA)-based hydrogel that can easily be manipulated to achieve appropriate hydrogel properties such as gelatin time, mechanical stiffness and degradation rate. Previous experiments with this hydrogel showed that gel weight decreases linearly over time [44], implying that the gel degrades linearly. Taking the forward finite difference approximation of the data revealed an approximately linear, increasing rate of change for the DCs released from the gel. As we chose not to model the gel-release mechanics, we assumed that a linear function would be sufficient to capture the gel-release rate.
We fixed d L , d S and f to the values obtained for Figure 2a and assumed u DC ptq would be where a D is the gradient of u DC ptq; b is the initial release rate from the gel (i.e., u DC p0q " b D ); and D I is the number of DCs in the gel, where dD I {dt "´u DC ptq and D I p0q " 2.5ˆ10 6 . Qualitatively, a D describes the increasing speed that treatment is released from the gel as a function of time, and b D describes the initial rate at which treatment is released when the gel is first put into the solution bath. See Figure 2b and Table 1 for the results. From the fitted parameter values, we obtain the approximate time that the gel degrades as t r " 5.945 days, which matched the observations of Oh et al. that by day 6 the gel had completely degraded.  Figure 2. Fitting the viable dendritic cell (DC) count with and without gel to parameters in Equations (1)-(9) with U " I " V " A I " A M " 0, where D S is the short-term DC population and D L is the long-term DC population. In panel (a), the decay rate of short-term and long-term DCs, d S and d L , and the fraction of the initial injection that consisted of short-term DCs, f , was fit to the viable DC count after a single injection, where u DC ptq " 0. In panel (b), the count of DCs released from the gel was fit to the release function u DC ptq fixing the value of d L , d S and f obtained in panel (a). In both figures, circles represent the number of released viable DCs, the black solid lines are the model's approximations to the data, the blue dashed lines are the long-term DCs in the dish, the yellow dotted lines are the short-term DCs in the dish and the purple dash-dot line is the number of DCs in the gel. The fitted parameter values are in Table 1. Table 1. Parameter values obtained from fitting the release rate of the DCs and the short-term and long-term decay rates of the DCs to the in vitro viable DC count, see Figure 2.

Tumour, Immune and Viral Parameters
To obtain remaining parameter estimates for the model, we sequentially optimised the tumour time-series measurements of Oh et al. [10] for PBS, gel, Ad/IL12/GMCSF, DC+Ad/IL12/GMCSF and DC+Ad/IL12/GMCSF+gel. To reduce the degrees of freedom in the model, the parameters  [36]. The remaining parameters were obtained by sequentially fitting sub-models of Equations (1)-(9), where fitted parameter values were fixed for higher level models in accordance with the different treatments investigated. Table 2 gives a summary of the fitting algorithm, with parameter values in Table 3. Table 2. Experiment-specific optimisation conditions for the measurements of Oh et al. [10]. Equations used to optimise each experiment are listed along with the state variables considered and parameters fitted or fixed. Some parameters were fixed to previous optimisation work in Jenner et al. [36].

Relevant equations
Equation (1) Equation (1) Equation (1) Equation (1) Equation (1) Equation (2) Equation (2) Equation (2) Equation (3) Equation (3) Equation (3) Equation (5) Equation (5) Equation (5) Equation (4) Equation (4) Equation (4) Equation (6) Equation (6) Equation (7) Equation (7) Equation (7) Equation (7) Equation (8) Equation (8) Equation (8) Equation (8) Equation (9) Equation (9) Equation (9) Equation (9) Variables to Tables 1 and 3 Table 3. Parameter estimates from the sequential optimisation of the model following the algorithm in Table 2 to the experimental measurements in [10]. Certain parameters were fixed to their value in [36], and this is indicated in the  As a control, Oh et al. [10] measured the tumour growth under a single injection of PBS and an empty gel. We assumed that the underlying tumour growth in both of these experiments should be identical and fit the growth rate, r; carrying capacity, L; and initial tumour size, U 0 , to both data sets, fixing I " V " A I " A M " H " K " 0 in Equations (1)  Following the control experiments, Oh et al. [10] measured the tumour size under a single injection of Ad/IL12/GMCSF. As there was no injection of immature DCs in this experiment, we assumed the immune stimulation would be driven solely by the injected virus with the effects of endogenous DCs negligible, i.e., we fixed D S " D L " A I " 0. Under this assumption, we removed any activation of APCs by uninfected tumour cells, and modelled APC stimulation as directly relating to infected cells I, i.e., s AI IpA I`DL`DS q " s AI I. The remaining parameters of the model describing the infection rate of the virus, β; the killing rate of the killer T cells, κ; and the initial tumour size, U 0 , were optimised, see Figure 3b.
Oh et al. [10] then measured tumour growth under a single injection of 2.5ˆ10 6 immature DCs. As there was no virus present, we fixed I " V " 0. The decay rate for immature APCs, d AI , was fixed to the short-term DC decay rate, d S , as we assumed that any recruited immature APCs would not reside at the tumour site for long. The parameters describing the stimulation rate of immature APCs, s AU ; the killing rate of killing cells, κ; and the initial tumour size were optimised and fixed for subsequent optimisations, see Figure 3c.
Following this, Oh et al. [10] combined a single injection of Ad/IL12/GMCSF with a single injection of immature DCs. Fixing the previous fit parameters, we optimised the kinetics for the activation of immature APCs and DCs using the full model in Equations (1)- (9). Optimising the rate that infected cells recruit inactivated APCs r AI , the stimulation of immature APCs s AI and initial tumour size U 0 gave the model curve in Figure 3d.
The final experiment of Oh et al. [10] measured the tumour volume after a gel loaded with Ad/IL12/GMCSF and immature DCs was injected adjacent to the tumour. Using these measurements, we obtained the release profile for the virus from the gel u V ptq. Fixing the release profile of the immature DCs from the gel, u DC ptq, to the profile determined in Section 3.1 and all parameters to the values for the DC+Ad/IL12/GMCSF experiment, we assumed the release rate for the virus would be linear, and the time of release t r would be equivalent to that obtained for the DCs, t r " 5.945 days, where a V ą 0 and b V ě 0. The value of b V was obtained by fixing the initial amount of virus as V 0 " 2ˆ10 10 and integrating u V ptq. Additionally, we obtained an upper and lower bounds of a V using the condition that u V ptq is an increasing function. Optimising the release function parameter a V was not sufficient, and a fit could only be obtained when the activation rate of APCs by uninfected cells was also allowed to vary. The resulting fit is in Figure 3e and the parameter values are in Table 3. The model dynamics for this experiment are plotted adjacent in Figure 3f. All parameter optimisations were undertaken using a least-squares nonlinear fitting algorithm "lsqnonlin" in MATLAB2017a. The termination tolerance was 1ˆ10´6. The maximum number of function evaluations was fixed as 100ˆnumber of parameters. The maximum number of interactions for each fit was 400. The solver ode45 was used to solve the model. The model was fit simultaneously to the set of individual mouse data for a particular experiment. When solving Equations (1)-(9) numerically, T was replaced by T` for " 0.001 to avoid a singularity as T Ñ 0. The goodness of fit statistics have been reported in Table 4.
It is evident visually in Figure 3 and from the goodness of fit measurements ( Table 4) that our model and the sequential fitting algorithm we used in Table 2 provide a reliable representation of the data. From the population dynamics in Figure 3f, it is clear that the immune response is driven by a large increase in helper T cells coinciding with mature APCs. This then prolongs the killer immune cell population's survival. From the model, it does appear though that it is the initial viral infection that drives the tumour population down significantly. Table 4. Goodness of fit statistics for the optimisation of the mathematical model Equations (1)- (9) with the experiments in Figure 3 following the optimisation algorithm in Table 2.

Optimising the Gel Release Profile
Oh et al. [10]'s gel-based medium improved the therapeutic efficacy of the DC+Ad/IL12/GMCSF treatment and reduced the tumour volume to 1005.9 mm 3 by day 20, see Figure 3e. It may be possible to further improve this treatment through the manipulation of the gel-release profile, u DC ptq and u V ptq. To investigate this, we considered gel-release profiles that were constant, linear and sigmoidal. We chose these functions as they had increasing complexity in their time dependence and were qualitatively similar for certain parameter regimes. Many therapeutics require a constant release rate for varying durations [45], and while this is difficult to achieve, there are several examples [46][47][48]. In comparison, variable drug release rates, such as linear and sigmoidal curves are more obtainable [49][50][51][52].
We, therefore, define three generalised release rate functions for u DC ptq and u V ptq: constant f c ptq, linear f l ptq and sigmoidal f s ptq. These have the forms where C is the constant release rate, a is the gradient of the linear release, k is the steepness of the sigmoidal release rate, x 0 is the midpoint of the sigmoidal function and L is the maximum sigmoidal release rate. We impose the constraint that the total virus and DCs released from the gel over t r days must be equal to V 0 and D 0 , i.e., Fixing all parameter values not related to the release curves to those in the DC+Ad/I/G+gel column in Table 3, we simulated the tumour size up to day 20 using Equations (1)-(9) for the different release profiles f c ptq, f l ptq or f s ptq for t ď t r , and 0 otherwise.

Constant Release Profile
Assuming both virus and DCs are released from the gel at a constant rate, f c ptq in Equation (12), and integrating and simplifying using Equation (15) gives We can then simulate the effect of different release periods, t r , for virus and DCs. The tumour volume on day 20 as a function of the release period is plotted in Figure 4a. The surface has a global minimum tumour size of 689.24 mm 3 , which corresponds to a 13-day release period of DCs and a single day for the virus.  To illustrate the tumour growth dynamics under different constant release profiles, model simulations are plotted in Figure 4b,c, corresponding to the red points on the surface plot in Figure 4a. There is a decrease in the initial gradient of the tumour growth in Figure 4b compared to Figure 4c, which in turn allows for a more gradual linear tumour growth over the period of 20 days. This suggests that releasing the virus rapidly and then slowly releasing the immature DCs over 13 days, reduces the growth rate of the tumour by providing continual stimulation for the immune system.

Linear Release Profile
Assuming both virus and DCs are released from the gel at a linear rate, f l ptq in Equation (13), which is increasing, applying the constraint in Equation (15) gives To investigate whether there was an optimal increasing linear release rate, we initially fixed t r as common for virus and DCs and varied the gradient of the release curve a, which was bounded by 0 ă a ă 2A 0 {t 2 r so as to be strictly positive and increasing. The resulting tumour size on day 20 is plotted in Figure 5a,b. As t r increases from 1 day to 5 days, the overall tumour size decreases irrespective of the choice of a. However, after a release period of t r " 6 days, the minimum achievable tumour size at day 20 increases and so does the dependence of the tumour size at day 20 on the steepness of the release rate a. A minimum tumour size of 887 mm 3 was achieved for t r " 11 days and t r " 13 days when a " 0. In other words, if the release time is common for DCs and virus, the optimal linear release curve is approximately constant. We then numerically simulated the increasing linear release rate with t r varying for virus and DCs and fixed a " 0.001, 3.3 and 7.4, see Figure 5c. Qualitatively, the dependence on the value of t r for u DC ptq and u V ptq is similar to that of the constant release profile in Figure 4a. A minimum tumour volume of 628.4 mm 3 was achieved for an interval of a between 3.3 and 4.6 for t r " 1 for the virus and t r between 11 and 13 for DCs. This suggests that it is possible to obtain an optimal efficacy with the linear release function, but only when the release time of the constituents are separated.
Although the gel-release mechanisms measured by Oh et al. [10] had an increasing gradient, it is worth considering how effective a treatment would be when released at a decreasing linear rate. A decreasing linear function crosses u DC ptq " 0 or u V ptq " 0 when t "´b{a. As this can occur before or after 20 days, we use this to constrain the release curve by fixing the total drug released between 0 and´b{a days equal to V 0 and D 0 . Integrating and simplifying gives where b is the initial release rate. Note that this means the gel releases until´b{a days, or 2V 0 {b for the virus and 2D 0 {b days for DCs. This can be outside the window of the evaluated time-frame and as such may lead to less virus and DCs released than V 0 and D 0 . Simulating a common value of b for virus and DCs gives the tumour volume on day 20 in Figure 5d. The tumour volume depends on the value of the initial release rate for the virus with larger values of b resulting in the lowest tumour volume. Note that when b increases, the release time t r decreases, so the optimal tumour burden for the decreasing linear release is achieved when b is large enough so that the total amount of D 0 and V 0 is released within 20 days. Interestingly, we do not see a dependence on the time frame in which the treatment is released (in contrast to the cases of constant release and linear increase). For the linearly decreasing release, the duration of the release is negligible, even considering separate release rates for virus and DCs. The numerical global minimum for the decreasing linear release rate is 815.5 mm 3 . The tumour dynamics for different linear release curve are plotted in Figure 5e,f.

Sigmoidal Release Profile
The final release curve considered was a sigmoidal function f s ptq, given by Equation (14). Integrating and simplifying using the constraint in Equation (15) gives We imposed additional constraints on the sigmoidal function so that the function was always positive, i.e., when k ą 0, the function within the "ln" would always be positive, and when k ă 0, the denominator was always negative.
Assuming initially that the characteristics of the release curve for virus and DCs are equivalent, gives the tumour size at day 20 in Figure 6a for a release period of 5 days, 9 days and 13 days. The sign of k determines whether the sigmoidal function is increasing or decreasing. This explains the clear dynamical switch at k " 0. For increasing sigmoidal release curves (k ą 0), the tumour size is minimised when the mid-point x 0 is larger. The overall minimum tumour size is achieved for increasing release curves when the midpoint is decreased below zero. The release period t r influences the minimum achieved with t r " 9 achieving a lower tumour volume than t r " 5, and the planes associated with t r " 9 and t r " 13 approaching the same limit for x 0 Ñ 0.
Investigating this further, the release time t r was fixed to 9 days and the value of k varied for virus and DCs for x 0 "´1, 1 and 3, see Figure 6b. Lower values of k for the DCs reduced the tumour volume when the midpoint x 0 was positive. As x 0 decreases, the variance in the plane reduces, and we see, overall, the minimum tumour burden is achieved irrespective of k for x 0 "´1. (f) Figure 6. Effectiveness of sigmoidal gel-release rate f s ptq (Equation (14)) for DCs, u DC ptq, and virus, u V ptq. In panel (a), the tumour size at day 20 as a function of k, x 0 and t r is plotted, where the parameters were common for virus and DC release functions. In panel (b), the tumour size at day 20 is plotted as a function of varying k for virus and DCs with t r " 9 and planes representing x 0 "´1, x 0 " 1 and x 0 " 3. In panel (c), the tumour size at day 20 is plotted as a function of varying x 0 for virus and DCs with t r " 9 and planes representing k "´1, 1 and k " 3. In panel (d), the tumour size at day 20 as a function of t r for virus and DCs where x 0 " 1 and k " 1. The red points in panels (c,d) correspond to the simulated release profiles in panels (e,f). In panel (e), k " 1, t r " 9 and x 0 " 0.5 for u V and x 0 " 3 for u DC . In panel (f), k " x 0 " 1 and t r " 1. Simulating varying x 0 ą 0 for virus and DCs and fixing k "´1, 1 and 3, we see that the value of the midpoint needed to minimise the tumour burden depends most heavily on the midpoint for u DC ptq, Figure 6c. If k ą 0, then a lower x 0 is needed, whereas the opposite is exhibited if k ă 0.
Unsurprisingly, simulating the sigmoidal function for fixed values of k " 1 and x 0 " 1 and varying t r for DCs and virus, we have similar dynamics to the linear and constant release curves, see Figure 6d. For an increasing sigmoidal function, the optimal tumour size is obtained when the virus is released quickly and the DCs have an extended release. Two model simulations are provided in Figure 6e,f to demonstrate the different sigmoidal release functions considered.

Implementing a Genetic Algorithm to Determine the True Optimal Release Curves
In the previous sections, the gel-release dynamics have been qualitatively linked with the tumour burden on day 20. To determine the optimal parameter sets for the constant, linear and sigmoidal dosage functions, the tumour volume on day 20 was minimised under each release profile using Matlab's genetic algorithm function ga [53]. Genetic algorithms are heuristic global optimisation routines inspired by natural selection that are frequently employed to estimate parameters in computational biology models. They have previously been applied to study optimal dosing routines in immunology [54,55] and combined oncolytic virotherapy and immunotherapy [29].
The tumour growth under the optimised constant, linear and sigmoidal release curves for virus and DCs are pictured in Figure 7. The parameter values from the genetic algorithm for each optimal release curve are given in Table 5. The constraints applied for the numerical simulations were also applied here for the different release curves. The optimal decreasing linear release curve has not been plotted as we know from the numerical simulations that it did not achieve a significant global minimum. Additionally, we found that the decreasing sigmoidal function's optimal parameters qualitatively matched that of the increasing sigmoidal function.  (14)) for Ad/IL12/GMCSF and immature DCs from a gel. The virus and DCs were allowed to have individualised curves. Figure 7 is a plot of the tumour volume under the five optimised curves.  Figure 7c shows that all three optimised release curves are able to reduce the tumour burden below the average size from Oh et al.'s experiments [10]. The dynamics of the immature APC population (D L`DS`AI ) are similar for the constant and sigmoidal release curves, see Figure 7a, whereas the optimal linear release profile delays the initiation of the APC accumulation. As the virus dynamics are the same for the three release curves (Figure 7b), this suggests that the delay on the immature APC population expansion causes the decrease in tumour size under the optimised linear release therapy in Figure 7c.

Constant Linear (Increasing) Sigmoidal
Overall, the parameters returned from the genetic algorithm matched the qualitative results of the numerical investigation for the constant and linear release functions. Both release curves required less than a day for the virus release and approximately 13 days for the DCs. The steepness of the linear curve a was much higher for the virus than the DCs, with the DC's value of a close to that plotted in Figure 5c. The sigmoidal release curve had the same dependence on the release curves for viruses and DCs as the constant and linear release period. Interestingly, the midpoint for both functions was significantly negative, meaning the curve would be similar to an increasing linear release profile. These results all suggest that the added complexity of the sigmoidal function with the time-varying gradient does not improve the efficacy of the gel-release mechanism.

Discussion
In this work, we developed a system of ordinary differential equations for the sustained release of virus and immature dendritic cells from a gel and the resulting interaction with a population of tumour cells and immune cells. Using our model, we recreated the in vitro and in vivo experiments of Oh et al. [10] (Figures 2 and 3). Using the calibrated model, we have been able to provide insight into how the gel-release profile may be optimised to improve treatment.
In the sequential fitting, attributes of the model were calibrated to the tumour response to combined Ad/IL12/GMCSF and immature DC therapy (Figure 3). To optimise the model to the Ad/IL12/GMCSF treatment of Lewis Lung Carcinoma (LLC) tumours by Oh et al. [10], we fixed most parameters to values obtained from our previous optimisation of the Ad/IL12/GMCSF treatment of melanoma [36], fitting only the virus infection rate, β, and killer T cell killing rate, κ (see Figure 3b). Both these parameters decreased for the treatment of LLC tumours compared to melanoma, suggesting that different tumour cells lines will ultimately result in different virus and immune characteristics.
From the optimisation of the single immature DC injection (Figure 3c) and the combined DCs+Ad/IL12/GMCSF injection (Figure 3d), we found that immature APCs are stimulated more rapidly by infected as opposed to uninfected cells, parameters s AI and s AU , respectively. Additionally, we note that the killing rate, κ, decreased for Ad/IL12/GMCSF alone, compared to DC+Ad/IL12/GMCSF. We hypothesise that the addition of immature DCs results in a heightened antitumour immune response, which results in a "division of resources" between the antitumour and antiviral immune response.
Oh et al. [10] reduced the tumour size by 50% on day 20 by using a GHPA hydrogel to release the DC+Ad/IL12/GMCSF treatment (Figure 3e). Using our calibrated model (Table 3), we numerically investigated alternative gel-release profiles to determine whether we could reduce the tumour burden further. Conserving the total amount of virus and DCs released from the gel, we first considered a constant gel-release profile. Simulating different release periods, t r , for the virus and DCs we found a minimum tumour size is attainable when the virus is released over a day, and the DCs are released for 13 days (Figure 4). This suggests that an initial burst of virus reduces the tumour volume through infection and stimulates an immune response that, under a prolonged constant release of DCs from the gel, results in a significantly reduced tumour volume.
A natural extension was then to examine the tumour volume under different linear releases. If the release period, t r , is fixed for both the DCs and the virus, then release periods between t r " 5 to 14 days reduce the tumour burden most significantly (Figure 5a,b). The global optimal of this regime occurs when a « 0, i.e., a constant release profile, suggesting that the optimal protocol is constant administration if release periods for the virus and DCs cannot be distinct. For unique release times, the optimal release period is similar to the constant release profiles: virus released rapidly and the DCs released slowly over time (Figure 5c). Overall, the gel designed by Oh et al. performed optimally and there are many increasing linear release profiles that would have performed significantly worse.
Considering the more complex sigmoidal gel-release profile, we found that the tumour burden was reduced when the midpoint was below 0 and the curve was an increasing function (Figure 6a,b). In other words, when the gradient of the release rate was reducing and tending towards a constant release profile. Furthermore, explicitly increasing the steepness, k, of the sigmoidal curve always increased the tumour size (Figure 6b,c). Unsurprisingly, the optimal release period for viruses and DCs was qualitatively similar to that of the constant and linear release periods (Figure 6d), implying that it is primarily the release period that optimises the treatment irrespective of the release curve.
Implementing a genetic algorithm, we were able to determine the true optimal constant, linear and sigmoidal release curves (Table 5 and Figure 7). Unsurprisingly, the optimal release curves each had a short release time for the virus and an extended release time of 13 days for the DCs. Comparing the impact of the different optimal release curves (Figure 7c), we see that the difference in the tumour size was minimal. Although the optimised linear release curve reduced the tumour size the most, all release curves improved the treatment by a similar amount. This suggests that extending the complexity of the release mechanisms may not be necessary, as a constant release curve performs almost as well as a linear release curve, and the optimised sigmoidal curve is more constant than sigmoidal. Additionally, this also suggests that for the constraints considered on the release curves it is not possible to optimise the gel so that the tumour is completely eradicated.
Previously, Barish et al. [38] optimised the scheduling of discrete doses of oncolytic virus and DCs. They showed that higher doses of DCs were needed to result in a robust optimised therapy. The fact that our results suggest the need for sustained release of DCs complements their findings and demonstrates that therapy can be optimised when the timing of DC administration is optimal. Our results also align with the results of Wares et al. [37], where they showed using a calibrated model that it is more effective to treat a tumour with immunostimulatory oncolytic viruses first and then follow up with a sequence of DCs than to alternate oncolytic virus and DC injections.
A limiting assumption of this work is that the complex gel mechanics are modelled phenomenologically. As the goal was to investigate the effect of changing the function u DC ptq, we believe it is not necessary to model the gel-release kinetics explicitly and, as qualitatively the model formulation for the measurements with and without the gel is able to capture the model dynamics (see Figure 2), this was a suitable simplification. Future work, however, would look to develop a more appropriate realisation of the gel kinetics and investigate mechanistically the optimal release profile. Another simplification in the model is that the rate at which immune cells are stimulated is independent of the antigen type. In reality, immune cells are antigen specific [18] with either tumour-antigen-or virus-antigen-specific cells. Future extensions of this work could consider two populations of immune cells similar to the work of Mahasa et al. [56] and investigate further why optimised gel-release profiles require an extended DC release.
Using a system of ODEs, we were able to replicate the experimental conditions of DC+Ad/ IL12/GMCSF therapy and determine possible improvements to the gel-release profile. This model can be applied to other viro-immunotherapies that aim to activate the adaptive antitumour immune response. As we only consider simple formulations for the gel-release rate, this model is easily translatable to other sustained-release delivery systems where the release functions u V ptq and u DC ptq may or may not be known. Our results on the optimal timing of virus and DC treatment provide insight into how future viro-immunotherapies using DCs need to be designed. From this work, we present a mathematical model that can be used to provide insight into other oncolytic virotherapies or viro-immunotherapies and predict the impact of sustained-delivery systems.

Conclusions
The mathematical model of tumour cell, virus and immune interactions derived in this study is a good representation of the response to viro-immunotherapy. We have shown how a sequential optimisation of experimental data provides insight into the underlying biology of the interaction. From our numerical optimisation of the gel developed by Oh et al. [10], we were able to show which key characteristics optimised or reduced the gel's efficacy. We found that if the virus is released instantaneously and the DCs are released over 13 days, then the tumour size under treatment can be improved. In turn, we were able to show that additional complexity in the gel-release profile does not add any improvement to the treatment outcome. Ultimately, to improve this treatment dramatically, it is not sufficient to solely manipulate the gel: more would need to be done to change the underlying efficacy of the virus and immune response to ultimately improve treatment.