Open access peer-reviewed chapter - ONLINE FIRST

# Mathematical Modeling and Dynamics of Oncolytic Virotherapy

By Abdullah Abu-Rqayiq

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

DOI: 10.5772/intechopen.96963

## Abstract

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.

### Keywords

• 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 [4]. 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. [19] and Wein et al. [14] proposed and analyzed some partial differential equations models to study some aspects of cancer virotherapy. For ordinary differential equations models, Wodarz in [15, 16], Komarova and Wodarz [7], Novozhilov et al. [10], and Bajzer et al. [3], Tian in [6, 13], and others. Wodarz and Komarova [18] 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 xand infected tumor cells, denoted y. The general model is based on the law of mass action and is given by

dxdt=xFxyβyGxydydt=βyGxyayE1

Where the function Fdescribes the growth properties of the uninfected tumor cells, denoted as x, and the function Gdescribes 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 ayrepresents the virus-infected cells die.

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

dxdt=rx1x+yCdxβxvdydt=βxvd+aydvdt=αyγv,E2

in which vstands for the free virus population and Cis maximal tumor size. The term αymodels the release of virions by infected tumor cells, and γvis 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 dxseems redundant, since it is included in the logistic model.

## 2. A basic model of oncolytic virotherapy

Tian [13] 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.

dxdt=λx1x+yCβxvdydt=βxvδydvdt=bδyβxvγvE3

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=Kŷ, v=Kv̂, and rename parameters r=λ,a=βK, and c=γ.

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

dxdt=rx1xyaxvdydt=axvydvdt=byaxvcvE4

It is assumed that all parameters are nonnegative.

Model (2.2) has three equilibrium points, E0=0,0,0, E1=1,0,0, and the positive equilibrium E+=xyv. The equilibrium E0is always unstable for all positive values of the burst size b. The equilibrium E1is 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 E1and E+. As the parameter bincreases, while μ1<b<μ2and bIp, E+is locally asymptotically stable. When b>b0and 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, μ1and μ2are thresholds, Ip=b>μ1:Hb>0, In=b>μ1:Hb<0, and Hbis defined next in formula (5).

Theorem 1.1

1. E0is unstable.

2. E1is globally asymptotically stable when b<1+ca; and unstable when b1+ca.

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

Two types of bifurcations occur in the system as the parameter bvaries. A transcritical bifurcation at b=μ1introduces the equilibrium point E+into the positive invariant domain D. The Hopf bifurcation at some value b>μ1gives 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

J=r2rxryavrxaxav1axavbaxcE5

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

JE0=r000100bcE6

The obtained eigenvalues of this lower triangular matrix are λ1=r, λ2=1, and λ3=c. Note that λ1is positive due the rbeing positive. The other eigenvalues are negative. Therefore, E0is 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 E1is given by

JE1=rra01a0bacE7

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 [10], show that for 1<b<1+ca, E1is a global asymptotically stable. This result can be obtained by applying appropriate Lyapunov functions such as V1xyv=y+vand V2xyv=12aba+cy2+a2byv+12a2v2. When the value of bexceeds the threshold μ1=1+ca, E1becomes unstable. The system exhibits a transcritical bifurcation with bifurcation value μ1and 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

JE+=r2rxryavrxaxav1axavbaxc.E8

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

H1=a1>0, H2=a1a31a2>0, and H3=a1a301a200a1a3>0. (9)

Computations based on the R-H stability criterion can lead to the following expression φb=ab1aba+rcaba+rc+abcbc+b1aba+rcb1abac, [13]. 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>μ1resulting in periodic solutions due to the pure imaginary eigenvalues ±ia2.

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

Hb=rcΦb1a2b13aba+rc,E10

where

Φx=aaxcax+rcx2+a+acx+rc+ac[c+1x+cax+rc+xaxcra]=a3x4+a23c+c2+raac+1x3+ac3rc+3a+rc2+3ac+r+r2a2x2+c23ar+2acr+r2c+2a2x+rc3r+a.

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

The analysis presented in [13] 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+rcis called the tumor load.

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

ParameterDescriptionValueDimentions
λ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 [11] developed the basic model by incorporating innate immune. The system is given by

dxdt=λx1x+yCβxvdydt=βxvμyzδydvdt=bδyβxvkvzγvdzdt=syzρzE11

where λis tumor growth rate, Cis 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, bis the burst size of oncolytic viruses (i.e., the number of new viruses coming out from a lysis of an infected cell), kis immune killing rate of viruses, γis clearance rate of viruses, sis 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

dxdt=rx1xyaxvdydt=axvcyzydvdt=byaxvdvzevdzdt=myznzE12

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 Krepresents 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;

E0=0,0,0,0,

E1=1,0,0,0,

E2=eab1reabaeab1aba+rerabaeaaba+re0,

E3=1NaArNAb1NeAcN+dA,

E4=1Nv2qNv2b1Nev2cN+dv2, and

E5=1Nv3qNv3b1Nev3cN+dv3.

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

The analysis of this model is similar but more complicated than the basic model proposed by Tian. For details, see [11].

Now we describe some results of Phan and Tian work, [11]. 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 E1with only tumor cells. This happens when the viruses are not powerful due to the burst size being smaller than some critical value. The equilibrium E2exists if the burst size is greater than the critical value (threshold), meaning that the viruses are powerful. With some conditions on the parameter Kand with another burst size critical value, two newly formed equilibrium points will be born for values of bexceed the critical point.

The dynamics of system (12) is largely similar to that of the basic model, but more complicated. For example, 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 E1and 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 bpasses 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/2derivative 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 [12]. 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 89hours, likely to be the sum of interferono-pharmacokinetics and pharmaco-dynamics as well as the intracellular delay of the ciral life cycle [8]. 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 [1] have formulated Tian’s model using Caputo derivative definition. The system can be formulated as

Dtαx=rαx1xyaαxvDtαy=aαxvyDtα=bαyaαxvcαvE13

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=L1abbe the class of Lebesgue integrable functions on ab, a<b<. The fractional integral of order νR+of the function ft,t>0f:RRis defined by

Iav=1Γvqttsv1fsds,t>0,E14

where Γ.is the Gamma function.

The fractional derivative of order αn1nof the function ftis 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 nthderivative

Dαaft=DαaInαa,

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

• Caputo’s fractional derivative: Start with a nthderivative 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.

## 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, [5].

The general form of the control function utis given by

Jut=ΨxT+01LxtutTtdt

where xtis the system state which evolves according to the state equation

ẋ=fxtuttx0=x0t0T

The Hamiltonian is defined as

HxΨut=ΨTtfxut+Lxut,

where Ψtis a vector of costate variables of the same dimension as the state variable xtsuch that, [9]

Ψ̇t=Hx.

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

dxdt=1utλx1x+yKβxvdydt=βxv1utδydvdt=bδ1utyβxvγvE15

The control is usually assumed to be bounded by maximim value less than 1and 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 uand it is given by

Jut=0Tyt+12Bu2dt.

Where Bis 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

Ω=uL0tf:0utumaxt[0T]

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

dxdt=1urx1xyaxvdydt=axv1uydvdt=b1uyaxvcvE16

with fixed initial conditions x0,y0,v0and 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, [2]

Theorem 1.2 System (16) along with the initial conditions and the final time Thas a unique optimal solution xyvassociated to an optimal control uon 0T.Moreover, there exists adjoint functions Ψ1,Ψ1,and Ψ3, such that with transversality conditions ΨiT=0,i=1,2,3.Furthermore,

ut=1xyΨ1yΨ2+bΨ3B.

The proof is given in [2].

## 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 [13] and summerized in Table 1. We also combing those results in the Figures 14 below. Note that α=1represents 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=5and μ2=27.766. When 5<b<27.766, E+is locally asymptotically stable while E1is unstable. The equilibrium point E0is always unstable. Figure 1 shows the treatment will eventually reach the equilibrium point E1that 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.98which 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 100days for burst size b=4, 200days for burst size b=9, and 1000days for the oscillation to capture that behavior as shown in the Figures 58.

Numerical results show that the existence of the control can improve the growth of the normal cells until approximately 60days 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.

chapter PDF

## More

© 2021 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## How to cite and reference

### Cite this chapter Copy to clipboard

Abdullah Abu-Rqayiq (April 23rd 2021). Mathematical Modeling and Dynamics of Oncolytic Virotherapy [Online First], IntechOpen, DOI: 10.5772/intechopen.96963. Available from: