# Stability and Bifurcation Analysis for a Delay Differential Equation of Hepatitis B Virus Infection

Authors: Xinchao Yang; Xiju Zong; Xingong Cheng; Zhenlai Han
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

The stability and bifurcation analysis for a delay differential equation of hepatitis B virus infection is investigated. We show the existence of nonnegative equilibria under some
appropriated conditions. The existence of the Hopf bifurcation with delayτ at the endemic equilibria is established by analyzing the distribution of the characteristic values. The explicit formulae which determine the direction of the bifurcations, stability, and the other properties of the bifurcating periodic solutions are given by using the normal form theory and the center manifold theorem. Numerical simulation verifies the theoretical results.


---

## Body

## 1. Introduction

Recently, a hepatitis B virus (HBV) model with time delay that was proposed and investigated in the literature [1–4] caught the attention of a lot of mathematicians. In practice, the HBV model has suffered time delay caused by the HBV incubation period, which varies from 45 to 180 days, and the delay in viral shedding which both suggest that viral production delay may significantly impact infection dynamics [1]. Precisely, the HBV model with time delay reads as the following:
(1)dx(t)dt=λ-dx(t)-βx(t)v(t)x(t)+e(t)+y(t),de(t)dt=-de(t)+βx(t)v(t)x(t)+e(t)+y(t)-βe-dτx(t-τ)v(t-τ)x(t-τ)+e(t-τ)+y(t-τ),dy(t)dt=βe-dτx(t-τ)v(t-τ)x(t-τ)+e(t-τ)+y(t-τ)-ay(t),dv(t)dt=ky(t)-μv(t),
where x(t) and y(t) represent the number of uninfected cells and infected cells, respectively. e(t) represents the number of exposed cells, that is, the cells that have acquired the virus but are not yet producing new virions. v(t) denotes the number of free virions. τ is the time delay for virion production. Here, the positive constant λ is the rate at which new uninfected live cells are generated. The positive constant d is the per capita death rate of uninfected live cells. Infected live cells are killed by immune cells at rate ay and produce free virions at rate ky, where k is what so-called “burst” constant. Free virions are cleared by lymphatic and other mechanisms at rate μv, where μ is a constant. β>0 is an incidence rate coefficient describing the infection process. The initial conditions for the system (1) are
(2)x(s)=x0(s),y(s)=y0(s),v(s)=v0(s),s∈[-τ,0],E0=β∫-τ0edsv0(s)x0(s)x0(s)+y0(s)+e0(s)ds,
where x0, y0, and v0 are nonnegative functions. Based on some observations of virus particles v, the system (1) is simplified in [1] as the following:
(3)dx(t)dt=λ-dx(t)-βkx(t)y(t)μ[x(t)+y(t)],de(t)dt=-de(t)+βkx(t)y(t)μ[x(t)+y(t)]-βke-dτx(t-τ)y(t-τ)μ[x(t-τ)+y(t-τ)],dy(t)dt=βke-dτx(t-τ)y(t-τ)μ[x(t-τ)+y(t-τ)]-ay(t).A direct computation shows that the basic infection reproduction number for the system (2) is
(4)R0=βke-dτaμ.For the sake of simplicity, letx-=x-(λ/d), and the system (3) is equivalent to the following system:
(5)dx-(t)dt=-dx-(t)-βk[x-(t)+(λ/d)]y(t)μ[x-(t)+(λ/d)+y(t)],de(t)dt=-de(t)+βk[x-(t)+(λ/d)]y(t)μ[x-(t)+(λ/d)+y(t)]-βke-dτ[x-(t-τ)+(λ/d)]y(t-τ)μ[x-(t-τ)+(λ/d)+y(t-τ)],dy(t)dt=βke-dτ[x-(t-τ)+(λ/d)]y(t-τ)μ[x-(t-τ)+(λ/d)+y(t-τ)]-ay(t),
which has two equilibria: the infection-free equilibrium Ef=(0,0,0)T and the infected equilibrium E*=(x*,e*,y*)T, where
(6)x*=λd+aedτ(R0-1)-λd,e*=a(edτ-1)dx*,y*=(R0-1)x*.The following results, Theorems1 and 2, come from [1].Theorem 1.
IfR0<1, the infection-free equilibrium Ef of the system (5) is globally asymptotically stable.Theorem 2.
IfR0>1, the chronic infected equilibrium E* of the system (5) is locally asymptotically stable.The initial conditions for the system (5) are
(7)(φ1(θ),φ2(θ),ψ(θ))∈C+=C([-τ,0],R+03),φk(0)>0,ψ(0)>0,k=1,2,
where R+03={(x1,x2,x3)∈R3:xk>0,k=1,2,3}.It is straightforward to show the following.Lemma 3.
The solution of (5) with an initial condition (6) is nonnegative for all t≥0.It is well known that the studies on the dynamical systems not only include the discussion of stabilities, attractivity, and persistence, but also include many dynamical behaviors such as periodic solutions, bifurcations, and chaos. Particularly, the properties of periodic solutions appearing through the Hopf bifurcation in delayed systems are of great interest [5–7]. In the present paper, our main objective is to investigate the bifurcation phenomena of the modified hepatitis B virus (HBV) model with time delay τ.This paper is organized as follows. In Section2, by analyzing the characteristic equation of the linearized system of the system (5) at the equilibria, we discuss the stability of the origin and the positive equilibrium and the existence of the Hopf bifurcations occurring at the chronic infected equilibrium. In Section 3, the formulae determining the direction of the Hopf bifurcations and the stability of bifurcating periodic solutions on the center manifold are obtained by using the normal form theory and the center manifold theorem due to Hassard et al. [8]. To verify the obtained results, some numerical simulations are included in Section 4. The paper ends with a brief discussion.


## 2. Stability of Equilibria and Existence of the Hopf Bifurcation

In this section, we will investigate the stability of the equilibria and the existence of the Hopf bifurcations occurring at the chronic infected equilibrium. Then, it is easy to check that the system (5) has an equilibrium Ef(0,0,0) for all nonnegative parameters. The characteristic equation of (5) at Ef is
(8)(z+d)2(z-βke-dτμ+a)=0.Hence,Ef is a saddle with dimWu(Ef)=1,dimWs(Ef)=2 for R0<1; Wu(Ef) and Ws(Ef) are the local unstable and stable manifolds of Ef, respectively. Ef is locally asymptotically stable for R0<1. In fact, Ef is globally asymptotically stable for R0<1, see [1, 4].Now, we will investigate the stability of the chronic infected equilibriumE*. Linearizing the system (5) at E* yields the following linear system:
(9)dx-(t)dt=-dx-(t)-{βky*(t)μ[x*(t)+(λ/d)+y*(t)]-βk[x*(t)+(λ/d)]y*(t)μ[x*(t)+(λ/d)+y*(t)]2}x-(t)-{βk[x*(t)+(λ/d)]μ[x*(t)+(λ/d)+y*(t)]-βk[x*(t)+(λ/d)]y*(t)μ[x*(t)+(λ/d)+y*(t)]2}y(t),de(t)dt=-de(t)+{βky*(t)μ[x*(t)+(λ/d)+y*(t)]-βk[x*(t)+(λ/d)]y*(t)μ[x*(t)+(λ/d)+y*(t)]2}x-(t)+{βk[x*(t)+(λ/d)]μ[x*(t)+(λ/d)+y*(t)]+βk[x*(t)+(λ/d)]y*(t)μ[x*(t)+(λ/d)+y*(t)]2}y(t)-{βke-dτy*(t)μ[x*(t)+(λ/d)+y*(t)]-βke-dτ[x*(t)+(λ/d)]y*(t)μ[x*(t)+(λ/d)+y*(t)]2}x-(t-τ)-{βke-dτ[x*(t)+(λ/d)]μ[x*(t)+(λ/d)+y*(t)]-βke-dτ[x*(t)+(λ/d)]y*(t)μ[x*(t)+(λ/d)+y*(t)]2}y(t-τ),dy(t)dt={βke-dτy*(t)μ[x*(t)+(λ/d)+y*(t)]-βke-dτ[x*(t)+(λ/d)]y*(t)μ[x*(t)+(λ/d)+y*(t)]2}x-(t-τ)-{βke-dτ[x*(t)+(λ/d)]μ[x*(t)+(λ/d)+y*(t)]-βke-dτ[x*(t)+(λ/d)]y*(t)μ[x*(t)+(λ/d)+y*(t)]2}y(t-τ)-ay(t),
whose characteristic equation reads as
(10)z3+(2d+A+a)z2+[d(2a+A+d)+aA]z+ad(d+A)=Be-zτe-dτ[z2+2dz+d2],
with
(11)A=βky*(t)μ[x*(t)+(λ/d)+y*(t)]-βk[x*(t)+(λ/d)]y*(t)μ[x*(t)+(λ/d)+y*(t)]2,B=βk[x*(t)+(λ/d)]μ[x*(t)+(λ/d)+y*(t)]-βk[x*(t)+(λ/d)]y*(t)μ[x*(t)+(λ/d)+y*(t)]2.Forτ=0, characteristic equation (10) reduces to
(12)z3+(2d+A+a)z2+[d(2a+A+d)+aA]z+ad(d+A)=B(z2+2dz+d2).By the Routh-Hurwitz criterion, we know that all the roots of (12) have negative real parts, that is, the chronic infected equilibrium E* is locally asymptotically stable for τ=0. We now give a definition, which can be found in [2, 9].Definition 4.
The equilibriumE* is called asymptotically stable if there is an ε>0 such that
(13)supτ≤θ≤∞[|φ1(θ)-x0*|+|φ2(θ)-x1*|+|ψ(θ)-y*|]<ε
implies that
(14)limt→∞(x0(t),x1(t),y(t))=(x0*,x1*,y*),
where (x0(t),x1(t),y(t)) is the solution of the system (5) with an initial condition (6).
Song et al. [10] investigated the distribution of roots of the following equation:
(15)z3+m2z2+m1z+m0+e-zτ(n2z2+n1z+n0)=0,
where mj,nj∈R(j=0,1,2) and ∑j=02‍nj2≠0. When τ=0, (15) reduces to
(16)z3+(m2+n2)z2+(m1+n1)z+n0=0.
Obviously,iω(ω>0) is a root of (16) if and only if ω satisfies
(17)-ω3i-m2ω2+m1ωi+m0+(-n2ω2+n1ωi+n0)(cosωτ-isinωτ)=0.
Separating the real and the imaginary parts, we have(18)m2ω2-m0=-n2ω2cosωτ+n1ωsinωτ+n0cosωτ,-ω3+m1ω=-n2ω2sinωτ-n1ωcosωτ+n0sinωτ.
Squaring the two equations and adding them give(19)ω6+pω4+qω2+r=0,
where p=m22-2m1-n22, q=m12-2m0m2+2n0n2-n12, and r=m02-n02. Song et al. [10] obtained the following results on the distribution of roots of (15) and (19).Lemma 5.
For the polynomial equation (19),(i)
ifr<0, then (19) has at least one positive root,(ii)
ifr≥0 and Δ=p2-3q≤0, then (19) has no positive roots,(iii)
ifr≥0 and Δ=p2-3q>0, then (19) has positive roots if and only if z1*=(-p+Δ)/3 and h(z1*)≤0, where
(20)h(z)=z3+pz2+qz+r.Lemma 6.
For the transcendental equation (19),(i)
ifr≥0 and Δ=p2-3q≤0, then all roots with positive real parts of (19) have the same sum as those of the polynomial equation (16), for all τ,(ii)
ifr<0 or r≥0, Δ=p2-3q>0,z1*=(-p+Δ)/3>0 and h(z1*)≤0, then all roots with positive real parts of (19) have the same sum as those of the polynomial equation (16), for τ∈[0,τ0).From Lemmas5 and 6, we can have the following lemma.Lemma 7.
(i) The chronic infected equilibriumE* of the system (5) is absolutely stable if and only if the equilibrium E* of the corresponding ordinary differential equation (ODE) system is asymptotically stable, and the characteristic equation (10) has no purely imaginary roots for any τ>0.
(ii) The chronic infected equilibriumE* of the system (5) is conditionally stable if and only if all roots of the characteristic equation (10) have negative real parts at τ=0 such that the characteristic equation (10) has a pair of purely imaginary roots iω0.
Then, one turns to an investigation of local stability of the chronic infected equilibriumE* in the case of R0<1.Theorem 8.
ForR0<1 holds, there exists a sequence of values for τ:
(21)τ0<τ1<⋯<τk<⋯,
such that (10) has a pair of purely imaginary roots iω0 when τ=τk, k=0,1,2,…. That is, the chronic infected equilibrium E* of the system (5) is conditionally stable.Proof.
From the above arguments, we know that all roots of characteristic equation (10) have negative real parts at τ=0. Next, we will show that there is a unique pair of purely imaginary roots iω0(ω0>0) for characteristic equation (10).
Assume that for someτ>0, iω(ω>0) is a root of (10), which implies that
(22)-iω3-(2d+A+a)ω2+[d(2a+A+d)+aA]iω+ad(d+A)=Be-iωτe-dτ[-ω2+2diω+d2].
Note thatA>0, B>0 because of the positivities of the parameters λ, μ, κ, β, d and the properties x*(t)+(λ/d)≤[x*(t)+(λ/d)+y*(t)], y*(t)≤[x*(t)+(λ/d)+y*(t)].
Separating the real and imaginary parts and using Euler’s formula give(23)-(2d+A+a)ω2+ad(d+A)=B(d2-ω2)e-dτcosωτ+2Bdωe-dτsinωτ,-ω3+[d(2a+A+d)+aA]ω=2Bdωe-dτcosωτ-B(d2-ω2)e-dτsinωτ,
which is equivalent to
(24)ω6+pω4+qω2+r=0,
where
(25)p=(2d+A+a)2-2[d(2a+A+d)+aA]=A2+2d2+a2+2dA>0,q=[d(2a+A+d)+aA]2-2(2d+A+a)ad(d+A),r=-e-2dτ[B2(d2-ω2)2+8B2d2]<0.
In order to solve (24), we first consider the following:
(26)ω6+pω4+qω2+r0=0,
where
(27)r0=-[B2(d2-ω2)2+8B2d2]<0.
By Lemma5, there is a unique positive ω0 satisfying (26). From (26), we get the corresponding τk′>0 such that (26) has a pair of purely imaginary roots
(28)τk′=1ω0arccos{{(B2(d2-ω2)2+8B2d2)-1}(ad(d+A)-(2d+A+a)ω02-ω03+[d(2a+A+d)+aA]ω0{(2d+A+a)ω02})×(B2(d2-ω2)2+8B2d2)-1}+2kπω0(k=0,1,2,…).
Therefore, by using Rouché’s theorem [3], there is a unique positive τk=τk′+o(1/k) satisfying (26), that is, the characteristic equation (10) has a pair of purely imaginary roots of the form ±iω0 as k→∞. By Lemma 7, we complete the proof of Theorem 8.
Next, we turn to show that(29)d(Reλ)dτ|τ=τk>0.
This will signify that there exists at least one eigenvalue with positive real partτ=τk. We first consider the following:
(30)z3+(2d+A+a)z2+[d(2a+A+d)+aA]z+ad(d+A)=Be-zτ[z2+2dz+d2].
Differentiating (30) with respect to τ, we have
(31)(dλdτ)-1=3z2+2(A+2d+a)z+[d(2a+A+d)+aA]-Bz(d+z)2e-zτ+τe-zτB(z+d)2-Bz(d+z)2e-zτ-2e-zτB(z+d)-Bz(d+z)2e-zτ=((d+z){3z2+2(A+2d+a)z+[d(2a+A+d)+aA]{(d+z){3z2+2(A+2d+a)z}})×(-Bz(d+z)3e-zτ)-1-2e-zτB(z+d)2-Bz(d+z)3e-zτ-τz=-(d{3z2+2(A+2d+a)z+[d(2a+A+d)+aA]{{3z2+2(A+2d+a)z}})×((d+z)[z3+(2d+A+a)z2+[d(2a+A+d)+aA]z+ad(d+A){[z3+(2d+A+a)z2}])-1-(z2-[d(2a+A+d)+aA]z-2ad(d+A))×((d+z)[z3+(2d+A+a)z2+[d(2a+A+d)+aA]z+ad(d+A){[z3+(2d+A+a)z2}])-1-τz.
For the sake of simplicity, leta1=A+2d+a, a2=d(2a+A+d)+aA, and a3=ad(d+A), and (28) can be written as follows:
(32)(dλdτ)-1=-d{3z2+2a1z+a2}(d+z)[z3+a1z2+a2z+a3]-z2-a2z-2a3(d+z)[z3+a1z2+a2z+a3]-τz=-(3d+1)z2-(2a1d-a2)z-(a2d-2a3)(d+z)[z3+a1z2+a2z+a3]-τz.
Therefore,(33)sign{d(Reλ)dτ|τ=τk}=sign{Re(dλdτ)-1}|τ=τk=signRe{-(3d+1)z2-(2a1d-a2)z-(a2d-2a3)(d+z)[z3+a1z2+a2z+a3]-τz{-(3d+1)z2-(2a1d-a2)z-(a2d-2a3)(d+z)[z3+a1z2+a2z+a3]}}=signRe{({×[(a32-a1ω02)2+(a2ω0-ω03)2])-1}([(3d+1)ω02-i(2a1d-a2)ω0-(a2d-2a3){[(3d+1)ω02-i(2a1d-a2)ω0}](d-iω0))×({×[(a32-a1ω02)2+(a2ω0-ω03)2])-1}(d2+ω02)×[(a32-a1ω02)2+(a2ω0-ω03)2])-1)×[(a32-a1ω02)-i(a2ω0-ω03)]}=sign{(3d+1)ω06+p2ω04+Bd2ω02(d2+ω02)[(a32-a1ω02)2+(a2ω0-ω03)2]}>0.
This root of characteristic equation (9) crosses the imaginary axis from the left to the right as τ continuously varies from a number less than τk to one greater than τk again by Rouché’s theorem [3]. Therefore, the transversality condition holds, and the conditions for Hopf bifurcation [11] are then satisfied at τ=τk. In conclusion, we have the following stability and bifurcation results to (5).Theorem 9.
Suppose thatR0<1 holds. Then, for each fixed τ>0, there exists a sequence of values for τ:
(34)τ0<τ1<⋯<τk<⋯,
such that the positive equilibrium E* is asymptotically stable when τ∈[0,τ0), and unstable when τ>τ0. Furthermore, (5) undergoes a Hopf bifurcation at E* when τ=τk=τk′+o(1/k),k=0,1,2,…, where τk′ is defined by (28).


## 3. Properties of the Hopf Bifurcations

In this section, we will study the properties of the Hopf bifurcations by using the normal theory and the center manifold theorem due to Hassard et al. [8]. Let x1=x--x*, x2=e-e*, x3=y-y*, x~i(t)=xi(τt), τ=τk+ϑ, and ϑ∈R=(-∞,∞). Then, ϑ=0 is the Hopf bifurcation value of system. We drop the tildes for simplification of notations, then the system (5) becomes an functional differential equation in C([-1,0],R3) as
(35)x˙(t)=Lϑ(xt)+f(ϑ,xt),
where C([-1,0],R3) is the Banach space of continuous functions mapping the interval [-τ,0] into R3, x(t)=(x1(t),x2(t),x3(t))T∈R3, xt(s)=x(t+s) for s∈[-τ,0] and Lϑ:C→R3, f:R×C→R3 are read, respectively, as
(36)Lϑ(ϕ)=(τk+ϑ)(-(d+A)0-BA-dB00-a)(ϕ1(0)ϕ2(0)ϕ3(0))Lϑ(ϕ)=+(τk+ϑ)(000-Ae-dτk0Be-dτkAe-dτk0-Be-dτk)(ϕ1(-1)ϕ2(-1)ϕ3(-1)),(37)f(ϑ,xt)=(τk+ϑ)×(-βk[ϕ1(0)+(λ/d)]ϕ3(0)μ[ϕ1(0)+(λ/d)+ϕ3(0)]βk[ϕ1(0)+(λ/d)]ϕ3(0)μ[ϕ1(0)+(λ/d)+ϕ3(0)]-βke-dτk[ϕ1(-1)+(λ/d)]ϕ3(-1)μ[ϕ1(-1)+(λ/d)+ϕ3(-1)]βke-dτk[ϕ1(-1)+(λ/d)]ϕ3(-1)μ[ϕ1(-1)+(λ/d)+ϕ3(-1)]).Obviously, Lϑ is a continuous linear function mapping C into R3, by the Riesz representation theorem, there exists a 3×3 matrix function η(s,μ) of bounded variation for s∈[-1,0], such that
(38)Lϑ(ϕ)=∫-10dη(s,μ)ϕ(s).In fact, we can choose(39)η(s,ϑ)=(τk+ϑ)(-(d+A)0-BA-dB00-a)δ(s)+(τk+ϑ)(000-Ae-dτk0Be-dτkAe-dτk0-Be-dτk)δ(s+1),
where δ(·) denotes the Dirac delta function.Ifϕ is any given function in C and x(ϕ) is the unique solution of the linearized equation x˙(t)=Lϑ(xt) of (35) with initial function ϕ at zero, then the solution operator T(t):C→C is defined by T(t)ϕ=xt(ϕ), for all t≥0. It is obvious that T(t), t≥0, is a strongly continuous semigroup of linear transformation on [0,+∞) and the infinitesimal generator Aϑ of T(t),t≥0, is
(40)Aϑ={dϕ(s)ds,s∈[-1,0)∫-10dη(s,μ)ϕ(s),s=0.
For ϕ∈C1([-1,0],R3), the space of functions mapping the interval [-1,0] into R3 which have a continuous first derivative also defines
(41)Rϑ(ϕ)={0,s∈[-1,0)f(ϑ,ϕ),s=0.Then, system (35) is equivalent to
(42)x˙(t)=Aϑxt+Rϑxt.Forψ∈C1([-1,0],(R3)*), the space of functions mapping interval [0,1] into the three-dimensional row vectors which have continuous first derivative defines
(43)Aϑ*={-dψ(s)ds,s∈(0,1]∫-10dηT(t,0)ϕ(-t),s=0,
and a bilinear inner product
(44)〈ψ(ζ),ϕ(s)〉=ψ-(0)ϕ(0)-∫-10∫ξ=0sψ-(ξ-s)dη(s)ϕ(ξ)dξ,
where η(s)=η(s,0). Then, A(0) and A* are adjoint operators. By the discussion in Section 2, we know that z are eigenvalues of A(0). Hence, they are also eigenvalues of A*. We first need to compute the eigenvectors of A(0) and A* corresponding to iω0τk and -iω0τk, respectively. Suppose that q(s)=(1,a,b)Teiω0τks is the eigenvectors of A(0) corresponding to iω0τk, then A(0)q(s)=iω0τkq(s). Then, from the definition of A(0) and (36), (38), and (39), we have
(45)τk(-(d+A)0-BA-dB00-a)q(0)+τk(000-Ae-dτk0Be-dτkAe-dτk0-Be-dτk)q(-1)=iω0τkq(0).Forq(-1)=q(0)e-iω0τk, then we obtain
(46)ρ=-1+1B+A(d-iω0)(a+iω0)B(d2+ω02),σ=-A+d+iω0B.Similarly, we can obtain the eigenvectorq*(s)=D(1,ρ*,σ*)eiω0τk of A* corresponding to -iω0τk, where
(47)ρ*=ρ=-1+1B+A(d-iω0)(a+iω0)B(d2+ω02),σ*=σ=-A+d+iω0B.We need to determine the value ofD to ensure that 〈q*(ς),q(s)〉=1. By (44), we have
(48)〈q*(ς),q(s)〉=D-(1,ρ*,σ*)(1,ρ,σ)T-∫-10∫ξ=0sD-(1,ρ*,σ*)e-iω0τk(ξ-s)dη(s)(1,ρ,σ)T×eiω0τkξdξ=D-[{∫-10‍}1+|ρ|2+|σ|2-∫-10‍(1,ρ*,σ*)seiω0τksdη(s)(1,ρ,σ)T]=D-[1+|ρ|2+|σ|2+e-dτke-iω0τk(A-Bσ)(ρ-*-σ-*)].Therefore, we can chooseD as
(49)D=11+|ρ|2+|σ|2+e-dτke-iω0τk(A-Bσ)(ρ-*-σ-*).Next, we will compute the coordinate to describe the center manifoldC0 at ϑ=0. Let xt be the solution of (42) when ϑ=0. Define(50)υ(t)=〈q*(ς),xt〉,W(t,s)=xt(s)-2Re{υ(t)q(s)}.On the center manifoldC0, we have
(51)W(t,s)=W(υ(t),υ-(t),s),
where
(52)W(υ(t),υ-(t),s)=W20(s)υ22+W10(s)υυ-+W02(s)υ-22+⋯.υ and υ- are local coordinates for the center manifold C0 in the direction of q* and q-*. Note that W is real if xt is real. We only consider real solutions. For solution xt∈C0 of (42), since ϑ=0, we have
(53)υ˙(t)=iω0τkυ(t)+q-*(0)f(0,W(υ,υ-,0)+2Re{υq(s)}):=iω0τkυ(t)+q-*(0)f0(υ,υ-).For convenience of calculation, let(54)q-*(0)f0(υ,υ-)=f20υ22+f11(s)υυ-+f02(s)υ-22+f21υ2υ-2+⋯.On the other hand, it follows from (50) and (52) that
(55)xt(s)=W(υ,υ-,0)+2Re{υq(s)}=W20(s)υ22+W10(s)υυ-+W02(s)υ-22+(1,a,b)Teiω0τksυ+(1,a-,b-)Te-iω0τksυ-+⋯.Notice that(56)ϕ1(0)=q1(0)v+q-1(0)v-+W20(1)(0)v22+W11(1)(0)vv-+W02(1)(0)v-22+⋯,(57)ϕ2(0)=q2(0)v+q-2(0)v-+W20(2)(0)v22+W11(2)(0)vv-+W02(2)(0)v-22+⋯,(58)ϕ3(0)=q3(0)v+q-3(0)v-+W20(3)(0)v22+W11(3)(0)vv-+W02(3)(0)v-22+⋯,(59)ϕ1(-1)=q1(0)ve-iω0*τ0*+q-1(0)v-eiω0*τ0*+W20(1)(-1)v22+W11(1)(-1)vv-+⋯,(60)ϕ2(-1)=q2(0)ve-iω0*τ0*+q-2(0)v-eiω0*τ0*+W20(2)(-1)v22+W11(2)(-1)vv-+⋯,(61)ϕ3(-1)=q3(0)ve-iω0*τ0*+q-3(0)v-eiω0*τ0*+W20(3)(-1)v22+W11(3)(-1)vv-+⋯,Taking (37) into account, we can obtain that(62)q-*(0)f0(υ,υ-)=q-*(0)f(0,xt)=τkD-(1,ρ-*,σ-*)(-βk[ϕ1(0)+(λ/d)]ϕ3(0)μ[ϕ1(0)+(λ/d)+ϕ3(0)]βk[ϕ1(0)+(λ/d)]ϕ3(0)μ[ϕ1(0)+(λ/d)+ϕ3(0)]-βke-dτk[ϕ1(-1)+(λ/d)]ϕ3(-1)μ[ϕ1(-1)+(λ/d)+ϕ3(-1)]βke-dτk[ϕ1(-1)+(λ/d)]ϕ3(-1)μ[ϕ1(-1)+(λ/d)+ϕ3(-1)])=τkD-{-βk[ϕ1(0)+(λ/d)]ϕ3(0)μ[ϕ1(0)+(λ/d)+ϕ3(0)]+ρ-*βk[ϕ1(0)+(λ/d)]ϕ3(0)μ[ϕ1(0)+(λ/d)+ϕ3(0)]-ρ-*βke-dτk[ϕ1(-1)+(λ/d)]ϕ3(-1)μ[ϕ1(-1)+(λ/d)+ϕ3(-1)]+σ-*βke-dτk[ϕ1(-1)+(λ/d)]ϕ3(-1)μ[ϕ1(-1)+(λ/d)+ϕ3(-1)]}=(ρ*-1)τkD-βkμΔ0[λd+q1(0)v+q-1(0)v-+W20(1)(0)v22+W11(1)(0)vv-+W02(1)(0)v-22+⋯]×[q3(0)v+q-3(0)v-+W20(3)(0)v22+W11(3)(0)vv-+W02(3)(0)v-22+⋯]×[λd+q1(0)v+q-1(0)v-+W20(1)(0)v22+W11(1)(0)vv-+W02(1)(0)v-22+⋯+q3(0)v+q-3(0)v-+W20(3)(0)v22+W11(3)(0)vv-+W02(3)(0)v-22+⋯]+(σ*-ρ*)τkD-βke-dτkμΔ-1×[λd+q1(0)ve-iω0*τ0*+q-1(0)v-eiω0*τ0*+W20(1)(0)v22+W11(1)(0)vv-+⋯]×[q3(0)ve-iω0*τ0*+q-3(0)v-eiω0*τ0*+W20(3)(0)v22+W11(3)(0)vv-+W02(3)(0)v-22+⋯]×[λd+q1(0)e-iω0*τ0*+q-1(0)v-eiω0*τ0*+W20(1)(0)v22+W11(1)(0)vv-+W02(1)(0)v-22+⋯+q3(0)ve-iω0*τ0*+q-3(0)v-eiω0*τ0*+W20(3)(0)v22+W11(3)(0)vv-+W02(3)(0)v-22+⋯],where(63)Δ0=|λd+q1(0)v+q-1(0)v-+W20(1)(0)v22+W11(1)(0)vv-+W02(1)(0)v-22+⋯+q3(0)v+q-3(0)v-+W20(3)(0)v22+W11(3)(0)vv-+W02(3)(0)v-22+⋯|2,Δ-1=|λd+q1(0)ve-iω0*τ0*+q-3(0)v-eiω0*τ0*+W20(1)(0)v22+W11(1)(0)vv-+W02(1)(0)v-22+⋯+q3(0)v+q-3(0)v-+W20(3)(0)v22+W11(3)(0)vv-+W02(3)(0)v-22+⋯|2.Straightforward calculation leads to(64)q-*(0)f0(υ,υ-)=(ρ*-1)βkμΔ0×[λd(q1(0)+q3(0))q3(0)v2+λdq1(0)q3(0)v2+(λd)2W20(3)(0)v22+λd(q1(0)+q3(0))q3(0)vv-+2Re{λdq1(0)q3(0)}vv-+(λd)2W11(3)(0)vv-+⋯]+(σ*-ρ*)τkD-βke-dτkμΔ-1×[λd(q1(0)+q3(0))q3(0)v2+λdq1(0)q3(0)v2+(λd)2W20(3)(-1)v22+λd(q1(0)+q3(0))q3(0)vv-+2Re{λdq1(0)q3(0)}vv-+(λd)2W11(3)(-1)vv-+⋯].Comparing the coefficients with (54), we have
(65)f20=(ρ*-1)βkμΔ0×[λd(q1(0)+q3(0))q3(0)+λdq1(0)q3(0)+12(λd)2W20(3)(0)]+(σ*-ρ*)τkD-βke-dτkμΔ-1×[λd(q1(0)+q3(0))q3(0)+λdq1(0)q3(0)+12(λd)2W20(3)(-1)],f11=λd(q1(0)+q3(0))q3(0)+2Re{λdq1(0)q3(0)}+(λd)2W11(3)(-1)+λd(q1(0)+q3(0))q3(0)+2Re{λdq1(0)q3(0)}+(λd)2W11(3)(-1),f02=(ρ*-1)βkμΔ0×[λd(q-1(0)+q-3(0))q-3(0)+λdq-1(0)q-3(0)+(λd)2W20(3)(0)v-22]+(σ*-ρ*)τkD-βke-dτkμΔ-1×[λd(q-1(0)+q-3(0))q-3(0)+λdq-1(0)q-3(0)+12(λd)2W20(3)(-1)],f21=(ρ*-1)βkμΔ0×[λd(q-1(0)+q-3(0))q-3(0)+λdq-1(0)q-3(0)+(λd)2W20(3)(0)v-22]+(σ*-ρ*)τkD-βke-dτkμΔ-1×[λd(q-1(0)+q-3(0))q-3(0)+λdq-1(0)q-3(0)+12(λd)2W20(3)(-1)].SinceW20(s) and W11(s) are in f21, we still need to compute them. From (40) and (44), we have
(66)W˙=x˙t-υ˙q-υ-˙q-={AW-2Re{q-*(0)f0q(s)},s∈[-1,0),AW-2Re{q-*(0)f0q(s)}+f0,s=0,
that is,
(67)W˙=AW+H(υ,υ-,s),
where
(68)H(υ,υ-,s)=H20(s)υ22+H10(s)υυ-+H02(s)υ-22+⋯.Substituting the corresponding series into (67) and comparing the coefficients, we have
(69)(A-2iω0τk)W20(s)=-H20(s),AW11(s)=-H11(s),….From (67), we know that for s∈[-1,0),
(70)H(υ,υ-,s)=-q-*(s)f0q(s)-q*(s)f-0q-(s)=-f(υ,υ-)q(s)-f-(υ,υ-)q-(s).Comparing the coefficients with (68) gives
(71)H20(s)=-f20q(s)-f-02q-(s),H11(s)=-f11q(s)-f-11q-(s).From the definition ofA and (69) and (71), we obtain
(72)W˙20(s)=2iω0τkW20(s)+f20q(s)+f-02q-(s).Forq(s)=(1,ρ,σ)Teiω0τks, hence
(73)W20(s)=if20ω0τkq(0)eiω0τks+if-023ω0τkq(0)e-iω0τksq-(s)+E1e2iω0τks,
where E1=(E1(1),E1(2),E1(3)) is a constant vector. Similarly, we know that
(74)W11(s)=-if11ω0τkq(0)eiω0τks+if-11ω0τkq(0)e-iω0τksq-(s)+E2,
where E2=(E2(1),E2(2),E2(3)) is a constant vector.In what follows, we will calculateE1 and E2. From the definition of A and (69), we have
(75)∫-10dηT(s)W20(s)=2iω0τkW20(0)-H20(0),(76)∫-10dηT(s)W11(s)=-H11(0),
where η(s)=η(s,0). When s=0, we obtain that
(77)H(υ,υ-,0)=-2Re{q-*(0)f0q(0)}+f0=-q-*(0)f0q(0)-q*(0)f-0q-(0)+f0,
that is,
(78)H20(s)υ22+H11(s)υυ-+H02(s)υ-22+⋯=-q(0)(f20υ22+f11(s)υυ-+f02(s)υ-22+f21υ2υ-2⋯)-q-(0)(f-20υ22+f-11(s)υυ-+f-02(s)υ-22+f21υ2υ-2⋯)+f0,
with
(79)f0=τk(-βk[ϕ1(0)+(λ/d)]ϕ3(0)μ[ϕ1(0)+(λ/d)+ϕ3(0)]βk[ϕ1(0)+(λ/d)]ϕ3(0)μ[ϕ1(0)+(λ/d)+ϕ3(0)]-βke-dτk[ϕ1(-1)+(λ/d)]ϕ3(-1)μ[ϕ1(-1)+(λ/d)+ϕ3(-1)]βke-dτk[ϕ1(-1)+(λ/d)]ϕ3(-1)μ[ϕ1(-1)+(λ/d)+ϕ3(-1)]).From (50), we have
(80)f0=(fv21fv22fv23)υ2+(fvv-1fvv-2fvv-3)vv-+⋯,
where
(81)fv21=-βkμΔ0[{if-023ω0τk}λd(q1(0)+q3(0))q3(0)+λdq1(0)q3(0)+12(λd)2×(if20ω0τkq(0)eiω0τks+if-023ω0τkq(0)e-iω0τksq-(s)+E1e2iω0τks{+if-023ω0τk})],(82)fv22=βkμΔ0[{if-023ω0τk}{if-023ω0τk}λd(q1(0)+q3(0))q3(0)+λdq1(0)q3(0)+12(λd)2×({if-023ω0τk}if20ω0τkq(0)eiω0τks+if-023ω0τkq(0)e-iω0τksq-(s)+E1e2iω0τks{if-023ω0τk})]+λd(q1(0)+q3(0))q3(0)+2Re{λdq1(0)q3(0)}+(λd)2W11(3)(-1)+λd(q1(0)+q3(0))q3(0)+2Re{λdq1(0)q3(0)}+(λd)2W11(3)(-1),(83)fv23=-βke-dτkμΔ-1[{if-023ω0τk}λd(q1(0)+q3(0))q3(0)+λdq1(0)q3(0)+12(λd)2×(if20ω0τkq(0)eiω0τks+if-023ω0τkq(0)e-iω0τksq-(s)+E1e2iω0τks{if-023ω0τk})].By (70) and (71), we have
(84)H20(s)=-f20q(s)-f-02q-(s)+fv2,H11(s)=-f11q(s)-f-11q-(s)+fvv-.Foriω0τk is the eigenvalue of A(0) and q(0) is the corresponding eigenvector, we obtain
(85)[iω0τkI-∫-10eiω0τksdη(s)]q(0)=0,[-iω0τkI-∫-10e-iω0τksdη(s)]q-(0)=0.So substituting (66) and (73) into (68), we obtain
(86)[2iω0τkI-∫-10e2iω0τksdη(s)]E1=(fv21fv22fv23),
that is, (87)(2iω0τk+(d+A)0BA(e-dτke-2iω0τk-1)2iω0τk+dB(e-dτke-2iω0τk-1)-Ae-dτke-2iω0τk02iω0τk+a+Be-dτke-2iω0τk)E1=(fv21fv22fv23).It follows that(88)E1(1)=1M1det(fv210Bfv222iω0τk+dB(e-dτke-2iω0τk-1)fv2302iω0τk+a+Be-dτke-2iω0τk),E1(2)=1M1det(2iω0τk+(d+A)fv21BA(e-dτke-2iω0τk-1)fv22B(e-dτke-2iω0τk-1)-Ae-dτke-2iω0τkfv232iω0τk+a+Be-dτke-2iω0τk),E1(3)=1M1det(2iω0τk+(d+A)0fv21A(e-dτke-2iω0τk-1)2iω0τk+dfv22-Ae-dτke-2iω0τk0fv23),where(89)M1=det(2iω0τk+(d+A)0BA(e-dτke-2iω0τk-1)2iω0τk+dB(e-dτke-2iω0τk-1)-Ae-dτke-2iω0τk02iω0τk+a+Be-dτke-2iω0τk).Similarly, substituting (67) and (74) into (69), we get(90)E2(1)=1M2det(fvv-10Bfvv-2dB(e-dτke-iω0τk-1)fvv-30a+Be-dτke-iω0τk),E2(2)=1M2det((d+A)fvv-1BA(e-dτke-iω0τk-1)fvv-2B(e-dτke-iω0τk-1)-Ae-dτke-iω0τkfvv-3a+Be-dτke-iω0τk),E2(3)=1M2det((d+A)0fvv-1A(e-dτke-iω0τk-1)dfvv-2-Ae-dτke-iω0τk0fvv-3),where(91)M2=det((d+A)0BA(e-dτke-iω0τk-1)dB(e-dτke-iω0τk-1)-Ae-dτke-iω0τk0a+Be-dτke-iω0τk).Thus, we can determineW20(s) and W11(s) from (66) and (67). Furthermore, we can compute f21 by (55). Thus, we can compute the following values:
(92)c(1)(0)=i2ω0τk(f20f11-2|f11|2-|f02|23)+f212,ϑ2=-Im{c1(0)}Re{λ′(τk)},σ2=2Re{c(1)(0)},T2=-Im{c1(0)}+ϑ2Re{λ′(τk)}ω0τk.Then from [8], we can give the properties of bifurcating periodic solution in the center manifold at the critical values τk. More specifically, we have the following theorem.Theorem 10.
The properties of the Hopf bifurcation are determined by the values in (75).(1)
ϑ2 determines the directions of the Hopf bifurcation: if ϑ2>0(ϑ2<0), then the Hopf bifurcation is supercritical (subcritical), and the bifurcating periodic solutions exist for τ>τk(τ<τk).(2)
ϑ2 determines the stability of the bifurcating periodic solutions: the bifurcating periodic solutions are stable (unstable) if ϑ2>0(ϑ2<0).(3)
T2 determines the period of the bifurcating periodic solutions: the period increases (decreases) if T2>0  (T2<0).


## 4. Numerical Simulations

As an example, we consider the system given in [1] with λ=1.5×105, d=0.011, β=9.8×10-7, k=150, a=0.1, and  μ=0.67, then, there is a chronic infected equilibrium E* which satisfies the conditions indicated in Theorem 9. When τ=0, the chronic infected equilibrium E*=(-1.3515×107,-1.0120×105,1.2090×105) is asymptotically stable (see Figure 1(a) which is the amplified part of Figure 1(b) during time interval [0,30]). It follows from the discussion in Section 2 that ω0=0.6029, τ0=8.1786 and λ(τ0)=0.2103-0.3942i. From the formulae (92) in Section 3, it follows that c1(0)=-0.1052-0.3061i, ϑ2=8.2963>0,  σ2=-0.2103<0, and T2=1.6737>0. Thus, the chronic infected equilibrium E* is asymptotically stable when 0≤τ<τ0 as illustrated by the computer simulations (see Figures 2 and 3).The response of the system (5) initial data (1.2e7,1.8e6,6e5)T, when τ=0.
(a)
(b)Figure 2
The response of the system (5) initial data (1.2e7,1.8e6,6e5)T, when τ=8.1786<τ0.Figure 3
3-phase portrait of the system (5) initial data (1.2e7,1.8e6,6e5)T, when τ=8.1786<τ0.Whenτ passes through the critical value τ0, E* loses its stability, and a Hopf bifurcation occurs, that is, a family of periodic solutions bifurcate from E*. Since ϑ2 and σ2=-0.2103<0, the Hopf bifurcation is supercritical, and the direction of the bifurcation is τ>τ0, and these bifurcating periodic solutions from E* at τ0 are stable, which are depicted in Figures 2 and 3.In Figure4, we try to reflect the changes of stability of chronic infected equilibrium E* as τ increases from 0 to 9. In Figure 4, each vertical blue strip corresponds to component of  x-(t),e(t), and y(t) at each τ, respectively. From Figure 3, we see that if τ∈(0,8.17864293), approximately, the vertical amplitudes of x-(t), e(t), and y(t) are as small as a point, suggesting that E* is asymptotically stable; if τ increases, the vertical amplitudes of x-(t), e(t), and y(t) will become larger and larger, showing that E* becomes more and more unstable. In particular, if τ=8.1787, points are well distributed around the positive equilibrium, and their amplitudes are all equivalent. This shows that periodic solution near positive equilibrium may occur.Figure 4
The response of the system (5) initial data (1.2e7,1.8e6,6e5)T, when τ=8.21>τ0.Figure 5
3-phase portrait of the system (5) initial data (1.2e7,1.8e6,6e5)T, when τ=8.21>τ0.


## 5. Conclusions

In this paper, a delay differential equation of hepatitis B virus infection is formulated. We analyzed the stability of the equilibria; a sufficient condition was given to guarantee the global stability of the origin. Local stability of the chronic infected equilibrium was considered. By choosing time delayτ as a bifurcated parameter, a sufficient condition has been presented for checking the existence of the Hopf bifurcation. The explicit formulae determining the direction, stability, and other properties of bifurcating periodic solutions were given by using the normal form theory and the center manifold theorem. Some numerical simulations were performed to support the analytical results found. Although bifurcations in a population dynamics with delay has been investigated by many researchers. However, to the best of our knowledge, there are few papers on the bifurcation of delay differential equation of hepatitis B virus infection dynamics with delay.



---

*Source: 875783-2013-03-06.xml*