Treating cancerous cells with viruses: insights from a minimal model for oncolytic virotherapy

Abstract In recent years, interest in the capability of virus particles as a treatment for cancer has increased. In this work, we present a mathematical model embodying the interaction between tumour cells and virus particles engineered to infect and destroy cancerous tissue. To quantify the effectiveness of oncolytic virotherapy, we conduct a local stability analysis and bifurcation analysis of our model. In the absence of tumour growth or viral decay, the model predicts that oncolytic virotherapy will successfully eliminate the tumour cell population for a large proportion of initial conditions. In comparison, for growing tumours and decaying viral particles there are no stable equilibria in the model; however, oscillations emerge for certain regions in our parameter space. We investigate how the period and amplitude of oscillations depend on tumour growth and viral decay. We find that higher tumour replication rates result in longer periods between oscillations and lower amplitudes for uninfected tumour cells. From our analysis, we conclude that oncolytic viruses can reduce growing tumours into a stable oscillatory state, but are insufficient to completely eradicate them. We propose that it is only with the addition of other anti-cancer agents that tumour eradication may be achieved by oncolytic virus.


Introduction
Oncolytic virotherapy is fast becoming a successful cancer treatment. Recent advances in genetically modified cancer-killing viruses have shown increasing promise both experimentally and clinically (see Aghi & Martuza, 2005;Jebar et al., 2015;Lawler, Speranza, Cho, & Chiocca, 2017;Prestwich et al., 2008;Russell, Peng, & Bell, 2012;Wang et al., 2016). Many naturally occurring viruses are being investigated for their use in cancer treatments, for example: Herpes Simplex Virus, adenovirus, measles, reovirus, Vesicular Stomatitis Virus (Russell et al., 2012). These viruses are currently being tested in clinical trials and are used to treat a range of cancer types such as Glioma, Ovarian cancer, Sarcoma, Pancreatic cancer, Prostate cancer and Bladder cancer (Prestwich et al., 2008;Russell et al., 2012). Recent success has been reported in treating metastatic melanoma with the herpes simplex virus (Andtbacka et al., 2015). However, oncolytic virotherapy is considered to be at a conceptual stage and consistent success in cancer eradication or control still remains elusive. Much is still unknown about the sensitivity of oncolytic virotherapy to tumour and viral heterogeneity. For example, protocols detailing doses, treatment lengths etc, are not yet universally established.
Interest in oncolytic virotherapy arises from its capacity to combine tumour-specific lysis with delivery properties for other anti-cancer drugs (Lawler et al., 2017). Oncolytic viruses are engineered to express genes that promote viral replication within tumour cells and inhibit replication within healthy cells. As such, the viral load remaining after the tumour cells have been eradicated will not affect the rest of the body and will decay over time. Once inside a tumour cell, virus particles will replicate until they reach a tumourcell-filling capacity, at which point the cell will burst, undergoing lysis, and release the viral progeny inside. By replicating within tumour cells, oncolytic viruses are able to rapidly increase their population at the tumour site. This makes oncolytic virotherapy, in theory, a successful cancer treatment.
For example, two of the first mathematical models for oncolytic virotherapy were developed by Wodarz (2001Wodarz ( , 2003. Wodarz et al. modelled a virus interacting with a population of uninfected and virus-infected tumour cells using a system of two ordinary differential equations (ODEs). Based on this work, Komarova and Wodarz developed a generalized system of two ordinary differential equations for oncolytic virotherapy (Komarova & Wodarz, 2010). All of these models provided valuable insight into the dynamics between an oncolytic virus and tumour cell population despite the viral population being modelled implicitly.
One of the first models that explicitly modelled the viral population was developed by Bajzer et al. (2008). Conducting a bifurcation analysis of their system, these authors discovered stable oscillations in the tumour cell population emerging from a Hopf bifurcation. While the concept of an oscillating tumour may seem unusual, such behaviour has been observed experimentally. Most recently, Titze et al. (2017) developed a mathematical model to characterize the relationship between tumour cell growth and different oncolytic viruses. Their model, analogous to that of Bajzer et al. (2008), used a system of three ordinary differential equations to embody glioblastoma cell growth in vitro under treatment from three different oncolytic adenoviruses. Optimizing their model to data, they obtained parameter estimates and predicted long-term tumour recurrence. Through the use of their models, both Bajzer et al. (2008) and Titze et al. (2017) were able to provide insight into the long-term behaviour of virus-tumour interactions.
Building on the work by Titze et al. (2017), we present a reduced system of ODEs that model an oncolytic virus interacting with a tumour cell population. While research has progressed and is still advancing, oncolytic virotherapy is at an early stage in its development and this work is primarily motivated by the study of model dynamics. We conduct a local stability analysis and bifurcation analysis of our system and find that stable equilibria only exist in the absence of tumour growth or viral decay. When considering a growing tumour with decaying viral treatment, we see the emergence of long period orbits in the model solution. Numerical evaluation of how the period and amplitude of the orbits depend on the parameter space is used to suggest improvements to oncolytic virotherapy. As we will see, these types of oscillations can biologically be interpreted as recursion or remission of the tumour. In this article, the focus is on which aspects of virus-tumour interactions drive the success of oncolytic virotherapy from a mathematical and biological point of view. Possible outcomes of oncolytic virotherapy predicted by our model are finally analyzed and discussed.

Model description
To model the interaction between an oncolytic virus and a growing tumour, we use a system of three ODEs. While ODE models do not address spatial spread, they do provide a mathematical framework within which the mean-field interactions between tumour cells and viral particles can be explored. The state variables in our model consist of where τ represents the number of days.
In this article we model an aggressive form of tumour, assuming that uninfected tumour cells replicate at a rate r proportional to their population. Biologically, this unbounded exponential tumour growth is infeasible due to nutrient and space limitations. However, exponential growth can still be seen as a sufficient approximation for tumour growth under treatment with an oncolytic virus, since the time scale of the interaction between virus particles and tumour cells is short.
To model the process of tumour cell infection, we use the law of mass action. Assuming that the rate of infection of the uninfected tumour cell population is proportional to the product of the virus and tumour cell populations  and that this occurs with rate constant β. While there are other more complex ways of modelling the viral infection of cells, e.g. a frequency-dependent infection rate or ratio-dependent model, see Hews, Eikenberry, Nagy, and Kuang (2010), Jenner et al. (in press), Novozhilov et al. (2006), we adopt a basic framework that will allow for analytical results, as it will be shown shortly. Further, it is assumed that, once infected, tumour cells are incapable of division as the virus particle within the cell takes control of the cellular machinery for self-replication. Virus-infected tumour cells will then burst due to lysis at a rate d I , releasing α new free virus particles. Figure 1 shows a schematic of the interaction between the uninfected tumour cell population u (or U in non-dimensional form), infected tumour cells i (or I in nondimensional form) and the virus population v (or V in non-dimensional form). The corresponding system of equations describing the interaction is given below: This model complements other oncolytic virotherapy models in the literature, see Bajzer et al. (2008), Berezovskaya et al. (2007), Dingli et al. (2006, 2009in press), Karev et al. (2006), Komarova and Wodarz (2010), Novozhilov et al. (2006), and is very close to the model by Titze et al. (2017). It differs to that of Titze et al. (2017), as tumour cell death due to factors unrelated to treatment are neglected. Unrelated tumour cell death is considered to be negligible in comparison to virus-induced tumour cell death.
Our equations also resemble some previous modelling work by Baccam, Beauchemin, Macken, Hayden, and Perelson (2006) on the kinetics of influenza in humans, i.e. the so called TIV model. Baccam, Beauchemin, Macken, Hayden, and Perelson (2006) derived a model for target-cell limited influenza infection, which is equivalent to Equations (1)-(3) when r = 0. Being our model essentially minimal and quite generic some of its results can be easily translated to influenza and infectious disease modelling.
With an appropriate change of variables, detailed in Appendix 1, the model in Equations (1)-(3) can be simplified to a dimensionless form: where ξ = r/d I and m = d V /d I and U, I, V and t represent the new dimensionless tumour and virus populations and time. This reduces the number of independent parameters and simplifies the mathematical analysis while preserving the essential properties of the model. The model dynamics are sensitive to initial conditions and there is a need to consider realistic initial tumour sizes and viral dosages. As we have assumed exponential growth in our model, if the initial tumour to virus ratio, i.e. U 0 /V 0 , is too high then the tumour growth will outcompete the virus and we will have unbounded growth in the tumour population. In this article, we present numerical simulations of Equations (4)-(6) for initial conditions close to the equilibria, and discuss the dynamical behaviour of the model.
Let us remark that this model pertains to an idealized situation of homogeneous tumour properties and virus spread. It is well documented that oncolytic virotherapy can fail due to intratumoural obstructions, such as the extracellular matrix, pressure and impermeable veins (Ariffin, Forde, Jahangeer, Soden, & Hinchion, 2014). Here, obstacles that may inhibit treatment efficacy are simply ignored.

Results
While parameter estimates for tumour cell replication, viral decay and viral infectivity are readily available in the literature (see Jenner et al., in press;Komarova & Wodarz, 2010;Titze et al., 2017), they represent only one adaptation of the tumour-virus interaction. In this section, we undertake a detailed local stability analysis and bifurcation analysis to quantify how our system behaves under various tumour and viral characteristics.

Equilibrium solutions
To understand long term behaviour, we calculate the equilibria for the non-dimensionalized system, Equations (4)-(6). We find one equilibrium at the origin and one non-zero equilibrium: For the specific case of m = 0 or ξ = 0, two more equilibria exist. For m to be equal to zero, we need d V = 0, i.e. viral particles are not decaying. This biologically represents the case when we have a virus that is not cleared by the immune system, for example because it is coated in an immunogenic polymer such as polyethylene glycol as in Kim et al. (2011). The resulting equilibrium for m = 0 will be U = 0, I = 0 and V ∈ IR.
Similarly, when ξ = 0, we have r = 0, i.e. tumour cells are not replicating. This can be thought of biologically as a stagnant or non-growing tumour. The oncolytic viruses will, therefore, only be removing existing tumour cells. The resulting equilibrium at ξ = 0 will be at I = 0, V = 0 and U ∈ IR. Therefore, we have four equilibria in total, two of which only exist for the specific cases m = 0 or ξ = 0.

Stability of the equilibrium at the origin: U
To achieve complete tumour eradication in our model, the equilibrium at the origin must be stable. Evaluating the Jacobian of the non-dimensionalized model, Equations (4)-(6), for the equilibrium at the origin we obtain the eigenvalues: Thus, the equilibrium is a stable node for m ≥ 0 and ξ ≤ 0 and a saddle point for all other regions in the parameter space as summarized in Figure 2.
Biologically, both ξ and m need to be non-negative real numbers. As such, the reasonable parameter values resulting in a stable node at the origin will be ξ = 0 and m ≥ 0. When ξ = 0, there is no tumour growth, i.e. r = 0, and we have a benign tumour or a malignant tumour growing at a negligible rate. Figure 3 shows a numerical simulation of Equations (4)-(6) for typical parameter values giving a stable node at the origin: all tumour cells and virus particles die out over time.

Stability of the non-zero equilibrium: U = m, I = mξ and V = ξ
As the equilibrium at the origin is unstable for ξ > 0 and m > 0, we now turn to the non-zero equilibrium to see whether our model has a stable fixed point for this parameter set. Evaluating the Jacobian for the non-dimensionalized model, Equations (4)- (6), at the non-zero equilibrium U = m, I = mξ and V = ξ , gives the characteristic equation: The eigenvalues corresponding to the non-zero equilibrium are the roots of the characteristic equation, Equation (9). To determine the stability of the non-zero equilibrium we calculate the position and nature of the stationary points of the characteristic equation and from this deduce the sign and number of real roots of Equation (9). See the diagram in Appendix 2 for a more detailed explanation.
Stationary points of the characteristic equation, Equation (9), occur for two values of λ: The first stationary point listed, λ * 1 , is fixed on the vertical axis λ = 0. The corresponding value of the characteristic equation at the stationary point λ * 1 is ρ(λ * 1 ) = −mξ . The second derivative at the stationary point λ * 1 is ρ (λ * 1 ) = −2(1 + m) and therefore λ * 1 will be a maximum for m > −1 and a minimum for m < −1. This is summarized in Figure 4  To determine whether there exists a (ξ , m)-combination resulting in stable solutions, i.e. all roots of the characteristic equation having negative real part, we use the Routh-Hurwitz criterion. The Routh-Hurwitz stability criterion is a necessary and sufficient condition for the stability of a linear time invariant control system. For the Routh-Hurwitz criterion to be satisfied, we need mξ < 0 and mξ > 0, which gives a contradiction. As such, there is no set of m and ξ that will result in all roots of the characteristic equation with negative real part and, therefore, the non-zero equilibrium will always be unstable.
The sign of the eigenvalues for the non-zero equilibrium, and hence the nature of the non-zero equilibrium, are determined by the position of the two stationary points λ * 1 and λ * 2 in the (λ, ρ(λ))-plane, refer to Appendix 2 for a summary of the eigenvalues based on the position and nature of the two turning points. The nature of the equilibrium for each region of the (ξ , m)-parameter space is plotted in Figure 5. We have three possible values of the non-zero equilibrium: an unstable focus node, a saddle focus and a saddle. For biologically reasonable parameters, ξ > 0 and m > 0, we have a saddle focus, which consists of one negative real eigenvalue and a pair of complex eigenvalues with positive real part.
To illustrate the behaviour of the saddle focus, the numerical solution to the model, Equations (4)- (6), is plotted in Figure 6 for initial conditions close to the non-zero equilibrium. From this, we see that for biologically reasonable parameters we have growing oscillations in all of the variables for the first 30 days of the interaction between the oncolytic virus and the tumour cells.

One-parameter bifurcation analysis
To determine how the solutions change, we conduct a bifurcation analysis of the nondimensionalized model. In Figure 7, the branches of equilibria are plotted for a typical value of m = 0.1. Given that ξ > 0, both equilibria are unstable, as previously illustrated  in the local stability analysis. As m approaches zero from above, the non-zero equilibrium value for U and I decreases and the value for V remains constant, while the stability of the equilibria stays the same. In Figure 7(a), at ξ = 0 we have a zero eigenvalue on the branch of equilibria labelled B2 (i.e. the axis U = 0). As such, between the two branch points (BP1 and BP2) and below the branch point at the origin (i.e. BP1) we have two eigenvalues with negative real part and one zero eigenvalue. Above BP2 we have one eigenvalue with negative real part, one eigenvalue with positive real part and one zero eigenvalue. Therefore branch B2 has a two dimensional manifold that is stable below BP2 and a one dimensional stable and one-dimensional stable manifold above BP2. To illustrate the behaviour of this branch, in Figure 8, we present a numerical simulation for Equations (4)-(6) when ξ = 0 and m = 0.1. Recall that this case represents a non-growing tumour so treatment will only  Figures 2 and 3, we know that for ξ = 0 the equilibrium at the origin is also stable. This means, in general, that we have a set of initial conditions that will tend to the origin and a set of initial conditions that will tend to a non-zero fixed point for ξ = 0 and m ≥ 0 (the case when the tumour cells are not replicating). In Figure 9, we have plotted the solution curves for the model for a fixed m = 0.1 and ξ = 0 and a range of initial conditions. We can see that for a subset of initial conditions in the parameter space, the resulting stable equilibrium will be non-zero for the uninfected tumour cells. This occurs for small vales of initially infected cells, I, and virus particles, V , and is a subset of the (larger) basin of attraction of the portion of branch B2 between BP1 and BP2. An interesting case occurs when m = 0, representing a decay rate of the virus of zero, i.e. the virus is not cleared from the tumour site. In Figure 10, we show the numerical model solution for m = 0 and ξ ≥ 0 and see that uninfected and infected tumour populations are quickly eradicated whilst the virus population tends to a non-zero fixed point. This corresponds to the branch of stable equilibria at U = I = 0 and V ∈ IR (see Section 3.1). While the virus population is still non zero in this situation, the tumour populations have been eradicated and, therefore, we have a positive outcome for oncolytic virotherapy. These results suggest that if we can eliminate viral decay by introducing an immunosuppressant material (such as polyethelene glycol), effective treatment, irrespective of the tumour growth rate, could be achieved.

Incomplete eradication and long period orbits
The goal of our analysis so far has been to determine a parameter regime for our model that will result in complete tumour eradication. As we have shown in the previous section, the equilibrium at the origin is unstable for m > 0, and ξ > 0 and tumour eradication cannot be achieved. Numerically simulating Equations (4)-(6) for a long time period shows the existence of stable long period orbits, see Figure 11. These orbits could be indicating that in the limit for m → 0 + , the system shows a homoclinic state. This investigation will be part of a forthcoming more theoretical communication. The meaning of these dynamics will be discussed in more detail below.
Simulating the non-dimensionalized system presented in Equations (4)-(6) for three biologically reasonable initial conditions, we produce time-series and phase portraits for a range of parameter values. Choosing m = 0.01 and ξ = 0.1 produces the plots in Figure 11(a) and (b). Similarly, fixing ξ = 0.06, we have plotted U as a function of time and the U vs I phase portrait for m = 0.01 and m = 0.001, see Figure 11(c)-(f). It is clear from Figure 11 that the lower the value of m or ξ , the closer the orbit gets to the long period orbit state.
To quantify the dependence of the orbits on the parameter values, the amplitude and period are numerically calculated as a function of m and ξ . In Figure 12(a), the period between oscillations has been calculated as a function of m (i.e.: d V /d I , the ratio of viral death to cell burst rate) and ξ (i.e.: r/d I , the ratio of tumour cell replication to cell burst rate). Decreasing m and increasing ξ results in a longer period of time between the oscillations. Therefore, a slower growing tumour relative to cell burst rate will produce longer intervals of no growth between its rapid burst-like growths. Equivalently, a more rapid clearance of the viral particles relative to cell burst rate would result in longer periods between the oscillations. This can happen when the immune system's response to a virus is heightened, for example through the expression of immunostimulatory cytokines. Figure 12(b)-(d) shows how the amplitude of the oscillation depends on m and ξ . We see an inverse relationship between the amplitude of the oscillation for uninfected tumour cells, U, and the amplitude of oscillation for infected tumour cells, I, and virus particles, V . Increasing ξ results in a lower amplitude for U, and a larger amplitude for I and V . Improving treatment correlates to obtaining the lowest possible tumour population and therefore the lower the amplitude of the uninfected tumour cells close to the long period orbit state, the more effective the treatment.

Discussion
In this article, we presented a mathematical model that captured some of the key dynamics of the interaction between an oncolytic virus and a population of tumour cells. To develop our model we drew on previous modelling in the literature, primarily the work of Titze et al. (2017). Titze et al. (2017) calibrated experimental viral titre and tumour measurements to obtain parameter estimates for a system of ODEs similar to the model we presented in Equations (1)-(3).
Before considering the specific effects of individual tumours and viral characteristics, an understanding of the main determinants of the general dynamic interaction between viruses and tumours is crucial. This was embodied in the parameters m and ξ representing the ratio or tumour replication and viral decay to cell burst rate. To discuss the possible long-term behaviour of the interaction between an oncolytic virus and tumour for a range of possible characteristics, we conducted a local stability analysis and bifurcation analysis of the reduced model, see Equations (4)-(6). We found that the equilibrium at the origin is unstable for all biologically reasonable parameter values, except the case where ξ = 0 and m = 0. When ξ = 0, there will be no tumour growth in the model: this can only occur for a specific type of tumour, i.e. a benign tumour or a tumour whose growth is extremely low. For this tumour type, the model predicts that complete tumour eradication can be obtained for a specific range of initial conditions. In Figure 9, we have plotted model solution curves as a function of the initial conditions, noting by the colour of the line the uninfected tumour cell number at the stable equilibrium. High enough initial viral dosages and initial tumour sizes result in complete tumour eradication: benign or slow growing tumours would do well under this treatment, given the right initial tumour sizes and viral dosages.
Considering the stability of the non-zero equilibrium, we found that this is stable in the absence of viral decay, i.e. for m = 0. In Figure 10 we have plotted the corresponding model simulation for m = 0 and ξ ≥ 0 and show how the system tends to equilibrium U = I = 0 and V ∈ R. Whilst developing a virus that rigorously does not decay is impossible, experimentalists are currently investigating ways to shield viral particles from immune detection and clearance (an example of this is the polymer polyethylene glycol discussed in Kim et al. (2011)). However, at the moment, we propose this as a purely hypothetical scenario, as the non-decaying virus, while harmless, will still need to be removed after it has eliminated the tumour. Our model predicts that if it were possible to genetically engineer an oncolytic virus to remain within the system indefinitely, complete tumour eradication can be obtained.
Analysis of equilibria for the case of an exponentially growing tumour and viral particles undergoing decay, i.e. m > 0 and ξ > 0, suggested that there is no way for treatment to eradicate the tumour, as both equilibria are unstable. As shown in Figures 5 and 6, the non-zero equilibrium is a saddle focus and the origin is a saddle, causing model solutions to spiral outwards with increasing amplitude. Further, numerically simulating the nondimensionalized model for a long time period showed the appearance of long period orbits, see Figure 11.
Oscillations in tumour cell population size have been seen in vivo in several other viral dynamic models (see Bajzer et al., 2008;Dingli et al., 2009;Wodarz, 2016). Komarova and Wodarz (2010), show that using mass action to model the viral infectivity leads to strong oscillations in the population of viruses and cancer cells. Titze et al. (2017) suggested that the rise and fall of tumour growth, seen in oscillations, could be due to the lack of bioincubators for viral replication. This is analogous to behaviours typical of predator-prey systems, where we see oscillations occurring due to a heavy dependence of each population on the other for survival.
There are two primary ways to interpret biologically the presence of long period orbits in oncolytic virotherapy: complete tumour eradication or tumour remission. A long period orbit can be considered as an example of complete tumour eradication: if the population of cells drops below certain levels, this could mean extinction. This in a more realistic setting could occur if we take into account increased likelihood of clearance or death due to nutrient deficiency. By definition of a long period orbit, the model solution spends a long period of time close to zero, and in this time frame we predict that other effects eradicate the negligible number of remaining cells.
Previously, Novozhilov et al. (2006) showed that for a model of oncolytic virotherapy there is a region of the parameter space where trajectories form a family of homoclinics to the origin. From the biological point of view, this occurrence implies that tumour cells can be eliminated with time and complete recovery is possible. Similarly, Berezovskaya et al. showed that the origin can have its own basin of attraction in the phase space (Berezovskaya, Karev, and Arditi, 2001), which corresponds to deterministic extinction of both species (Berezovskaya et al., 2007;Berezovsky, Karev, Song, & Castillo-Chavez, 2005;Hwang & Kuang, 2003;Jost & Arditi, 2000). Berezovskaya et al. (2007) showed how certain models possess a dynamical regime of deterministic extinction, through the presence of homoclinics.
Long period orbits can also be examples of tumour remission and recurrence. Analysis of the characteristics of the orbits as a function of the parameter space provides us with information about how to improve treatment effectiveness. Reducing the amplitude of the orbit and increasing the period between oscillations correspond to a more successful treatment. To quantify how the behaviour depends on parameter space, we compute the period of oscillations as a function of both m (the ratio of the decay of viral particles to cell burst rate) and ξ (the ratio of tumour replication to cell burst rate). In Figure 12(a) we see that decreasing the ratio of viral death to cell death, i.e. m, irrespective of the rate of tumour cell replication r and to cell burst rate d I , i.e. ξ , results in a longer period between oscillations, i.e. more time between tumour regrowth. Alternatively, decreasing the ratio of tumour cell replication to cell burst rate, i.e. ξ , the period increases. Therefore, one important insight the model provides is that decrease in both the ratio of the decay of viral particles to cell burst rate and the ratio of tumour replication rate to cell burst rate is a very effective strategy.
Another key feature of the oscillations in Figure 11 is their amplitude. In Figure 12(b)-(d) we examine the dependence of the amplitude of the oscillation on m (the ratio of the decay of viral particles to cell burst rate). We find that in all cases, the lower the value of m the lower the amplitude of the oscillation. However, if we look at ξ (the ratio of tumour replication to cell burst rate), we find that increasing ξ increases the amplitude of the oscillation for infected tumour cells and virus particles and decreases the amplitude for uninfected tumour cells. One of the primary objectives of a therapy is to reduce the number of uninfected tumour cells, and therefore reduce the amplitude of the oscillation for U. In that respect, we suggest a decrease in the ratio of the decay of viral particles to cell burst rate, m, and an increase in the ratio of tumour replication to cell burst rate, ξ . This will result in the lowest possible amplitude for the uninfected tumour cells. We note that this suggestion is also associated with a long period between oscillations, so it represents overall the most effective way of improving oncolytic virotherapy.

Conclusion
This paper has provided a detailed local stability and bifurcation analysis of a coupled ODE model that investigates the dynamics of an oncolytic virus and tumour cells. This model is a minimal embodiment of a complex intricate interaction and contains a number of biological limitations such as the absence of the immune interaction, the ability of the virus to infect already infected tumour cells and spatial heterogeneity. Despite this, we feel that our model still captures important aspects of the overall dynamics of the interaction between cancerous cells and oncolytic viruses.
From our analysis, we have determined that treatment of a benign or slowly growing tumour with oncolytic virotherapy will result in complete tumour eradication when the right viral dosage for a given initial tumour size is used. We show that modification of viral particles to remove the possibility of viral decay would also result in complete tumour eradication. For a growing tumour and decaying viral particles, long period orbits emerge indicating either eradication or remission. From our analysis, decreasing the ratio of the decay rate of viral particles to cell burst rate and increasing the ratio of tumour replication rate to cell burst rate results in low amplitudes for the uninfected tumour cells and long periods between oscillations.
Our results highlight potential difficulties in the treatment of aggressive tumours with virus therapy alone and suggest the need for combination therapies if complete tumour eradication is always to be achieved. Recent advancements in oncolytic virotherapy have looked at genetically modifying a virus to contain other cancer treatments. In this way, viral infection and drug delivery would occur simultaneously, increasing the treatment's antitumour potency. Combining the tumour-killing capability of oncolytic viruses with anti-tumour drugs, such as Cisplatin, Herceptin and chemotherapeutic agents Inoue et al. (2000), Toyoizumi et al. (1999), or with immunostimulatory cytokines such as IL-12, GM-CSF, 4-1BB ligand (Choi, Zhang, Choi, Kim, & Yun, 2012;Huang et al., 2010), allows for a dual attack on cancer. Increasing evidence suggests that, for a complete eradication, oncolytic virotherapy might need to be used in conjunction with other cancer therapies.