**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.

$$\begin{aligned} \frac{d\mathbf{x}}{dt} &= \lambda \mathbf{x} \left( 1 - \frac{\mathbf{x} + \mathbf{y}}{C} \right) - \beta \mathbf{x} v\\ \frac{d\mathbf{y}}{dt} &= \beta \mathbf{x} v - \delta \mathbf{y} \\ \frac{dv}{dt} &= b \delta \mathbf{y} - \beta \mathbf{x} v - \gamma v. \end{aligned} \tag{3}$$

*Mathematical Modeling and Dynamics of Oncolytic Virotherapy DOI: http://dx.doi.org/10.5772/intechopen.96963*

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*^*y*, *v* ¼ *K*^*v*, and rename parameters *r* ¼ *λ*, *a* ¼ *βK*, and *c* ¼ *γ:*

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

$$\begin{aligned} \frac{dx}{dt} &= r\mathbf{x}(1 - \mathbf{x} - \mathbf{y}) - a\mathbf{x}v \\\\ \frac{dy}{dt} &= a\mathbf{x}v - \mathbf{y} \\\\ \frac{dv}{dt} &= b\mathbf{y} - a\mathbf{x}v - cv \end{aligned} \tag{4}$$

It is assumed that all parameters are nonnegative.

Model (4) has three equilibrium points, *E*<sup>0</sup> ¼ ð Þ 0, 0, 0 , *E*<sup>1</sup> ¼ ð Þ 1, 0, 0 , and the positive equilibrium *<sup>E</sup>*<sup>þ</sup> <sup>¼</sup> *<sup>x</sup>*<sup>∗</sup> , *<sup>y</sup>* <sup>∗</sup> , *<sup>v</sup>* <sup>∗</sup> ð Þ. The equilibrium *<sup>E</sup>*<sup>0</sup> is always unstable for all positive values of the burst size *b*. The equilibrium *E*<sup>1</sup> is globally asymptotically stable when 0 <*b*<*μ*1, and it is unstable when *b*≥*μ*1.

At *b* ¼ *μ*1, the positive equilibrium *E*<sup>þ</sup> moves into the domain *D* ¼ f g ð Þ *x*, *y*, *v* : *x*≥0, *y*≥0, *v*≥0, 0 ≤*x* þ *y*≤1 , a type of transcritical bifurcation occurs with *E*<sup>1</sup> and *E*þ. As the parameter *b* increases, while *μ*<sup>1</sup> <*b*<*μ*<sup>2</sup> and *b*∈ *Ip*, *E*<sup>þ</sup> is locally asymptotically stable. When *b*>*μ*<sup>1</sup> and *b*∈*In*, *E*<sup>þ</sup> is unstable. Hopf bifurcations occur for some *b*≥ *μ*2, and these bifurcations give rise to one or three families of periodic solutions. Here, *μ*<sup>1</sup> and *μ*<sup>2</sup> are thresholds, *Ip* ¼ *b*>*μ*<sup>1</sup> f g : *H b*ð Þ> 0 , *In* ¼ *b*>*μ*<sup>1</sup> f g : *H b*ð Þ<0 , and *H b*ð Þ is defined next in formula (10).


Two types of bifurcations occur in the system as the parameter *b* varies. A transcritical bifurcation at *b* ¼ *μ*<sup>1</sup> introduces the equilibrium point *E*<sup>þ</sup> into the positive invariant domain *D*. The Hopf bifurcation at some value *b*>*μ*<sup>1</sup> 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

*<sup>E</sup>*<sup>0</sup> <sup>¼</sup> ð Þ 0, 0, 0 , *<sup>E</sup>*1ð Þ 1, 0, 0 , and the positive equilibrium *<sup>E</sup>*<sup>þ</sup> *<sup>x</sup>*<sup>∗</sup> , *<sup>y</sup>* <sup>∗</sup> , *<sup>v</sup>* <sup>∗</sup> ð Þ, where *<sup>x</sup>*<sup>∗</sup> <sup>¼</sup> *c a b*ð Þ �<sup>1</sup> , *<sup>y</sup>* <sup>∗</sup> <sup>¼</sup> *rc ab* ð Þ �*a*�*<sup>c</sup> a b*ð Þ �<sup>1</sup> ð Þ *ab*�*a*þ*rc* , and *<sup>v</sup>* <sup>∗</sup> <sup>¼</sup> *r ab* ð Þ �*a*�*<sup>c</sup> a ab* ð Þ �*a*þ*rc :*

The Jacobian matrix is given by

$$f = \begin{pmatrix} r - 2r\mathbf{x} - r\mathbf{y} - av & -r\mathbf{x} & -a\mathbf{x} \\ & av & -1 & ax \\ & -av & b & -a\mathbf{x} - c \end{pmatrix} \tag{5}$$

Proceeding in the linearization process and analyzing the eigenvalues, we find that the Jacobian matrix evaluated at *E*<sup>0</sup> is given as

$$J(E\_0) = \begin{pmatrix} r & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & -\mathbf{1} & \mathbf{0} \\ \mathbf{0} & b & -c \end{pmatrix} \tag{6}$$

The R-H criterion requires

*DOI: http://dx.doi.org/10.5772/intechopen.96963*

*H*<sup>1</sup> ¼ *a*<sup>1</sup> >0, *H*<sup>2</sup> ¼

expression *<sup>φ</sup>*ð Þ¼ *<sup>b</sup> a b*ð Þ �<sup>1</sup> ð Þ *ab*�*a*þ*rc*

imaginary eigenvalues �*i* ffiffiffiffiffi

where

is locally asymptotically stable if *ϕ*ð Þ *b* <*r* � *a*.

þ*x ax* ð Þ � *c* ð Þ� *r* � *a*

of the set *I*<sup>0</sup> ¼ *b*> *μ*<sup>1</sup> f g : *H b*ð Þ¼ 0 .

result of the Hope bifurcation (**Table 1**).

*<sup>x</sup>*<sup>4</sup> <sup>þ</sup> *<sup>a</sup>*<sup>2</sup> <sup>3</sup>*<sup>c</sup>* <sup>þ</sup> *<sup>c</sup>*

¼ �*a*<sup>3</sup>

þ*r* <sup>2</sup> � *<sup>a</sup>*<sup>2</sup> � *<sup>x</sup>*<sup>2</sup> <sup>þ</sup> *<sup>c</sup>*

the tumor load.

**Table 1.**

**143**

*Pameters' Description.*

*a*<sup>1</sup> *a*<sup>3</sup>

� � � � �

*ab*�*a*þ*rc*þ*abc* � ð Þ *bc*þ*b*�<sup>1</sup> ð Þ *ab*�*a*þ*rc*

>0, and *H*<sup>3</sup> ¼

Computations based on the R-H stability criterion can lead to the following

In addition to the transcritical bifurcation described above, the system develops a Hopf bifurcation for values of *b*>*μ*<sup>1</sup> resulting in periodic solutions due to the pure

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

<sup>Φ</sup>ð Þ¼� *<sup>x</sup> a ax* ð Þ � *<sup>c</sup>* ð Þ *ax* <sup>þ</sup> *rc <sup>x</sup>*<sup>2</sup> <sup>þ</sup> ðð Þ *<sup>a</sup>* <sup>þ</sup> *ac <sup>x</sup>* <sup>þ</sup> *rc* <sup>þ</sup> *ac*Þ½ð Þ ð Þ *<sup>c</sup>* <sup>þ</sup> <sup>1</sup> *<sup>x</sup>* <sup>þ</sup> *<sup>c</sup>* ð Þ *ax* <sup>þ</sup> *rc*

2 *<sup>c</sup>* <sup>þ</sup> <sup>2</sup>*a*<sup>2</sup> � �*<sup>x</sup>* <sup>þ</sup> *rc*<sup>3</sup>

Then, the another threshold, named *μ*<sup>2</sup> can be defined as the smallest number *b*

When the value of *b* satisfies *μ*<sup>1</sup> <*b*<*μ*2, *E*<sup>þ</sup> is locally asymptotically stable. For values of *b*<*μ*<sup>1</sup> the positive equilibrium *E*<sup>þ</sup> 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 *E*1. Finally, when *b*>*μ*2, periodic solutions will appear as a

**Parameter Description Value Dimentions** *<sup>λ</sup>* Tumor growth rate <sup>2</sup> � <sup>10</sup>�<sup>2</sup> <sup>1</sup>*=<sup>h</sup> δ* Death rate of infected tumor cells 1*=*18 1*=h <sup>β</sup>* Infection rate of the virus <sup>7</sup>*=*<sup>10</sup> � <sup>10</sup>�<sup>9</sup> *mm*3*h<sup>=</sup>* virusl *k* Immune killing rate of virus 10�<sup>8</sup> *mm*3*h=* immune cell *b* Burst size of free virus 50 viruses/cell *<sup>γ</sup>* Clearance rate of virus <sup>2</sup>*:*<sup>5</sup> � <sup>10</sup>�<sup>2</sup> <sup>1</sup>*=<sup>h</sup>*

<sup>2</sup> <sup>þ</sup> *<sup>r</sup>* � *<sup>a</sup>* � *ac* <sup>þ</sup> <sup>1</sup> � �*x*<sup>3</sup> <sup>þ</sup> *ac* <sup>3</sup>*rc* <sup>þ</sup> <sup>3</sup>*<sup>a</sup>* <sup>þ</sup> *rc*<sup>2</sup> <sup>þ</sup> <sup>3</sup>*ac* <sup>þ</sup> *<sup>r</sup>* �

*H b*ð Þ¼ *rc*Φð Þ *<sup>b</sup>* � <sup>1</sup> *<sup>a</sup>*<sup>2</sup>ð Þ *<sup>b</sup>* � <sup>1</sup> <sup>3</sup>

<sup>2</sup> <sup>3</sup>*ar* <sup>þ</sup> <sup>2</sup>*acr* <sup>þ</sup> *<sup>r</sup>*

partial success of virotherapy at a modest value of *b*. The expression *<sup>c</sup>*ð Þ <sup>1</sup>þ*<sup>r</sup>*

The analysis presented in [10] also shows that if *<sup>c</sup>*ð Þ <sup>1</sup>þ*<sup>r</sup>*

*a*<sup>1</sup> *a*<sup>3</sup> 0

� � � � � � � � �

>0*:* (9)

*:*

*ab*�*a*þ*rc* is called

� � � � � � � � �

1 *a*<sup>2</sup> 0

0 *a*<sup>1</sup> *a*<sup>3</sup>

ð Þ *ab* � *<sup>a</sup>* <sup>þ</sup> *rc* , (10)

ð Þ *r* þ *a*

*ab*�*a*þ*rc* <sup>&</sup>lt;1, then *<sup>E</sup>*<sup>þ</sup> represents a

ð Þ *<sup>b</sup>*�<sup>1</sup> ð Þ *ab*�*a*�*<sup>c</sup>* , [10]. the positive equilibrium *<sup>E</sup>*<sup>þ</sup>

� � � � �

*Mathematical Modeling and Dynamics of Oncolytic Virotherapy*

*a*2 p .

1 *a*<sup>2</sup>

The obtained eigenvalues of this lower triangular matrix are *λ*<sup>1</sup> ¼ *r*, *λ*<sup>2</sup> ¼ �1, and *λ*<sup>3</sup> ¼ �*c*. Note that *λ*<sup>1</sup> is positive due the *r* being positive. The other eigenvalues are negative. Therefore, *E*<sup>0</sup> is unstable. The local unstable invariant manifold lives in the *x*�axis. The stable invariant manifolds live in the *yv*�plane. 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 *E*0.

The Jacobian matrix evaluated at *E*<sup>1</sup> is given by

$$J(E\_1) = \begin{pmatrix} -r & -r & a \\ 0 & -1 & a \\ 0 & b & -a - c \end{pmatrix} \tag{7}$$

The eigenvalues are *λ*<sup>1</sup> ¼ �*r*, which is clearly negative, and *<sup>λ</sup>*2,3 <sup>¼</sup> <sup>1</sup> <sup>2</sup> �ð Þ� 1 þ *a* þ *c* ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð Þ 1 � *a* � *c* <sup>2</sup> <sup>þ</sup> <sup>4</sup>*ab* � � <sup>q</sup> . *λ*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<sup>&</sup>lt; *<sup>b</sup>*<sup>&</sup>lt; <sup>1</sup> <sup>þ</sup> *<sup>c</sup> <sup>a</sup>*, *E*<sup>1</sup> is a global asymptotically stable. This result can be obtained by applying appropriate Lyapunov functions such as *<sup>V</sup>*1ð Þ¼ *<sup>x</sup>*, *<sup>y</sup>*, *<sup>v</sup> <sup>y</sup>* <sup>þ</sup> *<sup>v</sup>* and *<sup>V</sup>*2ð Þ¼ *<sup>x</sup>*, *<sup>y</sup>*, *<sup>v</sup>* <sup>1</sup> <sup>2</sup> *ab a*ð Þ <sup>þ</sup> *<sup>c</sup> <sup>y</sup>*<sup>2</sup> <sup>þ</sup> *<sup>a</sup>*<sup>2</sup>*byv* <sup>þ</sup> <sup>1</sup> <sup>2</sup> *a*<sup>2</sup>*v*2. When the value of *<sup>b</sup>* exceeds the threshold *<sup>μ</sup>*<sup>1</sup> <sup>¼</sup> <sup>1</sup> <sup>þ</sup> *<sup>c</sup> <sup>a</sup>*, *E*<sup>1</sup> becomes unstable. The system exhibits a transcritical bifurcation with bifurcation value *μ*<sup>1</sup> and changes stability as this parameter varies near the bifurcation value. When *b*>*μ*1, the positive equilibrium point *E*<sup>þ</sup> appears. Proceeding with the linearization, the Jacobian materix at *E*<sup>þ</sup> is given by

$$f(E\_{+}) = \begin{pmatrix} r - 2r\mathbf{x} - r\mathbf{y} - av & -r\mathbf{x} & -a\mathbf{x} \\ av & -\mathbf{1} & a\mathbf{x} \\ & -av & b & -a\mathbf{x} - c \end{pmatrix}. \tag{8}$$

To analyze the more complicated eigenvalue expressions, we apply the so called Routh-Hurwitz criterion on the associated characteristic polynomial *<sup>P</sup>*ð Þ¼ *<sup>λ</sup> <sup>λ</sup>*<sup>3</sup> <sup>þ</sup> *<sup>a</sup>*1*λ*<sup>2</sup> <sup>þ</sup> *<sup>a</sup>*2*<sup>λ</sup>* <sup>þ</sup> *<sup>a</sup>*3, where *<sup>a</sup>*<sup>1</sup> <sup>¼</sup> *rc*þ*ab*�*a*þ*abc a b*ð Þ �<sup>1</sup> , *<sup>a</sup>*<sup>2</sup> <sup>¼</sup> *rc bc* ð Þ <sup>þ</sup>*b*�<sup>1</sup> *a b*ð Þ �<sup>1</sup> <sup>2</sup> <sup>þ</sup> *rc ab* ð Þ �*a*�*<sup>c</sup>* ð Þ *<sup>r</sup>*�*<sup>a</sup> a b*ð Þ �<sup>1</sup> ð Þ *ab*�*a*þ*rc* , and *<sup>a</sup>*<sup>3</sup> <sup>¼</sup> *rc ab* ð Þ �*a*�*<sup>c</sup> a b*ð Þ �<sup>1</sup> .

The R-H criterion requires

$$H\_1 = a\_1 > 0, \quad H\_2 = \begin{vmatrix} a\_1 & a\_3 \\ & \\ \mathbf{1} & a\_2 \end{vmatrix} > 0, \quad \text{and} \quad H\_3 = \begin{vmatrix} a\_1 & a\_3 & \mathbf{0} \\ \mathbf{1} & a\_2 & \mathbf{0} \\ \mathbf{0} & a\_1 & a\_3 \end{vmatrix} > 0. \tag{9}$$

Computations based on the R-H stability criterion can lead to the following expression *<sup>φ</sup>*ð Þ¼ *<sup>b</sup> a b*ð Þ �<sup>1</sup> ð Þ *ab*�*a*þ*rc ab*�*a*þ*rc*þ*abc* � ð Þ *bc*þ*b*�<sup>1</sup> ð Þ *ab*�*a*þ*rc* ð Þ *<sup>b</sup>*�<sup>1</sup> ð Þ *ab*�*a*�*<sup>c</sup>* , [10]. the positive equilibrium *<sup>E</sup>*<sup>þ</sup> is locally asymptotically stable if *ϕ*ð Þ *b* <*r* � *a*.

In addition to the transcritical bifurcation described above, the system develops a Hopf bifurcation for values of *b*>*μ*<sup>1</sup> resulting in periodic solutions due to the pure imaginary eigenvalues �*i* ffiffiffiffiffi *a*2 p .

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

$$H(b) = \frac{rc\Phi(b-1)}{a^2(b-1)^3(ab-a+rc)},\tag{10}$$

where

$$\begin{aligned} \Phi(\mathbf{x}) &= -a(a\mathbf{x} - c)(a\mathbf{x} + rc)\mathbf{x}^2 + ((a+ac)\mathbf{x} + rc + ac)[((c+1)\mathbf{x} + c)(a\mathbf{x} + rc) \\\\ &+ \mathbf{x}(a\mathbf{x} - c)(r - a)] \\\\ &= -a^3\mathbf{x}^4 + a^2(3c + c^2 + r - a - ac + 1)\mathbf{x}^3 + ac(3rc + 3a + rc^2 + 3ac + r) \\\\ &+ r^2 - a^2)\mathbf{x}^2 + c^2(3ar + 2acr + r^2c + 2a^2)\mathbf{x} + rc^3(r + a) \end{aligned}$$

Then, the another threshold, named *μ*<sup>2</sup> can be defined as the smallest number *b* of the set *I*<sup>0</sup> ¼ *b*> *μ*<sup>1</sup> f g : *H b*ð Þ¼ 0 .

The analysis presented in [10] also shows that if *<sup>c</sup>*ð Þ <sup>1</sup>þ*<sup>r</sup> ab*�*a*þ*rc* <sup>&</sup>lt;1, then *<sup>E</sup>*<sup>þ</sup> represents a partial success of virotherapy at a modest value of *b*. The expression *<sup>c</sup>*ð Þ <sup>1</sup>þ*<sup>r</sup> ab*�*a*þ*rc* is called the tumor load.

When the value of *b* satisfies *μ*<sup>1</sup> <*b*<*μ*2, *E*<sup>þ</sup> is locally asymptotically stable. For values of *b*<*μ*<sup>1</sup> the positive equilibrium *E*<sup>þ</sup> 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 *E*1. Finally, when *b*>*μ*2, periodic solutions will appear as a result of the Hope bifurcation (**Table 1**).


**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

$$\begin{aligned} \frac{d\mathbf{x}}{dt} &= \lambda \mathbf{x} \left( \mathbf{1} - \frac{\mathbf{x} + \mathbf{y}}{C} \right) - \beta \mathbf{x} v \\ \frac{d\mathbf{y}}{dt} &= \beta \mathbf{x} v - \mu y \mathbf{z} - \delta \mathbf{y} \\ \frac{dv}{dt} &= b \delta y - \beta \mathbf{x} v - kv \mathbf{z} - \gamma v \\ \frac{dz}{dt} &= syz - \rho \mathbf{z}, \end{aligned} \tag{11}$$

Where *<sup>N</sup>* <sup>¼</sup> *<sup>n</sup>=m*, *<sup>A</sup>* ¼ �*a*<sup>2</sup> <sup>þ</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

*Mathematical Modeling and Dynamics of Oncolytic Virotherapy*

*<sup>a</sup>*<sup>2</sup> <sup>¼</sup> *<sup>c</sup>*

*<sup>d</sup> <sup>N</sup>* <sup>þ</sup> *<sup>r</sup>*

*<sup>a</sup>* ð Þ *N* � 1 .

proposed by Tian. For details, see [13].

*DOI: http://dx.doi.org/10.5772/intechopen.96963*

values of *b* exceed the critical point.

different behavior over the dynamical landscape.

values of the burst size *b*.

undetectible level of tumor load.

**4. Fractional derivative approach**

definition. The system can be formulated as

**145**

*a*2 <sup>2</sup> � 3*a*<sup>1</sup> � � <sup>p</sup> *<sup>=</sup>*3, *<sup>a</sup>*<sup>1</sup> <sup>¼</sup> *rN*

The analysis of model (12) is similar but more complicated than the basic model

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 *E*0. If the effect of the immune system and the viruses is neglected, then the system has another equilibrium point *E*<sup>1</sup> 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 *E*<sup>2</sup> 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

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 *E*<sup>1</sup> and *E*2, and a Hopf bifurcation that occurs for larger

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 *E*3, *E*4, and *E*5, numerical simulations become a need in order to capture the

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

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 8 � 9 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

*<sup>a</sup>*<sup>2</sup> � *<sup>c</sup>*

*<sup>d</sup>* �*<sup>a</sup>* <sup>þ</sup> *aN* � *<sup>e</sup>* <sup>þ</sup> *<sup>d</sup>*

� �, and

*c*

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* ¼ *μC=δ*, *d* ¼ *kC=δ*,*e* ¼ *γ=δ*, *m* ¼ *sC=δ*, and *n* ¼ *ρ=δ:* Then system (3.1) becomes

$$\begin{aligned} \frac{dx}{dt} &= r\mathbf{x}(1-\mathbf{x}-\mathbf{y}) - a\mathbf{x}v\\ \frac{dy}{dt} &= a\mathbf{x}v - cy\mathbf{z} - \mathbf{y} \\ \frac{dv}{dt} &= by - a\mathbf{x}v - dv\mathbf{z} - ev \\ \frac{dz}{dt} &= my\mathbf{z} - v\mathbf{z} \end{aligned} \tag{12}$$

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;

$$\begin{aligned} E\_0 &= (\mathbf{0}, \mathbf{0}, \mathbf{0}, \mathbf{0}), \\ E\_1 &= (\mathbf{1}, \mathbf{0}, \mathbf{0}, \mathbf{0}), \\ E\_2 &= \left( \frac{\epsilon}{a(b-1)}, \frac{r(ab-a-\epsilon)}{a(b-1)(ab-a+\pi)}, \frac{r(ab-a-\epsilon)}{a(ab-a+\pi)}, \mathbf{0} \right), \\ E\_3 &= \left( \mathbf{1} - N - \frac{a\mathbf{A}}{r}, N, A, \frac{(b-1)N - \epsilon\mathbf{A}}{cN + d\mathbf{A}} \right), \\ E\_4 &= \left( \mathbf{1} - N - \frac{v\_2}{q}, N, v\_2, \frac{(b-1)N - \epsilon v\_2}{cN + dv\_2} \right), \text{ and} \\ E\_5 &= \left( \mathbf{1} - N - \frac{v\_3}{q}, N, v\_3, \frac{(b-1)N - \epsilon v\_3}{cN + dv\_3} \right). \end{aligned}$$

*Mathematical Modeling and Dynamics of Oncolytic Virotherapy DOI: http://dx.doi.org/10.5772/intechopen.96963*

Where *<sup>N</sup>* <sup>¼</sup> *<sup>n</sup>=m*, *<sup>A</sup>* ¼ �*a*<sup>2</sup> <sup>þ</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi *a*2 <sup>2</sup> � 3*a*<sup>1</sup> � � <sup>p</sup> *<sup>=</sup>*3, *<sup>a</sup>*<sup>1</sup> <sup>¼</sup> *rN <sup>a</sup>*<sup>2</sup> � *<sup>c</sup> <sup>d</sup>* �*<sup>a</sup>* <sup>þ</sup> *aN* � *<sup>e</sup>* <sup>þ</sup> *<sup>d</sup> c* � �, and *<sup>a</sup>*<sup>2</sup> <sup>¼</sup> *<sup>c</sup> <sup>d</sup> <sup>N</sup>* <sup>þ</sup> *<sup>r</sup> <sup>a</sup>* ð Þ *N* � 1 .

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 *E*0. If the effect of the immune system and the viruses is neglected, then the system has another equilibrium point *E*<sup>1</sup> 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 *E*<sup>2</sup> 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 *E*<sup>1</sup> and *E*2, 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 *E*3, *E*4, and *E*5, 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 8 � 9 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

$$\begin{aligned} \mathbf{D}\_t^a \mathbf{x} &= \mathbf{r}^a \mathbf{x} (\mathbf{1} - \mathbf{x} - \mathbf{y}) - \mathbf{a}^a \mathbf{x} v \\ \mathbf{D}\_t^a \mathbf{y} &= \mathbf{a}^a \mathbf{x} v - \mathbf{y} \\ \mathbf{D}\_t^a &= \mathbf{b}^a \mathbf{y} - \mathbf{a}^a \mathbf{x} v - \mathbf{c}^a v. \end{aligned} \tag{13}$$

where *D<sup>α</sup> <sup>t</sup>* 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 *<sup>L</sup>*<sup>1</sup> <sup>¼</sup> *<sup>L</sup>*<sup>1</sup> ½ � *a*, *b* be the class of Lebesgue integrable functions on ½ � *a*, *b* , *a*<*b*< ∞. The fractional integral of order *ν*∈ <sup>þ</sup> of the function *f t*ð Þ, *t* >0 ð Þ *f* : ! is defined by

$$I\_a^p = \frac{1}{\Gamma(v)} \int\_q^t (t - s)^{v - 1} f(s) ds, \ t > 0,\tag{14}$$

where Γð Þ*:* is the Gamma function.

The fractional derivative of order *α* ∈ð Þ *n* � 1, *n* of the function *f t*ð Þ is defined by several ways, the most common ones are:

1.Riemann-Liouville fractional derivative: Take the fractional integral of order ð Þ *<sup>n</sup>* � *<sup>α</sup>* and then apply the *<sup>n</sup>*th derivative

$$D\_a^a f(t) = D\_a^a I\_{n-a}^a,$$

where *D*<sup>∗</sup> *<sup>n</sup>* <sup>¼</sup> *dn dtn* , *n* ¼ 1, 2, … ;

2.Caputo's fractional derivative: Start with a *n*th derivative of the function, then take a fractional integral of order ð Þ *n* � *α*

$$D\_{\mathfrak{q}}^{a}f(t) = I\_{n-a}^{a}D\_{\mathfrak{q}}^{a}f(t), n = 1, 2, \dots$$

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 1**–**4**.

**5. Optimization by control theory**

*Dynamics of Tumor cells vs infected tumor cells of the model when α* ¼ *:*98*.*

**Figure 4.**

**147**

**Figure 2.**

**Figure 3.**

*α* ¼ 1*.*

*α* ¼ 1*.*

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

*Damped oscillators when b* ¼ 27 *and initial values x* ¼ 0*:*5*, y* ¼ 0*:*5*, and ν* ¼ 1*:*5*, for α* ¼ *:*96*, α* ¼ *:*98*, and*

*Mathematical Modeling and Dynamics of Oncolytic Virotherapy*

*DOI: http://dx.doi.org/10.5772/intechopen.96963*

*Dynamics of model when b* ¼ 28 *and initial values x* ¼ 0*:*5*, y* ¼ 0*:*5*, and ν* ¼ 1*:*5*, for α* ¼ *:*96*, α* ¼ *:*98*, and*

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

optimal ways that can minimize a tumor and the cost of the virotherapy.

**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*.*

*Mathematical Modeling and Dynamics of Oncolytic Virotherapy DOI: http://dx.doi.org/10.5772/intechopen.96963*

**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*.*
