# Stability and Hopf Bifurcation Analysis of an Oncolytic Virus Infection Model with Two Time Delays and Saturation Incidence

Authors: Xia Liu; Zhixing Hu
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In this paper, we study a model of oncolytic virus infection with two time delays, one of which is the time from the entry of viruses into tumor cells to start gene replication, and the other is the time from the entry of viruses into tumor cells to release new virus particles by infected tumor cells. In previous studies on oncolytic virus infection models, the infection rate was linear. Combined with the virus infection models, the saturated infection rate,βTV/1+qV is further considered to describe the dynamic evolution between viruses and tumor cells more objectively so as to further study the therapeutic effect of oncolytic viruses. This paper discusses the dynamics of the system under three conditions: (1) τ1=τ2=0, (2) τ1=0 and τ2>0, and (3) τ1>0 and τ2>0, and proves the global stability and local stability of the virusfree equilibrium, the stability of the infection equilibrium, and the existence of Hopf bifurcation. Finally, the conclusions of the paper are verified by MATLAB numerical simulations.


---

## Body

## 1. Introduction

An oncolytic virus is a natural or genetically modified virus that can specifically infect and kill tumor cells without damaging healthy cells (see [1]). Oncolytic virotherapy is a very promising cancer treatment method. Initially, The Lancet reported that influenza viruses can cause tumors in patients to regress, resulting in the emerging concept of oncolytic viruses. Since then, more than 160 different oncolytic viruses have been in preclinical studies and clinical trials (see [2]). According to current development trends, oncolytic viruses can be positioned as the next major breakthrough in cancer treatment after checkpoint inhibitors in immunotherapy.Many researchers in the fields of medicine, biology, and mathematics have been exploring and improving the related problems of oncolytic virotherapy, which has brought breakthroughs in the treatment of some malignant tumors, but there are still many problems to be solved at present (see [3]). For example, the proliferation of oncolytic viruses in host solid tumor cells is limited and the difference in tumor cell receptor expression limits the spread of the virus, which in turn affects the therapeutic effect. Immune-mediated pharmacokinetic therapy is relatively slow compared to drugs that directly kill tumor cells, so time delays need to be considered in clinical trial design and outcome evaluation. These problems have affected the antitumor effect of oncolytic viruses and failed to achieve the goal of systematic and reliable tumor elimination of oncolytic viruses. In the research process of oncolytic virotherapy, mathematical models can better understand the dynamic mechanism of tumor cells-virus-immune cells (Cancer-Virus-CTLs). It also provides an effective tool to predict the outcome of immunotherapy and optimize tumor treatment strategies.In order to describe the dynamic changes of the immune response in viral therapy, a series of ordinary differential equation models have been constructed in the past decades (see [4–11]. Dingli et al. published a model based on population dynamics in the mathematical biosciences (see [12]), which focused on the basic elements of radiation viroid therapy and discussed and proved the stability of the equilibrium of complete cure, local treatment, and treatment failure. On the basis of D. Dingli et al., Kim et al. proposed a virus model that considered both virus-specific immunity and tumor-specific immunity (see [13]). In addition, since the life cycle of viruses contains many intracellular processes, the time from the entry of viruses into tumor cells to the start of gene replication and the time of the release of new virus particles by the infected tumor cells are considered to further establish mathematical models with time delays. Wang et al. (see [14]) considered the cyclic period of viruses in cells and established a time-delay differential equation. Kim et al. (see [15]) proposed a delayed ordinary differential equation model. Li et al. (see [16]) further considered the time delay of virus self-replication based on the research of Kim et al. proposed an ordinary differential equation model with two time delays to analyze the influence of the time delay on the stability of the equilibrium.For most models of oncolytic virus infection, the infection rate of oncolytic viruses and uninfected tumor cells is bilinear, but the actual incidence is not strictly linear. Therefore, we further assume that the model has a saturated infection rateβTVa/1+qVba>0,b>0,q≥0. In this paper, we study a virus infection model with a saturated infection rate βTV/1+qVq>0.The rest of this paper is organized as follows: the second section builds the model. The third section solves the equilibrium of the system. The fourth section discusses the stability of the virusfree equilibrium and the virus infection equilibrium, as well as the existence conditions of Hopf bifurcation of the virus infection equilibrium, and explains the dynamic properties under three conditions: (1)τ1=τ2=0, (2) τ1=0 and τ2>0, and (3) τ1>0 and τ2>0. In the fifth section, numerical simulations are performed to verify the obtained results. Finally, the conclusions of this article are given.


## 2. Mathematical Model Formulation

On the basis of the research (see [16]), this chapter considers the impact of saturated infection rate and two time delays on the virus, and the following model is shown:(1)dTtdt=k1Tt1−bTt+It−βTt−τ1Vt−τ11+qVt−τ1−c1EtTt,dItdt=βTt−τ2Vt−τ21+qVt−τ2−c2EtIt−d1It,dVtdt=k2It−d2Vt,dEtdt=k3TtEt+k4ItEt−d3Et,where Tt,It,Vt, and Et represent the number of uninfected tumor cells, infected tumor cells, oncolytic viruses, and CTL immune cells, respectively. k1 represents the maximum growth rate of uninfected tumor cells. k2 represents the average rate at which the infected tumor cells release viruses. k3TtEt and k4ItEt represent the rate of CTL immune cells regeneration. 1/b represents the maximum carrying capacity of the human environment to tumor cells. τ1 represents the delay from viruses entering tumor cells to the beginning of gene replication. τ2 indicates the time from the entry of viruses into tumor cells to the release of new virus particles by infected tumor cells. Obviously, there is τ1≤τ2. c1 and c2 represent the rate of CTL immune cells clearing uninfected tumor cells and infected tumor cells, respectively. d1,d2, and d3 indicate the natural mortality of infected tumor cells, oncolytic viruses, and immune cells, respectively. The rupture of tumor cells releases tumor antigens, thereby inducing a systemic antitumor immune response, which may enhance the cytolytic activity of viruses.The initial condition of system (1) is(2)Tt=φ1θ,It=φ2θ,Vt=φ3θ,Et=φ4θ,θ∈−τ,0,where φiθ≥0,i=1,2,3,4. φ1θ,φ2θ,φ3θ,φ4θ∈C−τ,0,R+4. Note C+=C−τ,0,R+4, define the norm φ=supφ1θ,φ2θ,φ3θ,φ4θ, where φ:−τ,0⟶R+4 is a continuous function, C+ constitutes a Banach space.


## 3. The Existence of the Equilibrium

Through simple calculations,T≤1/b can be obtained from the first equation of system (1). The virusfree equilibrium of system (1) is Q0T0,0,0,E0, where T0=d3/k3 and E0=k11−bT0/c1>0.Let the infection equilibrium of system (1) be Q∗T∗,I∗,V∗,E∗, which is the positive solution of the following equations. Based on the model, there are reasonable assumptions as follows: k4>k3, k3>bd3 and k4/c2>k3/c1. All the above parameters are positive constants.(3)k1T∗1−bT∗+I∗−βT∗V∗1+qV∗−c1E∗T∗=0,βT∗V∗1+qV∗−c2E∗I∗−d1I∗=0,k2I∗−d2V∗=0,k3T∗E∗+k4I∗E∗−d3E∗=0.Through simple calculations, the following conclusions are obtained:(4)I∗=d2k2V∗,(5)T∗=c2E∗I∗+d1I∗1+qV∗βV∗,(6)E∗=bd3k1k21+qV∗+βV∗R1−1c11+qV∗,where(7)R1=k11+qV∗+bd2k1k4−k31+qV∗V∗bd3k1k21+qV∗+βV∗.WhenR1>1, we get E∗>0. Substituting (4)–(6) into the first equation of (3), we obtain the quadratic equation of one variable with respect to V∗.(8)bc2d22k1qk4−k3V∗2+w2−w1R0V∗−c2d2k1k3+c1d1d2k3−bc2d2d3k1k2R0−1=0,where(9)w1=k2qc2d2k1k3−bc2d2d3k1+cd1d2k3,w2=w1+βc1d3k22q+βc1d2k2k4−βc2d2k2k3+bc2d22k1k4−k3,R0=βk2T0d2c2E0+d1=βc1d3k2c2d2k1k3−bc2d2d3k1+c1d1d2k3.WhenR0>1, equation (8) has a unique positive root V∗. Therefore, R0 is called as the basic regeneration number of system (1), when R0>1 and R1>1, system (1) has a unique positive equilibrium Q∗T∗,I∗,V∗,E∗. R0 is the number of basic regeneration, and R1 is the number of immune response regeneration. R0 and R1 together ensure the existence of positive equilibrium Q∗.


## 4. Stability Analysis of the Equilibrium

In order to study the stability of the equilibrium, we consider the linear approximation equation of system (1) at the equilibrium T¯,I¯,V¯,E¯.(10)X′t=A˜Xt+B˜Xt−τ1+C˜Xt−τ2,where Xt=T,I,V,ET,(11)A˜=βV¯1+qV¯−bk1T¯−bk1T¯0−c1T¯0−c2E¯−d10−c2I¯0k2−d20k3E¯k4E¯00,B˜=−βV¯1+qV¯0−βT¯1+qV¯20000000000000,C˜=0000βV¯1+qV¯0βT¯1+qV¯2000000000.Let(12)detλJ−A˜−B˜e−λτ1−C˜e−λτ2=0,where J is a fourth-order identity matrix.


### 4.1. Global and Local Stability of Virusfree Equilibrium

Theorem 1.
IfR0≤βk2/βk2+k1bd2<1, the virusfree equilibrium Q0T0,0,0,E0 is globally asymptotically stable for any τ2≥0 and τ1=0.Proof.
Assumingc1/k3=c2/k4, we define the Lyapunov functional(13)Lt=Tt−T0−T0lnTtT0+It+d1+c2E0−k1bT0k2Vt+c1k3Et−E0−E0lnEtE0+∫t−τ2tβTsVs1+qVsds.
Further, calculating the derivative along with system (1)(14)L′t=T′t−T0TtT′t+I′t+d1+c2E0−k1bT0k2V′t+c1k3E′t−E0EtE′t+βTtVt1+qVt−βTt−τ2Vt−τ21+qVt−τ2=k1Tt1−bTt+It−βTtVt1+qVt−c1TtEt−k1T01−bTt+It+βT0Vt1+qVt+c1T0Et+βTt−τ2Vt−τ21+qVt−τ2−c2ItEt−d1It+d1+c2E0−k1bT0k2k2It−d2Vt+c1k3k3EtTt+k4ItEt−d3Et−c1k3k3E0Tt+k4ItE0−d3E0+βTtVt1+qVt−βTt−τ2Vt−τ21+qVt−τ2=−bk1Tt−T02−bk1TtIt+βT011+qVt+k1bd2βk2−1R0Vt≤−bk1Tt−T02−bk1TtIt−βT01R0−1−k1bd2βk2Vt.
WhenR0≤βk2/βk2+k1bd2<1, 1/R0−1−k1bd2/βk2≥0 is obtained, therefore L′t≤0. It is obvious to obtain that L′t=0 if and only if Tt=T0,It=0,Vt=0. Assuming that M̃ is the largest invariant set contained in T,I,V,E|L′t=0, it is easy to get M̃ =  T0,0,0,E0. According to the LaSalle’s invariant principle, the virusfree equilibrium Q0 is globally asymptotically stable.Theorem 2.
IfR0<1, the virusfree equilibrium Q0T0,0,0,E0 is locally asymptotically stable for any τ1≥0 and τ2≥0. If R0>1, Q0 is unstable for any τ1≥0 and τ2≥0.Proof.
According to equation (12), the characteristic equation of virusfree equilibrium Q0T0,0,0,E0 is as follows:(15)λ+bk1T000−k3E0bk1T0λ+c2E0+d1−k2−k4E0βe−λτ1T0−βe−λτ2T0λ+d20c1T000λ=0,
Through certain calculations, it follows that(16)L0λ,τ1,τ2≜λ+c2E0+d1λ+d2−βe−λτ2T0k2×λλ+bk1T0+c1k3E0T0=0.
Considering the following two equations separately(17)λλ+bk1T0+c1k3E0T0=0,(18)L1λ,τ1,τ2≜λ+c2E0+d1λ+d2−βe−λτ2T0k2=0.(i)
Equation (17) can be written as λ2+bk1T0λ+c1k3E0T0=0. Because bk1T0>0 and c1k3E0T0>0, the roots of equation (15) have negative real parts for all τ1≥0 and τ2≥0.(ii)
To discuss the case of the roots of equation (18), first of all, we consider τ1=τ2=0, so it is easy to get equation (18) as follows:(19)λ2+c2E0+d1+d2λ+d2c2E0+d11−R0=0.
WhenR0<1, it follows that d2c2E0+d11−R0>0, therefore, all the roots of (19) have negative real parts. When R0=1, λ=0 is a single root of (18), and the remaining roots have negative real parts. When R0>1, d2c2E0+d11−R0<0, therefore, (19) has a root of positive real parts.
Secondly, we considerτ1>0 and τ2>0. Assuming that for some τ1>0 and τ2>0, λ=ωiω>0 is the roots of (18), substituting λ=ωi into (18) and separating the real and imaginary parts to get that(20)cosωτ2=−c1ω2+d2c2k1−bc2k1T0+c1d1βc1k2T0,sinωτ2=−c1d2+c2k1−bc2k1T0+c1d1ωβc1k2T0.
Based onsin2ωτ2+cos2ωτ2=1, it follows that(21)c1ω4+H12−2c1J1ω2+J1J1+βc1k2T01−R0=0,where(22)H1=c1d2+c2k1−bc2k1T0+c1d1,J1=d2c2k1−bc2k1T0+c1d1>0.
IfR0<1, because of the following inequality(23)H12−2c1J1=c1d2+J1d22−2c1J1=c12d22+J12d22>0,J1J1+βc1k2T01−R0>0.
Equation (21) has no positive root. Thus, when R0<1, (18) has no pure imaginary root. Consequently, it can be seen from (i) and (ii) that when R0<1, Q0 is locally asymptotically stable for any τ1≥0 and τ2≥0. When R0>1, L1λ,τ1,τ2 is continuous with respect to λ,τ1, and τ2, respectively. For any τ1≥0,τ2≥0, we get L10,τ1,τ2=d2c2E0+d11−R0<0 and limλ⟶+∞L1λ,τ1,τ2=+∞. Therefore, there exists λ=λ¯τ1,τ2>0 such that L1λ¯τ1,τ2,τ1,τ2=0 holds. In summary, when R0>1, the equilibrium Q0 is unstable for any τ1≥0 and τ2≥0.


### 4.2. Local Stability Analysis of the Virus Infection Equilibrium

Suppose that the positive equilibrium of system (1) is Q∗T∗,I∗,V∗,E∗, and the characteristic equation corresponding to the positive equilibrium Q∗ of equation (12) is as follows:(24)λ+bk1T∗−βV∗1+qV∗+e−λτ1βV∗1+qV∗bk1T∗e−λτ1βT∗1+qV∗2c1T∗−e−λτ2βV∗1+qV∗λ+c2E∗+d1−e−λτ2βT∗1+qV∗2c2I∗0−k2λ+d20−k3E∗−k4E∗0λ=0,which is equivalent to the following equation:(25)Aλe−λτ1+Bλe−λτ2+Cλ=0,where(26)Aλ=A3λ3+A2λ2+A1λ+A0,Bλ=B2λ2+B1λ+B0,Cλ=λ4+C3λ3+C2λ2+C1λ+C0,A0=βc2d2k4E∗I∗V∗1+qV∗−βc2k2k3E∗I∗T∗1+qV∗2,A1=βV∗d1d2+c2d2E∗+c2k4E∗I∗1+qV∗,A2=βV∗c2E∗+d1+d21+qV∗,A3=βV∗1+qV∗,B0=βc1d2k4T∗V∗E∗1+qV∗−βc1k2k3E∗T∗21+qV∗2,B1=βbk1d2T∗V∗+βc1k4E∗T∗V∗1+qV∗−βbk1k2T∗21+qV∗2+β2k2T∗V∗1+qV∗3,B2=βbk1T∗V∗1+qV∗−βk2T∗1+qV∗2,C0=bk1T∗−βV∗1+qV∗c2d2k4E∗I∗−bc2d2k1k3T∗E∗I∗+c1d2k3E∗T∗c2E∗+d1,C1=bk1T∗−βV∗1+qV∗c2d2E∗+d1d2+c2k4E∗I∗+c2d2k4E∗I∗−bc2k1k3T∗E∗I∗+c1k3E∗T∗c2E∗+d1+d2,C2=bk1T∗−βV∗1+qV∗c2E∗+d1+d2+c2d2E∗+d1d2+c2k4E∗I∗+c1k3T∗E∗,C3=bk1T∗−βV∗1+qV∗+c2E∗+d1+d2.The caseτ1=0 and τ2=0.Whenτ1=0 and τ2=0, the characteristic equation (25) becomes(27)λ4+a3λ3+a2λ2+a1λ+a0=0,where(28)a0=βk2T∗I∗E∗1+qV∗c1k4−c2k3+βc2k2k3T∗I∗E∗1+qV∗−βc2k2k3T∗I∗E∗1+qV∗2+bc2d2k1k4−k3T∗I∗E∗+βc1k2k3T∗2E∗1+qV∗−βc1k2k3T∗2E∗1+qV∗2>0,a1=bk1T∗βd2V∗1+qV∗+βk2T∗1+qV∗−βk2T∗1+qV∗2+c2E∗I∗k4−k3+βc1k4E∗T∗V∗1+qV∗+β2k2T∗V∗1+qV∗3+c2d2k4I∗E∗+c1d1k3E∗T∗+c1d2k3E∗T∗+c1c2k3E∗2T∗>0,a2=bk1T∗d1+d2+c2E∗+βV∗1+qV∗+βk2T∗1+qV∗−βk2T∗1+qV∗2+c2k4I∗E∗+c1k3T∗E∗>0,a3=d1+d2+c2E∗+bk1T∗>0.(H1)a1a2a3−a12−a32a0>0.The following theorem is obtained by the Routh–Hurwitz criterion.Theorem 3.
Supposeτ1=0 and τ2=0, if R0>1,R1>1 and condition (H1) holds, then the equilibrium Q∗ is locally asymptotically stable.


#### 4.2.1. The Caseτ1=0 and τ2>0

Whenτ1=0 and τ2>0, the characteristic equation (25) is written as follows:(29)Dλ,0,τ2=λ4+M1λ3+M2λ2+M3λ+M4+M5λ2+M6λ+M7e−λτ2=0,where(30)M1=bk1T∗+c2E∗+d1+d2,M2=bk1T∗c2E∗+d1+d2+c2d2E∗+d1d2+c2k4I∗E∗+c1k3T∗E∗,M3=bk1T∗c2d2E∗+d1d2+c2I∗E∗k4−k3+c2d2k4I∗E∗+c1k3T∗E∗c2E∗+d1+d2,M4=bk1c2d2k4−k3T∗I∗E∗+c1d2k3T∗E∗c2E∗+d1−βc2k2k3T∗I∗E∗1+qV∗2,M5=βbk1T∗V∗1+qV∗−βk2T∗1+qV∗2,M6=βbk1d2T∗V∗+βc1k4E∗T∗V∗1+qV∗+β2k2T∗V∗1+qV∗3−βbk1k2T∗21+qV∗2,M7=βc1d2k4E∗T∗V∗1+qV∗−βc1k2k3T∗2E∗1+qV∗2.Assuming thatλ=ωiω>0 is a pure imaginary root of (29), substituting λ=ωi into (29) and separating the real and the imaginary parts to obtain that(31)sinωτ2=−M6ωω4−M2ω2+M4+M3ω−M1ω3M7−M5ω2M7−M5ω22+M6ω2≜P1,cosωτ2=−ω4−M2ω2+M4M7−M5ω2+M3ω−M1ω3M6ωM7−M5ω22+M6ω2≜P2.According tosin2ωτ2+cos2ωτ2=1, we can get that(32)ω8+l0ω6+l1ω4+l2ω2+l3=0,where(33)l0=M12−2M2,l1=M22+2M4−2M1M3−M52,l2=M32−2M2M4−M62+2M5M7,l3=M42−M72.Letω2=u, it follows from (32) that(34)fu≜u4+l0u3+l1u2+l2u+l3=0.Letukk=1,2,3,4 be the positive roots of (32), then (29) has pure imaginary roots ±ωki, where ωk=uk. From (31), the following solutions are obtained:(35)τ2kj=1ωkarc cosP2+2πj,P1>0,2π−arc cosP2+2πj,P1<0,where k=1,2,3,4, j=0,1,2,3,⋯. The distribution of the roots of (34) is discussed next, and the following conditions are given:(S1) (1)Γ>0,N1>0,N5>0,l2<0 and l3<0. (2) Γ<0 and l3<0.It is sufficient if any one of the above conditions is satisfied.(S2) (1)Γ>0,N1>0,N5>0,l1<0 and l3>0.(2)Γ<0 and −l0/4−sgnN2N6±N2/N6−N6+3N1>0.Meeting any one of the above conditions can be established.(S3)Γ>0,N1>0,N5>0,l2>0 and l3<0.(S4)Γ>0,N1>0,N5>0,l0<0,l1>0,l2<0 and l3>0.(S5) (1)Γ>0 and N5≤0.(2)Γ>0 and N1≤0.(3)Γ>0,N5>0,l0>0,l1>0,l2>0 and l3>0.(4)Γ<0 and −l0/4−sgnN2N6±N2/N6−N6+3N1<0.Any one of the above conditions can be satisfied, where(36)N1=l042−l16,N2=l24−3×l04×l16+2×l043,N3=l3−l0l24+3×l062,N4=4N13−N1N3−N22,N5=12N12−N3,Γ=N33−27N42,N6=12−N4+−Γ273+−N4−−Γ273+N1,sgnN2=1,ifN2>0,−1,ifN2<0.Next, we discuss the relationship between the roots of the characteristic equation and the coefficient. According to conclusions (see [17]), the following lemma is obtained.Lemma 1.
For equation (34), the following conclusions are drawn:(i)
If (S1) holds, then equation (34) has only one root of the positive real parts. Satisfying condition (1), the equation has one positive root and three negative roots. Satisfying condition (2), there is a positive root, a negative root, and a pair of conjugate virtual roots.(ii)
If (S2) holds, then equation (34) has two roots with positive real parts. Satisfying condition (1), the equation has two positive roots and two negative roots. Satisfying condition (2), there are two positive roots and a pair of conjugate virtual roots.(iii)
If (S3) holds, then equation (34) has three different positive real roots and one negative root.(iv)
If (S4) holds, then equation (34) has four different positive real roots.(v)
If (S5) holds, then equation (34) has no positive real root. That is, if (1) or (2) is satisfied, the equation has two pairs of different conjugate virtual roots. Satisfying condition (3), equation (34) has four different negative roots. Satisfying condition (4), there are two negative roots and a pair of conjugate virtual roots.According to Lemma1, the following lemma is obtained:Lemma 2.
For equation (29), the following conclusions are drawn:(i)
If (S1) holds andτ2=τ2kj, then equation (29) has a pair of pure imaginary roots ±ωki, where k=1 and j=0,1,2,⋯.(ii)
If (S2) holds andτ2=τ2kj, then equation (29) has two pairs of pure imaginary roots ±ωki, where k=2,3 and j=0,1,2,⋯.(iii)
If (S3) holds andτ2=τ2kj, then equation (29) has three pairs of pure imaginary roots ±ωki, where k=4,5,6 and j=0,1,2,⋯.(iv)
If (S4) holds andτ2=τ2kj, then equation (29) has four pairs of pure imaginary roots ±ωki, where k=7,8,9,10 and j=0,1,2,⋯.(v)
If (S5) holds andτ2=τ2kj, then equation (29) has no pure imaginary root.Lemma 3.
Ifωk is a solution to equation (32), then(37)signdReλdτ2|τ2=τ2kj=signdfdu|u=ωk2=sign4ωk6+3l0ωk4+2l1ωk2+l2.Proof.
Assuming thatλτ2=ξτ2+iωτ2 is the solution of equation (25), and both sides of the equation are differentiated from τ2 at the same time, it follows that(38)dλdτ2=λBe−λτ2Cλ′+Aλ′+e−λτ2Bλ′−τ2B,where Cλ′,Aλ′, and Bλ′ denote the derivation of λ, respectively. According to e−λτ2=−A+C/B, we obtain that(39)dλdτ2=λBλ′/B−Cλ′+Aλ′/C+A−τ2.
Through further calculations, we get that(40)signdReλdτ2|τ2=τ2kj=signReωkiBλ′/B−Cλ′+Aλ′/C+A−τ2kj=signRe1ωkiBλ′B−Cλ′+Aλ′C+A−τ2kj=signIm1ωkBλ′B−Cλ′+Aλ′C+A=signIm1ωkBλ′B¯−Cλ′+Aλ′C+A¯.When τ2=τ2kj, ξτ2=0, it follows that λτ2=iωkτ2, and then substitute λτ2=iωkτ2 into B2=A+C2 to get the following equation:(41)Im1ωkBλ′B¯−Cλ′+Aλ′C+A¯=dfdu|u=ωk2.
Therefore,(42)signIm1ωkBλ′B¯−Cλ′+Aλ′C+A¯=signdfdu|u=ωk2.The proof is completed.According to Lemmas1–3, the following theorem is obtained.Theorem 4.
For equation (29), the following conclusions are established:(i)
If (S1) is true, thensigndReλ/dτ2|τ2=τ21j>0.(ii)
If (S2) is true, thensigndReλ/dτ2|τ2=τ22j<0 and signdReλ/dτ2|τ2=τ23j>0.(iii)
If (S3) is true, thensigndReλ/dτ2|τ2=τ2kj>0k=4,6, dReλ/dτ2|τ2=τ25j<0.(iv)
If (S4) is true, thendReλ/dτ2|τ2=τ2kj<0k=7,9, dReλ/dτ2|τ2=τ2kj>0k=8,10.Proof.
(i)
If (S1) is true, then equation (34) has one positive root and three negative roots or one positive root, one negative root, and a pair of conjugate imaginary roots, and fu is monotonically increasing at the positive root u1, so we get f′u1>0. Therefore, the conclusion is true.(ii)
If (S2) is true, then equation (34) has two positive roots (u2,u3 and u2<u3) and two negative roots or has two positive roots and a pair of conjugate imaginary roots. Function fu decreases monotonically at the positive root u2 and increases monotonically at the positive root u3, so f′u2<0 and f′u3>0 are obtained. Therefore, the conclusion is true.(iii)
If (S3) is true, equation (34) has three different positive real roots (u4,u5,u6 and u4<u5<u6) and one negative root. Function fu increases monotonically at u4, decreases monotonically at u5, and increases monotonically at u6. So, we get f′u4>0, f′u5<0 and f′u6>0. Therefore, the conclusion is true.(iv)
If (S4) holds, equation (34) has four different positive real roots (ui,i=7,8,9,10 and u7<u8<u9<u10). Function fu is monotonically decreasing at u7 and u9, and monotonically increasing at u8 and u10. So we obtain f′u7<0, f′u8>0, f′u9<0 and f′u10>0. Therefore, the conclusion is true.According to the analysis in references (see [18–21]), Lemma 4 and Theorem 5 are obtained on the basis of Theorems 3 and 4.Lemma 4.
For equation (25), when R0>1 and R1>1, the following conclusions are obtained:(i)
If (H1) and (S1) are true, then forτ2∈0,τ21,0, all the roots have a negative real part. There is a pair of pure imaginary roots at τ2=τ21,0. There is at least one root with a positive real part at τ2>τ21,0.(ii)
If (H1) and (S2) are true, then there exists a positive integerθ that makes τ23,0<τ22,0<τ23,1<τ22,1<⋯<τ23,θ−2<τ22,θ−2<τ23,θ−1<τ23,θ<τ22,θ−1. This indicates that the system goes through θ stability switches from stable state to an unstable state. All the roots have negative real parts at τ2∈τ22,j,τ23,j+1τ22,−1=0,j=−1,0,…,θ−1. There are two pairs of pure imaginary roots at τ2=τ22,j or τ2=τ23,jj=0,1,…,θ−1, and the other roots have negative real parts. There is at least one root with a positive real part at τ2∈τ23,j,τ22,jj=−1,0,…,θ−1 or τ2>τ23,θ.(iii)
If (S2) is true and (H1) is not, then there exists a positive integerθ that makes τ22,0<τ23,0<τ22,1<τ23,1<⋯<τ22,θ−1<τ23,θ−1<τ23,θ<τ22,θ. This means that the system goes through θ stability switches from an unstable state to a stable state. There is at least one root with a positive real part at τ2∈τ23,j,τ22,j+1τ23,−1=0,j=−1,0,…,θ−2 or τ2>τ23,θ−1. There are two pairs of pure imaginary roots at τ2=τ22,j or τ2=τ23,jj=0,1,…,θ−1, and the other roots have negative real parts. All the roots have negative real parts at τ2∈τ22,j,τ23,jj=−1,0,…,θ−1.(iv)
If (H1) and (S3) are true, or (H1) and (S4) are true, then the system has at least one stability switch.(v)
If (H1) and (S5) are true, then for anyτ2>0, all roots have negative real parts.Theorem 5.
SupposeR0>1 and R1>1.(i)
If (H1) and (S1) are true, then the equilibriumQ∗ is locally asymptotically stable for any τ2∈0,τ21,0. When τ2=τ21,0, system (1) appears Hopf bifurcation at the equilibrium Q∗.(ii)
If (H1) and (S2) are true, whenτ2∈0∪τ22,j,τ23,j+1τ22,−1=0,j=−1,0,…,θ−2, the equilibrium Q∗ is locally asymptotically stable. When τ2∈τ23,j,τ22,jj=−1,0,…,θ−1 or τ2>τ23,θ, the equilibrium Q∗ is unstable. When τ2>τ23,m or τ2∈τ22,j∪τ23,j∪τ23,m−1j=−1,0,…,θ−2, system (1) appears Hopf bifurcation at Q∗.(iii)
If (S2) is true and (H1) is not. Whenτ2∈τ22,j,τ23,jj=−1,0,…,θ−1, the equilibrium Q∗ is locally asymptotically stable. When τ2∈0∪τ23,j,τ22,j+1j=−1,0,…,θ−1 or τ2>τ23,θ, Q∗ is unstable. When τ2>τ23,m or τ2∈τ22,j∪τ23,jj=−1,0,…,θ−1, Hopf bifurcation occurs at Q∗.(iv)
If (H1) and (S3), or (H1) and (S4) are true, then there is a stability switch atτ2=τ2∗. When τ2∗ satisfies ωτ2∗∈ωkj,k=4,5,6,7,8,9,10, Hopf bifurcation occurs at Q∗.(v)
If (H1) and (S5) are true, thenQ∗ is locally asymptotically stable for any τ2>0.


#### 4.2.2. The Caseτ1>0 and τ2>0

In order to discuss the possible Hopf bifurcation atτ1>0 and τ2>0, we regard τ2 as a parameter of τ1. According to the analyses and conclusions of (see [19,22]), system (1) is stable at τ1=0 and τ2>0, there is τ¯1τ2, and the system reaches stability at τ1∈0,τ¯1τ2. Therefore, the following theorem can be obtained:Theorem 6.
Let us assume thatR0>1 and R1>1.(i)
If (H1) and (S1) are true, then for anyτ2∈0,τ21,0, there exists τ¯1τ2, and the equilibrium Q∗ is locally asymptotically stable at τ1∈0,τ¯1τ2.(ii)
If (H1)and (S2) are true, there exists positive integerθ, when τ2∈0∪τ22,j,τ23,j+1τ22,−1=0,j=−1,0,…,θ−2, there exists τ¯1τ2, and the equilibrium Q∗ is locally asymptotically stable at τ1∈0,τ¯1τ2.(iii)
If (H1) and (S5) hold, then for anyτ2>0, there exists τ¯1τ2, and the equilibrium Q∗ is locally asymptotically stable at τ1∈0,τ¯1τ2.According to the analysis (see [23]), when τ2 is fixed, the characteristic equation of system (1) only depends on the change of time delay τ1. Therefore, we fix τ2 as a constant τ2=τ20 and study the stability of the system when τ1 increases from zero. In this case, the characteristic equation (25) becomes(43)D¯λ,τ1,τ20=Aλ,τ1,τ20e−λτ1+Bλ,τ1,τ20e−λτ20+Cλ,τ1,τ20=0,where(44)Aλ,τ1,τ20=R¯3τ1,τ20λ3+R¯2τ1,τ20λ2+R¯1τ1,τ20λ+R¯0τ1,τ20,Bλ,τ1,τ20=P¯2τ1,τ20λ2+P¯1τ1,τ20λ+P¯0τ1,τ20,Cλ,τ1,τ20=λ4+Q¯3τ1,τ20λ3+Q¯2τ1,τ20λ2+Q¯1τ1,τ20λ+Q¯0τ1,τ20.Assuming that the characteristic equation has pure imaginary rootλ=ωiω>0, we substitute it into the characteristic (43) and then separate the real parts and imaginary parts to obtain that(45)sinωτ1=R¯33ω3−R¯1ωP¯0−P¯2ω2cosωτ20+P¯1ωsinωτ20+ω4−Q¯2ω2+Q¯0R¯32ω6+R¯22−2R¯1R¯3ω4+R¯12−2R¯0R¯2ω2+R¯02+R¯0−R¯2ω2P¯2ω2−P¯0sinωτ20+P¯1ωcosωτ20−Q¯3ω3+Q¯1ωR¯32ω6+R¯22−2R¯1R¯3ω4+R¯12−2R¯0R¯2ω2+R¯02,cosωτ1=R¯33ω3−R¯1ωP¯2ω2−P¯0sinωτ20+P¯1ωcosωτ20−Q¯3ω3+Q¯1ωR¯32ω6+R¯22−2R¯1R¯3ω4+R¯12−2R¯0R¯2ω2+R¯02+R¯2ω2−R¯0P¯0−P¯2ω2cosωτ20+P¯1ωsinωτ20+ω4−Q¯2ω2+Q¯0R¯32ω6+R¯22−2R¯1R¯3ω4+R¯12−2R¯0R¯2ω2+R¯02.According toe−iωτ1=−Cωi,τ1,τ20+Bωi,τ1,τ20e−iωτ2/Aωi,τ1,τ20, it follows that(46)Cωi,τ1,τ20+Bωi,τ1,τ20e−iωτ2Aωi,τ1,τ202=1.Through further calculations, we get the following equation:(47)F¯ωi=Cωi,τ1,τ20+Aωi,τ1,τ20e−iωτ22−Bωi,τ1,τ202=ω8+h6τ1,τ20ω6+h5τ1,τ20ω5+h4τ1,τ20ω4+h3τ1,τ20ω3+h2τ1,τ20ω2+h1τ1,τ20ω+h0τ1,τ20=0,where(48)h6τ1,τ20=Q¯32−R¯32−2Q¯2−2P¯2cosωτ20,h5τ1,τ20=2P¯1sinωτ20−2P¯2Q¯3sinωτ20,h4τ1,τ20=P¯22+Q¯22+2Q¯0−R¯22+2P¯0cosωτ20−2P¯1Q¯3cosωτ20−2Q¯1Q¯3+2R¯1R¯3+2P¯2Q¯2cosωτ20,h3τ1,τ20=2P¯2Q¯1sinωτ20−2P¯1Q¯2sinωτ20+2P¯0Q¯3sinωτ20,h2τ1,τ20=P¯12+Q¯12−R¯12+2R¯0R¯2−2Q¯0Q¯2−2P¯0P¯2−2P¯2Q¯0cosωτ20+2P¯1Q¯1cosτ20−2P¯0Q¯2cosωτ20,h1τ1,τ20=2P¯1Q¯0sinωτ20−2P¯0Q¯1sinωτ20,h0τ1,τ20=Q¯02+P¯02−R¯02+2P¯0Q¯0cosωτ20.In order to study the influence of two time delays on the stability of the system, the following five criteria are given according to reference (see [23]):(U1) For anyτ1,τ20∈R02, there is A0,τ1,τ20+B0,τ1,τ20+C0,τ1,τ20≠0, that is, λ=0 is not the root of the characteristic equation (25).(U2) Ifλ=ωi, then for any τ1≥0, inequality Aωi,τ1,τ20+Bωi,τ1,τ20+Cωi,τ1,τ20≠0 can be obtained.(U3) For anyτ1, limλ⟶∞supAλ,τ1,τ20/e−λτ20Bλ,τ1,τ20+Cλ,τ1,τ20:Reλ≥0<1 holds.(U4) For anyτ1, the number of roots of the following function F¯ is limited.(49)F¯=e−ωiτ20Biω,τ1,τ20+Ciω,τ1,τ202−Aiω,τ1,τ202.(U5) As long asF¯=0 has positive roots ωτ1,τ20, then it is continuously differentiable with respect to τ1.Next, we prove that for anyτ1∈0,τmax, where τmax is the maximum value that enables the equilibrium to exist, each item in the characteristic (43) conforms to the above five criteria. Suppose that (U2) is true, and then we prove that the other four criteria are true.By calculations, we obtainF¯0,τ1,τ20=Q¯0+P¯02−R¯02, limω⟶+∞F¯ω,τ1,τ20=+∞. For the convenience of expression, we give conditions.(G1)Q¯0+P¯0<R¯0.If (G1) holds, then (U1) is correct andF¯ω,τ1,τ20=0 has at least one positive root. Next, when τ1∈0,τmax, obviously, (U3) holds. Then, on the basis of (U2), we can conclude that (U4) is true. Finally, we verify (U5). Suppose that ω0 is a positive real root of F¯ω,τ1,τ20=0 and ω0,τ1,τ20 is a point in the domain of F¯ω,τ1,τ20=0. Obviously, F¯ω′ and F¯τ1′ exist and they are continuous in some regions of ω0,τ1,τ20. In addition, there is F¯ω′ω0,τ1,τ20≠0. Therefore, according to the implicit function theorem, the verification of (U5) is true.IfI˜∈ℝ+ is the set of roots of F¯=0, for any ω∈I˜, we obtain F¯ω,τ1,τ20=0. Let us define θτ1∈0,2π, we get cosθτ1 and sinθτ1.Let us define the mapping(50)Sn:I˜⟶ℝ,Snτ1=τ1−θτ1+2nπωτ1,τ20,n∈ℕ0,τ1∈0,τmax.There exists a positive rootωτ1 of F¯ω,τ1,τ20=0 in I˜. Obviously, Sn is continuously differentiable with respect to τ1. Therefore, according to [23], the following theorem can be obtained:Theorem 7.
Suppose that one of the following conditions is correct.(i)
For any fixedτ20∈0,τ1,0, (H1) and (S1) hold.(ii)
For any fixedτ20∈0,τ3,0, (H1) and (S2) hold.
Then we obtain the following conclusions: ifR0>1 and R1>1, (H1) and (G1) are true, then F¯ω,τ1,τ20=0 has a pair of pure imaginary roots ±iωτ1∗,τ20, where ωτ1∗,τ20>0. For certain n∈N0, when τ1∗∈0,τmax, if Snτ1∗=0, then the root of the characteristic equation crosses the imaginary axis from left to right (from right to left) at Δτ1∗>0Δτ1∗<0, where Table 1(51)Δτ1∗=signdReλdτ|λ=iωτ1∗,τ20=signF¯ω′ωτ1∗,τ20,τ1∗,τ20signdSnτdτ|τ=τ1∗Table 1
Parameter of the system (1).
ParametersUnitsDefmitionsb−1mm3The maximum carrying capacity of tumor cellβDay−1 mm−3The infection rate of tumor cellsC1Day−1The clearance rate of CTLs to uninfected tumor cellsC2Day−1The clearance rate of CTLs to infected tumor cellsd1Day−1 mm−3The mortality of infected tumor cellsd2Day−1 mm−3The mortality of virusd3Day−1 mm−3The mortality of CTLsk1Day−1 mm−3The maximum growth rate ofuninfected tumor cellsk2Day−1 mm−3The average rate at which the infected tumor cells release the virusk3Day−1 mm−3Corresponding to the infection of uninfected tumor cells rate of CTLs simulation.k4Day−1 mm−3Corresponding to the infection of infected tumor cells rate of CTLs simu lationq—A constant greater than 0


## 4.1. Global and Local Stability of Virusfree Equilibrium

Theorem 1.
IfR0≤βk2/βk2+k1bd2<1, the virusfree equilibrium Q0T0,0,0,E0 is globally asymptotically stable for any τ2≥0 and τ1=0.Proof.
Assumingc1/k3=c2/k4, we define the Lyapunov functional(13)Lt=Tt−T0−T0lnTtT0+It+d1+c2E0−k1bT0k2Vt+c1k3Et−E0−E0lnEtE0+∫t−τ2tβTsVs1+qVsds.
Further, calculating the derivative along with system (1)(14)L′t=T′t−T0TtT′t+I′t+d1+c2E0−k1bT0k2V′t+c1k3E′t−E0EtE′t+βTtVt1+qVt−βTt−τ2Vt−τ21+qVt−τ2=k1Tt1−bTt+It−βTtVt1+qVt−c1TtEt−k1T01−bTt+It+βT0Vt1+qVt+c1T0Et+βTt−τ2Vt−τ21+qVt−τ2−c2ItEt−d1It+d1+c2E0−k1bT0k2k2It−d2Vt+c1k3k3EtTt+k4ItEt−d3Et−c1k3k3E0Tt+k4ItE0−d3E0+βTtVt1+qVt−βTt−τ2Vt−τ21+qVt−τ2=−bk1Tt−T02−bk1TtIt+βT011+qVt+k1bd2βk2−1R0Vt≤−bk1Tt−T02−bk1TtIt−βT01R0−1−k1bd2βk2Vt.
WhenR0≤βk2/βk2+k1bd2<1, 1/R0−1−k1bd2/βk2≥0 is obtained, therefore L′t≤0. It is obvious to obtain that L′t=0 if and only if Tt=T0,It=0,Vt=0. Assuming that M̃ is the largest invariant set contained in T,I,V,E|L′t=0, it is easy to get M̃ =  T0,0,0,E0. According to the LaSalle’s invariant principle, the virusfree equilibrium Q0 is globally asymptotically stable.Theorem 2.
IfR0<1, the virusfree equilibrium Q0T0,0,0,E0 is locally asymptotically stable for any τ1≥0 and τ2≥0. If R0>1, Q0 is unstable for any τ1≥0 and τ2≥0.Proof.
According to equation (12), the characteristic equation of virusfree equilibrium Q0T0,0,0,E0 is as follows:(15)λ+bk1T000−k3E0bk1T0λ+c2E0+d1−k2−k4E0βe−λτ1T0−βe−λτ2T0λ+d20c1T000λ=0,
Through certain calculations, it follows that(16)L0λ,τ1,τ2≜λ+c2E0+d1λ+d2−βe−λτ2T0k2×λλ+bk1T0+c1k3E0T0=0.
Considering the following two equations separately(17)λλ+bk1T0+c1k3E0T0=0,(18)L1λ,τ1,τ2≜λ+c2E0+d1λ+d2−βe−λτ2T0k2=0.(i)
Equation (17) can be written as λ2+bk1T0λ+c1k3E0T0=0. Because bk1T0>0 and c1k3E0T0>0, the roots of equation (15) have negative real parts for all τ1≥0 and τ2≥0.(ii)
To discuss the case of the roots of equation (18), first of all, we consider τ1=τ2=0, so it is easy to get equation (18) as follows:(19)λ2+c2E0+d1+d2λ+d2c2E0+d11−R0=0.
WhenR0<1, it follows that d2c2E0+d11−R0>0, therefore, all the roots of (19) have negative real parts. When R0=1, λ=0 is a single root of (18), and the remaining roots have negative real parts. When R0>1, d2c2E0+d11−R0<0, therefore, (19) has a root of positive real parts.
Secondly, we considerτ1>0 and τ2>0. Assuming that for some τ1>0 and τ2>0, λ=ωiω>0 is the roots of (18), substituting λ=ωi into (18) and separating the real and imaginary parts to get that(20)cosωτ2=−c1ω2+d2c2k1−bc2k1T0+c1d1βc1k2T0,sinωτ2=−c1d2+c2k1−bc2k1T0+c1d1ωβc1k2T0.
Based onsin2ωτ2+cos2ωτ2=1, it follows that(21)c1ω4+H12−2c1J1ω2+J1J1+βc1k2T01−R0=0,where(22)H1=c1d2+c2k1−bc2k1T0+c1d1,J1=d2c2k1−bc2k1T0+c1d1>0.
IfR0<1, because of the following inequality(23)H12−2c1J1=c1d2+J1d22−2c1J1=c12d22+J12d22>0,J1J1+βc1k2T01−R0>0.
Equation (21) has no positive root. Thus, when R0<1, (18) has no pure imaginary root. Consequently, it can be seen from (i) and (ii) that when R0<1, Q0 is locally asymptotically stable for any τ1≥0 and τ2≥0. When R0>1, L1λ,τ1,τ2 is continuous with respect to λ,τ1, and τ2, respectively. For any τ1≥0,τ2≥0, we get L10,τ1,τ2=d2c2E0+d11−R0<0 and limλ⟶+∞L1λ,τ1,τ2=+∞. Therefore, there exists λ=λ¯τ1,τ2>0 such that L1λ¯τ1,τ2,τ1,τ2=0 holds. In summary, when R0>1, the equilibrium Q0 is unstable for any τ1≥0 and τ2≥0.


## 4.2. Local Stability Analysis of the Virus Infection Equilibrium

Suppose that the positive equilibrium of system (1) is Q∗T∗,I∗,V∗,E∗, and the characteristic equation corresponding to the positive equilibrium Q∗ of equation (12) is as follows:(24)λ+bk1T∗−βV∗1+qV∗+e−λτ1βV∗1+qV∗bk1T∗e−λτ1βT∗1+qV∗2c1T∗−e−λτ2βV∗1+qV∗λ+c2E∗+d1−e−λτ2βT∗1+qV∗2c2I∗0−k2λ+d20−k3E∗−k4E∗0λ=0,which is equivalent to the following equation:(25)Aλe−λτ1+Bλe−λτ2+Cλ=0,where(26)Aλ=A3λ3+A2λ2+A1λ+A0,Bλ=B2λ2+B1λ+B0,Cλ=λ4+C3λ3+C2λ2+C1λ+C0,A0=βc2d2k4E∗I∗V∗1+qV∗−βc2k2k3E∗I∗T∗1+qV∗2,A1=βV∗d1d2+c2d2E∗+c2k4E∗I∗1+qV∗,A2=βV∗c2E∗+d1+d21+qV∗,A3=βV∗1+qV∗,B0=βc1d2k4T∗V∗E∗1+qV∗−βc1k2k3E∗T∗21+qV∗2,B1=βbk1d2T∗V∗+βc1k4E∗T∗V∗1+qV∗−βbk1k2T∗21+qV∗2+β2k2T∗V∗1+qV∗3,B2=βbk1T∗V∗1+qV∗−βk2T∗1+qV∗2,C0=bk1T∗−βV∗1+qV∗c2d2k4E∗I∗−bc2d2k1k3T∗E∗I∗+c1d2k3E∗T∗c2E∗+d1,C1=bk1T∗−βV∗1+qV∗c2d2E∗+d1d2+c2k4E∗I∗+c2d2k4E∗I∗−bc2k1k3T∗E∗I∗+c1k3E∗T∗c2E∗+d1+d2,C2=bk1T∗−βV∗1+qV∗c2E∗+d1+d2+c2d2E∗+d1d2+c2k4E∗I∗+c1k3T∗E∗,C3=bk1T∗−βV∗1+qV∗+c2E∗+d1+d2.The caseτ1=0 and τ2=0.Whenτ1=0 and τ2=0, the characteristic equation (25) becomes(27)λ4+a3λ3+a2λ2+a1λ+a0=0,where(28)a0=βk2T∗I∗E∗1+qV∗c1k4−c2k3+βc2k2k3T∗I∗E∗1+qV∗−βc2k2k3T∗I∗E∗1+qV∗2+bc2d2k1k4−k3T∗I∗E∗+βc1k2k3T∗2E∗1+qV∗−βc1k2k3T∗2E∗1+qV∗2>0,a1=bk1T∗βd2V∗1+qV∗+βk2T∗1+qV∗−βk2T∗1+qV∗2+c2E∗I∗k4−k3+βc1k4E∗T∗V∗1+qV∗+β2k2T∗V∗1+qV∗3+c2d2k4I∗E∗+c1d1k3E∗T∗+c1d2k3E∗T∗+c1c2k3E∗2T∗>0,a2=bk1T∗d1+d2+c2E∗+βV∗1+qV∗+βk2T∗1+qV∗−βk2T∗1+qV∗2+c2k4I∗E∗+c1k3T∗E∗>0,a3=d1+d2+c2E∗+bk1T∗>0.(H1)a1a2a3−a12−a32a0>0.The following theorem is obtained by the Routh–Hurwitz criterion.Theorem 3.
Supposeτ1=0 and τ2=0, if R0>1,R1>1 and condition (H1) holds, then the equilibrium Q∗ is locally asymptotically stable.


### 4.2.1. The Caseτ1=0 and τ2>0

Whenτ1=0 and τ2>0, the characteristic equation (25) is written as follows:(29)Dλ,0,τ2=λ4+M1λ3+M2λ2+M3λ+M4+M5λ2+M6λ+M7e−λτ2=0,where(30)M1=bk1T∗+c2E∗+d1+d2,M2=bk1T∗c2E∗+d1+d2+c2d2E∗+d1d2+c2k4I∗E∗+c1k3T∗E∗,M3=bk1T∗c2d2E∗+d1d2+c2I∗E∗k4−k3+c2d2k4I∗E∗+c1k3T∗E∗c2E∗+d1+d2,M4=bk1c2d2k4−k3T∗I∗E∗+c1d2k3T∗E∗c2E∗+d1−βc2k2k3T∗I∗E∗1+qV∗2,M5=βbk1T∗V∗1+qV∗−βk2T∗1+qV∗2,M6=βbk1d2T∗V∗+βc1k4E∗T∗V∗1+qV∗+β2k2T∗V∗1+qV∗3−βbk1k2T∗21+qV∗2,M7=βc1d2k4E∗T∗V∗1+qV∗−βc1k2k3T∗2E∗1+qV∗2.Assuming thatλ=ωiω>0 is a pure imaginary root of (29), substituting λ=ωi into (29) and separating the real and the imaginary parts to obtain that(31)sinωτ2=−M6ωω4−M2ω2+M4+M3ω−M1ω3M7−M5ω2M7−M5ω22+M6ω2≜P1,cosωτ2=−ω4−M2ω2+M4M7−M5ω2+M3ω−M1ω3M6ωM7−M5ω22+M6ω2≜P2.According tosin2ωτ2+cos2ωτ2=1, we can get that(32)ω8+l0ω6+l1ω4+l2ω2+l3=0,where(33)l0=M12−2M2,l1=M22+2M4−2M1M3−M52,l2=M32−2M2M4−M62+2M5M7,l3=M42−M72.Letω2=u, it follows from (32) that(34)fu≜u4+l0u3+l1u2+l2u+l3=0.Letukk=1,2,3,4 be the positive roots of (32), then (29) has pure imaginary roots ±ωki, where ωk=uk. From (31), the following solutions are obtained:(35)τ2kj=1ωkarc cosP2+2πj,P1>0,2π−arc cosP2+2πj,P1<0,where k=1,2,3,4, j=0,1,2,3,⋯. The distribution of the roots of (34) is discussed next, and the following conditions are given:(S1) (1)Γ>0,N1>0,N5>0,l2<0 and l3<0. (2) Γ<0 and l3<0.It is sufficient if any one of the above conditions is satisfied.(S2) (1)Γ>0,N1>0,N5>0,l1<0 and l3>0.(2)Γ<0 and −l0/4−sgnN2N6±N2/N6−N6+3N1>0.Meeting any one of the above conditions can be established.(S3)Γ>0,N1>0,N5>0,l2>0 and l3<0.(S4)Γ>0,N1>0,N5>0,l0<0,l1>0,l2<0 and l3>0.(S5) (1)Γ>0 and N5≤0.(2)Γ>0 and N1≤0.(3)Γ>0,N5>0,l0>0,l1>0,l2>0 and l3>0.(4)Γ<0 and −l0/4−sgnN2N6±N2/N6−N6+3N1<0.Any one of the above conditions can be satisfied, where(36)N1=l042−l16,N2=l24−3×l04×l16+2×l043,N3=l3−l0l24+3×l062,N4=4N13−N1N3−N22,N5=12N12−N3,Γ=N33−27N42,N6=12−N4+−Γ273+−N4−−Γ273+N1,sgnN2=1,ifN2>0,−1,ifN2<0.Next, we discuss the relationship between the roots of the characteristic equation and the coefficient. According to conclusions (see [17]), the following lemma is obtained.Lemma 1.
For equation (34), the following conclusions are drawn:(i)
If (S1) holds, then equation (34) has only one root of the positive real parts. Satisfying condition (1), the equation has one positive root and three negative roots. Satisfying condition (2), there is a positive root, a negative root, and a pair of conjugate virtual roots.(ii)
If (S2) holds, then equation (34) has two roots with positive real parts. Satisfying condition (1), the equation has two positive roots and two negative roots. Satisfying condition (2), there are two positive roots and a pair of conjugate virtual roots.(iii)
If (S3) holds, then equation (34) has three different positive real roots and one negative root.(iv)
If (S4) holds, then equation (34) has four different positive real roots.(v)
If (S5) holds, then equation (34) has no positive real root. That is, if (1) or (2) is satisfied, the equation has two pairs of different conjugate virtual roots. Satisfying condition (3), equation (34) has four different negative roots. Satisfying condition (4), there are two negative roots and a pair of conjugate virtual roots.According to Lemma1, the following lemma is obtained:Lemma 2.
For equation (29), the following conclusions are drawn:(i)
If (S1) holds andτ2=τ2kj, then equation (29) has a pair of pure imaginary roots ±ωki, where k=1 and j=0,1,2,⋯.(ii)
If (S2) holds andτ2=τ2kj, then equation (29) has two pairs of pure imaginary roots ±ωki, where k=2,3 and j=0,1,2,⋯.(iii)
If (S3) holds andτ2=τ2kj, then equation (29) has three pairs of pure imaginary roots ±ωki, where k=4,5,6 and j=0,1,2,⋯.(iv)
If (S4) holds andτ2=τ2kj, then equation (29) has four pairs of pure imaginary roots ±ωki, where k=7,8,9,10 and j=0,1,2,⋯.(v)
If (S5) holds andτ2=τ2kj, then equation (29) has no pure imaginary root.Lemma 3.
Ifωk is a solution to equation (32), then(37)signdReλdτ2|τ2=τ2kj=signdfdu|u=ωk2=sign4ωk6+3l0ωk4+2l1ωk2+l2.Proof.
Assuming thatλτ2=ξτ2+iωτ2 is the solution of equation (25), and both sides of the equation are differentiated from τ2 at the same time, it follows that(38)dλdτ2=λBe−λτ2Cλ′+Aλ′+e−λτ2Bλ′−τ2B,where Cλ′,Aλ′, and Bλ′ denote the derivation of λ, respectively. According to e−λτ2=−A+C/B, we obtain that(39)dλdτ2=λBλ′/B−Cλ′+Aλ′/C+A−τ2.
Through further calculations, we get that(40)signdReλdτ2|τ2=τ2kj=signReωkiBλ′/B−Cλ′+Aλ′/C+A−τ2kj=signRe1ωkiBλ′B−Cλ′+Aλ′C+A−τ2kj=signIm1ωkBλ′B−Cλ′+Aλ′C+A=signIm1ωkBλ′B¯−Cλ′+Aλ′C+A¯.When τ2=τ2kj, ξτ2=0, it follows that λτ2=iωkτ2, and then substitute λτ2=iωkτ2 into B2=A+C2 to get the following equation:(41)Im1ωkBλ′B¯−Cλ′+Aλ′C+A¯=dfdu|u=ωk2.
Therefore,(42)signIm1ωkBλ′B¯−Cλ′+Aλ′C+A¯=signdfdu|u=ωk2.The proof is completed.According to Lemmas1–3, the following theorem is obtained.Theorem 4.
For equation (29), the following conclusions are established:(i)
If (S1) is true, thensigndReλ/dτ2|τ2=τ21j>0.(ii)
If (S2) is true, thensigndReλ/dτ2|τ2=τ22j<0 and signdReλ/dτ2|τ2=τ23j>0.(iii)
If (S3) is true, thensigndReλ/dτ2|τ2=τ2kj>0k=4,6, dReλ/dτ2|τ2=τ25j<0.(iv)
If (S4) is true, thendReλ/dτ2|τ2=τ2kj<0k=7,9, dReλ/dτ2|τ2=τ2kj>0k=8,10.Proof.
(i)
If (S1) is true, then equation (34) has one positive root and three negative roots or one positive root, one negative root, and a pair of conjugate imaginary roots, and fu is monotonically increasing at the positive root u1, so we get f′u1>0. Therefore, the conclusion is true.(ii)
If (S2) is true, then equation (34) has two positive roots (u2,u3 and u2<u3) and two negative roots or has two positive roots and a pair of conjugate imaginary roots. Function fu decreases monotonically at the positive root u2 and increases monotonically at the positive root u3, so f′u2<0 and f′u3>0 are obtained. Therefore, the conclusion is true.(iii)
If (S3) is true, equation (34) has three different positive real roots (u4,u5,u6 and u4<u5<u6) and one negative root. Function fu increases monotonically at u4, decreases monotonically at u5, and increases monotonically at u6. So, we get f′u4>0, f′u5<0 and f′u6>0. Therefore, the conclusion is true.(iv)
If (S4) holds, equation (34) has four different positive real roots (ui,i=7,8,9,10 and u7<u8<u9<u10). Function fu is monotonically decreasing at u7 and u9, and monotonically increasing at u8 and u10. So we obtain f′u7<0, f′u8>0, f′u9<0 and f′u10>0. Therefore, the conclusion is true.According to the analysis in references (see [18–21]), Lemma 4 and Theorem 5 are obtained on the basis of Theorems 3 and 4.Lemma 4.
For equation (25), when R0>1 and R1>1, the following conclusions are obtained:(i)
If (H1) and (S1) are true, then forτ2∈0,τ21,0, all the roots have a negative real part. There is a pair of pure imaginary roots at τ2=τ21,0. There is at least one root with a positive real part at τ2>τ21,0.(ii)
If (H1) and (S2) are true, then there exists a positive integerθ that makes τ23,0<τ22,0<τ23,1<τ22,1<⋯<τ23,θ−2<τ22,θ−2<τ23,θ−1<τ23,θ<τ22,θ−1. This indicates that the system goes through θ stability switches from stable state to an unstable state. All the roots have negative real parts at τ2∈τ22,j,τ23,j+1τ22,−1=0,j=−1,0,…,θ−1. There are two pairs of pure imaginary roots at τ2=τ22,j or τ2=τ23,jj=0,1,…,θ−1, and the other roots have negative real parts. There is at least one root with a positive real part at τ2∈τ23,j,τ22,jj=−1,0,…,θ−1 or τ2>τ23,θ.(iii)
If (S2) is true and (H1) is not, then there exists a positive integerθ that makes τ22,0<τ23,0<τ22,1<τ23,1<⋯<τ22,θ−1<τ23,θ−1<τ23,θ<τ22,θ. This means that the system goes through θ stability switches from an unstable state to a stable state. There is at least one root with a positive real part at τ2∈τ23,j,τ22,j+1τ23,−1=0,j=−1,0,…,θ−2 or τ2>τ23,θ−1. There are two pairs of pure imaginary roots at τ2=τ22,j or τ2=τ23,jj=0,1,…,θ−1, and the other roots have negative real parts. All the roots have negative real parts at τ2∈τ22,j,τ23,jj=−1,0,…,θ−1.(iv)
If (H1) and (S3) are true, or (H1) and (S4) are true, then the system has at least one stability switch.(v)
If (H1) and (S5) are true, then for anyτ2>0, all roots have negative real parts.Theorem 5.
SupposeR0>1 and R1>1.(i)
If (H1) and (S1) are true, then the equilibriumQ∗ is locally asymptotically stable for any τ2∈0,τ21,0. When τ2=τ21,0, system (1) appears Hopf bifurcation at the equilibrium Q∗.(ii)
If (H1) and (S2) are true, whenτ2∈0∪τ22,j,τ23,j+1τ22,−1=0,j=−1,0,…,θ−2, the equilibrium Q∗ is locally asymptotically stable. When τ2∈τ23,j,τ22,jj=−1,0,…,θ−1 or τ2>τ23,θ, the equilibrium Q∗ is unstable. When τ2>τ23,m or τ2∈τ22,j∪τ23,j∪τ23,m−1j=−1,0,…,θ−2, system (1) appears Hopf bifurcation at Q∗.(iii)
If (S2) is true and (H1) is not. Whenτ2∈τ22,j,τ23,jj=−1,0,…,θ−1, the equilibrium Q∗ is locally asymptotically stable. When τ2∈0∪τ23,j,τ22,j+1j=−1,0,…,θ−1 or τ2>τ23,θ, Q∗ is unstable. When τ2>τ23,m or τ2∈τ22,j∪τ23,jj=−1,0,…,θ−1, Hopf bifurcation occurs at Q∗.(iv)
If (H1) and (S3), or (H1) and (S4) are true, then there is a stability switch atτ2=τ2∗. When τ2∗ satisfies ωτ2∗∈ωkj,k=4,5,6,7,8,9,10, Hopf bifurcation occurs at Q∗.(v)
If (H1) and (S5) are true, thenQ∗ is locally asymptotically stable for any τ2>0.


### 4.2.2. The Caseτ1>0 and τ2>0

In order to discuss the possible Hopf bifurcation atτ1>0 and τ2>0, we regard τ2 as a parameter of τ1. According to the analyses and conclusions of (see [19,22]), system (1) is stable at τ1=0 and τ2>0, there is τ¯1τ2, and the system reaches stability at τ1∈0,τ¯1τ2. Therefore, the following theorem can be obtained:Theorem 6.
Let us assume thatR0>1 and R1>1.(i)
If (H1) and (S1) are true, then for anyτ2∈0,τ21,0, there exists τ¯1τ2, and the equilibrium Q∗ is locally asymptotically stable at τ1∈0,τ¯1τ2.(ii)
If (H1)and (S2) are true, there exists positive integerθ, when τ2∈0∪τ22,j,τ23,j+1τ22,−1=0,j=−1,0,…,θ−2, there exists τ¯1τ2, and the equilibrium Q∗ is locally asymptotically stable at τ1∈0,τ¯1τ2.(iii)
If (H1) and (S5) hold, then for anyτ2>0, there exists τ¯1τ2, and the equilibrium Q∗ is locally asymptotically stable at τ1∈0,τ¯1τ2.According to the analysis (see [23]), when τ2 is fixed, the characteristic equation of system (1) only depends on the change of time delay τ1. Therefore, we fix τ2 as a constant τ2=τ20 and study the stability of the system when τ1 increases from zero. In this case, the characteristic equation (25) becomes(43)D¯λ,τ1,τ20=Aλ,τ1,τ20e−λτ1+Bλ,τ1,τ20e−λτ20+Cλ,τ1,τ20=0,where(44)Aλ,τ1,τ20=R¯3τ1,τ20λ3+R¯2τ1,τ20λ2+R¯1τ1,τ20λ+R¯0τ1,τ20,Bλ,τ1,τ20=P¯2τ1,τ20λ2+P¯1τ1,τ20λ+P¯0τ1,τ20,Cλ,τ1,τ20=λ4+Q¯3τ1,τ20λ3+Q¯2τ1,τ20λ2+Q¯1τ1,τ20λ+Q¯0τ1,τ20.Assuming that the characteristic equation has pure imaginary rootλ=ωiω>0, we substitute it into the characteristic (43) and then separate the real parts and imaginary parts to obtain that(45)sinωτ1=R¯33ω3−R¯1ωP¯0−P¯2ω2cosωτ20+P¯1ωsinωτ20+ω4−Q¯2ω2+Q¯0R¯32ω6+R¯22−2R¯1R¯3ω4+R¯12−2R¯0R¯2ω2+R¯02+R¯0−R¯2ω2P¯2ω2−P¯0sinωτ20+P¯1ωcosωτ20−Q¯3ω3+Q¯1ωR¯32ω6+R¯22−2R¯1R¯3ω4+R¯12−2R¯0R¯2ω2+R¯02,cosωτ1=R¯33ω3−R¯1ωP¯2ω2−P¯0sinωτ20+P¯1ωcosωτ20−Q¯3ω3+Q¯1ωR¯32ω6+R¯22−2R¯1R¯3ω4+R¯12−2R¯0R¯2ω2+R¯02+R¯2ω2−R¯0P¯0−P¯2ω2cosωτ20+P¯1ωsinωτ20+ω4−Q¯2ω2+Q¯0R¯32ω6+R¯22−2R¯1R¯3ω4+R¯12−2R¯0R¯2ω2+R¯02.According toe−iωτ1=−Cωi,τ1,τ20+Bωi,τ1,τ20e−iωτ2/Aωi,τ1,τ20, it follows that(46)Cωi,τ1,τ20+Bωi,τ1,τ20e−iωτ2Aωi,τ1,τ202=1.Through further calculations, we get the following equation:(47)F¯ωi=Cωi,τ1,τ20+Aωi,τ1,τ20e−iωτ22−Bωi,τ1,τ202=ω8+h6τ1,τ20ω6+h5τ1,τ20ω5+h4τ1,τ20ω4+h3τ1,τ20ω3+h2τ1,τ20ω2+h1τ1,τ20ω+h0τ1,τ20=0,where(48)h6τ1,τ20=Q¯32−R¯32−2Q¯2−2P¯2cosωτ20,h5τ1,τ20=2P¯1sinωτ20−2P¯2Q¯3sinωτ20,h4τ1,τ20=P¯22+Q¯22+2Q¯0−R¯22+2P¯0cosωτ20−2P¯1Q¯3cosωτ20−2Q¯1Q¯3+2R¯1R¯3+2P¯2Q¯2cosωτ20,h3τ1,τ20=2P¯2Q¯1sinωτ20−2P¯1Q¯2sinωτ20+2P¯0Q¯3sinωτ20,h2τ1,τ20=P¯12+Q¯12−R¯12+2R¯0R¯2−2Q¯0Q¯2−2P¯0P¯2−2P¯2Q¯0cosωτ20+2P¯1Q¯1cosτ20−2P¯0Q¯2cosωτ20,h1τ1,τ20=2P¯1Q¯0sinωτ20−2P¯0Q¯1sinωτ20,h0τ1,τ20=Q¯02+P¯02−R¯02+2P¯0Q¯0cosωτ20.In order to study the influence of two time delays on the stability of the system, the following five criteria are given according to reference (see [23]):(U1) For anyτ1,τ20∈R02, there is A0,τ1,τ20+B0,τ1,τ20+C0,τ1,τ20≠0, that is, λ=0 is not the root of the characteristic equation (25).(U2) Ifλ=ωi, then for any τ1≥0, inequality Aωi,τ1,τ20+Bωi,τ1,τ20+Cωi,τ1,τ20≠0 can be obtained.(U3) For anyτ1, limλ⟶∞supAλ,τ1,τ20/e−λτ20Bλ,τ1,τ20+Cλ,τ1,τ20:Reλ≥0<1 holds.(U4) For anyτ1, the number of roots of the following function F¯ is limited.(49)F¯=e−ωiτ20Biω,τ1,τ20+Ciω,τ1,τ202−Aiω,τ1,τ202.(U5) As long asF¯=0 has positive roots ωτ1,τ20, then it is continuously differentiable with respect to τ1.Next, we prove that for anyτ1∈0,τmax, where τmax is the maximum value that enables the equilibrium to exist, each item in the characteristic (43) conforms to the above five criteria. Suppose that (U2) is true, and then we prove that the other four criteria are true.By calculations, we obtainF¯0,τ1,τ20=Q¯0+P¯02−R¯02, limω⟶+∞F¯ω,τ1,τ20=+∞. For the convenience of expression, we give conditions.(G1)Q¯0+P¯0<R¯0.If (G1) holds, then (U1) is correct andF¯ω,τ1,τ20=0 has at least one positive root. Next, when τ1∈0,τmax, obviously, (U3) holds. Then, on the basis of (U2), we can conclude that (U4) is true. Finally, we verify (U5). Suppose that ω0 is a positive real root of F¯ω,τ1,τ20=0 and ω0,τ1,τ20 is a point in the domain of F¯ω,τ1,τ20=0. Obviously, F¯ω′ and F¯τ1′ exist and they are continuous in some regions of ω0,τ1,τ20. In addition, there is F¯ω′ω0,τ1,τ20≠0. Therefore, according to the implicit function theorem, the verification of (U5) is true.IfI˜∈ℝ+ is the set of roots of F¯=0, for any ω∈I˜, we obtain F¯ω,τ1,τ20=0. Let us define θτ1∈0,2π, we get cosθτ1 and sinθτ1.Let us define the mapping(50)Sn:I˜⟶ℝ,Snτ1=τ1−θτ1+2nπωτ1,τ20,n∈ℕ0,τ1∈0,τmax.There exists a positive rootωτ1 of F¯ω,τ1,τ20=0 in I˜. Obviously, Sn is continuously differentiable with respect to τ1. Therefore, according to [23], the following theorem can be obtained:Theorem 7.
Suppose that one of the following conditions is correct.(i)
For any fixedτ20∈0,τ1,0, (H1) and (S1) hold.(ii)
For any fixedτ20∈0,τ3,0, (H1) and (S2) hold.
Then we obtain the following conclusions: ifR0>1 and R1>1, (H1) and (G1) are true, then F¯ω,τ1,τ20=0 has a pair of pure imaginary roots ±iωτ1∗,τ20, where ωτ1∗,τ20>0. For certain n∈N0, when τ1∗∈0,τmax, if Snτ1∗=0, then the root of the characteristic equation crosses the imaginary axis from left to right (from right to left) at Δτ1∗>0Δτ1∗<0, where Table 1(51)Δτ1∗=signdReλdτ|λ=iωτ1∗,τ20=signF¯ω′ωτ1∗,τ20,τ1∗,τ20signdSnτdτ|τ=τ1∗Table 1
Parameter of the system (1).
ParametersUnitsDefmitionsb−1mm3The maximum carrying capacity of tumor cellβDay−1 mm−3The infection rate of tumor cellsC1Day−1The clearance rate of CTLs to uninfected tumor cellsC2Day−1The clearance rate of CTLs to infected tumor cellsd1Day−1 mm−3The mortality of infected tumor cellsd2Day−1 mm−3The mortality of virusd3Day−1 mm−3The mortality of CTLsk1Day−1 mm−3The maximum growth rate ofuninfected tumor cellsk2Day−1 mm−3The average rate at which the infected tumor cells release the virusk3Day−1 mm−3Corresponding to the infection of uninfected tumor cells rate of CTLs simulation.k4Day−1 mm−3Corresponding to the infection of infected tumor cells rate of CTLs simu lationq—A constant greater than 0


## 4.2.1. The Caseτ1=0 and τ2>0

Whenτ1=0 and τ2>0, the characteristic equation (25) is written as follows:(29)Dλ,0,τ2=λ4+M1λ3+M2λ2+M3λ+M4+M5λ2+M6λ+M7e−λτ2=0,where(30)M1=bk1T∗+c2E∗+d1+d2,M2=bk1T∗c2E∗+d1+d2+c2d2E∗+d1d2+c2k4I∗E∗+c1k3T∗E∗,M3=bk1T∗c2d2E∗+d1d2+c2I∗E∗k4−k3+c2d2k4I∗E∗+c1k3T∗E∗c2E∗+d1+d2,M4=bk1c2d2k4−k3T∗I∗E∗+c1d2k3T∗E∗c2E∗+d1−βc2k2k3T∗I∗E∗1+qV∗2,M5=βbk1T∗V∗1+qV∗−βk2T∗1+qV∗2,M6=βbk1d2T∗V∗+βc1k4E∗T∗V∗1+qV∗+β2k2T∗V∗1+qV∗3−βbk1k2T∗21+qV∗2,M7=βc1d2k4E∗T∗V∗1+qV∗−βc1k2k3T∗2E∗1+qV∗2.Assuming thatλ=ωiω>0 is a pure imaginary root of (29), substituting λ=ωi into (29) and separating the real and the imaginary parts to obtain that(31)sinωτ2=−M6ωω4−M2ω2+M4+M3ω−M1ω3M7−M5ω2M7−M5ω22+M6ω2≜P1,cosωτ2=−ω4−M2ω2+M4M7−M5ω2+M3ω−M1ω3M6ωM7−M5ω22+M6ω2≜P2.According tosin2ωτ2+cos2ωτ2=1, we can get that(32)ω8+l0ω6+l1ω4+l2ω2+l3=0,where(33)l0=M12−2M2,l1=M22+2M4−2M1M3−M52,l2=M32−2M2M4−M62+2M5M7,l3=M42−M72.Letω2=u, it follows from (32) that(34)fu≜u4+l0u3+l1u2+l2u+l3=0.Letukk=1,2,3,4 be the positive roots of (32), then (29) has pure imaginary roots ±ωki, where ωk=uk. From (31), the following solutions are obtained:(35)τ2kj=1ωkarc cosP2+2πj,P1>0,2π−arc cosP2+2πj,P1<0,where k=1,2,3,4, j=0,1,2,3,⋯. The distribution of the roots of (34) is discussed next, and the following conditions are given:(S1) (1)Γ>0,N1>0,N5>0,l2<0 and l3<0. (2) Γ<0 and l3<0.It is sufficient if any one of the above conditions is satisfied.(S2) (1)Γ>0,N1>0,N5>0,l1<0 and l3>0.(2)Γ<0 and −l0/4−sgnN2N6±N2/N6−N6+3N1>0.Meeting any one of the above conditions can be established.(S3)Γ>0,N1>0,N5>0,l2>0 and l3<0.(S4)Γ>0,N1>0,N5>0,l0<0,l1>0,l2<0 and l3>0.(S5) (1)Γ>0 and N5≤0.(2)Γ>0 and N1≤0.(3)Γ>0,N5>0,l0>0,l1>0,l2>0 and l3>0.(4)Γ<0 and −l0/4−sgnN2N6±N2/N6−N6+3N1<0.Any one of the above conditions can be satisfied, where(36)N1=l042−l16,N2=l24−3×l04×l16+2×l043,N3=l3−l0l24+3×l062,N4=4N13−N1N3−N22,N5=12N12−N3,Γ=N33−27N42,N6=12−N4+−Γ273+−N4−−Γ273+N1,sgnN2=1,ifN2>0,−1,ifN2<0.Next, we discuss the relationship between the roots of the characteristic equation and the coefficient. According to conclusions (see [17]), the following lemma is obtained.Lemma 1.
For equation (34), the following conclusions are drawn:(i)
If (S1) holds, then equation (34) has only one root of the positive real parts. Satisfying condition (1), the equation has one positive root and three negative roots. Satisfying condition (2), there is a positive root, a negative root, and a pair of conjugate virtual roots.(ii)
If (S2) holds, then equation (34) has two roots with positive real parts. Satisfying condition (1), the equation has two positive roots and two negative roots. Satisfying condition (2), there are two positive roots and a pair of conjugate virtual roots.(iii)
If (S3) holds, then equation (34) has three different positive real roots and one negative root.(iv)
If (S4) holds, then equation (34) has four different positive real roots.(v)
If (S5) holds, then equation (34) has no positive real root. That is, if (1) or (2) is satisfied, the equation has two pairs of different conjugate virtual roots. Satisfying condition (3), equation (34) has four different negative roots. Satisfying condition (4), there are two negative roots and a pair of conjugate virtual roots.According to Lemma1, the following lemma is obtained:Lemma 2.
For equation (29), the following conclusions are drawn:(i)
If (S1) holds andτ2=τ2kj, then equation (29) has a pair of pure imaginary roots ±ωki, where k=1 and j=0,1,2,⋯.(ii)
If (S2) holds andτ2=τ2kj, then equation (29) has two pairs of pure imaginary roots ±ωki, where k=2,3 and j=0,1,2,⋯.(iii)
If (S3) holds andτ2=τ2kj, then equation (29) has three pairs of pure imaginary roots ±ωki, where k=4,5,6 and j=0,1,2,⋯.(iv)
If (S4) holds andτ2=τ2kj, then equation (29) has four pairs of pure imaginary roots ±ωki, where k=7,8,9,10 and j=0,1,2,⋯.(v)
If (S5) holds andτ2=τ2kj, then equation (29) has no pure imaginary root.Lemma 3.
Ifωk is a solution to equation (32), then(37)signdReλdτ2|τ2=τ2kj=signdfdu|u=ωk2=sign4ωk6+3l0ωk4+2l1ωk2+l2.Proof.
Assuming thatλτ2=ξτ2+iωτ2 is the solution of equation (25), and both sides of the equation are differentiated from τ2 at the same time, it follows that(38)dλdτ2=λBe−λτ2Cλ′+Aλ′+e−λτ2Bλ′−τ2B,where Cλ′,Aλ′, and Bλ′ denote the derivation of λ, respectively. According to e−λτ2=−A+C/B, we obtain that(39)dλdτ2=λBλ′/B−Cλ′+Aλ′/C+A−τ2.
Through further calculations, we get that(40)signdReλdτ2|τ2=τ2kj=signReωkiBλ′/B−Cλ′+Aλ′/C+A−τ2kj=signRe1ωkiBλ′B−Cλ′+Aλ′C+A−τ2kj=signIm1ωkBλ′B−Cλ′+Aλ′C+A=signIm1ωkBλ′B¯−Cλ′+Aλ′C+A¯.When τ2=τ2kj, ξτ2=0, it follows that λτ2=iωkτ2, and then substitute λτ2=iωkτ2 into B2=A+C2 to get the following equation:(41)Im1ωkBλ′B¯−Cλ′+Aλ′C+A¯=dfdu|u=ωk2.
Therefore,(42)signIm1ωkBλ′B¯−Cλ′+Aλ′C+A¯=signdfdu|u=ωk2.The proof is completed.According to Lemmas1–3, the following theorem is obtained.Theorem 4.
For equation (29), the following conclusions are established:(i)
If (S1) is true, thensigndReλ/dτ2|τ2=τ21j>0.(ii)
If (S2) is true, thensigndReλ/dτ2|τ2=τ22j<0 and signdReλ/dτ2|τ2=τ23j>0.(iii)
If (S3) is true, thensigndReλ/dτ2|τ2=τ2kj>0k=4,6, dReλ/dτ2|τ2=τ25j<0.(iv)
If (S4) is true, thendReλ/dτ2|τ2=τ2kj<0k=7,9, dReλ/dτ2|τ2=τ2kj>0k=8,10.Proof.
(i)
If (S1) is true, then equation (34) has one positive root and three negative roots or one positive root, one negative root, and a pair of conjugate imaginary roots, and fu is monotonically increasing at the positive root u1, so we get f′u1>0. Therefore, the conclusion is true.(ii)
If (S2) is true, then equation (34) has two positive roots (u2,u3 and u2<u3) and two negative roots or has two positive roots and a pair of conjugate imaginary roots. Function fu decreases monotonically at the positive root u2 and increases monotonically at the positive root u3, so f′u2<0 and f′u3>0 are obtained. Therefore, the conclusion is true.(iii)
If (S3) is true, equation (34) has three different positive real roots (u4,u5,u6 and u4<u5<u6) and one negative root. Function fu increases monotonically at u4, decreases monotonically at u5, and increases monotonically at u6. So, we get f′u4>0, f′u5<0 and f′u6>0. Therefore, the conclusion is true.(iv)
If (S4) holds, equation (34) has four different positive real roots (ui,i=7,8,9,10 and u7<u8<u9<u10). Function fu is monotonically decreasing at u7 and u9, and monotonically increasing at u8 and u10. So we obtain f′u7<0, f′u8>0, f′u9<0 and f′u10>0. Therefore, the conclusion is true.According to the analysis in references (see [18–21]), Lemma 4 and Theorem 5 are obtained on the basis of Theorems 3 and 4.Lemma 4.
For equation (25), when R0>1 and R1>1, the following conclusions are obtained:(i)
If (H1) and (S1) are true, then forτ2∈0,τ21,0, all the roots have a negative real part. There is a pair of pure imaginary roots at τ2=τ21,0. There is at least one root with a positive real part at τ2>τ21,0.(ii)
If (H1) and (S2) are true, then there exists a positive integerθ that makes τ23,0<τ22,0<τ23,1<τ22,1<⋯<τ23,θ−2<τ22,θ−2<τ23,θ−1<τ23,θ<τ22,θ−1. This indicates that the system goes through θ stability switches from stable state to an unstable state. All the roots have negative real parts at τ2∈τ22,j,τ23,j+1τ22,−1=0,j=−1,0,…,θ−1. There are two pairs of pure imaginary roots at τ2=τ22,j or τ2=τ23,jj=0,1,…,θ−1, and the other roots have negative real parts. There is at least one root with a positive real part at τ2∈τ23,j,τ22,jj=−1,0,…,θ−1 or τ2>τ23,θ.(iii)
If (S2) is true and (H1) is not, then there exists a positive integerθ that makes τ22,0<τ23,0<τ22,1<τ23,1<⋯<τ22,θ−1<τ23,θ−1<τ23,θ<τ22,θ. This means that the system goes through θ stability switches from an unstable state to a stable state. There is at least one root with a positive real part at τ2∈τ23,j,τ22,j+1τ23,−1=0,j=−1,0,…,θ−2 or τ2>τ23,θ−1. There are two pairs of pure imaginary roots at τ2=τ22,j or τ2=τ23,jj=0,1,…,θ−1, and the other roots have negative real parts. All the roots have negative real parts at τ2∈τ22,j,τ23,jj=−1,0,…,θ−1.(iv)
If (H1) and (S3) are true, or (H1) and (S4) are true, then the system has at least one stability switch.(v)
If (H1) and (S5) are true, then for anyτ2>0, all roots have negative real parts.Theorem 5.
SupposeR0>1 and R1>1.(i)
If (H1) and (S1) are true, then the equilibriumQ∗ is locally asymptotically stable for any τ2∈0,τ21,0. When τ2=τ21,0, system (1) appears Hopf bifurcation at the equilibrium Q∗.(ii)
If (H1) and (S2) are true, whenτ2∈0∪τ22,j,τ23,j+1τ22,−1=0,j=−1,0,…,θ−2, the equilibrium Q∗ is locally asymptotically stable. When τ2∈τ23,j,τ22,jj=−1,0,…,θ−1 or τ2>τ23,θ, the equilibrium Q∗ is unstable. When τ2>τ23,m or τ2∈τ22,j∪τ23,j∪τ23,m−1j=−1,0,…,θ−2, system (1) appears Hopf bifurcation at Q∗.(iii)
If (S2) is true and (H1) is not. Whenτ2∈τ22,j,τ23,jj=−1,0,…,θ−1, the equilibrium Q∗ is locally asymptotically stable. When τ2∈0∪τ23,j,τ22,j+1j=−1,0,…,θ−1 or τ2>τ23,θ, Q∗ is unstable. When τ2>τ23,m or τ2∈τ22,j∪τ23,jj=−1,0,…,θ−1, Hopf bifurcation occurs at Q∗.(iv)
If (H1) and (S3), or (H1) and (S4) are true, then there is a stability switch atτ2=τ2∗. When τ2∗ satisfies ωτ2∗∈ωkj,k=4,5,6,7,8,9,10, Hopf bifurcation occurs at Q∗.(v)
If (H1) and (S5) are true, thenQ∗ is locally asymptotically stable for any τ2>0.


## 4.2.2. The Caseτ1>0 and τ2>0

In order to discuss the possible Hopf bifurcation atτ1>0 and τ2>0, we regard τ2 as a parameter of τ1. According to the analyses and conclusions of (see [19,22]), system (1) is stable at τ1=0 and τ2>0, there is τ¯1τ2, and the system reaches stability at τ1∈0,τ¯1τ2. Therefore, the following theorem can be obtained:Theorem 6.
Let us assume thatR0>1 and R1>1.(i)
If (H1) and (S1) are true, then for anyτ2∈0,τ21,0, there exists τ¯1τ2, and the equilibrium Q∗ is locally asymptotically stable at τ1∈0,τ¯1τ2.(ii)
If (H1)and (S2) are true, there exists positive integerθ, when τ2∈0∪τ22,j,τ23,j+1τ22,−1=0,j=−1,0,…,θ−2, there exists τ¯1τ2, and the equilibrium Q∗ is locally asymptotically stable at τ1∈0,τ¯1τ2.(iii)
If (H1) and (S5) hold, then for anyτ2>0, there exists τ¯1τ2, and the equilibrium Q∗ is locally asymptotically stable at τ1∈0,τ¯1τ2.According to the analysis (see [23]), when τ2 is fixed, the characteristic equation of system (1) only depends on the change of time delay τ1. Therefore, we fix τ2 as a constant τ2=τ20 and study the stability of the system when τ1 increases from zero. In this case, the characteristic equation (25) becomes(43)D¯λ,τ1,τ20=Aλ,τ1,τ20e−λτ1+Bλ,τ1,τ20e−λτ20+Cλ,τ1,τ20=0,where(44)Aλ,τ1,τ20=R¯3τ1,τ20λ3+R¯2τ1,τ20λ2+R¯1τ1,τ20λ+R¯0τ1,τ20,Bλ,τ1,τ20=P¯2τ1,τ20λ2+P¯1τ1,τ20λ+P¯0τ1,τ20,Cλ,τ1,τ20=λ4+Q¯3τ1,τ20λ3+Q¯2τ1,τ20λ2+Q¯1τ1,τ20λ+Q¯0τ1,τ20.Assuming that the characteristic equation has pure imaginary rootλ=ωiω>0, we substitute it into the characteristic (43) and then separate the real parts and imaginary parts to obtain that(45)sinωτ1=R¯33ω3−R¯1ωP¯0−P¯2ω2cosωτ20+P¯1ωsinωτ20+ω4−Q¯2ω2+Q¯0R¯32ω6+R¯22−2R¯1R¯3ω4+R¯12−2R¯0R¯2ω2+R¯02+R¯0−R¯2ω2P¯2ω2−P¯0sinωτ20+P¯1ωcosωτ20−Q¯3ω3+Q¯1ωR¯32ω6+R¯22−2R¯1R¯3ω4+R¯12−2R¯0R¯2ω2+R¯02,cosωτ1=R¯33ω3−R¯1ωP¯2ω2−P¯0sinωτ20+P¯1ωcosωτ20−Q¯3ω3+Q¯1ωR¯32ω6+R¯22−2R¯1R¯3ω4+R¯12−2R¯0R¯2ω2+R¯02+R¯2ω2−R¯0P¯0−P¯2ω2cosωτ20+P¯1ωsinωτ20+ω4−Q¯2ω2+Q¯0R¯32ω6+R¯22−2R¯1R¯3ω4+R¯12−2R¯0R¯2ω2+R¯02.According toe−iωτ1=−Cωi,τ1,τ20+Bωi,τ1,τ20e−iωτ2/Aωi,τ1,τ20, it follows that(46)Cωi,τ1,τ20+Bωi,τ1,τ20e−iωτ2Aωi,τ1,τ202=1.Through further calculations, we get the following equation:(47)F¯ωi=Cωi,τ1,τ20+Aωi,τ1,τ20e−iωτ22−Bωi,τ1,τ202=ω8+h6τ1,τ20ω6+h5τ1,τ20ω5+h4τ1,τ20ω4+h3τ1,τ20ω3+h2τ1,τ20ω2+h1τ1,τ20ω+h0τ1,τ20=0,where(48)h6τ1,τ20=Q¯32−R¯32−2Q¯2−2P¯2cosωτ20,h5τ1,τ20=2P¯1sinωτ20−2P¯2Q¯3sinωτ20,h4τ1,τ20=P¯22+Q¯22+2Q¯0−R¯22+2P¯0cosωτ20−2P¯1Q¯3cosωτ20−2Q¯1Q¯3+2R¯1R¯3+2P¯2Q¯2cosωτ20,h3τ1,τ20=2P¯2Q¯1sinωτ20−2P¯1Q¯2sinωτ20+2P¯0Q¯3sinωτ20,h2τ1,τ20=P¯12+Q¯12−R¯12+2R¯0R¯2−2Q¯0Q¯2−2P¯0P¯2−2P¯2Q¯0cosωτ20+2P¯1Q¯1cosτ20−2P¯0Q¯2cosωτ20,h1τ1,τ20=2P¯1Q¯0sinωτ20−2P¯0Q¯1sinωτ20,h0τ1,τ20=Q¯02+P¯02−R¯02+2P¯0Q¯0cosωτ20.In order to study the influence of two time delays on the stability of the system, the following five criteria are given according to reference (see [23]):(U1) For anyτ1,τ20∈R02, there is A0,τ1,τ20+B0,τ1,τ20+C0,τ1,τ20≠0, that is, λ=0 is not the root of the characteristic equation (25).(U2) Ifλ=ωi, then for any τ1≥0, inequality Aωi,τ1,τ20+Bωi,τ1,τ20+Cωi,τ1,τ20≠0 can be obtained.(U3) For anyτ1, limλ⟶∞supAλ,τ1,τ20/e−λτ20Bλ,τ1,τ20+Cλ,τ1,τ20:Reλ≥0<1 holds.(U4) For anyτ1, the number of roots of the following function F¯ is limited.(49)F¯=e−ωiτ20Biω,τ1,τ20+Ciω,τ1,τ202−Aiω,τ1,τ202.(U5) As long asF¯=0 has positive roots ωτ1,τ20, then it is continuously differentiable with respect to τ1.Next, we prove that for anyτ1∈0,τmax, where τmax is the maximum value that enables the equilibrium to exist, each item in the characteristic (43) conforms to the above five criteria. Suppose that (U2) is true, and then we prove that the other four criteria are true.By calculations, we obtainF¯0,τ1,τ20=Q¯0+P¯02−R¯02, limω⟶+∞F¯ω,τ1,τ20=+∞. For the convenience of expression, we give conditions.(G1)Q¯0+P¯0<R¯0.If (G1) holds, then (U1) is correct andF¯ω,τ1,τ20=0 has at least one positive root. Next, when τ1∈0,τmax, obviously, (U3) holds. Then, on the basis of (U2), we can conclude that (U4) is true. Finally, we verify (U5). Suppose that ω0 is a positive real root of F¯ω,τ1,τ20=0 and ω0,τ1,τ20 is a point in the domain of F¯ω,τ1,τ20=0. Obviously, F¯ω′ and F¯τ1′ exist and they are continuous in some regions of ω0,τ1,τ20. In addition, there is F¯ω′ω0,τ1,τ20≠0. Therefore, according to the implicit function theorem, the verification of (U5) is true.IfI˜∈ℝ+ is the set of roots of F¯=0, for any ω∈I˜, we obtain F¯ω,τ1,τ20=0. Let us define θτ1∈0,2π, we get cosθτ1 and sinθτ1.Let us define the mapping(50)Sn:I˜⟶ℝ,Snτ1=τ1−θτ1+2nπωτ1,τ20,n∈ℕ0,τ1∈0,τmax.There exists a positive rootωτ1 of F¯ω,τ1,τ20=0 in I˜. Obviously, Sn is continuously differentiable with respect to τ1. Therefore, according to [23], the following theorem can be obtained:Theorem 7.
Suppose that one of the following conditions is correct.(i)
For any fixedτ20∈0,τ1,0, (H1) and (S1) hold.(ii)
For any fixedτ20∈0,τ3,0, (H1) and (S2) hold.
Then we obtain the following conclusions: ifR0>1 and R1>1, (H1) and (G1) are true, then F¯ω,τ1,τ20=0 has a pair of pure imaginary roots ±iωτ1∗,τ20, where ωτ1∗,τ20>0. For certain n∈N0, when τ1∗∈0,τmax, if Snτ1∗=0, then the root of the characteristic equation crosses the imaginary axis from left to right (from right to left) at Δτ1∗>0Δτ1∗<0, where Table 1(51)Δτ1∗=signdReλdτ|λ=iωτ1∗,τ20=signF¯ω′ωτ1∗,τ20,τ1∗,τ20signdSnτdτ|τ=τ1∗Table 1
Parameter of the system (1).
ParametersUnitsDefmitionsb−1mm3The maximum carrying capacity of tumor cellβDay−1 mm−3The infection rate of tumor cellsC1Day−1The clearance rate of CTLs to uninfected tumor cellsC2Day−1The clearance rate of CTLs to infected tumor cellsd1Day−1 mm−3The mortality of infected tumor cellsd2Day−1 mm−3The mortality of virusd3Day−1 mm−3The mortality of CTLsk1Day−1 mm−3The maximum growth rate ofuninfected tumor cellsk2Day−1 mm−3The average rate at which the infected tumor cells release the virusk3Day−1 mm−3Corresponding to the infection of uninfected tumor cells rate of CTLs simulation.k4Day−1 mm−3Corresponding to the infection of infected tumor cells rate of CTLs simu lationq—A constant greater than 0


## 5. Numerical Simulations

In this section, the main conclusions drawn in the paper are verified through a series of numerical simulations, and the influence of two time delays on the stability of the system (1) is discussed. Some of the parameters of numerical simulations are selected with reference to [8, 16].Whenτ1=0andτ2=3, the following parameters are selected:(52)k1=0.206,k2=0.0001,k3=0.00016,k4=0.00016,β=0.39,b=1/2139.3,c1=0.0000673,c2=0.0000673,d1=0.703,d2=0.285,d3=0.3.Through calculation, the virusfree equilibrium Q01875,0,0,378.1618 and R0=0.3552<1 are obtained. Applying Theorem 1, we obtain that the virusfree equilibrium Q0 is globally asymptotically stable. Finally, Figure 1 verifies this conclusion.Figure 1
Global asymptotic stability of equilibriumQ0 at τ1=0andτ2=3.
(a)(b)Whenτ1=τ2=0, the parameters we selected are as follows:(53)k1=0.206,k2=0.0005,k3=0.00015,k4=3.2136,β=0.39,b=12139.3,c1=0.00673,c2=0.008749,d1=0.703,d2=0.285,d3=0.3,q=5.By simple calculations, we obtain that the positive equilibrium of system (1) is Q∗=1199.7,0.0374,6.5535×10−5,13.4394, R0=1.8994, and R1=1.4696. According to condition (i) of theorem 3, the equilibrium is locally asymptotically stable, and the numerical simulation is shown in Figure 2.Figure 2
Local asymptotic stability of equilibriumQ∗ at τ1=τ2=0.
(a)(b)Whenτ1=0 and τ2=3, (29) has no positive real part root. Moreover, H1=7.5036×10−4>0, Γ=1.8897×10−5>0, and N5=−0.0279<0. So, (H1) and (S5) hold. According to (v) in theorem 5, the equilibrium is locally asymptotically stable. As is shown in Figure 3, the number of uninfected tumor cells eventually tends to be stable.Figure 3
Local asymptotic stability ofQ∗ at τ1=0 and τ2=3.
(a)(b)We setτ1=4 and τ2=9, we still take the set of parameters above. By simple calculations, it can obtain Q∗1199.7,0.0374,6.5535×10−5,13.4394. We can compute H1=7.5096×10−4>0, Γ=1.8894×10−5>0, and N5=−0.0279<0. So, (H1) and (S5) hold, which satisfies the condition (iii) of theorem 6. Therefore, according to theorem 6, the equilibrium is locally asymptotically stable. As shown in Figure 4, the number of uninfected tumor cells gradually tends to stabilize Figure 5.Figure 4
Local asymptotic stability ofQ∗ at τ1=4 and τ2=9.
(a)(b)Figure 5
Local asymptotic stability ofQ∗ at τ1=0 and τ2=28.
(a)(b)Then we discuss the partial Hopf bifurcation with a time delay. Whenτ1=0 and τ2 is used as the bifurcation parameter, the selected parameters are k1=0.206, k2=0.0005, k3=0.00015, k4=3.2136, β=0.39, b=1/2139.3, c1=0.0673, c2=0.008749, d1=0.703, d2=0.285, d3=0.3, and q=5. When τ2=30 (τ2>τ2∗≈29), we obtain H1=−8.8217×10−4<0. According to condition (iii) of Theorem 6, the numerical simulation is shown in Figure 6. The number of uninfected tumor cells fluctuates periodically.Figure 6
Local asymptotic stability ofQ∗ at τ1=0 and τ2=30.
(a)(b)Next, we discuss the Hopf bifurcation with two time delays.τ2 is fixed and τ1 is taken as the bifurcation parameter. The parameters selected are k1=0.16, k2=0.9, k3=0.005, k4=0.39, β=0.59, b=1/2139.3, c1=0.00673, c2=0.9749, d1=0.703, d2=0.203, d3=0.3, and q=5. When τ2=20 is fixed and τ1=3.5 (τ1<τ1∗≈4.16), by calculations, it obtains Q∗=30.0303,0.3842,1.7035,7.7450, R0=6.7561, R1=19.3875, Γ=−2.5804×1011<0, and l3=−1.4581×10−5<0. As shown in Figure 7, the number of uninfected tumor cells tends to be ultimately stable. When τ2=20 is fixed and τ1=4.2 (τ1>τ1∗≈4.16), as shown in Figures 8(a) and 8(b), the number of uninfected tumor cells always fluctuates periodically. Figures 7 and8 demonstrate the condition (i) of Theorem 7, that there exists a Hopf bifurcation with two time delays.Figure 7
Whenτ2=20 is fixed and τ1=3.5 ,Q∗ is locally asymptotically stable.
(a)(b)Figure 8
Whenτ2=20 is fixed and τ1=4.2, there is a bifurcation periodic solution.
(a)(b)The parameters selected arek1=0.2226, k2=0.2524, k3=0.17, k4=9.2136, β=0.229, b=1/2139.3, c1=0.592, c2=0.58, d1=0.263,d2=0.25, d3=0.6, q=5, and τ1=0. It can obtain R0=1.6974 and R1=43.7981. When τ2=0, the positive equilibrium Q∗ is unstable. With the increase of τ2, the system switches from unstable to stable. The numerical simulations are shown in Figures 9 and 10, where Figure 10 is the bifurcation diagram showing the stability of positive equilibrium Q∗ changing with τ2. These images demonstrate Theorem 5.Figure 9
Whenτ1=0, the stability of equilibrium changes with the increase of τ2.
(a)(b)(c)Figure 10
Whenτ1=0, the bifurcation image with τ2.The parameters selected arek1=0.16, k2=0.05, k3=0.009, k4=3.2136, β=0.26, b=1/2139.3, c1=0.0673, c2=0.39, d1=0.503, d2=0.185, d3=0.3, q=5, and τ1=0. R0=1.6974 and R1=43.7981 are obtained. When τ2=0, the positive equilibrium Q∗ is stable. According to Theorem 5, with the increase of τ2, the system appears to switch phenomenons from stable to unstable. The numerical simulations are shown in Figures 11 and 12, where Figure 12 is a bifurcation diagram showing the change of the stability of Q∗ with τ2. The following images prove Theorem 5.Figure 11
The equilibrium changes from stable to unstable whenτ1=0 and τ2 increases.
(a)(b)(c)(d)Figure 12
The bifurcation image withτ2.


## 6. Conclusions

In this paper, we study a model with two time delays of oncolytic virus infection with saturation incidence. The first delay is the time from oncolytic viruses entering tumor cells to gene replication. The second delay is the time from viruses entering tumor cells to infected tumor cells releasing new virus particles. In the previous studies of oncolytic virus infection models, the infection rate was linear. In this paper, we further consider the nonlinear infection rate. The global asymptotic stability and local asymptotic stability of the uninfected equilibrium, as well as the conditions for the local asymptotic stability of the virus-infected equilibrium, are discussed, and the effect of time delay on the stability of the equilibrium is also discussed. Through theoretical analysis and numerical simulation, it is obtained that the local stability of the uninfected equilibrium is independent of the size of the time delay, but the change in time delay can promote or destroy the stability of the virus infection equilibrium. That is, the time delay will destroy the stable state of tumor treatment, and even produce Hopf bifurcation and lead to periodic oscillation.When only one time delay is considered, that isτ1=0 and τ2>0, taking τ2 as the bifurcation parameter, the conditions for generating Hopf bifurcation are obtained. It is verified that the time delay from the virus entering the tumor cells to the tumor cells cracking to produce new virus particles makes the stable equilibrium of the system become an unstable equilibrium, that is, there is a critical value of τ2∗. When τ2<τ2∗, the virus-infected equilibrium is locally asymptotically stable; when τ2>τ2∗, the virus-infected equilibrium is unstable. At τ2∗, system (1) generates Hopf bifurcation, which indicates that the time delay of oncolytic virus infection will affect the stability of tumor treatment. Therefore, reasonable control of time delay plays a key role in effectively curbing the growth and recurrence of tumor cells.When considering the two time delays of the system, that isτ1>0 and τ2>0, taking τ1 as the parameter of τ2 and τ1 as the bifurcation parameter, the conditions for generating Hopf bifurcation are discussed. When both delays exist, there are critical values τ1∗ and τ2∗. When τ2∈0,τ2∗, τ1∈0,τ1∗, the infection equilibrium is locally asymptotically stable, but when τ1>τ1∗, the infection equilibrium becomes unstable, and Hopf bifurcation occurs at the critical value τ1∗. When τ2 is fixed, the delay time of τ1 is small in a certain range, and the system is in a stable state. As the delay time increases, the system changes to an unstable state. In other words, when the time for tumor cell lysis to produce new virions is fixed, the faster the replication rate of oncolytic viruses is, the higher the concentration in the cell is, the more virions are released and the immune system is further activated. In this case, the system is in a stable state. When the amount of time delay increases, the replication rate of viruses is relatively slow, the number of tumor cells lysis and death due to infection with oncolytic virus is reduced, and the tumor immune system enters the escape stage in the tumor immune editing process. The tumor cells in the system are not effectively contained, and the system is in an unstable state. At the same time, the replication time of the oncolytic virus cannot be blindly accelerated, and the excessive intracellular concentration of the oncolytic virus cannot be pursued, because oncolytic viruses not only kill tumor cells, but also reduce the immune system response caused by tumor cell stimulation, which may also lead to the failure of tumor cell treatment. That is to say, too little concentration of oncolytic virus makes it difficult to effectively remove tumor cells, and a high concentration of oncolytic virus cannot fully enhance the immune system and may even have pathological effects. In the process of immune system interaction, different time delays can lead to different tumor states, so the selection of two time delays plays an important role in tumor cell therapy.In modern medical research, genetic engineering techniques such as DNA recombination technology can change or obtain more efficient oncolytic viruses and delay the time of infection. By controlling infection delay, the efficiency of viral infection can be further improved, which provides an effective method for the treatment of oncolytic viruses and has important application value.This paper explores the effect of discrete-time delays on the tumor treatment interaction system. In the future, we can extend the deterministic system to a stochastic system to more comprehensively and realistically study the stochastic dynamics of oncolytic virus infection models, which will enrich our understanding of the influence of biological delay in the process of virus infection and its interaction with immune response.



---

*Source: 3713439-2022-07-05.xml*