# A New Fractional Model for Cancer Therapy with M1 Oncolytic Virus

Authors: Majda El Younoussi; Zakaria Hajhouji; Khalid Hattaf; Noura Yousfi
Journal: Complexity
 (2021)

---

## Abstract

The aim of this work is to propose and analyze a new mathematical model formulated by fractional differential equations (FDEs) that describes the dynamics of oncolytic M1 virotherapy. The well-posedness of the proposed model is proved through existence, uniqueness, nonnegativity, and boundedness of solutions. Furthermore, we study all equilibrium points and conditions needed for their existence. We also analyze the global stability of these equilibrium points and investigate their instability conditions. Finally, we state some numerical simulations in order to exemplify our theoretical results.


---

## Body

## 1. Introduction

Cancer is a collection of related diseases where some of the body’s cells divide continuously and spread into surrounding tissues. Cancer is caused by certain changes to genes. It can start almost anywhere in the human body. Old or damaged cells survive when they should die; new cells form when the body does not need them. These extra cells can divide continuously and may form tumors. A tumor becomes dangerous when it begins to form extensions to neighboring areas (metastasis) [1]. This is why it is important to detect cancer as early as possible in order to avoid this migration. Cancer treatment is adapted according to each situation. There are different cancer treatments used alone or in combination, such as surgery, radiotherapy, chemotherapy, hormone therapy, immunotherapy, and virotherapy. Virotherapy is one of the new therapies; it consists in using a virus after having reprogrammed it. This virus is called oncolytic virus. Oncolytic viruses infect and destroy cancer cells; they use the cell’s genetic machinery to make copies of themselves and subsequently spread to surrounding uninfected cells [2].According to a medical experiment, in vitro, in vivo, and ex vivo studies showed potent oncolytic efficacy and high tumor tropism of alphavirus M1, which is a naturally occurring and a selective oncolytic virus targeting zinc-finger antiviral protein (ZAP) deficient cancer cells [3]. To model the role of the M1 virus in oncolytic virotherapy, Wang et al. [4] proposed a nonlinear system governed by ordinary differential equations (ODEs) that describe the growth of normal cells, tumor cells, and the M1 virus with limited nutrients. Elaiw et al. [5] extended the model presented in [4] by including spatial effects and anti-tumor immune response mediated by cytotoxic T lymphocyte (CTL) cells. The results in [5] indicated that the immune response has a negative impact on oncolytic M1 virotherapy, and it reduced its efficiency.On the other hand, all the above mathematical models neglected the memory effect by considering only integer-order derivatives. However, fractional-order derivative provides an excellent tool for describing memory and hereditary properties which exist in most biological systems. For instance, Cole [6] proved that the membranes of cells of the biological organism have fractional-order electrical conductance since the memory means that the system’s response is dependent not only on the current state but on its complete history. Therefore, the classical integer-order derivative does not reflect this memory effect because it is a local operator, unlike the fractional derivative.The main purpose of this study is to develop a mathematical model governed by fractional-order differential equations (FDEs) to study the effect of memory on the dynamics of oncolytic M1 virotherapy. So, the rest of the paper is outlined as follows: the next section is devoted to the formulation of the model, including the well-posedness and the existence of equilibria. Section3 focuses on stability analysis. Section 4 deals with numerical simulations in order to illustrate our main analytical results. Finally, a brief conclusion is given in Section 5.


## 2. Model Formulation and Preliminaries

In this section, we propose the following FDE model:(1)DαSt=A−dSt−β1StNt−β2StTt,DαNt=r1β1StNt−d+ε1Nt,DαTt=r2β2StTt−d+ε2Tt−β3TtVt,DαVt=B+r3β3TtVt−d+ε3Vt,where St, Nt, Tt, and Vt are the concentrations of nutrient, normal cells, tumor cells, and M1 virus at time t, respectively. The parameters A and B are the recruitment rates of nutrient and M1 virus, respectively. Also, B represents the minimum effective dosage of medication. The normal and tumor cells consume the nutrient at rates β1SN and β2ST, respectively. The growth rate of normal cells as a result of consuming the nutrient is given by r1β1SN, while the growth rate of tumor cells is given by r2β2ST. The virus infects and kills tumor cells at rate β3TV, and it replicates at rate r3β3TV. The parameter d is the washout constant rate of nutrient and bacteria. The parameters ε1, ε2, and ε3 are the natural death rates of normal cells, tumor cells, and M1 virus, respectively. The operator Dα denotes the Caputo fractional derivative with α∈0,1 that describes the memory effect.It is important to note that the ODE mathematical model verifying potent oncolytic efficacy of M1 virus [4] is a special case of our model presented by system (1), and it suffices to take α=1. Furthermore, to prove that our model is biologically well-posed, we assume that the initial conditions of (1) satisfy:(2)S0=ϕ10≥0,N0=ϕ20≥0,T0=ϕ30≥0,V0=ϕ40≥0.Theorem 1.
If the initial conditions (2) are given, then there exists a unique solution of system (1) defined on 0,+∞. Moreover, this solution remains nonnegative and bounded for all t≥0.Proof.
It is not hard to show that the vector function of system (1) satisfies the first condition of Lemma 4 in [7]. It remains to prove the second condition. Let(3)Xt=StNtTtVt,Y=A00B.
Then(4)DαXt=Y+A1Xt+A2StXt+A3TtXt,where(5)A1=−d0000−d+ε10000−d+ε20000−d+ε3,A2=0−β1−β200r1β10000r2β200000,A3=00000000000−β3000r3β3.
Thus,(6)DαXt≤Y+A1+A2S+A3TX.
This implies that the second condition of Lemma 4 in [7] is satisfied. Then system (1) has a unique solution on 0,+∞.
On the other hand and according to (1), we have(7)DαSt|S=0=A>0,DαNt|N=0=0≥0,DαTt|T=0=0≥0,DαVt|V=0=B>0.
By Lemmas 5 and 6 in [7], we deduce that the solution of (1) is nonnegative.
It remains to prove the boundedness of solutions. Then we consider the following function:(8)Ft=r1r2r3St+r2r3Nt+r1r3Tt+r1Vt.
Hence,(9)DαFt=r1r2r3A+r1B−r1r2r3dSt−r2r3d+ε1Nt−r1r3d+ε2Tt−r1d+ε3Vt≤r1r2r3A+B−dFt.
Then(10)Ft≤F0Eα−dtα+r1dr2r3A+B1−Eα−dtα.
Since0≤Eα−dtα≤1, we have(11)Ft≤F0+r1dr2r3A+B,which implies that S, N, T, and V are bounded. This completes the proof.
Now, we establish the equilibrium points of our model. It is obvious that any equilibrium point of system (1) satisfies the following algebraic equations:(12)A−dS−β1SN−β2ST=0,(13)r1β1SN−d+ε1N=0,(14)r2β2ST−d+ε2T−β3TV=0,(15)B+r3β3TV−d+ε3V=0.
From (13), we have N=0 or S=d+ε1/r1β1. Similarly, equation (14) leads to T=0 or r2β2S=d+ε2+β3V:(i)
ForN=0 and T=0, we have S=A/d and V=B/d+ε3. Then system (1) has an equilibrium point of the form E0S0,0,0,V0, where S0=A/d and V0=B/d+ε3.(ii)
ForN≠0 and T=0, we have S=d+ε1/r1β1, V=B/d+ε3 and N=d/β1A1−1, where(16)A1=Ar1β1dd+ε1.This number reflects the ability of absorbing nutrients by normal cells. It is called absorbing number [4]. When A1>1, system (1) has another equilibrium E1S1,N1,0,V1, where S1=d+ε1/r1β1, N1=d/β1A1−1, and V1=B/d+ε3.(iii)
ForN=0 and T≠0, we have S=β3V+d+ε2/r2β2, T=−d/β2+Ar2/β3V+d+ε2, and(17)a1V2+a2V+a3=0,where(18)a1=β3r3β3d+β2d+ε3,a2=a1β3d+ε2−β2β3B+r2r3A,a3=−Bβ2d+ε2.Sincea1>0 and a3<0, we have δ=a22−4a1a3≥0. Thus equation (17) has two roots given by(19)V±=−a2±δ2a1.Clearly,V+>0 and V−<0. As V>0, we have V=V+.It is obvious thatS>0. However, T>0 implies that A2>1+Bβ3/d+ε2d+ε3, where(20)A2=Ar2β2dd+ε2.This number reflects the ability of absorbing nutrients by tumor cells. It can be called the absorbing number of nutrients by tumor cells. Hence, system (1) has another equilibrium point when A2>1+Bβ3/d+ε2d+ε3. This equilibrium point is denoted by E2S2,0,T2,V2, where(21)V2=V+,S2=β3V2+d+ε2r2β2,T2=−dβ2+Ar2β3V2+d+ε2.(iv)
ForN≠0 and T≠0, we have S=d+ε1/r1β1 and V=d+ε2/β3A2/A1−1 as V>0 implies that A2>A1. From (15), we get T=−B+d+ε3V/r3β3V. Similarly, T>0 leads to A2>A1+ABr1β1β3/dd+ε3d+ε2d+ε1. Substituting S and T in (12), we obtain(22)N=d+ε2A2/A1−1Ar1r3β1β3−r3β3dd+ε1−β2d+ε1d+ε3r3β1β3d+ε1d+ε2A2/A1−1+Bβ2β3r3β1β3d+ε2A2/A1−1.
Furthermore,N>0 implies that(23)A1+β2Br3dd+ε2A2/A1−1>1+β2d+ε3r3β3d.
Thus, system (1) has another equilibrium point when A2>A1:(24)A1+β2Br3dd+ε2A2/A1−1>1+β2d+ε3r3β3d,A2>A1+ABr1β1β3dd+ε3d+ε2d+ε1.
This equilibrium point is denoted byE3S3,N3,T3,V3, where(25)S3=d+ε1r1β1,V3=d+ε2β3A2/A1−1,N3=d+ε2A2/A1−1Ar1r3β1β3−r3β3dd+ε1−β2d+ε1d+ε3r3β1β3d+ε1d+ε2A2/A1−1+Bβ2β3r3β1β3d+ε2A2/A1−1,T3=−B+d+ε3V3r3β3V3.
All the above cases are summerized in the following result.Theorem 2.
LetA1 and A2 be defined by (16) and (20). Then(i)
System (1) always has a competition-free equilibrium E0S0,0,0,V0.(ii)
System (1) has a tumor-free equilibrium E1S1,N1,0,V1 when A1>1.(iii)
System (1) has a treatment failure equilibrium E2S2,0,T2,V2 when A2>1+Bβ3/d+ε2d+ε3.(iv)
System (1) has a partial success equilibrium E3S3,N3,T3,V3 when(26)A2>A1+ABr1β1β3dd+ε3d+ε2d+ε1,A1+β2Br3dd+ε2A2/A1−1>1+β2d+ε3r3β3d.


## 3. Stability Analysis

In this section, we focus on the stability analysis of the equilibriaE0, E1, E2, and E3.Theorem 3.
The competition-free equilibriumE0 is globally asymptotically stable for A2≤1+Bβ3/d+ε2d+ε3 and A1≤1, and it is unstable if A2>1+Bβ3/d+ε2d+ε3 or A1>1.Proof.
In order to show the first part of this theorem, we consider the following Lyapunov functional:(27)L0t=S0ϕStS0+1r1Nt+1r2Tt+1r2r3V0ϕVtV0,where ϕx=x−lnx−1 for x>0.
Based on the property of fractional derivatives given in [8], we get(28)DαL0≤1−S0SDαS+1r1DαN+1r2DαT+1r2r31−V0VDαV=1−S0SA−dS−β1SN−β2ST+1r1r1β1SN−d+ε1N+1r2r2β2ST−d+ε2T−β3TV+1r2r31−V0VB+r3β3TV−d+ε3V.
ByS0=A/d and V0=B/d+ε3, we obtain(29)DαL0≤dS01−SS01−S0S+β1S0N+β2S0T−d+ε1r1N−d+ε2r2T+d+ε3r2r3V01−VV01−V0V−β3r2TV0=−dSS−S02+β1S0−d+ε1r1N+β2S0−d+ε2r2−β3V0r2T−d+ε3r2r3V−V02V=−dSS−S02+d+ε1r1A1−1N+d+ε2r2A2−1−Bβ3d+ε2d+ε3T−d+ε3r2r3V−V02V.
ThenDαL0≤0 when A1≤1 and A2≤1+Bβ3/d+ε2d+ε3. Clearly, DαL0=0 if and only if S=S0, N=0, T=0, and V=V0. Then the largest invariant set contained in S,N,T,V|DαL0t=0 is the singleton E0. By LaSalle’s invariance principale [9], we deduce that E0 is globally asymptotically stable for A1≤1 and A2≤1+Bβ3/d+ε2d+ε3.
It remains to investigate the dynamical property ofE0 in case when A1>1 or A2>1+Bβ3/d+ε2d+ε3. For this purpose, we compute the characteristic equation at E0 that is given by(30)λ−λ1λ−λ2λ−λ3λ−λ4=0,where(31)λ1=−d,λ2=−d−ε3,λ3=d+ε2A2−1−Bβ3d+ε2d+ε3,λ4=d+ε1A1−1.
We haveλ1<0, λ2<0, λ3>0 if A2>1+Bβ3/d+ε2d+ε3, and λ4>0 if A1>1. Consequently, E0 is unstable if A1>1 or A2>1+Bβ3/d+ε2d+ε3.Theorem 4.
Suppose thatA1>1. Then the tumor-free equilibrium E1 is globally asymptotically stable if(32)A2≤A1+ABr1β1β3dd+ε1d+ε2d+ε3,and it is unstable if(33)A2>A1+ABr1β1β3dd+ε1d+ε2d+ε3.Proof.
Consider the following Lyapunov functional:(34)L1t=S1ϕStS1+1r1N1ϕNtN1+1r2Tt+1r2r3V1ϕVtV1.
Then(35)DαL1≤1−S1SA−dS−β1SN−β2ST+1r11−N1Nr1β1SN−d+ε1N+1r2r2β2ST−d+ε2T−β3TV+1r2r31−V1VB+r3β3TV−d+ε3V.
ByV1=B/d+ε3 and S1=d+ε1/r1β1, we obtain(36)DαL1≤dS11−S1S1−SS1+β1S1N12−S1S−SS1+β1S1−d+ε1r1N+β2d+ε1r1β1−d+ε2r2−β3Br2d+ε3T+Br2r32−V1V−VV1=−d+β1N1S−S12S−Br2r3V−V12VV1+dd+ε1d+ε2Ar1r2β1A2−A1−ABr1β1β3dd+ε1d+ε2d+ε3T.
ThenDαL1≤0 when(37)A2≤A1+ABr1β1β3dd+ε1d+ε2d+ε3.
Obviously,DαL1=0 if and only if S=S1, N=N1, T=0, and V=V1. Then the largest invariant set contained in S,N,T,V|DαL1t=0 is the singleton E1. By LaSalle’s invariance principle, we deduce that E1 is globally asymptotically stable for(38)A2≤A1+ABr1β1β3dd+ε1d+ε2d+ε3.
On the contrary, the characteristic equation atE1 is given by(39)d+ε1+λr2β2S1−d−ε2−β3V1−λfλ=0,where fλ=d+λ+β1N1d+ε1+λ−r1β1S1+r1β12N1S1.
One of the eigenvalues of (39) is(40)λ1=r2β2S1−d−ε2−β3V1=dd+ε1d+ε2Ar1β1A2−A1−ABr1β1β3dd+ε1d+ε2d+ε3.
We observe thatλ1>0 if(41)A2>A1+ABr1β1β3dd+ε1d+ε2d+ε3.
Thus,E1 is unstable when(42)A2>A1+ABr1β1β3dd+ε1d+ε2d+ε3.Theorem 5.
Suppose thatA2>1+Bβ3/d+ε2d+ε3 and A2/A1>1. Then the treatment failure equilibrium E2 is globally asymptotically stable if(43)1+β2d+ε3r3β3d≥A1+Bβ2r3dd+ε2A2/A1−1,and becomes unstable if(44)1+β2d+ε3r3β3d<A1+Bβ2r3dd+ε2A2/A1−1.Proof.
Consider the following Lyapunov functional:(45)L2t=S2ϕStS2+1r1Nt+1r2T2ϕTtT2+1r2r3V2ϕVtV2.
Then(46)DαL2≤dS21−S2S1−SS2+β2S2T22−S2S−SS2+β1S2−d+ε1r1N+β2S2−d+ε2r2−β3r2V2T+Br2r32−V2V−VV2=−d+β2T2S−S22S+β1S2−S3N−Br2r3V−V22VV2.
By computation, we find(47)S2−S3=Ar1β1r3β3d+ε2A2A1−1+Bβ2β3d+ε1−β2d+ε1d+ε2d+ε3A2A1−1−r3β3dd+ε1d+ε2A2A1−1.
Thus,S2−S3≤0 implies that(48)A1+Bβ2r3dd+ε2A2/A1−1≤1+β2d+ε3r3β3d.
Consequently,DαL2≤0 when(49)A1+Bβ2r3dd+ε2A2/A1−1≤1+β2d+ε3r3β3d.
Clearly,DαL2=0 if and only if S=S2, N=0, T=T2, and V=V2. Then the largest invariant set contained in S,N,T,V|DαL2t=0 is the singleton E2. By LaSalle’s invariance principle, we deduce that E2 is globally asymptotically stable for(50)1+β2d+ε3r3β3d≥A1+Bβ2r3dd+ε2A2/A1−1.
On the other side, the characteristic equation atE2 is given by(51)r1β1S2−d−ε1−λgλ=0,where(52)gλ=d+β2T2+λd+ε2+λ+β3V2−r2β2S2r3β3T2−d−ε3−λ−r3β32T2V2+r2β22S2T2r3β3T2−d−ε3−λ.
One of the eigenvalues of (51) is(53)λ2=r1β1S2−d−ε1=r1β1S2−S3.
We can observe from the proof of part (a) thatλ2>0 if(54)1+β2d+ε3r3β3d<A1+Bβ2r3dd+ε2A2/A1−1.
Thus,E2 is unstable when(55)1+β2d+ε3r3β3d<A1+Bβ2r3dd+ε2A2/A1−1.Theorem 6.
The partial success equilibriumE3 is globally asymptotically stable if(56)A2>A1+ABr1β1β3dd+ε1d+ε2d+ε3,A1+Bβ2r3dd+ε2A2/A1−1>1+β2d+ε3r3β3d.Proof.
Consider the following Lyapunov functional(57)L3t=S3ϕStS3+1r1N3ϕNtN3+1r2T3ϕTtT3+1r2r3V3ϕVtV3.
Then(58)DαL3≤1−S3SA−dS−β1SN−β2ST+1r11−N3Nr1β1SN−d+ε1N+1r21−T3Tr2β2ST−d+ε2T−β3TV+1r2r31−V3VB+r3β3TV−d+ε3V.
ByA=dS3+β2S3T3+β1S3N3, β3/r2T3V3=d+ε3/r2r3V3−B/r2r3, β1S3N3=d+ε1/r1N3, and β2S3T3=d+ε2/r2T3+β3/r2T3V3, we get(59)DαL3≤dS31−S3S1−SS3+β2S3T32−S3S−SS3+β1S3−d+ε1r1N+β2S3−d+ε2r2−β3r2V3T+β1S3N32−S3S−SS3+Br2r32−V3V−VV3=−d+β2T3+β1N3S−S32S−Br2r3V−V32VV3.
ThereforeDαL3≤0, with equality if and only if S=S3 and V=V3. By a simple computation, we show that DαL3=0 if and only if S=S3, N=N3, T=T3 and V=V3. It follows from LaSalle’s invariance principale that E3 is globally asymptotically stable under the conditions that this point exists. ■


## 4. Numerical Simulations

In this section, we give some numerical simulations to illustrate and validate our theoretical results, and we present some biological interpretations. We choose the time interval fromt=0 to t=400 with a step size Δt=0.1. We take A=0.02, d=0.02, B=0.01, r1=0.8, and r3=0.5. The parameters β1, β2, β3, ε1, ε2, ε3, and r2 of model (1) are taken as free parameters.First, we takeβ1=0.03, β2=0.03, β3=0.1, ε1=0.04, ε2=0.01, ε3=0.008, and r2=0.8. These values give A1=0.4, A2=0.8, and 1+Bβ3/d+ε2d+ε3=2.1905. Thus, A1<1 and A2<1+Bβ3/d+ε2d+ε3. According to Theorem 1, the equilibrium E01,0,0,0.3571 is globally asymptotically stable which consists with our numerical simulation in Figure 1. This may reflect an extreme competition between normal and tumor cells, leading to the extinction of normal cells and eradicating tumor cells by the M1 virus, giving rise to patient death.Figure 1
Stability of the competition-free equilibriumE0.
(a)(b)(c)(d)Next, we takeβ1=0.1, β2=0.03, β3=0.1, ε1=0.008, ε2=0.01, ε3=0.006, and r2=0.8. For this case, we obtain(60)A1=2.8571,A2=0.8,A1+ABr1β1β3dd+ε1d+ε2d+ε3=6.5201.Thus we get(61)A1>1,A2≤A1+ABr1β1β3dd+ε1d+ε2d+ε3.In agreement with Theorem4, the tumor-free equilibrium E10.35,0.3714,0,0.3846 is globally asymptotically stable as exhibited in Figure 2. In this situation, we notice that the M1 virotherapy successfully eliminates tumor cells, and then normal cells have been restored. Consequently, the patient’s health will be improved.Figure 2
Stability of the tumor-free equilibriumE1.
(a)(b)(c)(d)In Figure3, we assume that β1=0.03, β2=0.1, β3=0.1, ε1=0.04, ε2=0.008, ε3=0.008, and r2=0.8. Thus,(62)A2=2.8571>2.2755=1+Bβ3d+ε2d+ε3,A2A1=7.1429>1,1+β2d+ε3r3β3d=3.8≥0.9814=A1+Bβ2r3dd+ε2A2/A1−1.Figure 3
Stability of the treatment failure equilibriumE2.
(a)(b)(c)(d)We can see that the equilibriumE20.8313,0,0.0406,0.385 is globally asymptotically stable that agrees with our result in Theorem 5. Biologically, our treatment fails in eliminating the tumor cells as that normal cells are lost. Hence, the patient’s health is in danger.Finally, we chooseβ1=0.15, β2=0.35, β3=0.1, ε1=0.008, ε2=0.008, ε3=0.008, and r2=0.9. We get(63)A2=11.25>4.2857=A1,A2>9.7522=A1+ABr1β1β3dd+ε1d+ε2d+ε3,1+β2d+ε3r3β3d=10.8<11.978=A1+Bβ2r3dd+ε2A2/A1−1.In agreement with Theorem6, the equilibrium E30.2333,0.1571,0.1204,0.455 is globally asymptotically stable as shown in Figure 4. Here, our treatment partially reduces tumor cells and increases normal cells’ levels. However, the treatment cannot wholly eliminate tumor cells, but it can prolong the patient’s life.Figure 4
Stability of the partial success equilibriumE3.
(a)(b)(c)(d)


## 5. Conclusion

In this paper, we have studied the dynamics of an oncolytic M1 virotherapy model, considering the memory effect denoted by the Caputo fractional derivative. The well-posedness of the proposed model was proved through nonnegativity and boundedness of solutions. We found that the model has four possible equilibrium points, namely, the competition-free equilibriumE0, the tumor-free equilibrium E1, the treatment failure equilibrium E2, and the partial success equilibrium E3. By constructing suitable Lyapunov functionals, the global stability of E0 is determined by two threshold parameters that are the absorbing number of nutrients by normal cells A1 and the absorbing number of nutrients by tumor cells A2, when A1≤1 and A2≤1+Bβ3/d+ε2d+ε3, E0 is globally asymptotically stable, and these conditions determine when normal and tumor cells are lost, which may not be useful to test the viability of treatment. The tumor-free equilibrium E1 exists and is globally asymptotically stable if A1>1 and(64)A2≤A1+ABr1β1β3dd+ε1d+ε2d+ε3,and these conditions show that the M1 virus succeeds to eliminate the tumor, which is helpful in improving virotherapy. The treatment failure equilibrium E2 exists and is globally asymptotically stable if(65)A2>A1,A2>1+Bβ3d+ε2d+ε3,1+β2d+ε3r3β3d≥A1+Bβ2r3dd+ε2A2/A1−1,and these conditions refer to the failure of the treatment, as indicated by his name. The partial success equilibrium E3 exists and is globally asymptotically stable if(66)A2>A1+ABr1β1β3dd+ε1d+ε2d+ε3,1+β2d+ε3r3β3d<A1+Bβ2r3dd+ε2A2/A1−1.These results indicate the partial success of M1 virus in decreasing tumor cells and increasing normal cells, which can reduce the tumor’s size and stabilize the disease progression.From the above analytical results, we remark that the Caputo fractional derivative’s memory does not affect the stability analysis of equilibria. Based on the numerical simulations, we observe that the fractional order affects the speed of convergence and the time for arriving to equilibria (Figures1–4).The results obtained in this study are based on the fractional derivative in sense of Caputo with singular kernel. It will be more interesting to model the dynamics of oncolytic M1 virotherapy by using the new generalized fractional derivative with nonsingular kernel [10]. Moreover, we will extend our model presented in (1) by taking into account other biological factors such as diffusion [11, 12] and immunity [13, 14].



---

*Source: 9934070-2021-06-07.xml*