Article Search
eISSN 0454-8124
pISSN 1225-6951

### Article

Kyungpook Mathematical Journal 2022; 62(1): 119-132

Published online March 31, 2022 https://doi.org/10.5666/KMJ.2022.62.1.119

### Hopf-bifurcation Analysis of a Delayed Model for the Treatment of Cancer using Virotherapy

Maharajan Rajalakshmi and Mini Ghosh*

Division of Mathematics, School of Advanced Sciences, Vellore Institute of Technology, Chennai Campus, Chennai-600127, India
e-mail : rajalakshmi.m2016@vitstudent.ac.in and minighosh@vit.ac.in

Received: February 23, 2020; Revised: October 7, 2020; Accepted: November 16, 2020

Virotherapy is an effective method for the treatment of cancer. The oncolytic virus specifically infects the lyse cancer cell without harming normal cells. There is a time delay between the time of interaction of the virus with the tumor cells and the time when the tumor cells become infectious and produce new virus particles. Several types of viruses are used in virotherapy and the delay varies with the type of virus. This delay can play an important role in the success of virotherapy. Our present study is to explore the impact of this delay in cancer virotherapy through a mathematical model based on delay differential equations. The partial success of virotherapy is guarenteed when one gets a stable non-trivial equilibrium with a low level of tumor cells. There exits Hopf-bifurcation by considering the delay as bifurcation parameter. We have estimated the length of delay which preserves the stability of the non-trivial equilibrium point. So when the delay is less than a threshold value, we can predict partial success of virotherapy for suitable sets of parameters. Here numerical simulations are also performed to support the analytical findings.

Keywords: Cancer, Virotherapy, Delay model, Hopf-bifurcation

Cancer is a disease which is caused by unusual cell growth. It is one of the leading cause of death in the world. There are several types of cancer, and they are generally described by the body part they originate in. Most cancers are curable if detected early enough. Cancer is a very complex disease to treat and treatment of cancer varies. The choice of treatment depends on the location of the tumor, size of the tumor, cell type and most importantly the overall health of the patient. For the small size malignant tumors, surgery is often recommended. Chemotherapy, radiotherapy, immunotherapy, gene therapy, virotherapy, hormone therapy etc. are treatment methods which are used alone or in combination depending upon the size of the cancer and patient's health. Mathematical modeling is an important tool which is being used in treatment of several diseases including cancer. As our human body is very complex, the success of any therapy not only depends upon the types of treatment but also at the duration of treatment, frequency of treatment and effectiveness of a specific therapy. Mathematical modeling and simulation help in designing suitable personalized therapy for a patient. There are several research findings based on mathematical models on cancer growth and treatment of cancer [19, 21, 9, 6, 2, 22, 18, 12, 11]. In [19], the authors formulated and analyzed a mathematical model for a malignant tumor system by considering prey-predator type interactions. Here the authors investigated both deterministic and stochastic models and obtained thresholds which play a key role in the control of malignant tumor growth. In [21], the author investigated a mathematical model for tumor growth with quiescence by incorporating delay. In [9], the authors considered a kinetic theory approach to model the growth of tumors and macrophages heterogeneity and plasticity by incorporating diffusion. A mathematical model for tumor-immune interaction with delay is formulated and analyzed in [2]. Here the authors have estimated the length of the delay parameter which preserve the stability of the system. In [22], the authors presented a competition model for the growth of a malignant tumor by incorporating immune response. Here the authors demonstrated the presence of a Hopf-bifurcation. A delayed model for tumor-immune response with chemoimmunotherapy and optimal control is investigated in [18]. Here the authors have emphasized the importance of combination therapy and optimal control. In [12], the authors described fundamentals of modeling of tumor growth and discussed different approaches to model tumor growth. Modeling tumor-immune dynamics with treatment is investigated in [11]. Mathematical modeling of the treatment of cancer using virotherapy by considering a system of ordinary differential equations is reported in [1, 4, 7, 20, 15, 17]. In [16], the authors have formulated and analyzed a delay differential equation model for the treatment of cancer using virotherapy. Here the authors have considered the delay in tumor cells to become infectious after getting infected with viruses. This work is an extension of the author's own work in [15] where the ordinary differential equation model was formulated and analyzed. In both papers, the authors considered logistic growth of the tumor cells whereas a generalized logistic growth model is more suitable to describe the dynamics of the tumor cells. Keeping this in mind in [17], the authors formulated and analyzeda mathematical model for the treatment of cancer using virotherapy by considering generalized logistic growth of the tumor cells. Some existing models are derived using the framework of population dynamics and others include space dynamics which is an important feature. With these two approaches of population dynamics and system dynamics, existing mathematical models can be further extended by including space dynamics through diffusion [10] as it shows a key role in delaying the dynamics, or through reaction-diffusion [3] as it relates to the dynamic linked with angiogenesis.

The present work is an extension of the work in [17] which does not incorporate latent delay. It is a proven fact that delay is capable of changing the whole dynamics of the system under consideration [21, 14, 13]. The remaining part of this paper is organized as follows: Section 2 describes the proposed mathematical model; Section 3 discusses the reproduction numbers and equilibria of the model; Section 4 demonstrates the stability analysis of the non-trivial equilibrium point of the model. Section 5 demonstrates the Hopf-bifurcation analysis. Section 6 exhibits the numerical simulation results. Finally, we summarize our results in Section 7.

As delay plays an important role in the dynamics of cancer growth and treatment, here we extend the model by Rajalakshmi and Ghosh [17] by incorporating latent delay. It is assumed that after the interaction of tumor cells with viruses, tumor cells get infected and produce new viruses. We refer to the time between interaction and production of new viruses as 'delay', The tumor cells are divided into two disjoint classes: uninfected tumor cells with population y(t), and infected tumor cells with population x(t). The population of free virus particles is denoted by v(t) and CTLs by z(t). Here we incorporate the delay (τ) that describes the time between virus attack on uninfected tumor cells and the production of new virus particles by the virus infected tumor cells. The growth of uninfected tumor cells follows generalized logistic growth. Keeping in view the above fact, we formulate our delay model as follows:

dydt=ry1 y+xK ϵβy(tτ)v(tτ)ρxyν1yz,dxdt=βy(tτ)v(tτ)δxν2xz,dvdt=αxωvβyv,dzdt=γ1yz+γ2xzμz.

Let C([τ,0],R+4) denote the Banach space of continuous functions mapping the interval [τ,0] into R+4 where R+4={(x1,x2,x3,x4):xi0, i=1,2,3,4}. The initial conditions for the system (2.1) are given by y(θ),x(θ),v(θ),z(θ)C([τ,0],R+4).

Here r is the intrinsic growth rate of uninfected tumor cells, K is the carrying capacity of uninfected tumor cells, ϵ is the parameter which determines the shape of the sigmoidal growth curve in the generalized logistic growth of uninfected tumor cells, and βxy is the rate at which uninfected tumor cells are becoming infected due to interaction with virus particles. The term ρxy corresponds to the rate at which uninfected tumor cells become infected through fusion with infected cells either single or in syncytia. The uninfected tumor cells and infected tumor cells are removed by CTLs at the rates ν1yz and ν2xz respectively. The CTL population increases at rate γ1yz+γ2xz and is removed at the rate μz. δx denotes the death rate of infected tumor cells and αx corresponds to the rate of release of virions by infected tumor cells. The term ωv corresponds to to the rate of elimination of virus particles due to various causes which includes non-specific bindings and defective particles. The term βyv in the third equation of this model is based on the assumption that by infecting tumor cells, the virus particles will be lost. We can note here that although the value of ϵ can be less than one or greater than one, in practice it is taken as the order of 1 i.e. 0. As K is the carrying capacity, x+y is at most K.

### 3. The Reproduction Number and Equilibria of the Model

As delay does not change the equilibria of the system, we get, for this delay model, the same equilibria as obtained for the model discussed in [17]. The basic reproduction number R0 and the immune response reproduction number R1 are computed using the next generation matrix method as described in [5, 8]. First we consider the basic reproduction R0 corresponding to virus infection. The virus-free equilibrium or a therapy failure equilibrium is given by E3y^,0,0,z^ with y^=μγ1 and z^=r1y^/Kϵ/ν1. Following the same notations as in [5, 8], we find the matrix F and V as follows:

F=βyv00 and V=δx+ν2xzαx+ωv+βyvγ1yz+γ2xzμz.

F is the Jacobian of F at therapy failure equilibrium E3=0βy¯0000000 and V is the Jacobian of V at therapy failure equilibrium E3=δ+ν2z¯00αω+βy¯0γ2z¯0μ.

The largest eigenvalue of FV1 is called the basic reproduction number {R0} and is obtained as follows:

R0=βαy¯(δ+ν2z¯)(ω+βy¯)=βαμγ1δ+ν2 r1μ γ1K ϵ ν1 ω+βμγ1 .

Next we compute the immune response reproduction number R1, (say) using the same next generation matrix method by considering all the variables y, x,v, z. Here it can be noted that the CTL-free equilibrium is given by E4y*,x*,v*,0. Here we find the matrix F and V as follows:

F=000γ1yz+γ2xz and V=ry1y+xKϵ+βyv+ρxy+ν1yzβyv+δx+ν2xzαx+ωv+βyvμz.

F is the Jacobian of F at E4=000000000000000γ1y*+γ2x* and V is the Jacobian of V at E4=v11v12v13v14βv*δβy*ν2x*βv*αω+βy*0000μ,

where

v11=r1y*+x*Kϵ+ry*ϵy* +x* Kϵ11K+βv*+ρx*,
v12=ry*ϵy*+x*Kϵ11K,v13=βy*, and v14=ν1y*.

Now the largest eigenvalue of FV-1 corresponds to immune response reproduction number R1 and is given by

R1=γ1y*+γ2x*μ=γ1ωδβ(αδ)+γ2x*μ,

where 0<x*<Kωδβ(αδ). The model system (2.1) has five equilibria as follows:

(i) E1(0,0,0,0), a trivial equilibrium,

(ii) E2(K,0,0,0), a boundary equilibrium

(iii) E3y^,0,0,z^ with y^=μ/γ1 and z^=r(1(y^/K)ϵ)/ν1, a therapy failure equilibrium,

(iv) E4y*,x*,v*,0 with y*=ωδ/(β(αδ)), a CTL-free equilibrium and (v) E5y¯,x¯,v¯,z¯, an endemic equilibrium which may corresponds to partial success of therapy. Here it is assumed that the parameter α is greater than the parameter α which is reasonable. Here v*=(αδ)x*/ω and x* is the positive root of the following equation:

g(x)=rωδβ(αδ)1ωδ β(αδ) +xK ϵδxρωδxβ(αδ)=0,

which exists for K>ωδβ(αδ),

x¯=μγ1y¯γ2,v¯=α(μγ1y¯)γ2(ω+βy¯),z¯=(αδ)(ω+βy¯)ωαν2(ω+βy¯).

It is easy to observe that for v¯ and z¯ to be positive we need the following condition:

ωδβ(αδ)<y¯<μγ1,

where y¯ is the positive root of the following equation:

h(y)=r1(γ2γ1)y+μγ2Kϵβαμγ1yγ2(ω+βy)ρμγ1yγ2  ν1(αδ)(ω+βy)αων2(ω+βy)=0.

For more details about the existence of these equilibria one can refer [17].

Our aim is to investigate the success or failure of virotherapy. Hence we need to concentrate on the stability of the partial success equilibrium point E5. One time injection of oncolytic viruses may not kill all tumor cells, so we need to see whether it can cause significant reduction in the tumor cells or not. Additionally, we should understand that whether the partial success equilibrium is stable or not. In the next section we shall investigate the stability of partial success equilibrium E5 in presence of delay.

We linearize the system (2.1) about the equilibrium point E5 as follows:

S(t)=AS(t)+BS(tτ),

where S(t)=(y(t),x(t),v(t),z(t))T and

A=a11¯a12¯0ν1y¯0(δ+ν2z¯)0ν2x¯βv¯α(ω+βv¯)0γ1z¯γ2z¯0γ1y¯+γ2x¯μ,
B=βv¯0βy¯0βv¯0βy¯000000000,

a11=r1x¯+y¯Kϵry¯ϵKx¯ +y¯ K ϵ1ρx¯ν1z¯,  a12=ry¯ϵKx¯ +y¯ K ϵ1ρy¯. Here the matrices A and B are computed at the equilibrium under consideration i.e. E5. The stability is determined by computing the roots of the following characteristic equation

det(A+BeλtλI)=0.

The characteristic equation can be written as,

P(λ)+eλτQ(λ)=0,
(λ4+λ3M1+λ2M2+λM3+M4)+eλτ(N1λ3+N2λ2+N3λ+N4)=0

where

P(λ)=λ4+λ3M1+λ2M2+λM3+M4,  Q(λ)=λ3N1+λ2N2+λN3+N4,

and

M1=a11¯+ν2z¯+ω+βy¯M2=δa11¯a11¯ν2z¯a11¯ω+δω+δβy¯+βy¯z¯ν2+βy¯+γ2ν2x¯z¯γ1ν1y¯z¯M3=δωa11¯δβa11¯y¯ων2a11¯z¯βy¯a11¯z¯ν2a11¯βy¯a11¯ν2x¯γ2z¯+ν2γ2ωx¯z¯  +ν2βy¯γ2z¯x¯+a12¯ν2γ1z¯x¯δγ1ν1y¯z¯γ1ν1ν2y¯ z ¯ 2ωγ1ν1y¯z¯βγ1ν1z¯ y ¯ 2M4=ν2a11¯z¯γ2ων2a11¯γ2βy¯x¯z¯+a12¯ωγ1ν2x¯z¯+βν2γ1z¯x¯y¯a12¯+δων1γ1y¯z¯  +δβν1γ1z¯ y ¯ 2+ωγ1ν1ν2y¯ y ¯ 2+βγ1ν1ν2 y ¯ 2 z ¯ 2N1=βv¯N2=δβv¯+ν2βz¯v¯+ωβv¯+βδv¯+αβy¯βv¯a12¯+β2¯yv¯N3=a11¯αβy¯+δωβv¯+δβ2y¯v¯+βων2v¯z¯+β2v¯y¯z¯ν2+β2y¯v¯+βγ2ν2v¯x¯z¯  βa12¯ωv¯β2y¯v¯a12¯a12¯β2v¯y¯+β2v¯δy¯+β2ν2z¯v¯y¯+βν1γ2z¯y¯N4=ν2βv¯x¯γ2z¯ω+ν2x¯z¯γ2β2y¯v¯γ2β2ν2x¯y¯v¯z¯αγ1βx¯y¯z¯βωγ2v¯z¯y¯ν1  β2ν1γ2z¯v¯ y ¯ 2+γ2ν1β2v¯ y ¯ 2z¯+αβν1γ1z¯ y ¯ 2.

When τ=0 we get

λ4+M1+N1λ3+M2+N2λ2+M3+N3λ+M4+N4=0.

Using Routh-Hurwitz criteria, all the roots of the above biquadratic equation will have negative real parts provided following conditions hold:

M1+N1>0,M4+N4>0,M1+N1M2+N2M3+N3>0,
M1+N1M2+N2M3+N3M1+N12M4+N4M3+N32>0.

Hence in absence of delay the equilibrium point E5 is locally asymptotically stable provided above mentioned conditions are satisfied.

Now let us discuss the case when delay is not zero. When τ0, let us assume λ=iω, (ω>0) is the purely imaginary root of the equation (4.3). Substituting λ=iω in the equation (4.3) and equating the real and imaginary parts we get the following two equations:

ω4ω2M2+M4=cosωτω2N2N4+sinωτω3N1ωN3
ωM3ω3M1=cosωτω3N1ωN3sinωτω2N2N4

Now squaring and adding equations (5.1) and (5.2), we get the following trancendental equation:

ω8+ω6M122M2N12+ω4M22+2M42M1M3+2N1N3N22        +ω2M322M2M4+2N2N4N32+M42N42=0

Assuming ω2=c we have

L(c)=c4+a1c3+a2c2+a3c+a4=0,

where

a1=M122M2N12,a2=M22+2M42M1M3+2N1N3N22,a3=M322M2M4+2N2N4N32,a4=M42N42.

Now we have the following theorems based on the roots of the equation L(c)=0.

Theorem 5.1. If the coefficients a1,a2,a3,a4 in L(c) satisfy the conditions of Routh-Hurwitz criterion (i.e. a1>0,a4>0,a1a2a3>0 and a1a2a3a12a4a32>0), then the interior equilibrium point E5(y¯,x¯,v¯,z¯) if it exists is asymptotically stable for all delay τ>0 provided it is stable in absence of delay.

Theorem 5.2. If the coefficients a1,a2,a3,a4 in L(c) satisfy Routh-Hurwitz criterion (i.e. a1>0,a4>0,a1a2a3>0 and a1a2a3a12a4a32>0), and the interior equilibrium point E5(y¯,x¯,v¯,z¯) is unstable at τ=0, then it will remain unstable for all τ0.

If L(c)=0 has a positive root, then we can have the following result.

Theorem 5.3. The endemic equilibrium point E5 of the system (2.1) is conditionally stable if and only if all the roots of the characteristic equation (4.3) have negative real parts at τ=0 and there exist some positive value of the delay τ such that the characteristic equation (4.3) has a pair of purely imaginary roots ±iω0 (say). The system will undergo a stability change for an infinite number of values of τ say τn*, where

τn*=1ω0cos1(ω0M3ω03M1P2)(ω03N1ω0N3)+(ω02N2N4)(ω04ω02M2+M4)(ω02N2N4)2+(ω03N1ω0N3)2        +2nπω0,  n=0,1,2,

Here it can be noted that if the equilibrium point E5 is locally asymptotically stable in absence of delay then it remains stable for delay τ<τ0* with n=0. At delay τ=τ0* Hopf-bifurcation occurs and system undergoes periodic oscillations.

### 5.1. Analysis of Hopf-bifurcation

Now, we shall investigate the Hopf bifurcation of the model system (2.1), for which we need to verify the transversality condition d(Reλ)dτ|τ=τ0>0 i.e. dξdτ|ξ=0>0 for λ(τ)=ξ(τ)+iω(τ). This will designates that there exists at least one eigenvalue with nonnegative real part for τ>τs. Also, the conditions for Hopf bifurcation are necessary to prove the existence of periodic solutions. Firstly, we are interested for purely complex roots λ=iω0 of (4.3). Equation (4.3) implies that |P(iω0)|=|Q(iω0)| and this determines the possible set of values of ω0. Now our goal is to observe the direction of motion of λ as τ is varied, for which we find, Ω=sign[d(Reλ)dτ]|τ=τn=sign[d(Reλ)dτ]1|τ=τn.

On differentiating (4.3) with respect to τ, we get

(3λ3+2λ2M1+2λM2+M3)+eλτ(3λ2N1+2λN2+N3)τeλτλ3N1+λ2N2          +λN3+N4dλdτλeλτ(λ3N1+λ2N2+λN3+N4)=0

dλdt1=eλτ(3λ3+2λ2M1+2λM2+M3)λ(λ3N1+λ2N2+λN3+N4)+(3λ2N1+2λN2+N3)λ(λ3N1+λ2N2+λN3+N4)τλ.

Since λ(τ0)=iω0 is a simple root of the characteristic equation (4.3), we can evaluate the

expressions involved in the above derivative at τ=τ0 as follows:

eλτ(3λ3+2λ2M1+2λM2+M3)|τ=τ0=δ1+iδ2,λ(λ3N1+λ2N2+λN3+N4)|τ=τ0=δ3+iδ4,(3λ2N1+2λN2+N3)|τ=τ0=δ5+iδ6,

where

δ1=(M32ω2M1)cosω0τ0+(3ω32ωM2)sinω0τ0,δ2=(2ωM23ω3)cosω0τ0+(M32ω2M1)sinω0τ0,δ3=ω02(ω2N1N3),δ4=ω(N4ω2N2)δ5=N33ω2N1δ6=2ωN2.

Now

dλdτ1|τ=τ0=ddτReλ(τ0)1=δ1δ3+δ2δ4+δ5δ3+δ4δ6δ32+δ42.

Using the equations in (5.1) and (5.2), we can rewrite above expression as follows:

dλdτ1|τ=τ0=ω02δ32+δ424ω06+3(M122M2N12)ω04+2 M22+2M42M1M3      +2N1N3N22 ω02+(M322M2M4+2N2N4N32)    =ω02δ32+δ42(4c3+3a1c2+2a2c+a3)|c=ω02=ω02δ32+δ42L(c)|c=ω02

Therefore

signddτReλ(τ0)=signddτReλ(τ0)1=signω02δ32+δ42L(c)|c=ω02.

As δ32+δ42>0, ω02>0 and L(c)|c=ω020, the signddτReλ(τ0) will be determined by the signL(c)|c=ω02.

We already have Re(λ(τ))=ξ(τ) and ξ(τ0)=0. Therefore if signL(c)|u=ω02<0, then there exists a ζ>0 such that ξ(τ) is decreasing in (τ0ζ,τ0) and ξ(τ)=0 at τ=τ0. Hence for all τ(τ0ζ,τ0), ξ(τ)>0, which contradicts the fact that roots of the characteristic equation (4.3) have negative real parts for all τ[0,τ0] and τ=τ0 is the minimum value of delay τ for which (4.3) will have purely imaginary roots. Hence signL(c)|c=ω02>0 which shows that there exists at least one λ(τ) with ξ(τ)>0 for τ>τ0. Therefore, the transversality condition holds successfully and the system undergoes Hopf bifurcation at ξ=ξ0,τ=τ0.

Here first we simulate our model system (2.1) for the following set of parameters without delay i.e. when τ=0:

r=0.2062,m=2139.258,ϵ=1.649,β=0.001305,ρ=0.005,ν1=0.00673,
ν2=0.10095,δ=0.25,α=4,ω=0.285,γ1=0.00015,γ2=0.027,μ=0.2.

For this set of parameters we get our non-trivial equilibrium E5 as (195.5, 6.3, 46.8, 16.2) and R0=1.78. And this equilibrium is locally asymptotically stable. This fact is shown in Figure 1. Next we consider our delay parameter τ as 1 and we find that the system tends to the same equilibrium point E5. This is shown in Figure 2. Further, we change the delay parameter τ as 2. In this case we observe that system shows damped oscillations but finally it converges to the equilibrium point E5. The time taken to converge the equilibrium point E5 increases with the increase in the delay parameter τ. This is demonstrated in Figure 3. Now when we take τ=2.0465, the model system (2.1) undergoes to Hopf-bifurcation and shows periodic oscillations. This is demonstrated in Figure 4. Here it can be noted that the critical value of this delay τ=2.0465 is in fact τ0* below which the non-trivial equilibrium point E5 is stable. This critical value of delay is very much dependent on other model parameters too. It is observed that increase in β causes the decrease in the threshold value of delay i.e. τ0*.

Figure 1. Variation of state variables with time with delay τ=0.

Figure 2. Variation of state variables with time with delay τ=1.0.

Figure 3. Variation of state variables with time with delay τ=2.0.

Figure 4. Variation of state variables with time with delay τ=2.0465.

This paper aims to study the impact of delay on the cancer virotherapy. Here we formulate a mathematical model by incorporating the time gap between the time of interaction of virus with tumor cells and the time when the tumor cells become infectious to produce new virus particles. We analyze the model by keeping in mind the partial success of virotherapy. There exists a threshold value of delay below which the equilibrium point corresponding to partial success of virotherapy is stable. This threshold value is derived analytically and computed numerically. It is also observed that this threshold value changes with change in the other parameters such as β and δ etc. We also verify the transversality condition for Hopf-bifurcation. Numerical simulation is performed to support the analytical findings.

1. Z. Bajzer, T. Carr, K. Josic, S. J. Russell and D. Dingli, Modeling of cancer virotherapy with recombinant measles viruses, J. Theoret. Biol., 252(2008), 109-122.
2. S. Banerjee and R. R. Sarkar, Delay-induced model for tumorimmune interaction and control of malignant tumor growth, Biosystems, 91(1)(2008), 268-288.
3. N. Bellomo, K. Painter, Y. Tao and M. Winkler, Occurrence vs. absence of taxis-driven instabilities in a May-Nowak model for virus infection, SIAM J. Appl. Math., 79(5)(2019), 1990-2010.
4. M. Biesecker, J. Kimn, H. Lu, D. Dingli and Z. Bajzer, Optimization of Virotherapy for Cancer, Bull. Math. Biol., 72(2010), 469-489.
5. C. Castillo-Chavez, Z. Feng and W. Huang, On the computation of R0 and its role on global stability, in: Mathematical Approaches for For Emerging and Reemerging Infectious Diseases, Springer-Verlag, (2002), 229-250.
6. M. A. J. Chaplain. Modelling aspects of cancer growth: insight from mathematical and numerical analysis and computational simulation, in Multiscale Problems in the Life Sciences. Vol. 1940 of Lecture Notes in Mathematics. Springer; 2008:147-200.
7. J. J. Crivellia, J. Földes, P. S. Kim and J. R. Wares, A mathematical model for cell cycle-specific cancer virotherapy, J. Biol. Dyn., 6(1)(2012), 104-120.
8. P. V. Driessche and J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission Mathematical Biosciences, Mathematical Biosciences, 180(2002), 29-48.
9. R. Eftimie and L. Gibelli, A kinetic theory approach for modelling tumour and macrophages heterogeneity and plasticity during cancer progression, Math. Mod. Meth. App. Sci., 30(2020), 659-683.
10. A. M. Elaiw and E. D. Al Agha, Analysis of a delayed and diffusive oncolytic M1 virotherapy model with immune response, Nonlinear Anal. Real World Appl., 55(2020), 103-116.
11. A. Eladdadi, L. D. Pillis and P. Kim, odelling tumour-immune dynamics, disease progression and treatment, Letters in Biomathematics, 5(2018), S1-S5.
12. H. Enderling, A. J. Mark and Chaplain, Mathematical Modeling of Tumor Growth and Treatment, Current Pharmaceutical Design, 20(30)(2014), 4934-4940.
13. S. Khajanchi, J. J. Nieto and Mathematical modeling of tumor-immune competitive system, considering the role of time delay, Appl. Math. Comput., 340(2019), 180-205.
14. S. Khajanchi and S. Banerjee, Stability and bifurcation analysis of delay induced tumor immune interaction model, Appl. Math. Comput., 248(2014), 652-671.
15. K. S. Kim, S. Kim and I. H. Jung, Dynamics of tumor virotherapy: A deterministic and stochastic model approach, Stoch. Anal. Appl., 34(3)(2016), 483-495.
16. K. S. Kim, S. Kim and I. H. Jung, Hopf bifurcation analysis and optimal control of Treatment in a delayed oncolytic virus dynamics, Math. Comput. Simulation, 149(2018), 1-16.
17. M. Rajalakshmi and M. Ghosh, Modeling treatment of cancer using virotherapy with generalized logistic growth of tumor cells, Stoch. Anal. Appl., 36(6)(2018), 1068-1086.
18. F. A. Rihan, D. H. Abdelrahman, F. Al-Maskari, F. Ibrahim and M. A. Abdeen, Delay differential model for tumour-immune response with chemoimmunotherapy and optimal control, Comput. Math. Methods Med., 2014(2014).
19. R. R. Sarkar and S. Banerjee, Cancer self remission and tumor stability-a stochastic approach, Math. Biosci., 196(2005), 65-81.
20. Z. Wang, Z. Guo and H. Peng, A mathematical model verifying potent oncolytic efficacy of M1 virus, Math. Biosci., 276(2016), 19-27.
21. R. Yafia, Dynamics analysis and limit cycle in a delayed model for tumor growth with quiescence, Nonlinear Anal. Model. Control, 11(1)(2006), 95-110.
22. R. Yafia, A study of differential equation modeling malignant tumor cells in competition with immune system, Int. J. Biomath., 4(2)(2011), 185-206.