# Hopf Bifurcation Analysis of a Two-Delay HIV-1 Virus Model with Delay-Dependent Parameters

Authors: Yu Xiao; Yunxian Dai; Jinde Cao
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this paper, a two-delay HIV-1 virus model with delay-dependent parameters is considered. The model includes both virus-to-cell and cell-to-cell transmissions. Firstly, immune-inactivated reproduction rateR0 and immune-activated reproduction rate R1 are deduced. When R1>1, the system has the unique positive equilibrium E∗. The local stability of the positive equilibrium and the existence of Hopf bifurcation are obtained by analyzing the characteristic equation at the positive equilibrium with the time delay as the bifurcation parameter and four different cases. Besides, we obtain the direction and stability of the Hopf bifurcation by using the center manifold theorem and the normal form theory. Finally, the theoretical results are validated by numerical simulation.


---

## Body

## 1. Introduction

AIDS is a very dangerous infectious disease caused by HIV-1 virus which attacks the human immune system. At present, it has been proved that there are two different mechanisms for the spread of HIV-1 virus in the host: virus-cell transmission and cell-cell transmission [1]. A great number of studies have considered the above two contagion mechanisms [1–5]. Wang and Zou [4] considered an HIV-I model with humoral immunity, which is a typical virus-to-cell transmissions. In 2017, Lin et al. [5] proposed a cell-to-cell transmission model and described global threshold dynamics of the model.In recent years, delay differential equation has attracted extensive attention worldwide. It has important applications in many fields such as physics, information, economy, and biomathematics. Depending on the different circumstances, many differential equation models with single delay or multiple delays have been proposed and studied deeply [4, 6–10]. Dong et al. [6] studied the dynamics of the tumor immune system interaction model and investigated the existence of Hopf bifurcation with two time delays as bifurcation parameters. In [7], the authors discussed the influence of awareness coverage and time delays on infectious diseases and found that the endemic equilibrium existed a Hopf bifurcation in both delayed and nondelayed system. A two-delay model with Holling II functional response and stage structure is considered in [8]. The authors have investigated a predator-prey model with a class of Beddington–DeAngelis functional response and two delays in [9]. Recently, in [10], a two-delay HIV-1 virus model with virus-to-cell and cell-to-cell transmission is considered as follows:(1)dxtdt=rxt1−xtxM−β1xtvt−β2xtyt,dytdt=β1αxt−τ1vt−τ1+β2αxt−τ1yt−τ1−ayt,dvtdt=kyt−uvt−pvtzt,dztdt=cvt−τ2zt−τ2−bzt,where xt,yt,vt, and zt represent the susceptible cells, infected cells, virus, and B cells, respectively. xM is the number of carrying target cells. β1 and β2 denote infection rates of virus-to-cell transmission and cell-to-cell transmission. r and b are the growth and death rate of susceptible cells. a and u stand for the mortality rates of infected cells and virus, respectively. k is the number of free virus particles produced by per infected cell. α is the surviving probability the time period from t−τ1 to t. pvtzt describes the virus killed by B cells, and cvtzt represents the new B cells produced when stimulated by antigen. The time delay τ1 stands for the time between virus entering into a cell and producing new virus or the time between infected cells spreading virus into uninfected cells and producing new virus. And τ2 is the time that the HIV-1 virus stimulates the production of B cells.In [10], the coefficients of equations are considered to be constants independent of the time delay. The authors investigated the stability of positive equilibrium and the existence of Hopf bifurcation. The explicit formula for determining the direction of Hopf bifurcation and the stability of bifurcating periodic solutions was derived. However, in practical situations, virus-cell transmission or cell-cell transmission is not instantaneous but takes some time defined by the average evaluation period τ. Moreover, the growth of virus at t is decided by the number of infected cells at t−τ and still alive at t. On this occasion, it is inappropriate to regard survival rate as a constant. Let m be death rate of infected cells, then the survival rate is e−mτ from t−τ to t. Adding time delays to parameters is helpful to further explore the complex dynamic behavior of delay differential equations and their generation mechanism. Based on Sun and Li’work, considering that the two contagion processes are not instantaneous and the probability of cell survival is proportional to the number of cells that can survive to time t at t−τ1, we propose the following delay model:(2)dxtdt=rxt1−xtxM−β1xtvt−β2xtyt,dytdt=β1e−mτ1xt−τ1vt−τ1+β2e−mτ1xt−τ1yt−τ1−ayt,dvtdt=kyt−uvt−pvtzt,dztdt=cvt−τ2zt−τ2−bzt,where e−mτ1 is the probability of surviving the time period from t−τ1 to t in case m is the death rate of infected but not yet virus-producing cells.The initial conditions for system (2) are given as(3)xθ=ϕ1θ,yθ=ϕ2θ,vθ=ϕ3θ,zθ=ϕ4θ,ϕiθ≥0,θ∈−τ,0,ϕi0>0,i=1,2,3,4.A number of scholars investigate the delay differential equations with delay-dependent parameters [11, 12]. Song et al. [11] studied a delayed viral infection model with lytic immune response. The characteristic equations of the system are like P0λ,τ+P1λ,τe−λτ=0. They discussed the influence of delay on the stability of the equilibrium and the local and global asymptotic stability of the disease-free equilibrium. Li and Ma introduced a method for determining the stability of the characteristic equations P0λ,τ+P1λ,τe−λτ+P2λ,τe−2λτ=0 with delay-dependent parameters in [13]. Reference [14] presented a practical geometric method to study the stability switching properties of the characteristic equation P0λ,τ1+P1λ,τ1e−λτ1+P2λ,τ1e−λτ1+τ2=0 which may result from a stability analysis of a model with two time delays and delay-dependent parameters that depend only on one of the time delay. Jiang and Guo [15] studied a model with double delays and a delay-dependent parameter considering the interaction between nutrients and plankton. The authors took the time delay as the parameter to carry out the dynamic analysis of the system, including the equilibrium stability and Hopf bifurcation existence by using the method in [14].This paper’s goal is to delve into the local stability and Hopf bifurcation analysis of a two-delay HIV-1 virus model with delay-dependent parameters. The luminescent spots are as follows: (1) in this paper, two different spreads of HIV-1 virus are studied, and we will introduce time delay into the coefficients to be more realistic; (2) the conditions of the local stability of positive equilibrium and the existence of Hopf bifurcation are discussed with the time delay as the bifurcation parameters; and (3) the characteristic equation isP0λ,τ+P1λ,τe−λτ+P2λ,τe−2λτ=0 when τ1=τ2=τ, and the stability is determined by the geometric stability switch criterion [13, 14]; and (4) the influences of different delay values on system stability are investigated.The rest of this paper is organized as follows. In Section2, the existence of equilibrium is given. By analyzing the characteristic equation, we discuss the stability of the positive equilibrium and existence of Hopf bifurcation in four different cases. In Section 3, we consider τ2 as a parameter and determine the direction and stability of Hopf bifurcation by using the center manifold theorem and the normal form theory. In Section 4, we perform numerical simulations to illustrate our results. Finally, the conclusions of the paper are given in Section 5.


## 2. Stability of the Positive Equilibrium and Existence of Hopf Bifurcation

Similar to [10], we can derive immune-inactivated reproduction rate R0=kβ1+uβ2/auxMe−mτ1 and immune-activated reproduction rate R1=ckrxMkβ1+uβ2/xMkβ1+uβ2kbβ1+ubβ2+ckrauemτ1=R0/1+bkβ1+uβ2R0/ckr, and then the equilibria of system (2) are as follows:(i)
R1<R0<1, then system (2) only has an infection-free equilibrium E0=xM,0,0,0.(ii)
R1<1<R0, then system (2) has a no-immune response equilibrium E1=x1,y1,v1,0 except for E0, where(4)x1=auemτ1kβ1+uβ2,y1=urkβ1+uβ21−1R0,v1=krkβ1+uβ21−1R0.(iii)
1<R1<R0, then system (2) has an immunity-activated equilibrium E∗=x∗,y∗,v∗,z∗ except for E0 and E1, where(5)x∗=aub+pbz∗emτ1kbβ1+β2ub+pbz∗,y∗=ub+pbz∗ck,v∗=bc,andz∗ is a positive, real root of the following quadratic equation:(6)p2β22bxMzt2+ackrpemτ1+β2pxM−ckr+2kbβ1+2ubβ2zt+xMkβ1+uβ2kbβ1+ubβ2+ckrauemτ11−R1=0.The characteristic equation of system (2) at Ex,y,v,z is(7)detλI−A−Be−λτ1−Ce−λτ2=0,where(8)A=r−2rxxM−β1v−β2y−β2x−β1x00−a000k−u−pz−pv000−b≜a11a12a1300a22000a32a33a34000a44,B=0000β1e−mτ1v+β2e−mτ1yβ2e−mτ1xβ1e−mτ1x000000000≜0000b21b22b23000000000,C=00000000000000czcv≜00000000000000c43c44,that is,(9)λ4+A1λ3+A2λ2+A3λ+A4+B1λ3+B2λ2+B3λ+B4e−λτ1+C1λ3+C2λ2+C3λ+C4e−λτ2+D1λ2+D2λ+D3e−λτ1+τ2=0,where(10)A1=−a11+a22+a33+a44,A2=a11a22+a11a33+a11a44+a22a33+a22a44+a33a44,A3=−a11a22a33+a11a22a44+a11a33a44+a22a33a44,A4=a11a22a33a44,B1=−b22,B2=a11b22−a12b21−a32b23+a33b22+a44b22,B3=a11a32b23−a11a33b22+a12a33b21−a13a32b21−a11a44b22+a12a44b21+a32a44b23−a33a44b22,B4=−a11a32a44b23+a11a33a44b22−a12a33a44b21+a13a32a44b21,C1=−c44,C2=a11c44+a22c44+a33c44−a34c43,C3=−a11a22c44−a11a33c44+a11a34c43−a22a33c44+a22a34c43,C4=a11a22a33c44−a11a22a34c43,D1=b22c44,D2=−a11b22c44+a12b21c44+a32b23c44−a33b22c44+a34b22c43,D3=−a11a32b23c44+a11a33b22c44−a11a34b22c43−a12a33b21c44+a12a34b21c43+a13a32b21c44.Then, the local stability of the positive equilibriumE∗ and the existence of local Hopf bifurcation are discussed in the following four cases.Case 1.
τ1=0,τ2=0.
Equation (9) at the equilibrium E∗x∗,y∗,v∗,z∗ reduces to(11)λ4+A1+B1+C1λ3+A2+B2+C2+D1λ2+A3+B3+C3+D2λ+A4+B4+C4+D3=0.
Let(12)H1R1>1,H2Δ1=A1+B1+C1>0,Δ2=A1+B1+C1A2+B2+C2+D1−A3+B3+C3+D2>0,Δ3=A1+B1+C1A2+B2+C2+D1A3+B3+C3+D2−A1+B1+C1A4+B4+C4+D3−A3+B3+C3+D22>0,Δ4=A4+B4+C4+D3>0.
From the Routh–Hurwitz criterion, we have the following theorem.Theorem 1.
If (H1) and (H2) hold, the positive equilibrium E∗ is locally asymptotically stable for τ1=0andτ2=0.Case 2.
τ1>0andτ2=0.
The characteristic equation of system (2) is(13)λ4+A1+C1λ3+A2+C2λ2+A3+C3λ+A4+C4+B1λ3+B2+D1λ2+B3+D2λ+B4+D3e−λτ1=0.
Denote(14)F1=A1+C1,F2=A2+C2,F3=A3+C3,F4=A4+C4,G1=B1,G2=B2+D1,G3=B3+D2,G4=B4+D3,then equation (13) is reduced to(15)Pλ,τ1+Qλ,τ1e−λτ1=0,where Pλ,τ1=λ4+F1λ3+F2λ2+F3λ+F4andQλ,τ1=G1λ3+G2λ2+G3λ+G4.
Denoteτmax=1/mlnxMkβ1+uβ2ckr−kbβ1−ubβ2/ckrau and I=0≤τ<τmax:there isωτ>0such thatFω,τ=0. According to the geometrical criterion established by Beretta and Kuang [16], we can easily verify the following conditions for ω>0,τ1∈I:(a)
P0,τ1+Q0,τ1=F4+G4≠0.(b)
Piω,τ1+Qiω,τ1=ω4−F2+G2ω2+F4+G4−F1+G1ω3−F3+G3iω+G3ωi≠0.(c)
limλ⟶∞Qλ,τ1/Pλ,τ1=limλ⟶∞G1λ3+G2λ2+G3λ+G4/λ4+F1λ3+F2λ2+F3λ+F4=k=1/λlimk⟶0G1k+G2k2+G3k3+G4k4/1+F1k+F2k2+F3k3+F4k4=0.(d)
Piω,τ12=ω8+F12−2F2ω6+F22+2F4−2F1F3ω4+F32−2F2F4ω2+F42,Qiω,τ12=G12ω6+G22−2G1G3ω4+G32−2G2G4ω2+G42,Fω,τ1=Piω,τ12−Qiω,τ12=ω8+F12−2F2−G12ω6+F22+2F4−2F1F3−G22+2G1G3ω4+F32−2F2F4−G32+2G2G4ω2+F42−G42.We can get easilyFω,τ1 has finite roots.(f)
Each positive rootωτ1 of Fω,τ1=0 is continuous and differentiable.
Letλ=iωω>0 be a root of equation (15) and separate the real and imaginary parts, then we get(16)ω4−F2ω2+F4=−G4−G2ω2cosωτ1+G1ω3−G3ωsinωτ1,−F1ω3+F3ω=G4−G2ω2sinωτ1+G1ω3−G3ωcosωτ1.
Then, we can obtain(17)cosωτ1=ω4−F2ω2+F4G2ω2−G4+F3ω−F1ω3G1ω3−G3ωG2ω2−G42+G1ω3−G3ω2,sinωτ1=ω4−F2ω2+F4G1ω3−G3ω−F3ω−F1ω3G2ω2−G4G2ω2−G42+G1ω3−G3ω2.
From (16), it follows that(18)ω4−F2ω2+F42+−F1ω3+F3ω2=G4−G2ω22+G1ω3−G3ω2,which is equivalent to Fω,τ1=0 for (15). Denote h=ω2, then we have Fω,τ1=Gh,τ1=0. Therefore, equation (13) has a pair of pure imaginary roots ±iω1∗ when τ1>0andτ2=0 if ω1∗ is the positive root of Fω,τ1. Then, we denote ωτ1τ1=θτ1+2nπ and Snτ1=τ1−θτ1+2nπ/ωτ1. Hence, iω1∗ is a pure imaginary root of equation (13) if and only if τ1∗ is a zero of Snτ1. We introduce the following theorem [16].Theorem 2.
Assume thatSnτ1∗=0 have some positive roots τ1∗∈I for some n∈N. Then, a pair of simple pure imaginary roots ±iω1τ1∗ exists which crosses the imaginary axis from left to right if δτ1∗>0 and crosses the imaginary axis from right to left if δτ1∗<0, where(19)δτ1∗=signFω1′ω1τ1∗,τ1∗signdSnτ1dτ1|τ1=τ1∗.Since∂Fω1,τ1/∂ω1=2ω1∂Gω12,τ1/∂ω1, equation (19) is equivalent to(20)signdReλdτ|λ=iω1τ1∗=sign2ω1∂G∂ω1ω12,τ1×signdSnτ1dτ1|τ1=τ1∗.Next, we can easily getSn0<0 and Snτ1n∈N is a monotonic increasing function for all τ1∈I. If S0 has no positive root in I, then Sn also does not have positive root in I. And if Snτ1 has a positive root τ1∈I for some n∈N, then there is at least one positive root which satisfies dSnτ1/dτ1τ1>0. We introduce the Hopf bifurcation theorem:Theorem 3.
A single parameter system of formx˙t=Fα,xt where Fα,ϕ has continuous first and second derivatives for α∈R,ϕ∈C. Define L:R×C⟶Rn as Lαψ=Fϕα,0ψ, where Fϕα,0 is the derivation of Fα,ϕ with respect to ϕ at ϕ=0, and define fα,ϕ=Fα,ϕ−Lαϕ. There are two hypotheses:(i)
The linear differential equation (L0) has a pure imaginary characteristic root λ0=iω0≠0 and none of the other roots are multiples of iω0(ii)
Reλ′0≠0
The hypotheses (i) and (ii) imply that there are nonconstant periodic solutions, and a Hopf bifurcation occurs.From all above analysis, we have the following theorem.Theorem 4.
Assume thatR1>1 holds, then(i)
IfS0τ1 has no positive root in I, E∗ is locally asymptotically stable for all τ1>0(ii)
IfS0τ1 has at least one positive root in I, there exists τ1∗ such that E∗ is locally asymptotically stable for all τ1∈0,τ1∗(iii)
If (ii) holds and∂F1/∂ω1ω12τ1∗,τ1∗>0, a Hopf bifurcation occurs at E∗ for τ1=τ1∗Case 3.
τ1=τ2=τ.
DenoteP0λ,τ=λ4+A1λ3+A2λ2+A3λ+A4,P1λ,τ=B1+C1λ3+B2+C2λ2+B3+C3λ+B4+C4,P2λ,τ=D1λ2+D2λ+D3, then equation (9) reduces to(21)P0λ,τ+P1λ,τe−λτ+P2λ,τe−2λτ=0.
Obviously,Plλ,τ,l=0,1,2, is a polynomial with λ, and its coefficients depend on τ. Next, we discuss the existence of pure imaginary roots of (21) by using the method introduced in [17].
Letτ∈0,τmax, we need to verify the following conditions:(a)
degP0λ,τ≥maxdegP1λ,τ,degP2λ,τ(b)
P00,τ+P10,τ+P20,τ≠0(c)
Ifλ=iω,ω∈R, then P0iω,τ+P1iω,τ+P2iω,τ≠0,τ∈R(d)
limλ⟶∞P1λ,τ/P0λ,τ+P2λ,τ/P0λ,τ<1(e)
Fω,τ=R2ω,τ−T2ω,τ−S2ω,τ=0 for each τ has at most a finite number of real zeros where(22)Rω,τ=P0iω,τ2−P2iω,τ2,Sω,τ=ReP1iω,τImP0iω,τ+ImP2iω,τ−ImP1iω,τReP0iω,τ+ReP2iω,τ,Tω,τ=ReP1iω,τReP2iω,τ−ReP0iω,τ+ImP1iω,τImP2iω,τ−ImP0iω,τ.(f)
Each positive rootωτ of Fω,τ=0 is continuous and differentiable in τ whenever it exists
Next, we verify the above six conditions.(a)
It is obviously satisfied.(b)
P00,τ+P10,τ+P20,τ=A4+B4+C4+D3 according to H2, P00,τ+P10,τ+P20,τ>0.(c)
P0iω,τ+P1iω,τ+P2iω,τ=ω4−A2+B2+C2+D1ω2+A4+B4+C4+D3+iω−A1+B1+C1ω2+A3+B3+C3+D2≠0.(d)
limλ⟶∞P1λ,τ/P0λ,τ+P2λ,τ/P0λ,τ=limλ⟶∞B1+C1λ3+B2+C2λ2+B3+C3λ+B4+C4/λ4+A1λ3+A2λ2+A3λ+A4+D1λ2+D2λ+D3/λ4+A1λ3+A2λ2+A3λ+A4=0.(e)
P0iω,τ=ω4−A2ω2+A4+A3ω−A1ω3i,P1iω,τ=−B2+C2ω2+B4+C4+−B1+C1ω3+B3+C3ωi,P2iω,τ=−D1ω2+D3+D2ωi,Rω,τ=ω4−A2ω2+A42+A3ω−A1ω32−D1ω2+D32+D22ω2,Sω,τ=−B2+C2ω2+B4+C4A3ω−A1ω3+D2ω−−B1+C1ω3+B3+C3ωω4−A2ω2+A4−D1ω2+D3,Tω,τ=−B2+C2ω2+B4+C4−D1ω2+D3−ω4−A2ω2+A4+A3ω−A1ω3D2ω−A3ω−A1ω3.SubstitutingRω,τ,Sω,τ,Tω,τ into Fω,τ=R2ω,τ−T2ω,τ−S2ω,τ=0, then we get this condition holds.(f)
The implicit function theorem shows that this condition is true.
Assuming thatλ=iωω>0 is a root of equation (13), then we have(23)ω4−A1iω3−A2ω2+A3iω+A4cosωτ+isinωτ−B1+C1iω3−B2+C2ω2+B3+C3iω+B4+C4+−D1ω2+D2iω+D3cosωτ−isinωτ=0.
Separating the real and imaginary parts, we have(24)cosωτω4−A2+D1ω2+A4+D3+sinωτA1ω3+D2−A3ω=B2+C2ω2−B4+C4,cosωτ−A1ω3+A3+D2ω+sinωτω4−A2−D1ω2+A4−D3=B1+C1ω3−B3+C3ω.
Therefore,ω=ωτ>0 needs to satisfy the following equations:(25)cosωτ=l5ω6+l6ω4+l7ω2+l8ω8+k1ω6+k2ω4+k3ω2+k4=Tω,τRω,τ,sinωτ=l1ω7+l2ω5+l3ω3+l4ωω8+k1ω6+k2ω4+k3ω2+k4=Sω,τRω,τ,where(26)k1=A12−2A2,k2=2A4+A22−D12−2A1A3,k3=A32−D22−2A2A4+2D1D3,k4=A42−D32,l1=B1+C1,l2=A1B2+C2−B3+C3−B1+C1A2+D1,l3=B1+C1A4+D3+B3+C3A2+D1−B2+C2A3+D2−A1B4,l4=B4+C4A3+D2−B3+C3A4+D3,l5=B2+C2−A1B1+C1,l6=−B4+C4+B1+C1A3−D2−B2+C2A2−D1+A1B3+C3,l7=B2+C2A4−D3−B3+C3A3−D2+B4A2−D1,l8=−B4+C4A4−D3.
Equation (25) is equivalent to the condition that ω=ωτ>0 is the root of Fω,τ=R2ω,τ−T2ω,τ−S2ω,τ=0.
LetI=τ:there  isωτ>0such thatFω,τ=0∪0,τmax, then we denote Snτ=τ−θτ+2nπ/ωτ,τ∈I,n∈N0, where θτ∈0,2π. ωττ=θτ+2nπ. Hence, iω∗ω∗=ωτ∗>0 is the pure imaginary root of (21) and is a necessary and sufficient condition for τ∗ being the root of Snτ=0. We need to introduce the following theorem [18].Theorem 5.
Assume thatSnτ=0 has some positive roots τ∗∈I for some n∈N0. Then, a pair of simple pure imaginary roots ±iωτ∗ exists when τ=τ∗ which crosses the imaginary axis from left to right if δτ∗>0 and crosses the imaginary axis from right to left if δτ∗<0, where(27)δτ∗=signdReλdτ|λ=iωτ∗=signRωτ,τFω′ωτ,τ|τ=τ∗×signSn′τ∗.From the above discussion, the following theorem can be obtained.Theorem 6.
Assume thatR1>1 holds, then(i)
IfI=ϕ or I≠ϕ, Snτ=0 has no positive root τ∗∈I and E∗ is locally asymptotically stable for all τ∈0,τmax(ii)
IfI≠ϕ, Snτ=0 has at least one positive root τ∗∈I and δτ∗≠0, E∗ is locally asymptotically stable when τ∈0,τ∗, and a Hopf bifurcation occurs at E∗ for τ=τ∗Case 4.
τ1∈0,τ1∗,τ2>0,τ1≠τ2.
In this case, we considerτ2 as a parameter and τ1=τ10∗∈0,τmax. Let λ=iωω>0 be a root of (9), then we get(28)ω4−A2ω2+A4+B4−B2ω2cosωτ1+B3ω−B1ω3sinωτ1=C2ω2−C4−D3−D1ω2cosωτ1−D2ωsinωτ1cosωτ2+C1ω3−C3ω+D3−D1ω2sinωτ1−D2ωcosωτ1sinωτ2,−A1ω3+A3ω−B4−B2ω2sinωτ1+B3ω−B1ω3cosωτ1=C4−C2ω2+D3−D1ω2cosωτ1+D2ωsinωτ1sinωτ2+C1ω3−C3ω+D3−D1ω2sinωτ1−D2ωcosωτ1cosωτ2.
From (29), we have(29)cosωτ2=h1ω6+h2ω5+h3ω4+h4ω3+h5ω2+h6ω+h7f1ω6+f2ω5+f3ω4+f4ω3+f5ω2+f6ω+f7,where(30)h1=C2−A1C1+D1−B1C1cosωτ1,h2=A1D1−D2−B1C2+B2C1sinωτ1,h3=A1C3−C4−A2C2+A3C1+B1D2−B2D1+−D3+A1D2−A2D1+B1C3−B2C2+B3C1cosωτ1,h4=A2D2−A1D3−A3D1+B1C4−B2C3+B3C2−B4C4sinωτ1,h5=A2C4−A3C3+A4C2+B2D3−B3D2+B4D1+A2D3−A3D2+A4D1+B2C4−B3C3+B4C2cosωτ1,h6=A3D3−A4D2−B3C4+B4C3sinωτ1,h7=−A4C4−B4D3−A4D3+B4C4cosωτ1,f1=C12,f2=−2C1D1sinωτ1,f3=C22−2C1C3+D12+2C2D1−2C1D2cosωτ1,f4=2C1D3−2C2D2+2C3D1sinωτ1,f5=−2C2C4+C32+D22−2D1D3+2C3D2−2C2D3−2C4D1cosωτ1,f6=2C4D2−2C3D3sinωτ1,f7=C42+D32+2C4D3cosωτ1.
We can also get(31)Fω≜ω8+I1ω6+I2ω4+I3ω2+I4+I5sinωτ1+I6cosωτ1=0,where(32)I1=−2A2+A12−C12+B12,I2=2A4+B22−D12+A22−C22−2A1A3−2B1B3+2C1C3,I3=B32−D22+A32−C32−2A2A4−2B2B4+2C2C4+2D1D3,I4=B42−D32+A42−C42,I5=−2B1ω7+2B3−2A1B2+2A2B1+2C1D1ω5+2A1B4−2A2B3+2A3B2−2A4B1+2C2D2−2C3D1−2C1D3ω3+2A4B3−2A3B4+2C3D3−2C4D2ω,I6=2B2+2A1B1ω6+2B4+2A2B2−2C2D1−2A1B3−2A3B1+2C1D2ω4+2A3B3−2C3D2−2A2B4−2A4B2+2C4D1+2C2D3ω2+2A4B4−2C4D3.
IfB42−D32+A42−C42−2A4B4−2C4D3<0, then F0<0 and limω⟶∞Fω⟶+∞. We obtain that (31) has at least one positive real root. Suppose that there are eight positive real roots of (31) which are ω1∗,ω2∗,…,ω8∗. According to (29), we have(33)τ2k∗j=1ωk∗arccosh1ω6+h2ω5+h3ω4+h4ω3+h5ω2+h6ω+h7f1ω6+f2ω5+f3ω4+f4ω3+f5ω2+f6ω+f7+2jπ,where k=1,…,8,j=0,1,2,…. Define(34)τ20∗=τ2k0∗0=mink∈1,…,8τ2k0∗0,ω20∗=ω2k0∗.
Taking the derivative of (9) with respect to τ2, we can get(35)dλdτ2−1=1λQ1+Q2e−λτ1+Q3e−λτ2+Q4e−λτ1+τ2Q5e−λτ2+Q6e−λτ1+τ2−τ2λ,where(36)Q1=4λ3+3A1λ2+2A2λ+A3,Q2=3B1λ2+2B2λ+B3−τ1B1λ3+B2λ2+B3λ+B4,Q3=3C1λ2+2C2λ+C3,Q4=3C1λ2+2C2λ+C3−τ1D1λ2+D2λ+D3,Q5=C1λ3+C2λ2+C3λ+C4,Q6=D1λ2+D2λ+D3.
From (35), it follows that(37)Redλτ2k∗jdτ2λ=iω2k∗−1=R1R3+R2R4R32+R42,where(38)R1=−3A1ω3+A3+−3B1ω2+B3+B2τ1ω2−B4τ1cosωτ1+2B2ω+B1τ1ω3−B3τ1ωsinωτ1+−3C1ω2+C3cosωτ2+2C2ωsinωτ2+−3C1ω2+C3+D1τ1ω2−D3τ1cosωτ1+τ2+2C2ω−D2τ1ωsinωτ1+τ2,R2=−4ω3+2A2ω−−3B1ω2+B3+B2τ1ω2+B4τ1sinωτ1+2B2ω+B1τ1ω3−B3τ1ωcosωτ1−−3C1ω2+C3sinωτ2+2C2ωcosωτ2−−3C1ω2+C3+D1τ1ω2−D3τ1sinωτ1+τ2+2C2ω−D2τ1ωcosωτ1+τ2,R3=C4ω−C2ω3sinωτ2−C3ω2−C1ω4cosωτ2+D3ω−D1ω3sinωτ1+τ2−D2ω2cosωτ1+τ2,R4=C4ω−C2ω3cosωτ2+C3ω2−C1ω4sinωτ2+D3ω−D1ω3cosωτ1+τ2+D2ω2sinωτ1+τ2.
IfH5R1R3+R2R4≠0 holds, we can obtain the existence of a Hopf bifurcation as follows.Theorem 7.
Forτ1∈0,τ1∗, if H1,H2,H5 holds and B42−D32+A42−C42−2A4B4−2C4D3<0, then E∗ is locally asymptotically stable when τ2∈0,τ20∗. System (2) undergoes Hopf bifurcations at E∗ for τ2=τ2k∗jk=1,…,8,j=0,1,2,….


## 3. Direction and Stability of the Hopf Bifurcation

In this section, we considerτ2 as a parameter and discuss the direction and stability of the Hopf bifurcation of (2) when τ1=τ10∗∈0,τ1∗ by using the normal form and the center manifold theorem. Denote τ2=τ2k∗jk=1,2,…,8,j=0,1,2,… as τ20∗. We can obtain the existence of a Hopf bifurcation at E∗=x∗,y∗,v∗,z∗. Let u1=x−x∗,u2=y−y∗,u3=v−v∗,u4=z−z∗,t=t/τ2,τ2=τ20∗+μ,andμ∈R4. System (2) can be written as an FDE in C=C−1,0,R4 as the following form:(39)u˙t=Lμut+fμ,ut,where ut=xt,yt,vt,ztT∈R4, and Lμ:C⟶R4,f:R×C⟶R4 are shown as follows:(40)Lμϕ=τ20∗+μAϕ10ϕ20ϕ30ϕ40+τ20∗+μBϕ1−τ1τ20∗ϕ2−τ1τ20∗ϕ3−τ1τ20∗ϕ4−τ1τ20∗+τ20∗+μCϕ1−1ϕ2−1ϕ3−1ϕ4−1,fμ,ϕ=τ20∗+μf1f2f3f4,where(41)f1=−rϕ102xM−β1ϕ10ϕ30−β2ϕ10ϕ20+…,f2=β1e−mτ1ϕ1−τ1τ20∗ϕ3−τ1τ20∗−β2e−mτ1ϕ1−τ1τ20∗ϕ2−τ1τ20∗+…,f3=−pϕ30ϕ40+…,f4=cϕ3−1ϕ4−1+…,and ϕ=ϕ1,ϕ2,ϕ3,ϕ4∈R4.By the Riesz representation theorem, there exists a functionηθ,μ of bounded variation for θ∈−1,0, such that(42)Lμϕ=∫−10dηθ,μϕθ,ϕ∈−1,0,R4.Next, we choose(43)ηθ,μ=τ20∗+μAδθ+Cδθ+1+Bδθ+τ1τ20∗,where δ is the Dirac delta function. Define(44)Aμϕ=dϕθdθ,θ∈−1,0,∫−10dηs,μϕs,θ=0,Rμϕ=0,θ∈−1,0,fμ,ϕ,θ=0,where ϕ∈C1−1,0,R4. When θ=0, equation (39) becomes u˙t=Aμut+Rμut, where utθ=ut+θ=xt+θ,yt+θ,vt+θ,zt+θ,θ∈−1,0. Define(45)A∗ψs=−dψsds,s∈0,1,∫−10dηTt,0ψ−t,s=0,(46)ψs,ϕθ=ψ¯0ϕ0−∫−10∫ξ=0θψ¯ξ−θdηθϕξdξ,where ηθ=ηθ,0. A and A∗ are adjoint operators. Let A=A0, we know that ±iω20∗τ20∗ are eigenvalues of both A and A∗. Assume that qθ=1,α,β,γTeiθω20∗τ20∗ is the eigenvector of A corresponding to iω20∗τ20∗. Then, we have(47)τ20∗iω20∗−a11−a12−a130−b21e−iω20∗τ1iω20∗−a22−b22e−iω20∗τ1−b23e−iω20∗τ100−a32iω20∗−a33−a3400−c43e−iω20∗τ20∗iω20∗−a44−c44e−iω20∗τ20∗1αβγ=0000,which yields(48)q0=1,b23e−iω20∗τ1iω20∗−a11+b21a13e−iω20∗τ1a12b23e−iω20∗τ1+a13iω20∗−a22−b22e−iω20∗τ1,iω20∗−a11−a12αa13,c43βe−iω20∗τ20∗iω20∗−a44−c44e−iω20∗τ20∗T.Similarly, it can be verified thatq∗s=D1,α∗,β∗,γ∗eisω20∗τ20∗ is the eigenvector of A∗ corresponding to −iω20∗τ20∗, where(49)α∗=−iω20∗+a11b21eiω20∗τ1,β∗=−a12+iω20∗+a22+b22eiω20∗τ1α∗a32,γ∗=−a34β∗iω20∗+a44+c44eiω20∗τ20∗.According to equation (46), we obtain that(50)q∗s,qθ=D¯1,α∗¯,β∗¯,γ∗¯1,α,β,γT−∫−10∫ξ=0θD¯1,α∗¯,β∗¯,γ∗¯e−iξ−θω20∗τ20∗dηθ1,α,β,γTeiξω20∗τ20∗dξ=D¯1+αα∗¯+ββ∗¯+γγ∗¯−∫−101,α∗¯,β∗¯,γ∗¯θeiθω20∗τ20∗dηθ1,α,β,γT=D¯1+αα∗¯+ββ∗¯+γγ∗¯−τ20∗e−iω20∗τ1e−mτ1α¯∗β1v∗+β2y∗+αα¯∗x∗+βα¯∗x∗−τ20∗e−iω20∗τ20∗βγ¯∗cz∗+γγ¯∗cv∗.Thus, we have(51)D¯=1+αα∗¯+ββ∗¯+γγ∗¯−τ20∗e−iω20∗τ1e−mτ1α¯∗β1v∗+β2y∗+αα¯∗x∗+βα¯∗x∗−τ20∗e−iω20∗τ20∗βγ¯∗cz∗+γγ¯∗cv∗−1,such that q∗s,qθ=1.Then, we compute the coordinates that describe the center manifoldC0 at μ=0. Let ut be the solution of equation (39) at μ=0 and define(52)zt=q∗,ut,Wt,θ=utθ−2Reztqθ.On the center manifoldC0, we have(53)Wt,θ=Wzt,z¯t,θ=W200z22+W11θzz¯+W02θz¯22+W30θz36+⋯,wherez and z¯ are local coordinates for center manifold C0 in the direction of q and q¯.Whenθ=0, we have(54)z˙t=iω20∗τ20∗z+q¯∗0f0z,z¯≜iω20∗τ20∗z+gz,z¯,where(55)gz,z¯=q∗¯0f0z,z¯=g20z22+g11zz¯+g02z¯22+g21z¯z22+⋯.Then,(56)gz,z¯=q∗¯0f0z,z¯=D¯τ20∗1,α¯∗,β¯∗,γ¯∗f10f20f30f40=D¯τ20∗−ru1t20xM−β1u1t0u3t0−β2u1t0u2t0+α¯∗β1e−mτ1u1t−τ1τ20∗u3t−τ1τ20∗+β2e−mτ1u1t−τ1τ20∗u2t−τ1τ20∗−β¯∗pu3t0u4t0+γ¯∗cu4t−1u3t−1.From (52), we have utθ=Wt,θ+zqθ+z¯q¯θ. Comparing the coefficients with (54), we get(57)g20=2D¯τ20∗−rxM−β1β−β2α+α∗¯β1e−mτ1βe−2iω20∗τ1+β2e−mτ1αe−2iω20∗τ1−β∗¯pβγ+γ∗¯cβγe−2iω20∗τ20∗,g11=D¯τ20∗−2rxM−β1β+β¯−β2α+α¯+α∗¯β1e−mτ1β¯+β+β2e−mτ1α¯+α−β∗¯pβ¯γ+βγ¯+γ∗¯cβ¯γ+βγ¯,g02=2D¯τ20∗−rxM−β1β¯−β2α¯+α∗¯β1e−mτ1β¯e2iω20∗τ20∗+β2e−mτ1α¯e2iω20∗τ20∗−β∗¯pβ¯γ¯+γ∗¯cβ¯γ¯e2iω20∗τ20∗,g21=2D¯τ20∗−rxM2W1110+W2010−β1βW1110+12β¯W2010+W130+12W2030−β2αW1110+12α¯W2010+W120+12W2020+α∗¯β1e−mτ1β1e−iω20∗τ20∗W111−τ1τ20∗+12β¯eiω20∗τ20∗W201−τ1τ20∗+e−iω20∗τ20∗W113−τ1τ20∗+12eiω20∗τ20∗W203−τ1τ20∗+β2e−mτ1αe−iω20∗τ20∗W111−τ1τ20∗+12α¯eiω20∗τ20∗W201−τ1τ20∗+e−iω20∗τ20∗W112−τ1τ20∗+12eiω20∗τ20∗W202−τ1τ20∗−β∗¯pβW1140+12β¯W2040+γW1130+12γ¯W2030+γ∗¯cβe−iω20∗τ20∗W114−1+12β¯eiω20∗τ20∗W204−1+γe−iω20∗τ20∗W113−1+12γ¯eiω20∗τ20∗W203−1,where(58)W20θ=ig20ω20τ20∗q0eiω20τ20∗θ+ig¯023ω20τ20∗q¯0e−iω20τ20∗θ+E1e2iω20τ20∗θ,W11θ=−ig11ω20τ20∗q0eiω20τ20∗θ+ig¯11ω20τ20∗q¯0e−iω20τ20∗θ+E2,E1=22iω20∗−a11−a12−a130−b21e−2iω20∗τ12iω20∗−a22−b22e−2iω20∗τ1−b23e−2iω20∗τ100−a322iω20∗−a33−a3400−c43e−2iω20∗τ12iω20∗−a44−c44e−2iω20∗τ1−1M,E2=−a11−a12−a130−b21−a22−b22−b2300−a32−a33−a3400−c43−a44−c44−1N,M=−rxM−β1β−β2αββ1e−mτ1−2iω20∗τ1+αβ2e−mτ1−2iω20∗τ1−pβγcβγe−2iω20∗τ20∗,N=−2rxM−2β1Reβ−2β2Reα2β1e−mτ1Reβ+2β2e−mτ1Reα−2pReβ¯γ2cReβ¯γ.Next, we can calculate the following values:(59)C10=i2ω20∗τ20∗g11g20−2g112−g0223+g212,μ2=−ReC10Reλ′τ0∗,T2=−ImC10+μ2Imλ′τ0∗ω20∗τ20∗,β2=2ReC10.Then, for system (2), when τ1∈0,τ1∗,τ2>0,τ1≠τ2, the direction and the stability of periodic solution of Hopf bifurcation are determined by (59).


## 4. Numerical Simulation

In this section, we consider the following system:(60)dxtdt=rxt1−xt10000−4.8×10−5xtvt−4.7×10−5xtyt,dytdt=4.8×10−5e−0.05τ1xt−τ1vt−τ1+4.7×10−5e−0.05τ1xt−τ1yt−τ1−ayt,dvtdt=kyt−13vt−0.06vtzt,dztdt=0.01vt−τ2zt−τ2−zt,and the numerical simulation is carried out under the following four cases: (i)
τ1=τ2=0.Leta=1, k=200, and r=0.01, then we have R1=1.91>1, and H1 is satisfied. Meanwhile, Δ1=179.53>0, Δ2=2.97×104>0, Δ3=2.66×106>0, Δ4=0.79>0, and H2 holds. Then, E∗ is locally asymptotically stable from Theorem 1 (see Figure 1).Leta=0.9, k=170, and r=60; by simple calculation, we get R1=7.49>1, and H1 is satisfied. And Δ1=250.12>0, Δ2=2.89×106>0, Δ3=3.05×1010>0, and Δ4=4.56×103>0, we know that H2 holds. Then, E∗ is locally asymptotically stable from Theorem 1 (see Figure 2).(ii)
τ1>0andτ2=0.Choosinga=1, k=200, and r=0.01, E∗ is locally asymptotically stable when τ1=30(see Figure 3), and Hopf bifurcation occurs when τ1=72.98(see Figure 4). It supports the results of Theorem 4.(iii)
τ1=τ2=τ>0.Consideringa=0.9, k=170, and r=60, we obtain E∗ is locally asymptotically stable when τ=20(see Figure 5), and Hopf bifurcation occurs when τ=40.1(see Figure 6). The results of Theorem 6 are verified.(iv)
τ1∈0,τmaxandτ2>0Consideringa=0.9, k=170, and r=60, by simple calculation, we get τmax=40.29 and we choose τ2 as a parameter and let τ1=40.2. When τ2=0(see Figure 7) and τ2=5(see Figure 8), E∗ is locally asymptotically stable. And when τ2=48, Hopf bifurcation occurs (see Figure 9).Figure 1
E∗ is asymptotically stable when τ1=τ2=0. The initial value is 9930, 85, 95, and 2760.Figure 2
E∗ is asymptotically stable when τ1=τ2=0. The initial value is 9998, 110, 90, and 2940.Figure 3
E∗ is asymptotically stable when τ1=30 and τ2=0. The initial value is 29980, 45, 95, and 1340.Figure 4
The trajectories of system (6) with τ1=72.98andτ2=0. The initial value is 48930, 6, 95, and 1.Figure 5
E∗ is asymptotically stable when τ1=τ2=20. The initial value is 9998, 20, 90, and 470.Figure 6
The trajectories of system (6) with τ1=τ2=40.1. The initial value is 9998, 7, 90, and 3.Figure 7
The trajectories of system (6) with τ1=40.2andτ2=0. The initial value is 9995, 7, 95, and 1.Figure 8
E∗ is asymptotically stable when τ1=40.2andτ2=5. The initial value is 9995, 7, 95, and 1.Figure 9
The trajectories of system (6) with τ1=40.2andτ2=48. The initial value is 9995, 7, 95, and 1.


## 5. Conclusion

On the basis of [10], considering that e−mτ1 is more practical that stands for the probability of surviving the time period from t−τ1 to t, this paper establishes and investigates an HIV-1 virus model with two delay and delay-dependent parameters to describe both virus-to-cell and cell-to-cell transmissions. We find out basic reproductive rate R1, and system (2) has the unique positive equilibrium E∗ when R1>1. Without time delays, E∗ is asymptotically stable (see Figures 1 and 2). In other words, system (6) does not have any excitable nature. When τ2=0, there exists a threshold limit beyond which the stability exchange takes place and system shows periodic orbits around E∗(Figures 3 and 4). It shows that the number of cells and virus can fluctuate in time. Similarly, we find out that as τ1 and τ2 vary, E∗ loses its stability and Hopf bifurcations occur when τ1=τ2=τ and τ1∈0,τmax,τ2>0(see Figures 5–9). Moreover, we obtain the direction and stability of the Hopf bifurcation when τ1∈0,τ1∗,τ2>0 using the normal form theory and center manifold theorem. Finally, numerical simulations help us illustrate the main results of model (2). Clearly, we can show that delay value is responsible for changing the stable system to system with periodic cycles. In the discussion of the characteristic equations, when neither τ1 nor τ2 is zero, the characteristic equation of system (2) is P0λ,τ1+P1λ,τ1e−λτ1+P2λ,τ1e−λτ2+P3λ,τ1e−λτ1+τ2=0. We only discuss about the case when τ1=τ2=τ>0 by applying the geometric stability switch. When τ1≠τ2, it is very difficult to discuss the root of such characteristic equation. We leave the analysis of the more complicated bifurcations as the future work.



---

*Source: 2521082-2021-11-08.xml*