Open access peer-reviewed chapter

Mathematical Modeling and Dynamics of Oncolytic Virotherapy

Written By

Abdullah Abu-Rqayiq

Submitted: October 28th, 2020 Reviewed: March 2nd, 2021 Published: April 23rd, 2021

DOI: 10.5772/intechopen.96963

Chapter metrics overview

265 Chapter Downloads

View Full Metrics


Oncolytic virotherapy is a cancer treatment that uses competent replicating viruses to destroy cancer cells. This field progressed from earlier observations of accidental viral infections causing remission in many malignancies to virus drugs targeting and killing cancer cells. In this chapter, we study some basic models of the oncolytic virotherapy and their dynamics. We show how the dynamical system’s theory can capture the behavior of the solutions of those models and provide different approaches to studying the models. We study the thresholds that enable us to classify asymptotic dynamics of the solutions. Fractional-derivative approach tells us about the memory of the derivative and related solutions of the models. We also study the affect of introducing control parameters on the cost of the therapy.


  • Oncolytic Virotherapy
  • Stability
  • Bifurcation
  • Burst Size
  • Optimization
  • Immunotherapy

1. Introduction

Oncolytic viruses are a form of immunotherapy that uses viruses to infect and destroy cancer cells. These viruses can selectively replicate in cancer cells but leave healthy normal cells largely intact. In oncolytic virotherapy, the free viruses infect tumor cells and replicate themselves in tumor cells; upon analysis of infected tumor cells, new virion particles burst out and proceed to infect additional tumor cells. This idea was initially tested in the middle of the last century and merged with renewed ones over the last three decades due to the technological advances in virology and in the use of viruses as vectors for gene transfer. Over the last decade, great efforts have been made for understanding dynamics and molecular mechanics of viral cytotoxicity of oncolytic viruses. Those efforts provided an interesting possible alternative therapeutic approach to help cure cancer patients. However, the outcomes of virotherapy depends in a complex way on interactions between viruses and tumor cells [1]. One of the main advantages of applying the oncolytic virotherapy is that it can selectively damage cancerous tissues leaving normal cells unharmed. In addition, oncolytic viruses can mediate the killing of the normal cells by indirect mechanisms such as the destruction of tumor blood vessels, the amplification of specific anticancer immune responses or through the specific activities of transgene-encoded proteins expressed from engineered viruses.

During the last two decades, several mathematical models have been applied to understanding oncolytic virotherapy. For example, Wu et al. [2] and Wein et al. [3] proposed and analyzed some partial differential equations models to study some aspects of cancer virotherapy. For ordinary differential equations models, Wodarz in [4, 5], Komarova and Wodarz [6], Novozhilov et al. [7], and Bajzer et al. [8], Tian in [9, 10], and others. Wodarz and Komarova [11] have modeled the dynamics of the oncolytic virus replication by ordinary differential equations that describe the development of the average population sizes of cells and viruses over time. For this purpose, they used a generalized approach and considered a class of models instead of a specific model and took into account two populations: uninfected tumor cells, denoted by x and infected tumor cells, denoted y. The general model is based on the law of mass action and is given by


where the function F describes the growth properties of the uninfected tumor cells, and the function G describes the rate at which tumor cells become infected by the virus. The two functions can take several forms depending on the biological content and meaning that we may want to incorporate into the model. The parameter β represents the infectivity of the virus, and the death rate ay represents the virus-infected cells die.

Then, a three populations model was introduced by Wodarz [12] as


in which v stands for the free virus population and C is maximal tumor size. The term αy models the release of virions by infected tumor cells, and γv is the clearance rate of free virus particles by various causes including non-specific binding and generation of defective interfering particles. The death rate of tumor cells dx seems redundant, since it is included in the logistic model.


2. A basic model of oncolytic virotherapy

Tian [10] has proposed a modified model where the burst size was incorporated. The burst size of a virus is the number of new viruses released from a lysis of an infected cell. It us known that different types of viruses have different burst sizes. Viruses of the same type have almost the same burst size. The burst size is an important parameter of virus replicability.


Considering this last model as a starting point of our discussion of mathematical models of oncolytic virotherapy, we first show some analytical results. In this model, there are two threshold values for the burst size. When the burst size is smaller than the first threshold value, virotherapy always fails. When the burst size is in the between of the two threshold values, we have a partial success of virotherapy represented by the stable positive equilibrium solution. Since the tumor load is a decreasing function of the burst size, the minimum tumor load can be reached by genetically increasing the burst size of the virus up to the second threshold value. If the set in which the positive equilibrium solution is stable has more than one open intervals, we can increase the burst size up to the supreme value of this set, and still have stable partial therapeutic success with even lower tumor load. Once the burst size is greater than the second threshold value, there are one or three families of stable periodic solutions to the system of virotherapy dynamics.

For simplicity, the above system can be non-dimensionalized by setting τ=δt, x=Kx̂, y=Kŷ, v=Kv̂, and rename parameters r=λ,a=βK, and c=γ.

Then dropping all hats over the variables and write τ as t, we have


It is assumed that all parameters are nonnegative.

Model (4) has three equilibrium points, E0=0,0,0, E1=1,0,0, and the positive equilibrium E+=xyv. The equilibrium E0 is always unstable for all positive values of the burst size b. The equilibrium E1 is globally asymptotically stable when 0<b<μ1, and it is unstable when bμ1.

At b=μ1, the positive equilibrium E+ moves into the domain D=xyv:x0y0v00x+y1, a type of transcritical bifurcation occurs with E1 and E+. As the parameter b increases, while μ1<b<μ2 and bIp, E+ is locally asymptotically stable. When b>μ1 and bIn, E+ is unstable. Hopf bifurcations occur for some bμ2, and these bifurcations give rise to one or three families of periodic solutions. Here, μ1 and μ2 are thresholds, Ip=b>μ1:Hb>0, In=b>μ1:Hb<0, and Hb is defined next in formula (10).

  • E0 is unstable.

  • E1 is globally asymptotically stable when b<1+ca ; and unstable when b1+ca .

  • When μ1<b<μ2, the equilibrium solution E+ is locally asymptotically stable

Two types of bifurcations occur in the system as the parameter b varies. A transcritical bifurcation at b=μ1 introduces the equilibrium point E+ into the positive invariant domain D. The Hopf bifurcation at some value b>μ1 gives rise to the periodic solutions. The system (4) is a basic model of the oncplytic virotherapy. Three equilibrium points can be found: the trivial equilibrium E0=0,0,0, E11,0,0, and the positive equilibrium E+xyv, where x=cab1, y=rcabacab1aba+rc, and v=rabacaaba+rc.

The Jacobian matrix is given by


Proceeding in the linearization process and analyzing the eigenvalues, we find that the Jacobian matrix evaluated at E0 is given as


The obtained eigenvalues of this lower triangular matrix are λ1=r, λ2=1, and λ3=c. Note that λ1 is positive due the r being positive. The other eigenvalues are negative. Therefore, E0 is unstable. The local unstable invariant manifold lives in the xaxis. The stable invariant manifolds live in the yvplane. This observation might be interpreted by saying that in the absence of viruses and the infected tumor cells, the tumor cells will grow away from E0.

The Jacobian matrix evaluated at E1 is given by


The eigenvalues are λ1=r, which is clearly negative, and

λ2,3=121+a+c±1ac2+4ab. λ2, the eigenvalue with negative square root, is negative as well. However, for λ3, it is not so clear when it is negative. Computations in [7], show that for 1<b<1+ca, E1 is a global asymptotically stable. This result can be obtained by applying appropriate Lyapunov functions such as V1xyv=y+v and V2xyv=12aba+cy2+a2byv+12a2v2. When the value of b exceeds the threshold μ1=1+ca, E1 becomes unstable. The system exhibits a transcritical bifurcation with bifurcation value μ1 and changes stability as this parameter varies near the bifurcation value. When b>μ1, the positive equilibrium point E+ appears. Proceeding with the linearization, the Jacobian materix at E+ is given by


To analyze the more complicated eigenvalue expressions, we apply the so called Routh-Hurwitz criterion on the associated characteristic polynomial Pλ=λ3+a1λ2+a2λ+a3, where a1=rc+aba+abcab1, a2=rcbc+b1ab12+rcabacraab1aba+rc, and a3=rcabacab1.

The R-H criterion requires


Computations based on the R-H stability criterion can lead to the following expression φb=ab1aba+rcaba+rc+abcbc+b1aba+rcb1abac, [10]. the positive equilibrium E+ is locally asymptotically stable if ϕb<ra.

In addition to the transcritical bifurcation described above, the system develops a Hopf bifurcation for values of b>μ1 resulting in periodic solutions due to the pure imaginary eigenvalues ±ia2.

In [10], bifurcation analysis was accomplished by the aid of defining a function




Then, the another threshold, named μ2 can be defined as the smallest number b of the set I0=b>μ1:Hb=0.

The analysis presented in [10] also shows that if c1+raba+rc<1, then E+ represents a partial success of virotherapy at a modest value of b. The expression c1+raba+rc is called the tumor load.

When the value of b satisfies μ1<b<μ2, E+ is locally asymptotically stable. For values of b<μ1 the positive equilibrium E+ does not live in the positive invarient domain D, and as b increases to μ1, the equilibrium point moves into the domain D and it coalecses with E1. Finally, when b>μ2, periodic solutions will appear as a result of the Hope bifurcation (Table 1).

λTumor growth rate2×1021/h
δDeath rate of infected tumor cells1/181/h
βInfection rate of the virus7/10×109mm3h/ virusl
kImmune killing rate of virus108mm3h/ immune cell
bBurst size of free virus50viruses/cell
γClearance rate of virus2.5×1021/h

Table 1.

Pameters’ Description.


3. Model with innate immune response

Phan and Tian [13] developed the basic model by incorporating innate immune. The system is given by


where λ is tumor growth rate, C is the carrying capacity of tumor cell population, β is the infected rate of the virus, μ is immune killing rate of infected tumor cells, δ is death rate of infected tumor cells, b is the burst size of oncolytic viruses (i.e., the number of new viruses coming out from a lysis of an infected cell), k is immune killing rate of viruses, γ is clearance rate of viruses, s is the stimulation rate of the innate immune system, and ρ is immune clearance rate.

We non-dimensionalize the system by setting τ=δt,x= Cx¯,y=Cy¯,v=Cv¯,z=Cz¯ and rename parameters r=λ/δ a=/δ,c=μC/δ,d=kC/δ,e=γ/δ,m=sC/δ, and n=ρ/δ. Then system (3.1) becomes


All parameters are assumed to be nonnegative. The effects of the innate immune system on the virotherapy in the model are encoded in the parameters c,d, and m. To understand how the innate immune system affects the dynamics of oncolytic virotherapy, they use three combined parameters, the viral burst size b, the relative immune killing rate K=c/d, and the relative immune response decay rate N=n/m, to describe the solution behaviors of the model. Note that K represents the ratio of the rate of immune killing infected tumor cells over the rate of immune killing viruses, which can be considered as a relative immune killing rate of viral therapy since it describes the ability of the innate immune system destroying infection versus destroying viruses.

The system (12) have the following equilibrium points;





E4=1Nv2qNv2b1Nev2cN+dv2, and


Where N=n/m, A=a2+a223a1/3, a1=rNa2cda+aNe+dc, and a2=cdN+raN1.

The analysis of model (12) is similar but more complicated than the basic model proposed by Tian. For details, see [13].

Now we describe some results of Phan and Tian work, [13]. It is clear that if no tumor cells existed, then there exists the equilibrium point E0. If the effect of the immune system and the viruses is neglected, then the system has another equilibrium point E1 with only tumor cells. This happens when the viruses are not powerful due to the burst size being smaller than some critical value. The equilibrium E2 exists if the burst size is greater than the critical value (threshold), meaning that the viruses are powerful. With some conditions on the parameter K and with another burst size critical value, two newly formed equilibrium points will be born for values of b exceed the critical point.

Analysis and numerical simulation can show the existence of two types of bifurcations around the threshold values; the transcritical bifurcation which occurs with the equilibrium points E1 and E2, and a Hopf bifurcation that occurs for larger values of the burst size b.

Due to the complexity of expressions and knowing that it is impossible to find closed forms of bifurcation parameters, especially for the positive equilibrium points E3, E4, and E5, numerical simulations become a need in order to capture the different behavior over the dynamical landscape.

Below the first threshold value of the burst size, the tumor always grows to its maximum size (carrying capacity), then as the bfurcation parameter b passes the first threshold value, the first locally stable positive equilibrium is born through the transcritical bifurcation. When the parameter value is at or exceeds the second threshold, families of periodic solutions arise from the Hopf bifurcation leading to undetectible level of tumor load.


4. Fractional derivative approach

Fractional derivative is a generalization of the usual derivative to include all orders of derivations. It can be traced back to the times of the invention of the calculus itself. The question about the 1/2 derivative was first asked by L’hopital as a reply to Leibniz letter which introduces the notation of the nth derivative. Most of biological systems have long-range temporal memory. Modeling such systems by fractional models provides the systems with a long-time memory and extra degrees of freedom. Despite of the fact that differential equations with integer-orders have long been used in modeling cancer, the fractional-order differential equations (FODEs) have been recently used to model many biological phenomena. One of the advantages of using FODEs to model such phenomena is that models become more consistent with the biological model. This is due to the fact that fractional order derivatives can capture the memory and hereditary properties of those models [14]. The classical mathematical models with integer-orders ignore the intermediate cellular interactions and memory effects. For example, the kinetic of the viral decline in patients responding to interferon is characterized by bi-phase shape following a delay about 89 hours, likely to be the sum of interferono-pharmacokinetics and pharmaco-dynamics as well as the intracellular delay of the ciral life cycle [15]. Therefore, modeling of the biological systems by fractional order differential equations has more advantages than classical integer-order mathematical modeling, in which such effects are neglected. Abu-Rqayiq and Zannon [16] have formulated Tian’s model using Caputo derivative definition. The system can be formulated as


where Dtα is the Caputo fractional derivative and 0<α1. We assume that all parameters are nonnegative.

The fractional order integration and fractional order can be defined as: the definition of fractional order integration and fractional order. Let L1=L1ab be the class of Lebesgue integrable functions on ab, a<b<. The fractional integral of order νR+ of the function ft,t>0 f:RR is defined by


where Γ. is the Gamma function.

The fractional derivative of order αn1n of the function ft is defined by several ways, the most common ones are:

  1. Riemann-Liouville fractional derivative: Take the fractional integral of order nα and then apply the nth derivative


    where Dn=dndtn,n=1,2,;

  2. Caputo’s fractional derivative: Start with a nth derivative of the function, then take a fractional integral of order nα


Since fractional-order models possess memory, FODE gives us a more realistic way to model oncolytic virotherapy and study their dynamics. The presence of a fractional derivative in a differential equation can lead to an increase in the complexity of the observed behavior. On the other hand, it can show how the solution is continuously dependent on all the previous states. The numerical results of applying the fractional approach will be show in Figures 14.

Figure 1.

Dynamics of virotherapy when b=4 and initial values x=0.5, y=0.5, and ν=1.5, for α=1, α=.8, and α=.9.

Figure 2.

Damped oscillators when b=27 and initial values x=0.5, y=0.5, and ν=1.5, for α=.96, α=.98, and α=1.

Figure 3.

Dynamics of model when b=28 and initial values x=0.5, y=0.5, and ν=1.5, for α=.96, α=.98, and α=1.

Figure 4.

Dynamics of Tumor cells vs infected tumor cells of the model when α=.98.


5. Optimization by control theory

In this section, we develop a model for the controlled infected brain tumor cells. optimal control theory is applied to the cost functional and is supposed to achieve the ultimate goal of optimizing that functional and find a best strategy for minimizing the cost of the virotherapy. The goal here is to model, analyze, and explore optimal ways that can minimize a tumor and the cost of the virotherapy.

Optimal control theory is a branch of the applied mathematics that deals with finding the best possible control that can take a dynamical system from one state to another. The Hamiltonian of optimal control theory was developed by the Russian mathematician Lev Pontryagin as a part of his investigation into the maximum principle. Pontryagin proved that the necessary condition for solving certain optimal control problems is that the control should be chosen in such a way that minimizes the Hamiltonian, [17].

The general form of the control function ut is given by


where xt is the system state which evolves according to the state equation


The Hamiltonian is defined as


where Ψt is a vector of costate variables of the same dimension as the state variable xt such that, [18]


Applying the control theory approach, we reformulate the basic model by introducing a control function ut which represents efforts on damaging the tumor cells AND 1ut represents the growth rate of the infected cells. After incorporating the control u into the basic model, we obtain the following model with control


The control is usually assumed to be bounded by maximim value less than 1 and greater than 0. For our current model, we assume the maximum value is 0.9, a choice that make our proposed model more realistic from a medical view point.

The objective function will be the function that will host our optimal value u and it is given by


Where B is a measure of the relative cost of interventions associated to the control ut. Our goal is to minimize the number of the infected tumor cells by choosing an appropriate strategy that can lower the number of free viruses as well. As a result of that, the cost of treatment will be lowered.

The admissible set of control functions is defined as


Rescaling the system to ease mathematical treatment, we use same parameter rescaling that was used for the previous models and we get


with fixed initial conditions x0,y0,v0 and a fixed final time T.

According to The Pontryagin’s Maximum Principle, if u.Ω is optimal for the problem under consideration, the minimizer with the initial conditions and fixed final time T, then there exists a nontrivial absolutely continuous mapping Ψ:01R3.

Now we come to the main result in this section, [19]

System (16) along with the initial conditions and the final time T has a unique optimal solution xyv associated to an optimal control u on 0T. Moreover, there exists adjoint functions Ψ1,Ψ1, and Ψ3, such that with transversality conditions ΨiT=0,i=1,2,3. Furthermore,


The proof is given in [19].


6. Numerical simulation and discussion

For the simulations and numerical results of the basic model and those of the fractional approach, we use the same parameter values used in [10] and summarized in Table 1. We also combing those results in the Figures 14 below. Note that α=1 represents the simulation of the basic model (2.2), whereas the other values of α describe the memory of the derivatives of the basic model. The parameter values are r=0.36, a=0.11, and c=0.2. By considering b=9, the following equilibrium points can be obtained E0=0,0,0,0, E1=1,0,0, E2=0.6,0.0730,2.5729. Here the bifurcation parameter values are μ1=5 and μ2=27.766. When 5<b<27.766, E+ is locally asymptotically stable while E1 is unstable. The equilibrium point E0 is always unstable. Figure 1 shows the treatment will eventually reach the equilibrium point E1 that is locally asymptotically stable. Figures 2 and 3 show periodic solutions rising from Hopf bifurcation, and Figure 4 shows the dynamics of tumor cells vs infected tumor cells when α=0.98 which is almost the same as the result of the usual derivative α=1. E+ is locally asymptotically stable when 5<b<27.766.

The numerical results of the optimal control system (5.2) can be obtained by implementing forward fourth-order Runge-Kutta method for state system and the backward one for the adjoint system. The method depends on the choice of an initial guess for the value of the control u. The optimal control system is estimated to predict the evolution of the tumors cells relative to specific choices of virus bust size. Simulation shows the results in time scale of 100 days for burst size b=4, 200 days for burst size b=9, and 1000 days for the oscillation to capture that behavior as shown in the Figures 58.

Figure 5.

Optimal state variables for the controlled and the uncontrolled systems subject to the initial values x=0.5,y=0.5, and v=1.5, b=4, and the admissible control set versus trajectories without control measures.

Figure 6.

Optimal state variables for the controlled and the uncontrolled systems subject to the initial values x=0.5,y=0.5, and v=1.5, b=9, and the admissible control set versus trajectories without control measures.

Figure 7.

Damped oscillators appear for the controlled and the uncontrolled systems subject to the initial values x=0.5,y=0.5, and v=1.5, b=26, and the admissible control set versus trajectories without control measures.

Figure 8.

The optimal control u for the Oncolytic virotherapy model subject to the initial values x=0.5,y=0.5, and v=1.5, b=9, and the admissible control.

Numerical results show that the existence of the control can improve the growth of the normal cells until approximately 60 days of the therapy and will be stabilized after then. Whereas the number of the infected cells will be dropped significantly after the fifth day of the treatment until they are completely terminated in the day 50. The dynamics is hugely determined by the burst size in addition to the other control parameter values. The numerical results clearly show that the virotherapy can reduce the tumor load within days of the therapy and reduces number of the free viruses that are needed in the therapy. As a result, the cost of the therapy is minimized. See Figures 58.


  1. 1. E. A., Chiocca, Oncolytic Viruses, Nature Reviews Cancer, 2, (2002), 938–950
  2. 2. J. T. Wu, H. M. Byrne, D. H. Kirn, and L. M. Wein, Modeling and analysis of a virus that replicates selectively in tumor cells, Bulletin of Mathematical Biology, 63, no. 4,(2001), 731–768
  3. 3. L.M. Wein, J.T.Wu, and D.H.Kirn, Validation and analysis of a mathematical model of a replication-competent oncolytic virus for cancer treatment: Implications for virus design and delivery, Cancer Research, 63, no. 6, (2003), 1317–1324
  4. 4. D. Wodarz, Computational approaches to study oncolytic virutherapy: insights and challenges, Gene Ther Mol Biol, 8(2004), 137–146
  5. 5. D. Wodarz, Viruses as antitumor weapons: defining conditions for tumor remission, Cancer Research, 61, no 8, 2001, 3501–3507
  6. 6. L. Komarova and D. Wodarz, ODE models for oncolytic virus dynamics, Journal of Theoretical Biology, 263(2010), 530–543
  7. 7. A. S. Novozhilov, F., S. Berezovskaya, E. V. Koonin, and G.P. Karev, Mathematical modeling of tumor therapy with oncolytic viruses: regimes with complete tumor elimination within the framework of deterministic models, Biology Direct, 1, article 6,(2006), 1–18
  8. 8. Z. Bajzer, T. Carr, K. Josic, S. J. Russell, and D. Dingli, Modeling of cancer virotherapy with recombinant measles viruses, Journal of Theoretical Biology, 252 (2008), 109–122
  9. 9. A. Friedman, P. Tian, G. Fulci, A. Chiocca, and J. Wang, Glioma Virotherapy: Effect of Innate Immune Suppression and Increased Viral Replication Capacity, Cancer research, vol. 66, no. 4, (2006), 2314–2319
  10. 10. P. Tian, The Replicability of Oncolytic Virus: Defining Conditions in Tumor Virotherapy, Mathematical Biosciences and Engineering, vol 8, n0. 3, (2011), 841–860
  11. 11. D. Wodarz and N. Komarova, Towards predictive computational models of oncolytic virus therapy: Basis for experimental validation and model selection, PloS ONE, 4 (2009), e4271
  12. 12. D. Wodarz, Gene Therapy for Killing p53-Negative Cancer Cells: Use of Replicating Versus Nonreplicating Agents, HUM. Gean Ther,, 14, (2003), 143–159
  13. 13. T. Phan and P. Tian, The Role of the Innate Immune System in Oncolytic Virotherapy, Comp. & Math. Methods in Medicine, 6, (2017), 1–17
  14. 14. I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, (1999)
  15. 15. J. A. Machado, Entropy analysis of integer and fractional dynamical systems, Nonlinear Dynam., 62, (2010), 371–378
  16. 16. A. Abu-Rqayiq and M. Zannon, On Dynamics of Fractional Order Oncolytic Virotherapy Models, J. Math. Computer Sci., 20 (2020), 79–87
  17. 17. W. H. Fleming and R. W. Rishel, Deterministic and Stochastic Optimal Control, Springer Verlag, New York, (1975)
  18. 18. R. Neilan, and S. Lenhart, An Introduction to Optimal Control with an Application in Disease Modeling, DIMACS Series in Discrete Mathematics and Theoretical Computer Science. Vol. 75, (2010)
  19. 19. A. Abu-Rqayiq, H. Alayed, and M. Zannon, Optimal Control of a Basic Oncolytic Virotherapy Model, J. Math. Computer Sci., In Press

Written By

Abdullah Abu-Rqayiq

Submitted: October 28th, 2020 Reviewed: March 2nd, 2021 Published: April 23rd, 2021