Mathematical Modeling and Dynamics of Oncolytic Virotherapy

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.

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

$$\begin{array}{ll}\frac{\mathit{dx}}{\mathit{dt}}& =\mathit{xF}\left(x,y\right)-\mathit{\beta yG}\left(x,y\right)\\ \frac{\mathit{dy}}{\mathit{dt}}& =\mathit{\beta yG}\left(x,y\right)-\mathit{ay},\end{array}$$E1

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 $\beta$ represents the infectivity of the virus, and the death rate $\mathit{ay}$ represents the virus-infected cells die.

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

$$\begin{array}{ll}\frac{\mathit{dx}}{\mathit{dt}}& =\mathit{rx}\left(1-\frac{x+y}{C}\right)-\mathit{dx}-\mathit{\beta xv}\\ \frac{\mathit{dy}}{\mathit{dt}}& =\mathit{\beta xv}-\left(d+a\right)y\\ \frac{\mathit{dv}}{\mathit{dt}}& =\mathit{\alpha y}-\mathit{\gamma v},\end{array}$$E2

in which $v$ stands for the free virus population and $C$ is maximal tumor size. The term $\mathit{\alpha y}$ models the release of virions by infected tumor cells, and $\mathit{\gamma 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 $\mathit{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.

$$\begin{array}{ll}\frac{\mathit{dx}}{\mathit{dt}}& =\mathit{\lambda x}\left(1-\frac{x+y}{C}\right)-\mathit{\beta xv}\\ \frac{\mathit{dy}}{\mathit{dt}}& =\mathit{\beta xv}-\mathit{\delta y}\\ \frac{\mathit{dv}}{\mathit{dt}}& =\mathit{b\delta y}-\mathit{\beta xv}-\mathit{\gamma v}.\end{array}$$E3

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 $\tau =\mathit{\delta t}$, $x=K\widehat{x}$, $y=K\widehat{y}$, $v=K\widehat{v}$, and rename parameters $r=\lambda ,a=\mathit{\beta K}$, and $c=\gamma .$

Then dropping all hats over the variables and write $\tau$ as $t$, we have

$$\begin{array}{ll}\frac{\mathit{dx}}{\mathit{dt}}& =\mathit{rx}\left(1-x-y\right)-\mathit{axv}\\ \frac{\mathit{dy}}{\mathit{dt}}& =\mathit{axv}-y\\ \frac{\mathit{dv}}{\mathit{dt}}& =\mathit{by}-\mathit{axv}-\mathit{cv}\end{array}$$E4

It is assumed that all parameters are nonnegative.

Model (4) has three equilibrium points, $E_0=\left(\mathrm{0,0,0}\right)$, $E_1=\left(\mathrm{1,0,0}\right)$, and the positive equilibrium $E_{+}=\left(x^{\ast },y^{\ast },v^{\ast }\right)$. The equilibrium $E_0$ is always unstable for all positive values of the burst size $b$. The equilibrium $E_1$ is globally asymptotically stable when $0<b<{\mu }_1$, and it is unstable when $b\ge {\mu }_1$.

At $b={\mu }_1$, the positive equilibrium $E_{+}$ moves into the domain $D=\left\{\left(x,y,v\right):x\ge 0,y\ge 0,v\ge 0,0\le x+y\le 1\right\}$, a type of transcritical bifurcation occurs with $E_1$ and $E_{+}$. As the parameter $b$ increases, while ${\mu }_1<b<{\mu }_2$ and $b\in I_p$, $E_{+}$ is locally asymptotically stable. When $b>{\mu }_1$ and $b\in I_n$, $E_{+}$ is unstable. Hopf bifurcations occur for some $b\ge {\mu }_2$, and these bifurcations give rise to one or three families of periodic solutions. Here, ${\mu }_1$ and ${\mu }_2$ are thresholds, $I_p=\left\{b>{\mu }_1:H\left(b\right)>0\right\}$, $I_n=\left\{b>{\mu }_1:H\left(b\right)<0\right\}$, and $H\left(b\right)$ is defined next in formula (10).

• $E_0$ is unstable.

• $E_1$ is globally asymptotically stable when $b<1+\frac{c}{a}$ ; and unstable when $b\ge 1+\frac{c}{a}$ .

• When ${\mu }_1<b<{\mu }_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={\mu }_1$ introduces the equilibrium point $E_{+}$ into the positive invariant domain $D$. The Hopf bifurcation at some value $b>{\mu }_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 $E_0=\left(\mathrm{0,0,0}\right)$, $E_1\left(\mathrm{1,0,0}\right)$, and the positive equilibrium $E_{+}\left(x^{\ast },y^{\ast },v^{\ast }\right)$, where $x^{\ast }=\frac{c}{a\left(b-1\right)}$, $y^{\ast }=\frac{\mathit{rc}\left(\mathit{ab}-a-c\right)}{a\left(b-1\right)\left(\mathit{ab}-a+\mathit{rc}\right)}$, and $v^{\ast }=\frac{r\left(\mathit{ab}-a-c\right)}{a\left(\mathit{ab}-a+\mathit{rc}\right)}.$

The Jacobian matrix is given by

$$J=\left(\begin{array}{ccc}r-2\mathit{rx}-\mathit{ry}-\mathit{av}& -\mathit{rx}& -\mathit{ax}\\ \mathit{av}& -1& \mathit{ax}\\ -\mathit{av}& b& -\mathit{ax}-c\end{array}\right)$$E5

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

$$J\left(E_0\right)=\left(\begin{array}{ccc}r& 0& 0\\ 0& -1& 0\\ 0& b& -c\end{array}\right)$$E6

The obtained eigenvalues of this lower triangular matrix are ${\lambda }_1=r$, ${\lambda }_2=-1$, and ${\lambda }_3=-c$. Note that ${\lambda }_1$ is positive due the $r$ being positive. The other eigenvalues are negative. Therefore, $E_0$ is unstable. The local unstable invariant manifold lives in the $x-$axis. The stable invariant manifolds live in the $\mathit{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_1$ is given by

$$J\left(E_1\right)=\left(\begin{array}{ccc}-r& -r& a\\ 0& -1& a\\ 0& b& -a-c\end{array}\right)$$E7

The eigenvalues are ${\lambda }_1=-r$, which is clearly negative, and

${\lambda }_{2,3}=\frac{1}{2}\left(-\left(1+a+c\right)\pm \sqrt{{\left(1-a-c\right)}^2+4\mathit{ab}}\right)$. ${\lambda }_2$, the eigenvalue with negative square root, is negative as well. However, for ${\lambda }_3$, it is not so clear when it is negative. Computations in [7], show that for $1<b<1+\frac{c}{a}$, $E_1$ is a global asymptotically stable. This result can be obtained by applying appropriate Lyapunov functions such as $V_1\left(x,y,v\right)=y+v$ and $V_2\left(x,y,v\right)=\frac{1}{2}\mathit{ab}\left(a+c\right)y^2+a^2\mathit{byv}+\frac{1}{2}{a}^2{v}^2$. When the value of $b$ exceeds the threshold ${\mu }_1=1+\frac{c}{a}$, $E_1$ becomes unstable. The system exhibits a transcritical bifurcation with bifurcation value ${\mu }_1$ and changes stability as this parameter varies near the bifurcation value. When $b>{\mu }_1$, the positive equilibrium point $E_{+}$ appears. Proceeding with the linearization, the Jacobian materix at $E_{+}$ is given by

$$J\left(E_{+}\right)=\left(\begin{array}{ccc}r-2\mathit{rx}-\mathit{ry}-\mathit{av}& -\mathit{rx}& -\mathit{ax}\\ \mathit{av}& -1& \mathit{ax}\\ -\mathit{av}& b& -\mathit{ax}-c\end{array}\right).$$E8

To analyze the more complicated eigenvalue expressions, we apply the so called Routh-Hurwitz criterion on the associated characteristic polynomial $P\left(\lambda \right)={\lambda }^3+a_1{\lambda }^2+a_2\lambda +a_3$, where $a_1=\frac{\mathit{rc}+\mathit{ab}-a+\mathit{abc}}{a\left(b-1\right)}$, $a_2=\frac{\mathit{rc}\left(\mathit{bc}+b-1\right)}{a{\left(b-1\right)}^2}+\frac{\mathit{rc}\left(\mathit{ab}-a-c\right)\left(r-a\right)}{a\left(b-1\right)\left(\mathit{ab}-a+\mathit{rc}\right)}$, and $a_3=\frac{\mathit{rc}\left(\mathit{ab}-a-c\right)}{a\left(b-1\right)}$.

The R-H criterion requires

$$H_1=a_1>0,H_2=\left|\begin{array}{cc}{a}_1& a_3\\ 1& a_2\end{array}\right|>0,\text{and}{H}_3=\left|\begin{array}{ccc}{a}_1& a_3& 0\\ 1& a_2& 0\\ 0& a_1& a_3\end{array}\right|>0.$$E9

Computations based on the R-H stability criterion can lead to the following expression $\phi \left(b\right)=\frac{a\left(b-1\right)\left(\mathit{ab}-a+\mathit{rc}\right)}{\mathit{ab}-a+\mathit{rc}+\mathit{abc}}-\frac{\left(\mathit{bc}+b-1\right)\left(\mathit{ab}-a+\mathit{rc}\right)}{\left(b-1\right)\left(\mathit{ab}-a-c\right)}$, [10]. the positive equilibrium $E_{+}$ is locally asymptotically stable if $\varphi \left(b\right)<r-a$.

In addition to the transcritical bifurcation described above, the system develops a Hopf bifurcation for values of $b>{\mu }_1$ resulting in periodic solutions due to the pure imaginary eigenvalues $\pm i\sqrt{a_2}$.

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

$$H\left(b\right)=\frac{\mathit{rc}\mathrm{\Phi }\left(b-1\right)}{a^2{\left(b-1\right)}^3\left(\mathit{ab}-a+\mathit{rc}\right)},$$E10

where

$$\begin{array}{ll}\mathrm{\Phi }\left(x\right)& =-a\left(\mathit{ax}-c\right)\left(\mathit{ax}+\mathit{rc}\right)x^2+\left(\left(a+\mathit{ac}\right)x+\mathit{rc}+\mathit{ac}\right)[\left(\left(c+1\right)x+c\right)\left(\mathit{ax}+\mathit{rc}\right)\\ & +x\left(\mathit{ax}-c\right)\left(r-a\right)]\\ & =-a^3{x}^4+a^2\left(3c+c^2+r-a-\mathit{ac}+1\right)x^3+\mathit{ac}\left(3\mathit{rc}+3a+{\mathit{rc}}^2+3\mathit{ac}+r\right.\\ & \left.+r^2-a^2\right)x^2+c^2\left(3\mathit{ar}+2\mathit{acr}+r^{2}c+2a^2\right)x+{\mathit{rc}}^3\left(r+a\right)\end{array}.$$

Then, the another threshold, named ${\mu }_2$ can be defined as the smallest number $b$ of the set $I_0=\left\{b>{\mu }_1:H\left(b\right)=0\right\}$.

The analysis presented in [10] also shows that if $\frac{c\left(1+r\right)}{\mathit{ab}-a+\mathit{rc}}<1$, then $E_{+}$ represents a partial success of virotherapy at a modest value of $b$. The expression $\frac{c\left(1+r\right)}{\mathit{ab}-a+\mathit{rc}}$ is called the tumor load.

When the value of $b$ satisfies ${\mu }_1<b<{\mu }_2$, $E_{+}$ is locally asymptotically stable. For values of $b<{\mu }_1$ the positive equilibrium $E_{+}$ does not live in the positive invarient domain $D$, and as $b$ increases to ${\mu }_1$, the equilibrium point moves into the domain $D$ and it coalecses with $E_1$. Finally, when $b>{\mu }_2$, periodic solutions will appear as a result of the Hope bifurcation (Table 1).

Parameter	Description	Value	Dimentions
$\lambda$	Tumor growth rate	$2\times {10}^{-2}$	$1/h$
$\delta$	Death rate of infected tumor cells	$1/18$	$1/h$
$\beta$	Infection rate of the virus	$7/10\times {10}^{-9}$	${\mathit{mm}}^{3}h/$ virusl
$k$	Immune killing rate of virus	${10}^{-8}$	${\mathit{mm}}^{3}h/$ immune cell
$b$	Burst size of free virus	$50$	viruses/cell
$\gamma$	Clearance rate of virus	$2.5\times {10}^{-2}$	$1/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

$$\begin{array}{l}\frac{\mathit{dx}}{\mathit{dt}}=\mathit{\lambda x}\left(1-\frac{x+y}{C}\right)-\mathit{\beta xv}\\ \frac{\mathit{dy}}{\mathit{dt}}=\mathit{\beta xv}-\mathit{\mu yz}-\mathit{\delta y}\\ \frac{\mathit{dv}}{\mathit{dt}}=\mathit{b\delta y}-\mathit{\beta xv}-\mathit{kvz}-\mathit{\gamma v}\\ \frac{\mathit{dz}}{\mathit{dt}}=\mathit{syz}-\mathit{\rho z},\end{array}$$E11

where $\lambda$ is tumor growth rate, $C$ is the carrying capacity of tumor cell population, $\beta$ is the infected rate of the virus, $\mu$ is immune killing rate of infected tumor cells, $\delta$ 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, $\gamma$ is clearance rate of viruses, $s$ is the stimulation rate of the innate immune system, and $\rho$ is immune clearance rate.

We non-dimensionalize the system by setting $\tau =\mathit{\delta t},x=$ $C\overline{x},y=C\overline{y},v=C\overline{v},z=C\overline{z}$ and rename parameters $r=\lambda /\delta$ $a=\mathit{C\beta }/\delta ,c=\mathit{\mu C}/\delta ,d=\mathit{kC}/\delta ,e=\gamma /\delta ,m=\mathit{sC}/\delta ,$ and $n=\rho /\delta .$ Then system (3.1) becomes

$$\begin{array}{l}\frac{\mathit{dx}}{\mathit{dt}}=\mathit{rx}\left(1-x-y\right)-\mathit{axv}\\ \frac{\mathit{dy}}{\mathit{dt}}=\mathit{axv}-\mathit{cyz}-y\\ \frac{\mathit{dv}}{\mathit{dt}}=\mathit{by}-\mathit{axv}-\mathit{dvz}-\mathit{ev}\\ \frac{\mathit{dz}}{\mathit{dt}}=\mathit{myz}-\mathit{nz}.\end{array}$$E12

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;

$E_0=\left(\mathrm{0,0,0,0}\right)$,

$E_1=\left(\mathrm{1,0,0,0}\right)$,

$E_2=\left(\frac{e}{a\left(b-1\right)},\frac{\mathit{re}\left(\mathit{ab}-a-e\right)}{a\left(b-1\right)\left(\mathit{ab}-a+\mathit{re}\right)},\frac{r\left(\mathit{ab}-a-e\right)}{a\left(\mathit{ab}-a+\mathit{re}\right)},0\right)$,

$E_3=\left(1-N-\frac{\mathit{aA}}{r},N,A,\frac{\left(b-1\right)N-\mathit{eA}}{\mathit{cN}+\mathit{dA}}\right)$,

$E_4=\left(1-N-\frac{v_2}{q},N,v_2,\frac{\left(b-1\right)N-{\mathit{ev}}_2}{\mathit{cN}+{\mathit{dv}}_2}\right)$, and

$E_5=\left(1-N-\frac{v_3}{q},N,v_3,\frac{\left(b-1\right)N-{\mathit{ev}}_3}{\mathit{cN}+{\mathit{dv}}_3}\right)$.

Where $N=n/m$, $A=\left(-a_2+\sqrt{a_2^2-3a_1}\right)/3$, $a_1=\frac{\mathit{rN}}{a^2}\cdot \frac{c}{d}\left(-a+\mathit{\text{aN}}-e+\frac{d}{c}\right)$, and $a_2=\frac{c}{d}N+\frac{r}{a}\left(N-1\right)$.

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_1$ 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_2$ 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_1$ 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{array}{ll}{\mathrm{D}}_{\mathrm{t}}^{\alpha }x={\mathrm{r}}^{\alpha }x\left(1-x-\mathrm{y}\right)-{\mathrm{a}}^{\alpha }\mathrm{x}v& \\ {\mathrm{D}}_{\mathrm{t}}^{\alpha }y={\mathrm{a}}^{\alpha }\mathit{xv}-y& \\ {\mathrm{D}}_{\mathrm{t}}^{\alpha }={\mathrm{b}}^{\alpha }\mathrm{y}-{\mathrm{a}}^{\alpha }\mathrm{x}v-{\mathrm{c}}^{\alpha }v.& \end{array}$$E13

where $D_t^{\alpha }$ is the Caputo fractional derivative and $0<\alpha \le 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 $L^1=L^1\left[a,b\right]$ be the class of Lebesgue integrable functions on $\left[a,b\right]$, $a<b<\infty$. The fractional integral of order $\nu \in {\mathbb{R}}^{+}$ of the function $f\left(t\right),t>0$ $\left(f:\mathbb{R}\to \mathbb{R}\right)$ is defined by

$$I_a^v=\frac{1}{\mathrm{\Gamma }\left(v\right)}{\int }_q^t{\left(t-s\right)}^{v-1}f\left(s\right)\mathit{ds},t>0,$$E14

where $\mathrm{\Gamma }\left(.\right)$ is the Gamma function.

The fractional derivative of order $\alpha \in \left(n-1,n\right)$ of the function $f\left(t\right)$ is defined by several ways, the most common ones are:

1. Riemann-Liouville fractional derivative: Take the fractional integral of order $\left(n-\alpha \right)$ and then apply the $n^{\mathrm{th}}$ derivative

   $$D_{\alpha }^{a}f\left(t\right)=D_{\alpha }^a{I}_{n-\alpha }^a,$$

   where $D_n^{\ast }=\frac{d^n}{{\mathit{dt}}^n},n=1,2,\dots$;

2. Caputo’s fractional derivative: Start with a $n^{\mathrm{th}}$ derivative of the function, then take a fractional integral of order $\left(n-\alpha \right)$

$$D_{\alpha }^{a}f\left(t\right)=I_{n-\alpha }^a{D}_n^{a}f\left(t\right),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.

Figure 1.

Dynamics of virotherapy when $b=4$ and initial values $x=0.5$, $y=0.5$, and $\nu =1.5$, for $\alpha =1$, $\alpha =.8$, and $\alpha =.9$.

Figure 2.

Damped oscillators when $b=27$ and initial values $x=0.5$, $y=0.5$, and $\nu =1.5$, for $\alpha =.96$, $\alpha =.98$, and $\alpha =1$.

Figure 3.

Dynamics of model when $b=28$ and initial values $x=0.5$, $y=0.5$, and $\nu =1.5$, for $\alpha =.96$, $\alpha =.98$, and $\alpha =1$.

Figure 4.

Dynamics of Tumor cells vs infected tumor cells of the model when $\alpha =.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 $u\left(t\right)$ is given by

$$\mathrm{J}\left(\mathfrak{u}\left(\mathrm{t}\right)\right)=\mathrm{\Psi }\left(x\left(\mathrm{T}\right)\right)+{\int }_0^1\mathrm{L}\left(x\left(\mathrm{t}\right),\mathrm{u}\left(\mathrm{t}\right),\mathrm{T}\left(\mathrm{t}\right)\right)\mathrm{dt}$$

where $x\left(t\right)$ is the system state which evolves according to the state equation

$$\dot{x}=\mathbf{f}\left(x\left(t\right),u\left(t\right),t\right)x\left(0\right)=x_{0}t\in \left[0,T\right]$$

The Hamiltonian is defined as

$$\mathrm{H}\left(x,\mathrm{\Psi },u,\mathrm{t}\right)={\mathrm{\Psi }}^{\mathrm{T}}\left(\mathrm{t}\right)\mathrm{f}\left(x,u,\mathrm{t}\right)+\mathrm{L}\left(x,u,\mathrm{t}\right),$$

where $\mathrm{\Psi }\left(\mathrm{t}\right)$ is a vector of costate variables of the same dimension as the state variable $x\left(\mathrm{t}\right)$ such that, [18]

$$\dot{\mathrm{\Psi }}\left(t\right)=-\frac{\partial H}{\partial x}.$$

Applying the control theory approach, we reformulate the basic model by introducing a control function $u\left(t\right)$ which represents efforts on damaging the tumor cells AND $\left(1-u\left(t\right)\right)$ represents the growth rate of the infected cells. After incorporating the control $u$ into the basic model, we obtain the following model with control

$$\begin{array}{l}\frac{\mathrm{dx}}{\mathrm{dt}}=\left(1-\mathrm{u}\left(\mathrm{t}\right)\right)\lambda \mathrm{x}\left(1-\frac{\mathrm{x}+\mathrm{y}}{\mathrm{K}}\right)-\beta \mathrm{x}v\\ \frac{\mathrm{dy}}{\mathrm{d}t}=\mathit{\beta xv}-\left(1-\mathrm{u}\left(\mathrm{t}\right)\right)\delta \mathrm{y}\\ \frac{\mathrm{d}v}{\mathrm{dt}}=\mathrm{b}\delta \left(1-\mathrm{u}\left(\mathrm{t}\right)\right)\mathrm{y}-\mathit{\beta xv}-\mathit{\gamma v}.\end{array}$$E15

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^{\ast }$ and it is given by

$$\mathrm{J}\left(u\left(\mathrm{t}\right)\right)={\int }_0^{\mathrm{T}}\mathrm{y}\left(\mathrm{t}\right)+\frac{1}{2}{\mathrm{Bu}}^2\mathrm{dt}.$$

Where $B$ is a measure of the relative cost of interventions associated to the control $u\left(t\right)$. 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

$$\mathrm{\Omega }=\left\{u\left(\cdot \right)\in {\mathrm{L}}^{\infty }\left(0,t_f\right):0⩽u\left(t\right)⩽u_{\max },\forall t\in [0,T]\right\}$$

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

$$\begin{array}{l}\frac{\mathrm{d}x}{\mathrm{d}t}=\left(1-u\right)\mathit{rx}\left(1-x-y\right)-\mathit{axv}\\ \frac{\mathrm{d}y}{\mathrm{d}t}=\mathit{axv}-\left(1-u\right)y\\ \frac{\mathrm{d}v}{\mathrm{d}t}=b\left(1-u\right)y-\mathit{axv}-\mathit{cv}.\end{array}$$E16

with fixed initial conditions $x\left(0\right),y\left(0\right),v\left(0\right)$ and a fixed final time $\mathrm{T}$.

According to The Pontryagin’s Maximum Principle, if $u\left(.\right)\in \mathrm{\Omega }$ 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 $\mathrm{\Psi }:\left[0,1\right]\to {\mathbb{R}}^3$.

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 $\left(x^{\ast }\left(\cdot \right),y^{\ast }\left(\cdot \right),v^{\ast }\left(\cdot \right)\right)$ associated to an optimal control $\left(u^{\ast }\right)$ on $\left[0,T\right].$ Moreover, there exists adjoint functions ${\mathrm{\Psi }}_1^{\ast },{\mathrm{\Psi }}_1^{\ast },$ and ${\mathrm{\Psi }}_3^{\ast }$, such that with transversality conditions ${\mathrm{\Psi }}_i^{\ast }\left(\mathrm{T}\right)=0,i=\mathrm{1,2,3}.$ Furthermore,

$$u^{\ast }\left(t\right)=\frac{\mathit{r}x\left(1-x-y\right){\mathrm{\Psi }}_1-y\left({\mathrm{\Psi }}_2+b{\mathrm{\Psi }}_3\right)}{B}.$$

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 1–4 below. Note that $\alpha =1$ represents the simulation of the basic model (2.2), whereas the other values of $\alpha$ 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 $E_0=\left(\mathrm{0,0,0,0}\right)$, $E_1=\left(\mathrm{1,0,0}\right)$, $E_2=\left(\mathrm{0.6,0.0730,2.5729}\right)$. Here the bifurcation parameter values are ${\mu }_1=5$ and ${\mu }_2=27.766$. When $5<b<27.766$, $E_{+}$ is locally asymptotically stable while $E_1$ is unstable. The equilibrium point $E_0$ is always unstable. Figure 1 shows the treatment will eventually reach the equilibrium point $E_1$ 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 $\alpha =0.98$ which is almost the same as the result of the usual derivative $\alpha =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 5–8.

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^{\ast }$ 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 5–8.