Application of control theory in a delayed-infection and immune-evading oncolytic virotherapy

1.   Introduction

Oncolytic viruses are genetically engineered viruses that preferentially target and destroy cancer cells ^[1]. This emerging anticancer strategy exploits the lytic nature of viral replication to enhance the killing of malignant cells. Over the past decade, many oncolytic viruses have been tested in clinical trials ^[2], with one of these approved by the FDA for treatment of metastatic melanoma ^[3]. With the advancement of a growing number of oncolytic viruses to clinical development, treatment protocol is becoming a crucial factor for achieving optimal therapeutic efficacy.

Host immunity and the tumor microenvironment contribute significantly to inefficient virus delivery ^[4]. To overcome these obstacles, a range of novel mechanisms have been developed: gel-based mediums, nanoparticles, immunomodulatory agents and molecules that can manipulate the tumor microenvironment ^[4]. Gel-based mediums have been used successfully to improve virotherapy, chemotherapy and immunotherapy ^{[5,6,7,8,9,10,11]}. In these therapies, therapeutic agents are either loaded onto gels or coated to provide sustained therapy release and therapeutic efficacy. Nanoparticles have been investigated as a viral DNA and RNA delivery system as they can be engineered to have a decreased immune response ^[12], and their physical properties can be used to provide controlled viral release and diminish infectivity, maintaining an elevated local concentration ^[4,13,14]. Alternatively, polymer-based nanomaterials, such as polyethyleneglycol (PEG), have been shown to be effective at shielding particles from the extracellular environment and preventing clearance ^[12].

Each of these mechanisms looks to extend the treatment activity by controlling delivery and modulating immune system involvement. Using these therapeutic devices it could be possible to control viral diffusion through the tumour bulk prior to initial infection and avoid activation of immune clearance. Previously, Jenner et al. ^[15] developed a Voronoi cell-based model that predicted that modifying viruses to delay their initial infection of tumour cells and avoid immune clearance could improve virotherapy. They showed that the length of the delay before initial infection has an impact on the efficacy of this modified treatment. While insightful, due to the complexity and computation time of their model they did not investigate what proportion of delayed and non-delayed virus would be optimal. Additionally, they were not able to consider whether a distribution of delayed-infection times throughout the injected virus population could improve treatment. In this work, we develop a deterministic formulation for their Voronoi Cell-Based model to be able to analyse in-depth their treatment predictions using an optimal control approach.

Mathematical modelling has been used to improve the understanding of oncolytic viruses. Systems of ordinary differential equations (ODEs) that consider the basic interaction between uninfected cells, infected cells and viruses have been used to suggest significant improvements to treatment protocols ^{[16,17,18,19]}. More complex systems of ODEs that consider the heterogeneity in the cancer cell cycle have also been used to determine effective dosage protocols using in silico clinical trials ^[20]. Partial differential equations (PDEs) have also been successful at understanding treatment pitfalls and suggesting improved therapies by considering the spatial aspect of this therapy ^[21,22,23].

Optimal control theory is a useful approach used to understand optimal viral characteristics and dosage protocols ^{[17,18,20,24,25]}. Finite-horizon optimal control provides a useful tool to determine ways of minimising a variable in an ODE system. Zurakowski and Wodarz ^[25] used this method to find the constant concentration of drug that minimised the lowest excursion of the total tumour size during the optimisation horizon. Genetic algorithms are heuristic global optimization routines frequently employed to estimate parameters in computational biology models. Cassidy and Craig ^[20] used this method to generate personalised optimal schedules for patients in an in silico virotherapy trial.

Coating an oncolytic virus to delay infection and evade the immune system using either gel-based mediums or nanoparticles could be instrumental in improving oncolytic therapy. Drawing on previous modelling of the delayed-infection and immune-evading coated oncolytic virus by Jenner et al. ^[15], we derive a PDE system for this therapy and investigate its potential. Using optimal control theory, we determine a dosage combination for the modified (coated) virus that can minimise the tumour size. Simultaneously, we investigate the impact of the viral clearance on the outcome of therapy by investigating the optimal dosage regime for different viral-clearance models.

2.   Model formulation

As previously demonstrated in the work by Jenner et al. ^[15], the efficacy of oncolytic virotherapy can be improved by altering the virus's delivery mechanism to delay infection and avoid immune clearance. To investigate the applicability of this suggested treatment improvement and to optimise the protocol, we have developed a system of PDEs based on their model, that incorporates the interaction between an oncolytic virus and a population of tumour cells. Virotherapy is predominantly administered intratumourally or intravenously, resulting in individual infection sites and a non-homogeneous distribution of the virus throughout the tumour. We simplify this aspect to consider a single injection of virus that initially coats the tumour periphery uniformly, and consider the dynamics of the virus-tumour interaction are radially symmetric. To consider the impact of clearance on the virus's effectiveness, we develop four different modelling assumptions. Parameters for the model were then taken from previous model optimisations to data.

2.1.   Model description and assumptions

Consider a density of susceptible tumour cells $S(x, t)$, growing under a logistic growth rate with a normalised carrying capacity, where $x$ is the radial distance from the centre of a circular tumour. Two types of viruses are considered in the model: coated virus particles $V_{C}(x, t)$ and uncoated virus particles $V_N(x, t)$. We assume that the coating inhibits the virus from infecting susceptible cancer cells. Uncoated virus infects susceptible cells with rate constant $\beta$, creating a density of infected tumour cells $I(x, t)$. Infected cells die due to lysis at rate $d_I$ and produce $\alpha$ uncoated viruses. Infected cells grow at the same logistic rate as the susceptible tumour cells. In the model developed by Jenner et al. ^[15] they assumed virus-infected cells do not replicate; however, there is evidence that this can occur for different modified viruses ^[26,27]. Both susceptible and infected tumour cells have a diffusion coefficient $D_T$.

The coating on the coated virus particles degrades at rate $d_C$. We assume coated viruses are not detectable by the immune system and as such are removed from the system, through loss into the tissue, lymphatic system or vasculature at a slow rate $\delta_C$. Uncoated virus, however, is detectable by the immune system and as such is cleared more rapidly at rate $Imm(S, I, V_N)$. This function represents the immune-driven clearance of the virus and will be described in the following subsection. Realistically, virus diffusion would be influenced by the fluid flow through the tumour and obstructions in the microenvironment, such as densely packed cells and vasculatures ^[15]. Since, we do not model additional factors such as these, we assume one-dimensional radially-symmetric diffusion is able to account for the average dynamics of virus diffusion. Coated and uncoated viruses, therefore, diffuse through a tumour with the same diffusion coefficient $D_V$.

To model the decay of coating from the virus in more detail, we consider $K$ levels of viral coating. The $K$th level representing fully coated virus, which degrades until there is no coating (i.e., to an uncoated virus). We model coated virus as having discrete levels of coating, and the density of virus particles with an $i$th level of coating, $1\le i\le K$, is given by $V_{C, i}(x, t)$.

A summary of the dynamics for coated and uncoated virus treatments of a growing tumour is presented in Figure 1. The corresponding partial differential equation (PDE) system is provided in Equation 2.1–Equation 2.7, with corresponding boundary conditions and initial conditions in Equation 2.8.

In this system, $t$ represents days and $x$ is the radial distance from the tumour centre in $\mu$m. The tumour grows within the fixed domain $[0, L]$ and the initial radius of the tumour, $P$, must satisfy $P < < L$. We model a single injection of virus at the periphery $P$ of the tumour and consider an initial combination of both uncoated and coated virus with different levels of coating. The initial injection of virus is administered at the periphery of a circular tumour and, as such, the area of the injection is $A_I = \pi[(P+\epsilon)^2-(P-\epsilon)^2]$, and as a result the total virus injected at the periphery of the one-dimensional problem considered is $V_{T} = \frac{V_0}{A_I}$. The Heaviside function $\mathcal{H}(x)$ is used to model the initial cell densities and virus injection. The proportion of the initial injection of virus, $V_0$, that contains each level of coated virus is given by $p_i$, where $i$ is the coating level. Note that $0\leq p_i\leq 1$ and $\sum_{i = 1}^K p_i\leq 1$.

No-flux boundary conditions are used for the tumour populations at $x = L$ to simulate that tumour cells do not leave the edge of the domain. Since virus particles are much smaller than cells, we assume they can diffuse past $x = L$ and be lost, whether into the surrounding tissue or the vasculature system, and we impose that $V_{C, i}(L, t) = V_N(L, t) = 0$. As we are modelling a radially symmetric tumour, no-flux boundary conditions are also applied for cells and virus at $x = 0$.

The biological mechanisms modelled in Equation 2.1–Equation 2.7 are further described as follows:

● In Equation 2.1, virus particles at the maximum coating level $V_{C, K}$ diffuse through the tumour with coefficient $D_V$. The coating decays at rate $d_C$, and as such, viruses leave the $K$th coating level at rate $d_CV_{C, K}$ and enter the subsequent lower-coating-level population $V_{C, K-1}$. Coated viruses are cleared at rate $\delta_C$.

● In Equation 2.2, virus particles at the next maximum coating level $V_{C, K-1}$ diffuse through the tumour with coefficient $D_V$. Viruses from the $K$th coating level join the $(K-1)$th coating level at rate $d_CV_{C, K}$, and viruses leave the $(K-1)$th coating level at rate $d_CV_{C, K-1}$ and enter the subsequent lower coating level $V_{C, K-2}$. Coated viruses are cleared at rate $\delta_C$.

● Equation 2.3 is the general formulation for coated viruses. Virus particles at the $i$th coating level diffuse through the tumour with coefficient $D_V$. The coating decays at rate $d_C$, meaning viruses from the coating levels above enter at rate $d_CV_{C, i+1}$, and viruses at the currect coating level leave at rate $d_CV_{C, i}$. Coated viruses are cleared at rate $\delta_C$.

● In Equation 2.4, virus particles at the last coating level diffuse and are cleared at the same rates as the preceding coating levels. As this is the last coating level, virus lose the remaining coating and enter the uncoated virus population $V_N$ at rate $d_C V_{C, 1}$.

● In Equation 2.5, uncoated virus, $V_N$, diffuses through the tumour with the same diffusion coefficient $D_V$. Viruses are cleared by an immune response described by the function $Imm(S, I, V_N)$, which can be dependent on the susceptible and infected tumour cells, $S$ and $I$. Uncoated viruses infect susceptible and infected tumour cells at rate $\beta(S+I)V_N$, and new uncoated viruses are created through cell lysis which occurs at rate $d_I I$ and produces $\alpha$ new viruses.

● In Equation 2.6, susceptible tumour cells diffuse with coefficient $D_T$. They undergo proliferation at a logistic growth rate $r$, normalised to a carrying capacity of 1. Susceptible cells are infected by uncoated-virus at rate $\beta S V_N$.

● In Equation 2.7, infected tumour cells also diffuse with coefficient $D_T$. They similarly undergo logistic growth. Infected cells arise through infection of the susceptible tumour cell population at rate $\beta S V_N$. These cells undergo lysis at rate $d_I I.$

2.2.   Modelling the immune clearance of uncoated virus

The immune response to an oncolytic virus is complex and depends heavily on the tumour type, level of heterogeneity and the genetic modification of the virus. In this section, we develop four hierarchically-related formulations for the clearance rate of uncoated virus by the immune system $Imm(S, I, V_N)$. Each of the models considers a different aspect of the immune response to virotherapy. Our primary goal is not to prove which model is correct, but to investigate how different assumptions play a role in the outcome of therapy.

Virus particles that infect tumour cells can instigate both an anti-viral and anti-cancer immune response. In general, an anti-viral response is one that targets either extracellular virus particles or infected cells and an anti-cancer response is one that elicits an immune response that can apoptose both infected and susceptible tumour cells. The clearance of viruses by the immune system has been shown to be a major hindrance to oncolytic virotherapy ^[28]. In this work, we choose to model the clearance of extracellular virus particles by the immune system with varying degrees of dependency on the uncoated virus density, susceptible tumour cell density and infected tumour cell density.

The simplest way to model clearance of the uncoated virus is at a rate proportional to the number of particles:

where $\delta_N$ is the clearance rate constant. This model is regularly used in deterministic models ^[16]. It assumes that the immune system activation is not dependent on the activities of the virus or the presence of the tumour, but on the density of the virus population.

In reality, the anti-viral immune response would be proportional to the number of infected cells as the immune system becomes stimulated by the presentation of viral antigen on the surface of infected cells ^[32]. As the number of infected cells increases and more immune cells are activated, the strength of the immune response will also increase. As such, we also consider the clearance of virus particles as proportional to the number of infected tumour cells:

In this model, we assume that the proportion of the tumour that is uninfected does not play a role in the activation of an immune response.

In addition, the presence of tumour cells can elicit an immune response resulting in an influx of immune cells at the tumour site. These cells can then become aware of the virus particles and contribute to viral clearance. Assuming, therefore, that the immune system can be stimulated by the presence of both susceptible and infected tumour cells provides a third viral clearance model:

where the activation of clearance is equally proportional to both types of tumour cells.

In Equation 2.11 we assume that the rate at which susceptible and infected cells stimulate the clearance of virus particles by immune cells is equivalent; however, this may not be the case. In the adaptive immune system, cells such as dendritic cells and macrophages are stimulated by the presentation of antigen ^[32]. Susceptible cells will only present tumour antigen (as they do not contain any viruses) whereas infected cells will present virus antigen (as they have been infected by viruses) and tumour antigen. As such there may be two different immune responses: an anti-cancer response and an anti-viral response. These will affect the clearance of the virus in different ways. As expected the virus immune response will rapidly clear virus particles, as this immune response is specific to viruses. The tumour immune response will increase immune cell presence at the tumour site, which we assume increases the likelihood of virus immune stimulation. As such, the viral clearance can be modelled with different dependencies on the tumour cell populations:

where $\delta_{NS}$ and $\delta_{NI}$ are the respective clearance rates of the anti-tumour and anti-viral responses.

In this work, we investigate the dynamics of Equation 2.1–Equation 2.7 with each of the different models of the immune response to virotherapy. To date, no one has investigated how different viral-clearance models may impact the efficacy of therapy, and we aim to provide motivation for this discussion.

2.3.   Parameter values and model simulation

The majority of the parameter values used in the following analysis were taken from Jenner et al.'s investigation of a PEG-modified adenovirus conjugated with herceptin ^[30]. In this work, a reduced ODE system similar to the PDE system in this work was used to optimise the infectivity rate, replication rate, lysis burst rate and size (see Table 1). These parameter values were also similar to the ones used in the VCBM model that conducted the preliminary investigation into the efficacy of this viral coating ^[15].

The remaining parameters were then estimated at values relative to our biological understanding of the system. For the diffusion coefficients, we assumed the virus would diffuse faster than the tumour, as such we chose $D_V = 0.001$ and $D_T = 0.0001$. Additionally, to simulate the immune evasion of the coated virus we chose $\delta_C < < \delta_N$, and as such set $\delta_C = 0.0001$.

Lastly, we wanted a reasonable estimate for the degradation rate of viral coating. Oh et al. ^[9] reported that gelatin-hydroxyphenyl propionic acid based hydrogel loaded with oncolytic adenovirus degraded after 6 days. Alternatively, Croyle et al. ^[6] found that administration of polyethylene glycol-coated viral vectors extended viral-gene expression by 38 days. Degradation of viral coating may also be a function of particular molecular concentrations, with Tseng et al. ^[14] designed a hypoxia-responsive vector carrier that oxidises with lactate oxidase. As we are not modelling a specific coating device, but instead the theoretical concept of a coating that instigates a delay in viral infection, we chose $d_C = 0.1$, which is equivalent to the rate of cell lysis and a half-life for the coating of approximately 7 days. In this work, we examine a range of different initial proportions and levels of coating, which implicitly determines the impact of varying the coating's degradation rate.

The domain was fixed to $[0, 2000]$, i.e., the maximum radial growth of the tumour was $L = 2000\mu$m. The initial radius of the tumour was $P = 1000\mu$m and the width of the initial virus injection rim was $2\epsilon = 0.1\mu$m. To demonstrate the effectiveness of the coated therapy, model simulations for different levels of coating are plotted in Figure 2. Three coating situations are considered: no coated virus, one level of coated virus $K = 1$, and five levels of coated virus $K = 5$. It is clear that treatment does significantly better by 50 and 100 days when coated virus is administered. Additionally, when the coating level is increased to $K = 5$, it is clear the treatment is significantly more effective. This is most likely a result of the coating providing further intratumoural diffusion before infection and clearance can commence. All parameter values in the simulations were taken from Table 1 and the initial proportion of coated virus at all coating levels was equal. The clearance model for the uncoated virus was Equation 2.9.

2.4.   Parameter sensitivity

As a few of the parameters were estimated, we conducted a sensitivity analysis to determine the variability in our model. Fixing all parameters to the values in Table 1 and using the clearance model in Equation 2.9 and initial proportion of uncoated virus $p = 0.5$, we measured the change in the total tumour cells at 100 days from perturbations of $\pm10\%$ for the parameters noted in Figure 3. The lysis burst size $\alpha$ most significantly affected the tumour size at 100 days. This is not surprising, as increasing or decreasing the amount of virus created through lysis would have a significant effect on the infection and subsequent death of uninfected tumour cells. The next highest variability in the tumour volume was exhibited from perturbations in the clearance rate of coated virus, $\delta_C$, and the degradation rate of the coating, $d_C$. This suggests, that the dynamics of the coating mechanism will influence the tumour volume reached under treatment. Overall, since the magnitude of the change in tumour volume is not large, fixing these parameters and investigating an optimised scheduling for a delayed-infection coating can provide a reliable approximation to the potential of this treatment.

3.   Optimal maximum level of coating

To develop an optimised therapeutic protocol for the delayed-infection and immune-evading viral treatment we first investigated whether coating the virus was optimal when compared to uncoated virus. From this, we determined for the simple viral-clearance model Equation 2.9 the optimal proportion of coated to uncoated virus in the initial injection. Furthermore, we then investigated whether there is an optimal maximum level of coating $K$ for each of the viral-clearance models in Equation 2.9–Equation 2.12.

3.1.   Preliminaries

We consider two mathematical formulations to measure the efficacy of treatment. The first aims to minimise the susceptible tumour density at a final assessment time $T_f$ days after the initial treatment injection. In this metric, we ignore the infected cell population as once a cell becomes infected it will eventually die through lysis. Additionally, we wish to suppress tumour growth throughout treatment, as a treatment that allows for an excessive amount of tumour growth, irrespective of the size at a particular point in time, would be unrealistic. These two metrics of treatment efficacy result in the following mathematical formalisms:

where $S(x, {\bf{p}}, t)$ is the density of susceptible tumour cells $S(x, t)$ defined in Equation 2.1–Equation 2.7 for a given vector ${\bf{p}}$, which is the initial proportion of coated virus in each coating level.

Initially, we investigated the effects of a single coating level (i.e., $K = 1$), allowing for the proportion of initial coated virus $p = p_1$ to vary. In Figure 4(a), we plot the value of the objective function Equation 3.1 with respect to $T_f$ and $p$. In Figure 4(b), the optimal proportion for $p$ for the constraint in Equation 3.1 is plotted as a function of $T_f$. After initial injection, there is a short period of days where the treatment is optimal when there is no coated virus. As the final assessment time $T_f$ increases, we see that this optimal proportion reaches $p = 1$. This is crucial for defining $T_f$ in the objective functions above (Equation 3.1–Equation 3.2), as we need to set $T_f$ so that a non-trivial amount of coated virus is necessary to minimise the amount of susceptible cancer cells. Interestingly, there appears to be a critical threshold at 5 days, before which coated virus treatment is ineffective.

After 5 days, there is a non-zero proportion ($p > 0$) of initial coated virus that is optimal. From this we can deduce that coating the virus improves the efficacy of treatment. In this simulation, the viral clearance was modelled using the simple model in Equation 2.9. To investigate whether the viral-clearance model affects the effectiveness of the coated virus, we now investigate the optimal maximum number of coating levels, $K$, for each viral-clearance model presented in Section 2.2.

As in Figure 4, in Figure 5 and Figure 6, we plot our metric (Equation 3.1) of tumour growth under treatment with varying maximum coating levels $K \ge 1$. To simplify the study, we set ${\bf{p}}$ by fixing the virus's initial condition so that 50% of the virus is uncoated and the remaining 50% is equally distributed between the $K$ levels of coating, i.e., $p_i = 0.5/K$ for $i = 1, \ldots, K$. We model the efficacy of this therapy with each of the viral-clearance models defined in Equation 2.9–Equation 2.12.

For the viral-clearance model in Equation 2.9, dynamics of the metric (Equation 3.1 with ${\bf{p}}$ set as above) are plotted in Figure 5(a) and (b). This case considers the viral clearance to depend only on the density of virus. We observe that the model eventually reaches a steady state irrespective of the maximum coating level. Following the line in Figure 5(a) or (b), we see that the optimal number of coating levels increases with time $T_f$, similar to what was already observed in Figure 4. Assuming, instead, that viral clearance is proportional to the amount of infected cancer cells (Equation 2.10) gives the optimal coating levels over time $T_f$ plotted in Figure 5(c) and (d). For smaller time horizons ($T_f < 100$), a faster-acting virus is more optimal, so the optimal coating level remains low, but as $T_f$ increased beyond 200 days, the optimal coating level increases to around $K = 20$.

We wished to then examine whether modelling viral clearance as dependent on the susceptible tumour population would result in similar dynamics to Figure 5. Modelling the viral-clearance rate to be $\delta_N V_N \int_0^NS(\tilde{x}, t)d\tilde{x}$, results in almost complete tumour eradication, irrespective of the maximum coating, see Figure 6(a). This may be a result of the rapid infection of susceptible cancer cells reducing the viral clearance rate. In comparison, when both susceptible and infected cells contribute at different rates to the viral clearance term (Equation 2.12), we see the observed dynamics in Figure 6(b). Under this assumption, the tumour size over time depends more significantly on the maximum coating level. When the density of susceptible and infected tumour cells contribute equally to the clearance rate (Equation 2.11), we obtain the dynamics in Figure 6(c) and (d).

It is clear from Figure 5 and Figure 6 that the different viral-clearance models result in significantly different treatment efficacy. Additionally, the resulting impact on the optimal maximum coating level is noticeable. Based on this, we chose to determine the optimal proportion of virus in each coating level using the Equation 2.11. The reason for this is that Figure 6(d) illustrated a change in the optimal maximum coating as time horizon $T_f$ increased. Additionally, we feel that allowing both susceptible and infected cells to elicit and immune response that results in direct (or indirect) viral clearance is more biologically reasonable. To set our final time $T_f$ and maximum number of coats $K$ to consider, we observe the minimum point of the objective function in Figure 6(c). It is clear that the cancer regrows to its initial size after about 300 to 400 days and achieves its minimum point for $K\leq30$ and $T\leq300$, so we fix these as our final time $T_f$ and maximum coating level $K$ in the following optimisation.

3.2.   Objectives

Based on the observations in the previous section, we define the following optimal control problems:

and

where $T_f = 365$ and $K = 30$ have been chosen large enough so that a reasonable optimal point may be obtained. Both Equation 3.3 and Equation 3.4 are based on our original optimisation formulations in Equation 3.1 and Equation 3.2. Each objective function aims to optimise the proportion ${\bf{p}}$ of coated virus in each coating level, where there is a maximum of $K = 30$ coating levels. The solution to Equation 3.3 optimises ${\bf{p}}$ at the final time $T_f = 365$, while the solution to Equation 3.4 optimises ${\bf{p}}$ for the total cancer cell density over 365 days. The objective given by Equation 3.4 does not guarantee a minimum tumour size at the final time; therefore, the resulting optimal ${\bf{p}}$ may be different.

4.   Optimal fraction of the initial injection in each coating level

As demonstrated in the previous section, there is a non-trivial optimal maximum viral coating that reduces the total tumour size, see Figure 6(c) and (d). To investigate whether the treatment efficacy could be improved, we relaxed the condition on the initial proportion of virus in each coating level. Below we describe the method undertaken to obtain the optimal proportion of initial virus in each coating level. We considered the objective function set-up described in Section 3.2.

Fixing the optimal control problems to be defined by Equation 3.3 and Equation 3.4, we first generated a vector for the initial proportion in each coating layer ${\bf{p}} = [p_1, p_2, ....p_K]$ where values for $p_i$ were generated from the uniform distribution sets constraining $\sum_i p_i\leq1$ and $p_i\geq0$. This was done several times, and the fmincon optimisation method in MATLAB was then applied. This optimisation method uses the interior-point algorithm ^[29] for solving the optimal control problem. We investigated different objective final times $T_f = 91,182,274$ and $365$ for Equation 3.3 and observed the changes of optimal vector ${\bf{p}}$ as the objective time increased. The results for the optimisation of for Equation 3.3 and Equation 3.4 are presented in Figure 7 and Figure 8 respectively.

The results of the optimisation of both Equation 3.3 and Equation 3.4 imply that a mixture of high-coating levels and low-coating levels are necessary to achieve a minimum tumour size. If we consider the physical meaning of high and low-coating levels, this may correspond to thick and thin coating. As such, there is some non-trivial dependence of the treatment efficacy on the decay rate from the coated to uncoated virus states.

Moreover, the final-time optimisation results in Figure 7 suggest that the maximum coating increases as time increases. When we classify classes of coating levels as either thin or thick (i.e., $V_\text{thin}$ and $V_\text{thick}$) the optimal initial proportion of virus can be summarised as

which allows us to clearly stratify how the proportion of thickly versus thinly coated virus changes with time. As time increases, the optimal initial injection requires thicker coated virus. Overall, the optimal ratio of uncoated to thin and thick virus is similar irrespective of the final time considered, with $p_0:p_{\text{thin}}:p_\text{thick} \simeq 1:5:4$. This does not include the case where $T_f = 91$, the shortest objective final time. On the other hand, the results for optimisation of the integrated tumour size, Figure 8, is slightly different with the ratio of uncoated to thin and thick coated virus as $1:6:3$.

5.   Discussion

Oncolytic virotherapy is a promising cancer treatment due to its ability to induce tumour-specific cell lysis and as a delivery vector for other tumour-targeting agents. One major pitfall of these viruses is their inhomogeneous diffusion throughout the tumour, a result of immune clearance ^[1] or high multiplicity of infection of a single tumour cell. In this work, we look at overcoming these obstacles by simulating a delayed-infection and immune-evading coated virus. We optimised the initial dosage of this virus and demonstrated that modifying a virus in this manner could improve treatment.

Experimentalists have developed ways to avoid the immune clearance of viruses through modifications with gel-based mediums, coating with polymers such as polyethylene glycol (PEG) ^[31], immunogenic-engineered viral capsid and additional manipulations of the tumour microenvironment with hormones such as relaxin. In previous work, Jenner et al. modelled a PEG-coated oncolytic adenovirus with a system of ODEs and analysed the virus's ability to avoid immune detection. In other previous work by Jenner et al. ^[15], a Voronoi cell-based model was developed to show the effectiveness of a virus that was able to avoid immune detection and delay its infection. Extending both of the modelling frameworks in ^[30] and ^[15], we developed a system of PDEs that considered the spatial interaction of an oncolytic virus with a growing population of susceptible tumour cells. We investigated the efficacy of an undefined coating mechanism that delayed the infection of cells and avoided immune clearance. We also developed four viral-clearance models (Equation 2.9–Equation 2.12), that modelled different aspects of the clearance process.

Coating the virus to delay infection and avoid immune clearance improves the efficacy of virotherapy (see Figure 2). Simulating the model with and without coated viruses we found that coated virus diffused farther intratumourally before the onset of infection. This resulted in a significant increase in the infected tumour cell number and reduced the overall tumour size (see Figure 2 last column). To investigate the efficacy of the coated virus treatment, we considered one level of coated virus and varying proportions of initial virus that was coated (Figure 4). For the first 5 days, the optimal proportion of initially coated virus was negligible and the coated therapy was not as effective. But as time increased, the optimal proportion of coated virus in the initial injection increased until all virus needed to be coated to achieve a tumour minimum. This suggests that coating the virus improves virotherapy in a time-dependent manner.

As the coating mechanism we are modelling in this work is theoretical, we decided to extend the model to consider $K$ levels of discrete coating. We then determined the optimal maximum coating level that would minimise the tumour size, given equal proportions of coated virus in the initial injection (Figure 5 and Figure 6). The optimal maximum coating was simulated for a range of different viral-clearance models. The goal was to quantify how the efficacy of treatment depended on both the maximum coating level and the viral-clearance assumptions.

Initially, we considered viral clearance to be proportional to the amount of virus (see Figure 4(a) and (b), and Figure 5(a) and (b)). Simulating Equation 3.1, we found that only a small maximum coating level ($K\le 5$) was necessary for the first 250 days, after which the optimal maximum coating increased linearly with time (Figure 5(b)). This suggests that the effectiveness of the coating is limited by the initial rapid clearance of the uncoated virus. As such, the coated virus is needed to instigate the infection and lysis of tumour cells. As time goes on, we see a tipping point at 250 days, after which the ability of the coated virus to diffuse is sufficient to warrant a larger maximum coating.

When the rate of clearance depends on the number of either susceptible or infected cells, the maximum coating for the virus appears to approximately stabilise for $20 < K < 25$ (Figure 5(d) and Figure 6(d))). In addition, when the viral clearance is equally proportional to the number of infected and susceptible tumour cells, the tumour is able to evade elimination and reaches a steady state value (see Figure 6(c)). The dynamics with this model are similar to those where the clearance is proportional to the number of virus particles (Figure 5(b)); however, the optimal maximum coating is completely different. This difference in the optimal maximum coating trace is due to a larger initial clearance in the case of Figure 6(c). This requires a larger delay in the initial viral infection and therefore a larger maximum coating of the virus.

Realistically, the clearance of viral particles will not be equally proportional to the susceptible and infected cells (Equation 2.12), see Figure 6(b). Simulating this viral-clearance model results in a linear relationship between the optimal number of coats and time (figure not included). As time increases the number of coats required to achieve a minimum tumour burden increased with a gradient of 0.25. In contrast, when the clearance of viruses is equally proportional to the number of susceptible and infected cells, the optimal number of coats roughly stabilises after 300 days. In other words, to reduce the tumour volume over a long period of time approximately the same number of coats is needed when we consider immune clearance to be equally proportional to both susceptible and infected tumour cells. Whereas, if the clearance is not equally proportional to both cell types, the relationship between the number of optimal coating levels and time is linear. Overall, this suggests that the type of immune response exhibited in response to virotherapy can be crucial in determining, long-term, how many coating levels is optimal.

The different viral-clearance models considered results in significantly polarised dynamics: tumour eradication or tumour stabilisation. For example, the treatments' effectiveness is qualitatively similar when the virus is cleared at a rate proportional to the density of infected or susceptible tumour cells (Figure 5(c) and Figure 6(a)). This is possibly explained by the fact that the uncoated virus would rapidly infect susceptible cells resulting in an approximately equivalent rate of viral clearance that is slow enough to allow for tumour eradication. In comparison, when viral clearance is proportional to both susceptible and infected cells, the tumour manages to stabilise (see Figure 6(c)). This suggests that there is a rapid switching between effective and less effective treatment based on the viral clearance rate.

In summary, modelling clearance as dependent on either the total susceptible or infected population results in quasi-eradication of the tumour, irrespective of the number of coating levels (see Figure 5(c) and Figure 6(a)). In comparison, assuming that the virus is cleared at a rate equally proportional to both susceptible and infected tumour cells results in tumour stabilisation (see Figure 6(c) and Equation 2.11). Interestingly, modulating the clearance rate's dependence on this population provides a significantly non-uniform long-term result to the coating level of the virus (see Figure 6(b)). This suggests that the clearance rate is significant and should be considered more carefully in future when modelling virotherapy.

Following these simulations, we wished to investigate the effects of relaxing the condition on the initial proportion of coated virus at each coating level. From the preliminary simulations (Figure 4–Figure 6), it was possible to define the final objective time and maximum number of coatings ($T_f = 365$ and $K = 30$), so that all observable possible behaviours would occur. We then defined two objective functions (Equation 3.3 and Equation 3.4): the first was to minimise the tumour size at a specific point in time $T_f$ days and the second was to minimise the overall tumour size throughout the entire time period.

If we consider then the optimal initial composition for specific times ($T_f = 91,182,274$ and 365), we see that the proportion of thicker coated virus ($K > 10$) increases as the final time of the experiment increases. Overall, solutions for $T_f > 91$ have very similar ratios of uncoated, thinly coated ($0 < i < 10$) and thickly coated ($i\ge10$) virus, see Equation 4.1. This suggests that the thicker coated virus takes a longer time to become effective, which is why in the 91 day window we do not see any thicker coated virus ($i > 10$) being necessary for achieving an optimal tumour minimum.

Overall, for a long time period, a significant proportion of the virus needed to have a thicker coating around $i \approx 13$, see Figure 7(d) and Figure 8. Minimising the tumour size on day 365 requires a different initial composition than minimising the tumour concentration over the entire 365 day period. This is not surprising as minimising the tumour volume over a period of time requires a different dynamical approach, as opposed to simply minimising the tumour volume at 365 days. What is surprising is the similarity between the two optimal initial conditions. Both require a significant amount of thinly coated virus ($i \in (0, 5]$), demonstrating that for long-term treatment effectiveness, a short-term delayed-infected immune-evading treatment is key.

A future simplification of our model would be to consider a continuous coating of the virus, as opposed to the discrete coating levels we have modelled. To achieve this, the PDE system could be transformed into an age-structured model where the level of coating would become a continuous variable. We decided not to model it in this way as experimentally it would be easier to create a discretely coated virus that has specific characteristics.

Normally, oncolytic virotherapy is administered in multiple dosages (typically three dosages every second day). This work does not consider the effects of additional dosages on the tumour size. Future work could investigate whether there are different optimal injection configurations for subsequent injections. Alternatively, we could consider injections of solely uncoated or coated virus. For the current preliminary investigation we consider a single dosage to illustrate the potential of this therapy and the existence of an optimal solution.

Lastly, because we did not prove the uniqueness of the optimal solution, our solution does not guarantee the global minimum. In other words, it is possible for another solution to attain a smaller value of the objective function. As we optimised by choosing randomized ${\bf{p}}$ several times, we cannot establish whether our solution is a global minimum or not. Therefore, proving uniqueness and existence is another open problem.

A natural extension to this model, would be to consider the dynamics in a non-radially symmetric, heterogeneous tumour microenvironment as non-uniform spatial interactions between virus particles and tumour cells have been shown to significantly impact the outcome of virotherapy ^[15,21]. Additionally, since we have demonstrated the importance of immune clearance, it would be worth explicitly modelling the immune cells in a spatial environment with restrictions on their point of entry to the tumour, to determine whether this may impact the optimisation results.

The next stage of this work would be to investigate experimentally whether designing an injection with proportions of uncoated and coated virus defined by Figure 7 result in a reduced tumour size. From this, it would be interesting to then investigate why this combination of uncoated, thinly coated and thickly coated virus produces an optimal reduction in the tumour by investigating intratumoural dissemination as a function of time.

6.   Conclusion

Overall, our results suggest that it is possible to improve virotherapy by creating a delayed-infection and immune-evading virus through coating or another mechanism. We have demonstrated that this treatment can be further optimised by considering multiple coating levels. We have determined the composition of the initial injection that should be coated at thin and thick levels. We also simulated our system under different viral clearance assumptions and found that the rate of viral clearance is significant in determining whether treatment results in tumour eradication or tumour stabilisation. While exciting, it is still a theoretical principle and would need significant further experimental investigation to verify its promise as an improvement on the current treatment approach.

Acknowledgments

TL, ALJ, PSK and JL gratefully acknowledge the support for this work through the University of Sydney and Yonsei University Joint Research Funding Scheme. Furthermore, ALJ and PSK gratefully acknowledge support from the Australian Research Council (DP180101512). The work of JL was also supported by NRF-2015R1A5A1009350 and NRF-2016R1A2B4014178. The work of ALJ was also supported by the Australian Mathematical Society Lift-Off Fellowship and Australian Federation of Graduate Women Tempe Mann Travelling Scholarship.

Conflict of interest

The authors acknowledge that there were no conflicts of interest in the development of this work.