# Bifurcation Analysis of a Delayed Infection Model with General Incidence Function

Authors: Suxia Zhang; Hongsen Dong; Jinhu Xu
Journal: Computational and Mathematical Methods in Medicine
 (2019)

---

## Abstract

In this paper, an infection model with delay and general incidence function is formulated and analyzed. Theoretical results reveal that positive equilibrium may lose its stability, and Hopf bifurcation occurs when choosing delay as the bifurcation parameter. The direction of Hopf bifurcation and the stability of the periodic solutions are also discussed. Furthermore, to illustrate the numerous changes in the local stability and instability of the positive equilibrium, we conduct numerical simulations by using four different types of functional incidence, i.e., bilinear incidence, saturation incidence, Beddington–DeAngelis response, and Hattaf–Yousfi response. Rich dynamics of the model, such as Hopf bifurcations and chaotic solutions, are presented numerically.


---

## Body

## 1. Introduction

Mathematical modeling has been proven to be valuable in exploring mechanisms and dynamical behaviors of the viral infection process. The analysis of these models can provide insights into developing effective control strategies for infections and evaluating antiviral therapies. Based on the basic virus infection model introduced by Nowak et al. [1], various models have been formulated by many authors to describe the dynamics of virus population qualitatively and quantitatively, and lots of interesting phenomena have been discussed, such as the global stability of the equilibria, bifurcations, periodic oscillations, limit cycles, and the effects of time delays.For models of virus infections, it is observed that the time delays introduced into them can predict the viral infection process better when compared to models without delays, such as the time between viral entry into a target cell and the production of subsequent viral particles, the time necessary for the newly produced virus to become mature and then infectious particles, and that needed to activate the immune response, all of which cannot be ignored when describing the interaction of them [2, 3]. In the field of virus dynamics, many models with discrete or distributed delays have been studied, showing that rich dynamics including stability switches, Hopf bifurcation, and chaotic oscillations can be generated by the mechanisms of time delays. Specifically, stability switches can occur due to a delayed immune response and exponentially decayed delay-dependent parameters [4–9], which leads to the switches in the stability of equilibria because of the change in the value of parameters related with these delays.It is known that the function forms of the incidence rate of the infection have a crucial role in the modeling of the virus dynamics, which are important in determining qualitative behavior of the proposed models and in giving a reasonable description of the dynamics. Huang et al. [10] and Zhang and Xu [11] formulated a virus dynamical model with Beddington–DeAngelis (BD) infection rate βxυ/1+a1x+a2υ, where x represents the concentration of uninfected cells and υ represents that of virus, with a1,a2≥0 being constants. This function is similar to the well-known Holling type II functional response, and the term a2υ in the denominator reflects the mutual interference between viruses. In [12], Zhou and Cui used a Crowley–Martin (CM) function response of the form βxυ/1+a1x+a2υ+a1a2xυ, with a1,a2,a3≥0, which are constants. It can be seen that in BD and CM functional responses, when a1>0,a2=0, these two functions are simplified to the Holling type II functional response. And when a1=0,a2>0, they express a saturation response. Moreover, when a1=0,a2=0, they are the mass action process (or Holling type I functional response). Other more general incidence rates have also been proposed. McCluskey and Yang used function fx,υ as incidence rate under some biologically motivated assumptions to study the global stability of a virus model [13]. Hattaf et al. in [14, 15] used an incidence rate fx,y,υ that covers a variety of incidence functions such as BD response when a3=0 and CM response when a3=a1a2 if fx,y,υ=βx/1+a1x+a2υ+a3xυ, which will be called Hattaf–Yousfi response in the following, with y being infected cells.In this paper, we will study the influence of general incidence function and time delay on the dynamical behaviors of an infection model. Motivated by the works of [7, 13], we will focus on the dynamics of a delayed model that incorporates the immune response in antiviral defense, time delay in activating the immune response for the virus, and two general incidence functions for the transmission of virus-to-cell and cell-to-cell, respectively. The main aim of this paper is to study stability switches of the positive equilibrium, which are generated by the mechanism of time delay, implying the existence of changes in local stability and the instability of this equilibrium. By using the characteristic equation with delay-dependent parameters, normal theory, and center manifold theorem, we can indicate the existence of pure imaginary roots and then Hopf bifurcation. Furthermore, to verify the theoretical results, numerical simulations are performed for four different forms of incidence functions, including bilinear incidence, saturation incidence, Beddington–DeAngelis response, and Hattaf–Yousfi functional response, respectively.The organization of the paper is as follows. In the next section, we introduce the virus dynamic model. The boundedness and nonnegativity of solutions of the model as well as the existence and uniqueness of equilibria are discussed in Section3. In Section 4, both the stabilities of two equilibria and the conditions for the existence of Hopf bifurcation are studied. Furthermore, numerical simulations are presented in Section 5. A brief summary is given in Section 6, and the properties of the Hopf bifurcation solutions have been provided in Appendix.


## 2. Model Formulation

To account for the possible effect of the latent period in viral infection, we introduce a constant time delayτ to represent the time needed to activate the immune response. In this model, two transmissions of virus-to-cell and cell-to-cell have been both taken into consideration, and we hope to employ the function of incidence rate given by more general form. Consequently, the recruitment of infected cells is given by two functions, i.e., f1u,vv and f2u,ww, which were used in [16–18], here u, w, and v denote the concentrations of uninfected cells, infected cells, and free virus particles. z is used to denote the concentration of CTL cells. Assume that the uninfected cells are produced at a constant rate s. d1, d2, d3, and d4 represent the death rates of uninfected cells, infected cells, free virus, and CTL cells. p denotes the killing rate of infected cells by CTL cells. k and c are the production rates of free viral particles and the effector cells, respectively. Then, we construct the following model:(1)dudt=s−f1u,vvt−f2u,wwt−d1ut,dwdt=f1u,vvt+f2u,wwt−d2wt−pwtzt,dvdt=kwt−d3vt,dzdt=cwt−τ−d4zt,with initial conditions(2)uθ=φ1θ,wθ=φ2θ,vθ=φ3θ,zθ=φ4θ,θ∈−τ,0,where φ=φ1,φ2,φ3,φ4∈C−τ,0,ℝ+4 is the Banach space of continuous functions mapping the interval −τ,0 into ℝ+4, with φiθ≥0 (θ∈−τ,0,i=1,2,3,4), and all constant parameters are positive.Functionfi:ℝ+2⟶ℝ+ (i=1,2) is continuously differentiable and is assumed to satisfy the following conditions:(3)f10,v=f20,w=0,for allv≥0andw≥0,(4)∂f1∂u>0,∂f2∂u>0,and∂f1∂v≤0,∂f2∂w≤0,for allu≥0,v≥0andw≥0,(5)∂f1v∂v=f1+∂f1∂v,v≥0,∂f2w∂w=f2+∂f2∂w,w≥0.These fundamental assumptions are biologically motivated, and it is easy to check that class of functionsf1 (and f2) satisfying these hypotheses include incidence functions such as f1=βu (mass incidence), βu/1+av (saturation incidence), βu/1+a1u+a2v (Beddington–DeAngelis response), and βu/1+a1u+a2v+a3uv (Hattaf–Yousfi response).


## 3. Preliminaries and Equilibria

We denote the state space of (1) by X=C−τ,0,ℝ+4, equipped with the sup-norm φ=sup−τ≤θ≤0φθ. For any φ∈X, the existence and uniqueness of the solution(6)Yt,φ=ut,φ,wt,φ,vt,φ,zt,φ,of model (1) with initial condition (2) follow from the standard theory of functional differential equations [19]. Furthermore, X is positively invariant for (1). In the following, we discuss the boundedness [20] and nonnegativity of solutions of model (1) as well as the existence of equilibria.Theorem 1.
The solutionsYt,φ of model (1) with initial condition (2) are nonnegative and ultimately uniformly bounded for all t≥0.Proof.
Forut, if there exists t0 such that ut>0 for t∈0,t0 and ut0=0, then it is obvious that u′t0=s>0, which implies that for sufficiently small ε>0, ut<0 for t∈t0−ε,t0, contradicting with ut>0 for all t∈0,t0. Thus, it is valid that ut>0 for all t≥0. In order to show wt>0, vt>0, and zt>0 for all t≥0, assume that there exists t1>0 such that minwt1,vt1,zt1 = 0 for the first time. Firstly, if wt1=0, then from the second equation of (1); we have(7)dwt1dt=f1ut1,vt1,vt1>0,and then there exists ε1 small enough such that dwt/dt>0 for t∈t1−ε1,t1, which contradicts with the facts that wt>0 for t∈t1−ε1,t1 and wt1=0. Thus wt>0 for all t≥0. For vt and zt, since(8)vt=v0e−d3t+∫0tkwθe−d3t−θdθ,zt=z0e−d4t+∫0tcwθ−τe−d4t−θdθ,it is obvious that vt1=0 or zt1=0 contradicts with the above expression about vt and zt at t=t1. Therefore, the solution Yt,φ with initial condition (2) is nonnegative for all t>0.
To prove the boundedness of the solutions, lettingN=u+w, we have(9)dNdt=s−d1u−d2w−pwz≤s−dN,where d=mind1,d2, then limsupt⟶∞N≤s/d, implying that ut and wt are ultimately bounded. For vt, it follows from the third equation of (1) that(10)dvdt≤ksd−d3v,which leads to limsupt⟶∞v≤ks/dd3. Similarly, we have limsupt⟶∞z≤cs/dd4. Therefore, the solutions of system (1) are ultimately uniformly bounded.
System (1) always has an infection-free equilibrium E0=u0,0,0,0, where u0=s/d1. In the following, we analyze the existence and uniqueness of the positive equilibrium E∗=u∗,w∗,v∗,z∗.
For convenience, we denote(11)f10=f1u0,0,f20=f2u0,0,f1∗=f1u∗,v∗,f2∗=f2u∗,w∗.
Define the basic reproductive number as(12)ℛ0=kf10d2d3+f20d2≔ℛ01+ℛ02.
From biological viewpoint,ℛ0 can be divided into two parts, with ℛ01 measuring the average number of secondary infected generation caused by an existing free virus and ℛ02 being that caused by an infected cell, which give the basic reproduction number corresponding to virus-to-cell and cell-to-cell infections, respectively.
It is easy to obtainv∗=k/d3w∗ and z∗=c/d4w∗. By the first and second equations in (1), we have(13)s−d1u∗−d2w∗−pw∗z∗=0,which gives u∗=1/d1s−w∗d2+pc/d4w∗. In order to get the infected steady state, it is necessary that s−w∗d2+pc/d4w∗≥0, and then, we must have 0<w≤w¯, with w¯=−d2d4+d4Δ/2pc and Δ=d22+4spc/d4.
From the equation ofw, it follows that(14)f1∗v∗+f2∗w∗−d2+pcd4w∗w∗=0.
Denote the left side of (11) as Gw∗, and then the positive equilibrium of model (1) are given by Gw∗=0 for w∗∈0,w¯. Noting that G0=0 and(15)Gw¯=−d2+pcd4w¯,w¯<0,for u=0 when w=w¯, and it is sufficient to show that G′0>0 if there is a positive equilibrium. In fact, when ℛ0>1, we have(16)G′0=f1sd1,0kd3+f2sd1,0−d2=d2ℛ0−1>0,and then Gw∗ is positive for sufficiently small w∗; therefore, there exists an infected steady state E∗.
Now we examine the derivative ofGw∗ at E∗. Noting that f1∗v∗+f2∗w∗=d2w∗+pw∗z∗, and then(17)f1∗v∗+f2∗w∗w∗=kf1∗v∗+f2∗w∗d3v∗=d2+pcd4w∗,and thus, we have(18)G′w∗=f1∗kd3+∂f1∗∂u−d2−2pcd4w∗kd3w∗+∂f1∗∂vkd3kd3w∗+f2∗+∂f2∗∂u−d2−2pcd4w∗w∗+∂f2∗∂ww∗−d2+2pcd4w∗=f1∗kd3+f2∗−d2+pcd4w∗+T−=f1∗kd3+f2∗−f1∗kd3−f2∗kw∗d3v∗+T−=T−,where(19)T−=∂f1∗∂u−d2−2pcd4w∗kd3w∗+∂f1∗∂vkd3kd3w∗+∂f2∗∂u−d2−2pcd4w∗w∗+∂f2∗∂ww∗−pcd4w∗,denoting the sum of negative terms, and thus function G is strictly decreasing at each of its zeros. Suppose there exists more than one infected steady state, then there must exist an equilibrium E¯∗=u¯∗,w¯∗,v¯∗,z¯∗ such that Gw¯∗≥0, which contradicts with the above fact. Therefore, there is only one positive equilibrium E∗ when the basic reproductive number ℛ0>1, and it is easy to verify that there exists no infection equilibrium when ℛ0<1. Then the following result is obtained.Theorem 2.
For system (1), when ℛ0<1, there exists only the infection-free equilibrium E0=s/d1,0,0,0. If ℛ0>1, there is a unique infection equilibrium E∗=u∗,w∗,v∗,z∗ as well as E0.
Whenℛ0>1, by using the persistence theory in [21], it can be shown that system (1) is uniformly persistent.


## 4. Stability and Hopf Bifurcation

In this section, we will study the stability of the two equilibria and the conditions for the existence of Hopf bifurcation.


### 4.1. Stability ofE0

Firstly, we linearize model (1) at the steady states to discuss its local stability. The Jacobian matrix leads to the following characteristic equation:(20)λ+d1+∂f1∂uv+∂f2∂uwf2+∂f2∂wwf1+∂f1∂vv0−∂f1∂uv+∂f2∂uwλ+d2+pz−f2−∂f2∂ww−f1+∂f1∂vvpw0−kλ+d300−ce−λτ0λ+d4=0,and thus, the linearization of system (1) at E0 can be expressed by(21)λ+d1λ+d4λ2+d2+d3−f20λ+d2d31−ℛ0=0.It is clear that (21) has two negative real roots λ1=−d1 and λ2=−d4. The stability of E0 is determined completely by using the following equation:(22)λ2+d2+d3−f20λ+d2d31−ℛ0=0.For the caseℛ0<1, the roots of (22) have only negative real parts. Then the infection-free equilibrium E0 is locally asymptotically stable. Otherwise, when ℛ0>1, equation (22) has at least one root with positive real part, which implies that E0 is unstable.Theorem 3.
Ifℛ0<1, then the infection-free equilibrium E0 of model (1) is globally asymptotically stable, and it is unstable when ℛ0>1.Proof.
For a continuous and bounded functionft, we define(23)f∞=limsupt⟶∞ft,f∞=liminft⟶∞ft.
Then for any solution of (1), we have(24)0≤u∞≤u∞<∞,0≤w∞≤w∞<∞,0≤v∞≤v∞<∞,0≤z∞≤z∞<∞.
By using the fluctuation lemma in [22], we know that there is a sequence tn with tn⟶∞ such that(25)wtn⟶w∞,andw′tn⟶0,n⟶∞.
Substituting the sequencetn into the first equation of (1) and taking the limit gives(26)limn⟶∞u′tn≤s−d1limn⟶∞utn,and then we have u∞≤λ/d1=u0. A similar argument to the second and third equations of model (1) yields(27)d2w∞≤f1u∞,v∞v∞+f2u∞,w∞w∞,and(28)v∞=kd3w∞.
By the assumption of (4), combining equality (28) into (27) gives(29)d2w∞≤f1u0,0kd3w∞+f2u0,0w∞,which leads to d2w∞≤d2ℛ0w∞. Noticing that w∞ is nonnegative since it is the supremum of the function wt, then there are possible cases of w∞>0 or w∞=0. If w∞>0, then we have ℛ0≥1, which gives a contradiction. Therefore, w∞=0 is valid. From (28), we have v∞=0. When considering the fourth equation about zt in model (1), we obtain z∞=0. As for ut, we can find that limt⟶∞ut=s/d1 by applying the limit theory to the first equation in system (1), and thus, the proof is completed.


### 4.2. Stability ofE∗

For convenience, we introduce the following notations:(30)∂f1∗∂uv∗+∂f2∗∂uw∗=Q0,∂f1∗∂vv∗+f1∗=Q1,∂f2∗∂ww∗+f2∗=Q2.From the assumption (5), it follows that Qi≥0i=0,1,2.In the following calculation, we will use the following equilibrium equations:(31)s−f1∗v∗−f2∗w∗−d1u∗=0,f1∗v∗+f2∗w∗=d2w∗+pw∗z∗,kw∗=d3v∗,cw∗=d4z∗.Then, we have(32)d2+pz∗−Q2=d2+pz∗−f2∗−∂f2∗∂ww∗=f1∗v∗+f2∗w∗w∗−f2∗−∂f2∗∂ww∗=f1∗kd3−∂f2∗∂w,w∗>0,−kf1∗=−kd2w∗+pw∗z∗−f2∗w∗v∗=−kd2d3k+pz∗d3k−f2∗d3k=−d2d3−d3pz∗+d3f2∗,d3f2∗−d3Q2=−d3∂f2∗∂w,w∗>0.Linearized model (1) at the infection equilibrium E∗ and the characteristic equation can be expressed as(33)λ4+A3λ3+A2λ2+A1λ+A0+B2λ2+B1λ+B0e−λτ=0,where Ai=Aiτi=0,1,2,3 and Bj=Bjτj=0,1,2 with(34)A3=d1+d3+d4+Q0+d2+pz∗−Q2,A2=d1+Q0d3+d4+d2+pz∗−Q2+d3d4+d4d2+pz∗−Q2+Q0Q2−kv∗∂f1∗∂v−d3w∗∂f2∗∂w,A1=d3d4d1+Q0+Q0Q2d3+d4+d4d1+Q0d2+pz∗−Q2+d1+d4+Q0−kv∗∂f1∗∂v−d3w∗∂f2∗∂w,A0=d4d1+Q0−kv∗∂f1∗∂v−d3w∗∂f2∗∂w+d3d4Q0Q2+kQ1,andB2=pw∗c,B1=pw∗cd1+d3+Q0,B0=pw∗cd3d1+Q0.It can be seen thatAii=0,1,2,3 and Bjj=0,1,2 are all positive. When the delay τ=0, Aii=0,1,2,3 and Bjj=0,1,2 are independent with τ and the characteristic equation (33) can be reduced to(35)λ4+A3λ3+A2+B2λ2+A1+B1λ+A0+B0=0.By using the Routh–Hurwitz criterion, we know that all solutions of (35) have negative real parts if and only if the following conditions are satisfied:(36)H1=A3A2+B2−A1+B1>0,H2=A3A2+B2A1+B1−A32A0+B0−A1+B12>0.Letf=f1v+f2w,f∗=f1∗v∗+f2∗w∗,fu=fu,w∗,v∗ and assume(37)1−ffuffu−vv∗≥0,for allu,w,v>0.Theorem 4.
If (37) holds and τ=0, then the infection equilibrium E∗ is globally asymptotically stable when ℛ0>1.Proof.
Construct a Lyapunov functional as(38)Lt=u−u∗−∫u∗uf∗fsds+w−w∗lnww∗+d2+pz∗kv−v∗lnvv∗+pz∗cz−z∗lnzz∗,and then calculating the time derivation of Lt along system (1) leads to(39)dLtdt=d1u∗1−f∗fu1−uu∗+f∗ffu−vv∗+pz∗−wz∗z−wzz∗+2pwz∗+3f∗−f∗wv∗w∗v−f∗f∗fu−f∗fw∗f∗w=d1u∗1−f∗fu1−uu∗−pwz∗z∗z+zz∗−2+f∗gffu−gvv∗−f∗gwv∗w∗v+gf∗fu+gw∗fwf∗,with gx=x−1−lnx. By using (4), (37), and the property of function gx, we have 1−f∗/fu1−u/u∗≤0 and gf/fu≤gv/v∗, which implies that dLt/dt≤0 and it can verified that E∗ is the largest compact invariant set where dLt/dt=0. Therefore, E∗ is globally asymptotically stable by using LaSalle’s invariance principle.
However, if the delayτ≠0, it may destabilize the infected steady state and lead to Hopf bifurcation, which will be discussed in the following subsection.


### 4.3. Hopf Bifurcation

It is shown that all roots of (33) locate in the left side of the imaginary axis if τ=0 and (36) is satisfied, which leads to local stability of the infected equilibrium E∗. When τ increases from 0 to variant positive delays, it is possible that the roots of (35) pass through the imaginary axis and enter the right side in the complex plane. Therefore, the stability switch may occur, and it is necessary to study the transcendental equation (33) when λ=iω, which is a critical value under small perturbation. In the following analysis, we will study the occurrence of any possible stability switching in this case.Substitutingλ=iω into (33) and separating the real and imaginary parts give(40)−ω4+A2ω2−A0=B0−B2ω2cosωτ+B1ωsinωτ,−A3ω3+A1ω=B0−B2ω2sinωτ−B1ωcosωτ.Squaring and adding the two equalities lead to(41)Fω=ω8+M3ω6+M2ω4+M1ω2+M0=0,where(42)M3=A32−2A2,M2=A22+2A0−2A1A3−B22,M1=A12−2A2A0+2B2B0−B12,M0=A02−B02.Note that equationFω=0 has at least one positive root when M0<0, i.e., A0<B0. And let z=ω2, then the characteristic equation (33) has a purely imaginary root iω which is equivalent to Fz=0 has a positive real root z. Take the transformation y=z+M3/4; then for F′z=0, it is equivalent to y3+m1y+m0=0 with m1=M2/2−3M32/16, m0=M33/32−M3M2/8+M1/4, and it has roots of yi, correspondingly, zi=yi−M3/4 (i=1,2,3); thus we have the following result about the roots of Fz=0 [23].Lemma 1.
DenotingΔ=m0/22+m1/33, for the roots of Fz=0,(i)
IfM0<0, then there exists at least one positive root(ii)
IfM0≥0 and Δ≥0, then there are positive roots if and only if z1>0 and Fz1≤0(iii)
IfM0≥0 and Δ<0, then there are positive roots if and only if there exists at least a positive z∗∈z1,z2,z3 such that Fz∗≤0
Now assume thatFz=0 has four roots zi∗i=1,2,3,4 with zi∗>0; then by the relation of z=ω2, we have the fact of ωi=zi∗. From the equalities in (40), we can express cosωτ and sinωτ as(43)cosωτ=ω4−A2ω2+A0B2ω2−B0+B1ωA3ω3−A1ωB12ω2+B0−B2ω22≜Fc,sinωτ=B1ω−ω4+A2ω2−A0+B0−B2ω2−A3ω3+A1ωB12ω2+B0−B2ω22≜Fs.
Then forj=0,1,2,3,…, we can obtain the following expression of delay τ:(44)τji=arccosFc+2πjωi,Fs≥0,2π−arccosFc+2πjωi,Fs<0.Lemma 2.
For the roots of characteristic equation (35), suppose conditions in (36) hold, then there are two possibilities:(i)
Equation (35) has only roots with negative real parts for τ∈0,mini,jτji if any one of the following conditions (c1)∼(c3) holds(ii)
Equation (35) has only roots with negative real parts for τ≥0 if the conditions (c1)∼(c3) are not satisfied, with(c1) M0<0(c2) M0≥0, Δ≥0, z1>0 and Fz1<0(c3) M0≥0, Δ<0, there is a positive z∗∈z1,z2,z3 such that Fz∗≤0
From this lemma, we can see that the infection equilibrium is asymptotically stable for allτ≥0 if the conditions (c1)∼(c3) are not satisfied. Otherwise, if one of these three conditions is satisfied, then the infection equilibrium is asymptotically stable for τ∈0,mini,jτji, and a Hopf bifurcation can occur at this equilibrium when τ=τ∗ with τ∗ being the critical value of τji.Theorem 5.
Suppose (36) holds, when τ∗=τji(correspondingly ω∗=ωi for some i=1,2,3,4), then the characteristic equation (35) admits a pair of simple conjugate pure imaginary roots λ=iω∗ and λ=−iω∗, which crosses the imaginary axis from left to right (from right to left) if δ<0>0, where(45)δ=signdReλdττ=τ∗=signF′ω∗2.Proof.
Differentiating the characteristic equation (33) with respect to delay τ and arranging it can give the expression of dλ/dτ−1 as(46)dλdτ−1=4λ3+3A3λ2+2A2λ+A1−λλ4+A3λ3+A2λ2+A1λ+A0+2B2λ+B1λB2λ2+B1λ+B0−τλ.
Meanwhile, noticing the fact(47)signdReλdττ=τ∗=signRedλdτ−1τ=τ∗,and then some computation leads to(48)signdReλdττ=τ∗=signRe4λ3+3A3λ2+2A2λ+A1−λλ4+A3λ3+A2λ2+A1λ+A0τ=τ∗+Re2B2λ+B1λB2λ2+B1λ+B0τ=τ∗=sign4ω∗6+3M3ω∗4+2M2ω∗2+M1B2ω∗2−B02+B12ω∗2=signF′ω∗2B2ω∗2−B02+B12ω∗2=signF′ω∗2,which gives the desired result.
When model (1) undergoes stability switch at E∗ for τ=τ∗, the conditions and theorem for the direction and stability of Hopf bifurcation can be discussed and proven by the normal theory and center manifold theorem in [24], which are presented in Appendix.


## 4.1. Stability ofE0

Firstly, we linearize model (1) at the steady states to discuss its local stability. The Jacobian matrix leads to the following characteristic equation:(20)λ+d1+∂f1∂uv+∂f2∂uwf2+∂f2∂wwf1+∂f1∂vv0−∂f1∂uv+∂f2∂uwλ+d2+pz−f2−∂f2∂ww−f1+∂f1∂vvpw0−kλ+d300−ce−λτ0λ+d4=0,and thus, the linearization of system (1) at E0 can be expressed by(21)λ+d1λ+d4λ2+d2+d3−f20λ+d2d31−ℛ0=0.It is clear that (21) has two negative real roots λ1=−d1 and λ2=−d4. The stability of E0 is determined completely by using the following equation:(22)λ2+d2+d3−f20λ+d2d31−ℛ0=0.For the caseℛ0<1, the roots of (22) have only negative real parts. Then the infection-free equilibrium E0 is locally asymptotically stable. Otherwise, when ℛ0>1, equation (22) has at least one root with positive real part, which implies that E0 is unstable.Theorem 3.
Ifℛ0<1, then the infection-free equilibrium E0 of model (1) is globally asymptotically stable, and it is unstable when ℛ0>1.Proof.
For a continuous and bounded functionft, we define(23)f∞=limsupt⟶∞ft,f∞=liminft⟶∞ft.
Then for any solution of (1), we have(24)0≤u∞≤u∞<∞,0≤w∞≤w∞<∞,0≤v∞≤v∞<∞,0≤z∞≤z∞<∞.
By using the fluctuation lemma in [22], we know that there is a sequence tn with tn⟶∞ such that(25)wtn⟶w∞,andw′tn⟶0,n⟶∞.
Substituting the sequencetn into the first equation of (1) and taking the limit gives(26)limn⟶∞u′tn≤s−d1limn⟶∞utn,and then we have u∞≤λ/d1=u0. A similar argument to the second and third equations of model (1) yields(27)d2w∞≤f1u∞,v∞v∞+f2u∞,w∞w∞,and(28)v∞=kd3w∞.
By the assumption of (4), combining equality (28) into (27) gives(29)d2w∞≤f1u0,0kd3w∞+f2u0,0w∞,which leads to d2w∞≤d2ℛ0w∞. Noticing that w∞ is nonnegative since it is the supremum of the function wt, then there are possible cases of w∞>0 or w∞=0. If w∞>0, then we have ℛ0≥1, which gives a contradiction. Therefore, w∞=0 is valid. From (28), we have v∞=0. When considering the fourth equation about zt in model (1), we obtain z∞=0. As for ut, we can find that limt⟶∞ut=s/d1 by applying the limit theory to the first equation in system (1), and thus, the proof is completed.


## 4.2. Stability ofE∗

For convenience, we introduce the following notations:(30)∂f1∗∂uv∗+∂f2∗∂uw∗=Q0,∂f1∗∂vv∗+f1∗=Q1,∂f2∗∂ww∗+f2∗=Q2.From the assumption (5), it follows that Qi≥0i=0,1,2.In the following calculation, we will use the following equilibrium equations:(31)s−f1∗v∗−f2∗w∗−d1u∗=0,f1∗v∗+f2∗w∗=d2w∗+pw∗z∗,kw∗=d3v∗,cw∗=d4z∗.Then, we have(32)d2+pz∗−Q2=d2+pz∗−f2∗−∂f2∗∂ww∗=f1∗v∗+f2∗w∗w∗−f2∗−∂f2∗∂ww∗=f1∗kd3−∂f2∗∂w,w∗>0,−kf1∗=−kd2w∗+pw∗z∗−f2∗w∗v∗=−kd2d3k+pz∗d3k−f2∗d3k=−d2d3−d3pz∗+d3f2∗,d3f2∗−d3Q2=−d3∂f2∗∂w,w∗>0.Linearized model (1) at the infection equilibrium E∗ and the characteristic equation can be expressed as(33)λ4+A3λ3+A2λ2+A1λ+A0+B2λ2+B1λ+B0e−λτ=0,where Ai=Aiτi=0,1,2,3 and Bj=Bjτj=0,1,2 with(34)A3=d1+d3+d4+Q0+d2+pz∗−Q2,A2=d1+Q0d3+d4+d2+pz∗−Q2+d3d4+d4d2+pz∗−Q2+Q0Q2−kv∗∂f1∗∂v−d3w∗∂f2∗∂w,A1=d3d4d1+Q0+Q0Q2d3+d4+d4d1+Q0d2+pz∗−Q2+d1+d4+Q0−kv∗∂f1∗∂v−d3w∗∂f2∗∂w,A0=d4d1+Q0−kv∗∂f1∗∂v−d3w∗∂f2∗∂w+d3d4Q0Q2+kQ1,andB2=pw∗c,B1=pw∗cd1+d3+Q0,B0=pw∗cd3d1+Q0.It can be seen thatAii=0,1,2,3 and Bjj=0,1,2 are all positive. When the delay τ=0, Aii=0,1,2,3 and Bjj=0,1,2 are independent with τ and the characteristic equation (33) can be reduced to(35)λ4+A3λ3+A2+B2λ2+A1+B1λ+A0+B0=0.By using the Routh–Hurwitz criterion, we know that all solutions of (35) have negative real parts if and only if the following conditions are satisfied:(36)H1=A3A2+B2−A1+B1>0,H2=A3A2+B2A1+B1−A32A0+B0−A1+B12>0.Letf=f1v+f2w,f∗=f1∗v∗+f2∗w∗,fu=fu,w∗,v∗ and assume(37)1−ffuffu−vv∗≥0,for allu,w,v>0.Theorem 4.
If (37) holds and τ=0, then the infection equilibrium E∗ is globally asymptotically stable when ℛ0>1.Proof.
Construct a Lyapunov functional as(38)Lt=u−u∗−∫u∗uf∗fsds+w−w∗lnww∗+d2+pz∗kv−v∗lnvv∗+pz∗cz−z∗lnzz∗,and then calculating the time derivation of Lt along system (1) leads to(39)dLtdt=d1u∗1−f∗fu1−uu∗+f∗ffu−vv∗+pz∗−wz∗z−wzz∗+2pwz∗+3f∗−f∗wv∗w∗v−f∗f∗fu−f∗fw∗f∗w=d1u∗1−f∗fu1−uu∗−pwz∗z∗z+zz∗−2+f∗gffu−gvv∗−f∗gwv∗w∗v+gf∗fu+gw∗fwf∗,with gx=x−1−lnx. By using (4), (37), and the property of function gx, we have 1−f∗/fu1−u/u∗≤0 and gf/fu≤gv/v∗, which implies that dLt/dt≤0 and it can verified that E∗ is the largest compact invariant set where dLt/dt=0. Therefore, E∗ is globally asymptotically stable by using LaSalle’s invariance principle.
However, if the delayτ≠0, it may destabilize the infected steady state and lead to Hopf bifurcation, which will be discussed in the following subsection.


## 4.3. Hopf Bifurcation

It is shown that all roots of (33) locate in the left side of the imaginary axis if τ=0 and (36) is satisfied, which leads to local stability of the infected equilibrium E∗. When τ increases from 0 to variant positive delays, it is possible that the roots of (35) pass through the imaginary axis and enter the right side in the complex plane. Therefore, the stability switch may occur, and it is necessary to study the transcendental equation (33) when λ=iω, which is a critical value under small perturbation. In the following analysis, we will study the occurrence of any possible stability switching in this case.Substitutingλ=iω into (33) and separating the real and imaginary parts give(40)−ω4+A2ω2−A0=B0−B2ω2cosωτ+B1ωsinωτ,−A3ω3+A1ω=B0−B2ω2sinωτ−B1ωcosωτ.Squaring and adding the two equalities lead to(41)Fω=ω8+M3ω6+M2ω4+M1ω2+M0=0,where(42)M3=A32−2A2,M2=A22+2A0−2A1A3−B22,M1=A12−2A2A0+2B2B0−B12,M0=A02−B02.Note that equationFω=0 has at least one positive root when M0<0, i.e., A0<B0. And let z=ω2, then the characteristic equation (33) has a purely imaginary root iω which is equivalent to Fz=0 has a positive real root z. Take the transformation y=z+M3/4; then for F′z=0, it is equivalent to y3+m1y+m0=0 with m1=M2/2−3M32/16, m0=M33/32−M3M2/8+M1/4, and it has roots of yi, correspondingly, zi=yi−M3/4 (i=1,2,3); thus we have the following result about the roots of Fz=0 [23].Lemma 1.
DenotingΔ=m0/22+m1/33, for the roots of Fz=0,(i)
IfM0<0, then there exists at least one positive root(ii)
IfM0≥0 and Δ≥0, then there are positive roots if and only if z1>0 and Fz1≤0(iii)
IfM0≥0 and Δ<0, then there are positive roots if and only if there exists at least a positive z∗∈z1,z2,z3 such that Fz∗≤0
Now assume thatFz=0 has four roots zi∗i=1,2,3,4 with zi∗>0; then by the relation of z=ω2, we have the fact of ωi=zi∗. From the equalities in (40), we can express cosωτ and sinωτ as(43)cosωτ=ω4−A2ω2+A0B2ω2−B0+B1ωA3ω3−A1ωB12ω2+B0−B2ω22≜Fc,sinωτ=B1ω−ω4+A2ω2−A0+B0−B2ω2−A3ω3+A1ωB12ω2+B0−B2ω22≜Fs.
Then forj=0,1,2,3,…, we can obtain the following expression of delay τ:(44)τji=arccosFc+2πjωi,Fs≥0,2π−arccosFc+2πjωi,Fs<0.Lemma 2.
For the roots of characteristic equation (35), suppose conditions in (36) hold, then there are two possibilities:(i)
Equation (35) has only roots with negative real parts for τ∈0,mini,jτji if any one of the following conditions (c1)∼(c3) holds(ii)
Equation (35) has only roots with negative real parts for τ≥0 if the conditions (c1)∼(c3) are not satisfied, with(c1) M0<0(c2) M0≥0, Δ≥0, z1>0 and Fz1<0(c3) M0≥0, Δ<0, there is a positive z∗∈z1,z2,z3 such that Fz∗≤0
From this lemma, we can see that the infection equilibrium is asymptotically stable for allτ≥0 if the conditions (c1)∼(c3) are not satisfied. Otherwise, if one of these three conditions is satisfied, then the infection equilibrium is asymptotically stable for τ∈0,mini,jτji, and a Hopf bifurcation can occur at this equilibrium when τ=τ∗ with τ∗ being the critical value of τji.Theorem 5.
Suppose (36) holds, when τ∗=τji(correspondingly ω∗=ωi for some i=1,2,3,4), then the characteristic equation (35) admits a pair of simple conjugate pure imaginary roots λ=iω∗ and λ=−iω∗, which crosses the imaginary axis from left to right (from right to left) if δ<0>0, where(45)δ=signdReλdττ=τ∗=signF′ω∗2.Proof.
Differentiating the characteristic equation (33) with respect to delay τ and arranging it can give the expression of dλ/dτ−1 as(46)dλdτ−1=4λ3+3A3λ2+2A2λ+A1−λλ4+A3λ3+A2λ2+A1λ+A0+2B2λ+B1λB2λ2+B1λ+B0−τλ.
Meanwhile, noticing the fact(47)signdReλdττ=τ∗=signRedλdτ−1τ=τ∗,and then some computation leads to(48)signdReλdττ=τ∗=signRe4λ3+3A3λ2+2A2λ+A1−λλ4+A3λ3+A2λ2+A1λ+A0τ=τ∗+Re2B2λ+B1λB2λ2+B1λ+B0τ=τ∗=sign4ω∗6+3M3ω∗4+2M2ω∗2+M1B2ω∗2−B02+B12ω∗2=signF′ω∗2B2ω∗2−B02+B12ω∗2=signF′ω∗2,which gives the desired result.
When model (1) undergoes stability switch at E∗ for τ=τ∗, the conditions and theorem for the direction and stability of Hopf bifurcation can be discussed and proven by the normal theory and center manifold theorem in [24], which are presented in Appendix.


## 5. Numerical Simulation

In model (1), the incidence rates for transmission of virus-to-cell and cell-to-cell are taken as the general form. By choosing four specific types of functions, i.e., bilinear incidence, saturation incidence, Beddington–DeAngelis response, and Hattaf–Yousfi response, we conduct numerical simulations to investigate the complex dynamics that model (1) can have, when taking time delay τ as the bifurcation parameter. The values of parameters, d3=2.4, d4=1.618, p=0.812, and k=200, remain the same in the following different cases.Case 1. Consider the functions of f1u,v=β1u and f2u,w=β2u, then the incidence rates have the bilinear forms of β1uv and β2uw. We take the other parameters as c=0.05, d1=0.01, d2=0.4, β1=0.00025, and β2=0.00065. Figure 1 presents the different results for varying producing rate of uninfected cells, plotted with s=2 (Figure 1(a)) and s=10 (Figure 1(b)), which shows the effect of parameter s on dynamical behavior of the model in this case. In this figure, we can see that there exists periodic solutions bifurcated from the infection equilibrium. The maximal and minimal values of vt are denoted by the two curves, and the line implies the local stability of E∗. It is clearly shown that the value of vt increases with that of s. Figure 2 shows that the infection equilibrium E∗ is asymptotically stable when τ=15 and a periodic solution exists when τ=20, which corresponds to the bifurcation diagram in Figure 1(b).Case 2. The type of saturation incidence is used in this case, i.e., f1u,v=β1u/1+av and f2u,w=β2u/1+bw; then the dynamical behaviors of the model are simulated in Figures 3 and 4, with the parameters taken as c=1, d1=0.01, and d2=0.03. In Figures 3(a), we set s=10, β1=0.002, β2=0.003, a=1, and b=1, and in Figure 3(b), we set s=20, β1=0.02, β2=0.03, a=0.1, and b=0.4. By the set of parameters, Figure 3(a) shows that the time delay does cause a bifurcation at τ=4.7403. The line is given by v∗=198.0481, and the infected steady state E∗ is locally asymptotically stable when τ<4.7403. More complicated dynamics for model (1) are caused by the delay perturbation in Figure 3(b). A period-three solution occurs for τ=35 and τ=45, which are displayed clearly by time series of vt for the time interval 5000,5500 and their corresponding phase portraits of vt and ut plotted in Figure 4, respectively.Case 3. Beddington–DeAngelis response is used in this case; here, f1u,v=β1u/1+a1u+a2v and f2u,w=β2u/1+a1u+a2w. By this type of incidence function, chaotic motions also occur when increasing τ from 0 to 50 with s=20, c=1, d1=0.1, d2=0.3, β1=0.02, and β2=0.03. When altering parameter values of a1 and a2, it can lead to the different bifurcation diagrams, which can be seen in Figures 5 and 6. We take a1=0.08, a2=0.5 and a1=0.01, a2=0.2, respectively, and then the simulation results in Figure 5 show the effect of these two parameters on the system. From Figure 5, we can see that model (1) with these parameter values can experience complicated dynamics. Furthermore, for τ=10,22,and45 in Figure 5(b), the solution of the model is presented by its orbit, plot of time series, and phase portrait correspondingly in Figure 6. For τ=10, there exists a Hopf bifurcation, and when it increases to τ=22, chaotic motion occurs, shown in Figures 6(a) and 6(b). The plot of time series of vt and phase portrait of vt and ut are given to illustrate the dynamical behaviors when τ=45 in Figures 6(c) and 6(d), showing that there exists a period-five solution.Case 4. When Hattaf–Yousfi response is used in this case, the functions take the form of f1u,v=β1u/1+a1u+a2v+a3uv and f2u,w=β2u/1+a1u+a2w+a3uw, respectively. If the parameter values are chosen as s=15, c=1, d1=0.01, d2=0.03, β1=0.5, β2=0.3, a1=1, a2=1, and a3=1, only Hopf bifurcation exists and no chaotic motions occur when τ∈0,50, shown in Figure 7(a). For a different set of values, let s=20, c=1, d1=0.1, d2=0.3, β1=0.02, β2=0.03, a1=0.01, a2=0.3, and a3=0.001, then the resulting bifurcation diagram in Figure 7(b) shows that system (1) can have solution of period three. Figure 8 gives the time series of vt and phase portrait of vt and ut for τ=49 to illustrate the solution of system (1), which corresponds to Figure 7(b).Figure 1
Bifurcation diagrams for model (1) with bilinear incidence rate. (a) s=2. (b) s=10. Other parameters are the same.
(a)
(b)Figure 2
Solutions of the model forτ=15 (a, b) and τ=20 (c, d), corresponding to Figure 1(b).
(a)
(b)
(c)
(d)Figure 3
Bifurcation diagrams for model (1) with saturation incidence rate. (a) s=10, β1=0.002, β2=0.003, a=1, and b=1. (b) s=20, β1=0.02, β2=0.03, a=0.1, and b=0.4. Other parameters are the same.
(a)
(b)Figure 4
Solutions of the model forτ=35 (a, b) and τ=45 (c, d), corresponding to Figure 3(b).
(a)
(b)
(c)
(d)Figure 5
Bifurcation diagrams for model (1) with Beddington–DeAngelis response. (a) a1=0.08 and a2=0.5. (b) a1=0.01 and a2=0.2. Other parameters are the same.
(a)
(b)Figure 6
Solutions of the model forτ=10 (a), τ=22 (b), and τ=45 (c, d), corresponding to Figure 5(b).
(a)
(b)
(c)
(d)Figure 7
Bifurcation diagrams for model (1) with Hattaf–Yousfi response. (a) s=15, c=1, d1=0.01, d2=0.03, β1=0.5, β2=0.3, a1=1, a2=1, and a3=1. (b) s=20, c=1, d1=0.1, d2=0.3, β1=0.02, β2=0.03, a1=0.01, a2=0.3, and a3=0.001.
(a)
(b)Figure 8
Solutions of the model forτ = 49, corresponding to Figure 7(b).
(a)
(b)


## 6. Discussion

To model the mechanisms of infectious diseases mathematically and explore dynamical behaviors of infection processes, some types of functional incidence have been used in [10, 12–15], which play an important role in determining qualitative behaviors of the proposed models and in giving reasonable descriptions of the dynamics. General incidence rate has been introduced into many epidemic models with an aim to include different situations as much as possible. In our model, we incorporate the effect of time delay needed to activate the immune response for the virus and two general incidence functions for the transmission of virus-to-cell and cell-to-cell. The basic reproductive number, which works as a crucial threshold and determines the dynamics of the model, has been defined as the sum of two parts related with infections of virus-to-cell and cell-to-cell, respectively.In this paper, we are concerned with the stabilities of two equilibria and the existence of Hopf bifurcation through which the positive equilibrium loses its stability and periodic solutions occur. In the analysis, the time delay is chosen as the bifurcation parameter, which can destabilize the positive equilibrium when it increases. By using the characteristic equation with delay-dependent parameters, normal theory, and center manifold theorem, the existence of pure imaginary roots is verified and then the system experiences Hopf bifurcation. To illustrate the dynamical behavior of stability switches, simulations are conducted to show the process numerically when taking general incidence function as four specific types of bilinear incidence, saturation incidence, Beddington–DeAngelis response, and Hattaf–Yousfi response. In each case, the dynamical behavior is simulated when increasingτ from 0 to 50, and two sets of parameters are taken to compare the bifurcation results. Furthermore, for some fixed values of τ, we also present the plot of time series, phase portrait, or solution orbit to show the complicated dynamics. The effect of delay perturbation and different forms of incidence can be seen in the figures.It should be noted that when other factors are taken into consideration for more realistic mechanism, such as the delays describing the intracellular latency for virus-to-cell infection and cell-to-cell infection, or varying producing rate of the uninfected cells instead of constant, the model can be extended reasonably and the dynamical analysis will become more challenging than the present one, which we may discuss in future work.



---

*Source: 1989651-2019-07-09.xml*