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

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.

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: d𝑥(𝑡)(𝑡)−(𝑡)𝑣(𝑡)d𝑡=𝜆−𝑑𝑥𝛽𝑥𝑥(𝑡)(𝑡)(𝑡),(𝑡)+𝑒+𝑦d𝑒(𝑡)+(𝑡)𝑣(𝑡)d𝑡=−𝑑𝑒𝛽𝑥𝑥(𝑡)(𝑡)(𝑡)−+𝑒+𝑦𝛽𝑒^−𝑑𝜏𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)𝑥()()(),(𝑡)𝑡−𝜏+𝑒𝑡−𝜏+𝑦𝑡−𝜏d𝑦=d𝑡𝛽𝑒^−𝑑𝜏𝑥(𝑡−𝜏)𝑣(𝑡−𝜏)𝑥()()()(𝑡),(𝑡)𝑡−𝜏+𝑒𝑡−𝜏+𝑦𝑡−𝜏−𝑎𝑦d𝑣(𝑡)(𝑡),d𝑡=𝑘𝑦−𝜇𝑣 (1) where 𝑥(𝑡) and 𝑦(𝑡) represent the number of uninfected cells and infected cells, respectively. 𝑒(𝑡) represents the number of exposed cells, that is, the cells that have acquired the virus but are not yet producing new virions. 𝑣(𝑡) 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 𝑑 is the per capita death rate of uninfected live cells. Infected live cells are killed by immune cells at rate 𝑎𝑦 and produce free virions at rate 𝑘𝑦, where 𝑘 is what so-called “burst” constant. Free virions are cleared by lymphatic and other mechanisms at rate 𝜇𝑣, where 𝜇 is a constant. 𝛽>0 is an incidence rate coefficient describing the infection process. The initial conditions for the system (1) are 𝑥(𝑠)=𝑥₀(𝑠),𝑦(𝑠)=𝑦₀(𝑠),𝑣(𝑠)=𝑣₀(𝑠)[],𝐸,𝑠∈−𝜏,0₀=𝛽⁰_−𝜏𝑒^𝑑𝑠𝑣₀(𝑠)𝑥₀(𝑠)𝑥₀(𝑠)+𝑦₀(𝑠)+𝑒₀(s)d𝑠, (2) where 𝑥₀, 𝑦₀, and 𝑣₀ are nonnegative functions. Based on some observations of virus particles 𝑣, the system (1) is simplified in [1] as the following: d𝑥(𝑡)(𝑡)−(𝑡)𝑦(𝑡)d𝑡=𝜆−𝑑𝑥𝛽𝑘𝑥𝜇𝑥(𝑡)(𝑡,(𝑡)+𝑦)d𝑒(𝑡)+(𝑡)𝑦(𝑡)d𝑡=−𝑑𝑒𝛽𝑘𝑥𝜇𝑥(𝑡)(𝑡−+𝑦)𝛽𝑘𝑒^−𝑑𝜏𝑥(𝑡−𝜏)𝑦(𝑡−𝜏)𝜇𝑥()(,(𝑡)𝑡−𝜏+𝑦𝑡−𝜏)d𝑦=d𝑡𝛽𝑘𝑒^−𝑑𝜏𝑥(𝑡−𝜏)𝑦(𝑡−𝜏)𝜇𝑥()((𝑡).𝑡−𝜏+𝑦𝑡−𝜏)−𝑎𝑦 (3)

A direct computation shows that the basic infection reproduction number for the system (2) is 𝑅₀=𝛽𝑘𝑒^−𝑑𝜏.𝑎𝜇 (4)

For the sake of simplicity, let 𝑥=𝑥−(𝜆/𝑑), and the system (3) is equivalent to the following system: d𝑥(𝑡)d𝑡=−𝑑𝑥(𝑡)−𝛽𝑘[𝑥(𝑡)+(𝜆/𝑑)]𝑦(𝑡)𝜇𝑥(𝑡)+()(𝑡,(𝑡)𝜆/𝑑+𝑦)d𝑒(𝑡)+[d𝑡=−𝑑𝑒𝛽𝑘𝑥(𝑡)+(𝜆/𝑑)]𝑦(𝑡)𝜇𝑥(𝑡)+()(𝑡−𝜆/𝑑+𝑦)𝛽𝑘𝑒^−𝑑𝜏[𝑥(𝑡−𝜏)+(𝜆/𝑑)]𝑦(𝑡−𝜏)𝜇𝑥()+()(,(𝑡)𝑡−𝜏𝜆/𝑑+𝑦𝑡−𝜏)d𝑦=d𝑡𝛽𝑘𝑒^−𝑑𝜏[𝑥(𝑡−𝜏)+(𝜆/𝑑)]𝑦(𝑡−𝜏)𝜇𝑥()+()((𝑡),𝑡−𝜏𝜆/𝑑+𝑦𝑡−𝜏)−𝑎𝑦 (5) which has two equilibria: the infection-free equilibrium 𝐸_𝑓=(0,0,0)^𝑇 and the infected equilibrium 𝐸^∗=(𝑥^∗,𝑒^∗,𝑦^∗)^𝑇, where 𝑥^∗=𝜆𝑑+𝑎𝑒^𝑑𝜏𝑅₀−𝜆−1𝑑,𝑒^∗=𝑎𝑒^𝑑𝜏−1𝑑𝑥^∗,𝑦^∗=𝑅₀−1𝑥^∗. (6)

The following results, Theorems 1 and 2, come from [1].

Theorem 1. If 𝑅₀<1, the infection-free equilibrium 𝐸_𝑓 of the system (5) is globally asymptotically stable.

Theorem 2. If 𝑅₀>1, the chronic infected equilibrium 𝐸^∗ of the system (5) is locally asymptotically stable.

The initial conditions for the system (5) are 𝜑₁(𝜃),𝜑₂(𝜃),𝜓(𝜃)∈𝐶₊=𝐶[−𝜏,0],𝑅³₊₀,𝜑_𝑘(0)>0,𝜓(0)>0,𝑘=1,2, (7) where 𝑅³₊₀={(𝑥₁,𝑥₂,𝑥₃)∈𝑅³∶𝑥_𝑘>0,𝑘=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 𝑡≥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 Section 2, 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 𝐸_𝑓(0,0,0) for all nonnegative parameters. The characteristic equation of (5) at 𝐸_𝑓 is (𝑧+𝑑)²𝑧−𝛽𝑘𝑒^−𝑑𝜏𝜇+𝑎=0. (8)

Hence, 𝐸_𝑓 is a saddle with dim𝑊^𝑢(𝐸_𝑓)=1,dim𝑊^𝑠(𝐸_𝑓)=2 for 𝑅₀<1; 𝑊^𝑢(𝐸_𝑓) and 𝑊^𝑠(𝐸_𝑓) are the local unstable and stable manifolds of 𝐸_𝑓, respectively. 𝐸_𝑓 is locally asymptotically stable for 𝑅₀<1. In fact, 𝐸_𝑓 is globally asymptotically stable for 𝑅₀<1, see [1, 4].

Now, we will investigate the stability of the chronic infected equilibrium 𝐸^∗. Linearizing the system (5) at 𝐸^∗ yields the following linear system: d𝑥(𝑡)d𝑡=−𝑑𝑥(𝑡)−𝛽𝑘𝑦^∗(𝑡)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡−𝑥)𝛽𝑘^∗(𝑡)+(𝜆/𝑑)𝑦^∗(𝑡)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡)²𝑥(𝑡)−𝛽𝑘𝑥^∗(𝑡)+(𝜆/𝑑)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡−𝑥)𝛽𝑘^∗(𝑡)+(𝜆/𝑑)𝑦^∗(𝑡)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡)²𝑦(𝑡),d𝑒(𝑡)(𝑡)+d𝑡=−d𝑒𝛽𝑘𝑦^∗(𝑡)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡−𝑥)𝛽𝑘^∗(𝑡)+(𝜆/𝑑)𝑦^∗(𝑡)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡)²𝑥(𝑡)+𝛽𝑘𝑥^∗(𝑡)+(𝜆/𝑑)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡+𝑥)𝛽𝑘^∗(𝑡)+(𝜆/𝑑)𝑦^∗(𝑡)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡)²𝑦(𝑡)−𝛽𝑘𝑒^−𝑑𝜏𝑦^∗(𝑡)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡−)𝛽𝑘𝑒^−𝑑𝜏𝑥^∗(𝑡)+(𝜆/𝑑)𝑦^∗(𝑡)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡)²𝑥()−𝑡−𝜏𝛽𝑘𝑒^−𝑑𝜏𝑥^∗(𝑡)+(𝜆/𝑑)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡−)𝛽𝑘𝑒^−𝑑𝜏𝑥^∗(𝑡)+(𝜆/𝑑)𝑦^∗(𝑡)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡)²𝑦(),(𝑡)𝑡−𝜏d𝑦=d𝑡𝛽𝑘𝑒^−𝑑𝜏𝑦^∗(𝑡)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡−)𝛽𝑘𝑒^−𝑑𝜏𝑥^∗(𝑡)+(𝜆/𝑑)𝑦^∗(𝑡)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡)²𝑥()−𝑡−𝜏𝛽𝑘𝑒^−𝑑𝜏𝑥^∗(𝑡)+(𝜆/𝑑)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡−)𝛽𝑘e^−𝑑𝜏𝑥^∗(𝑡)+(𝜆/𝑑)𝑦^∗(𝑡)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡)²𝑦()(𝑡),𝑡−𝜏−𝑎𝑦 (9) whose characteristic equation reads as 𝑧³+(2𝑑+𝐴+𝑎)𝑧²+[𝑑(2𝑎+𝐴+𝑑)+𝑎𝐴]𝑧+𝑎𝑑(𝑑+𝐴)=𝐵𝑒^−𝑧𝜏𝑒^−𝑑𝜏𝑧²+2𝑑𝑧+𝑑², (10) with 𝐴=𝛽𝑘𝑦^∗(𝑡)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡−𝑥)𝛽𝑘^∗(𝑡)+(𝜆/𝑑)𝑦^∗(𝑡)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡)²,𝑥𝐵=𝛽𝑘^∗(𝑡)+(𝜆/𝑑)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡−𝑥)𝛽𝑘^∗(𝑡)+(𝜆/𝑑)𝑦^∗(𝑡)𝜇𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡)². (11)

For 𝜏=0, characteristic equation (10) reduces to 𝑧³+(2𝑑+𝐴+𝑎)𝑧²+[𝑑()]()𝑧2𝑎+𝐴+𝑑+𝑎𝐴𝑧+𝑎𝑑𝑑+𝐴=𝐵²+2𝑑𝑧+𝑑². (12)

By the Routh-Hurwitz criterion, we know that all the roots of (12) have negative real parts, that is, the chronic infected equilibrium 𝐸^∗ is locally asymptotically stable for 𝜏=0. We now give a definition, which can be found in [2, 9].

Definition 4. The equilibrium 𝐸^∗ is called asymptotically stable if there is an 𝜀>0 such that sup_{𝜏≤𝜃≤∞}𝜑₁(𝜃)−𝑥^∗₀+𝜑₂(𝜃)−𝑥^∗₁+𝜓(𝜃)−𝑦^∗<𝜀 (13) implies that lim_𝑡→∞𝑥₀(𝑡),𝑥₁(𝑡),𝑦(𝑡)=𝑥^∗₀,𝑥^∗₁,𝑦^∗, (14) where (𝑥₀(𝑡),𝑥₁(𝑡),𝑦(𝑡)) 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: 𝑧³+𝑚₂𝑧²+𝑚₁𝑧+𝑚₀+𝑒^−𝑧𝜏𝑛₂𝑧²+𝑛₁𝑧+𝑛₀=0, (15) where 𝑚_𝑗,𝑛_𝑗∈𝑅(𝑗=0,1,2) and ∑²_𝑗=0𝑛²_𝑗≠0. When 𝜏=0, (15) reduces to 𝑧³+𝑚₂+𝑛₂𝑧²+𝑚₁+𝑛₁𝑧+𝑛₀=0. (16)
Obviously, 𝑖𝜔(𝜔>0) is a root of (16) if and only if 𝜔 satisfies −𝜔³𝑖−𝑚₂𝜔²+𝑚₁𝜔𝑖+𝑚₀+−𝑛₂𝜔²+𝑛₁𝜔𝑖+𝑛₀(cos𝜔𝜏−𝑖sin𝜔𝜏)=0. (17)
Separating the real and the imaginary parts, we have 𝑚₂𝜔²−𝑚₀=−𝑛₂𝜔²cos𝜔𝜏+𝑛₁𝜔sin𝜔𝜏+𝑛₀cos𝜔𝜏,−𝜔³+𝑚₁𝜔=−𝑛₂𝜔²sin𝜔𝜏−𝑛₁𝜔cos𝜔𝜏+𝑛₀sin𝜔𝜏. (18)
Squaring the two equations and adding them give 𝜔⁶+𝑝𝜔⁴+𝑞𝜔²+𝑟=0, (19) where 𝑝=𝑚²₂−2𝑚₁−𝑛²₂, 𝑞=𝑚²₁−2𝑚₀𝑚₂+2𝑛₀𝑛₂−𝑛²₁, and 𝑟=𝑚²₀−𝑛²₀. 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)if 𝑟<0, then (19) has at least one positive root,(ii)if 𝑟≥0 and Δ=𝑝²−3𝑞≤0, then (19) has no positive roots,(iii)if 𝑟≥0 and Δ=𝑝²−3𝑞>0, then (19) has positive roots if and only if 𝑧^∗₁=(−𝑝+√Δ)/3 and ℎ(𝑧^∗₁)≤0, where ℎ(𝑧)=𝑧³+𝑝𝑧²+𝑞𝑧+𝑟. (20)

Lemma 6. For the transcendental equation (19),(i)if 𝑟≥0 and Δ=𝑝²−3𝑞≤0, then all roots with positive real parts of (19) have the same sum as those of the polynomial equation (16), for all 𝜏,(ii)if 𝑟<0 or 𝑟≥0, Δ=𝑝²−3𝑞>0,𝑧^∗₁=(−𝑝+√Δ)/3>0 and ℎ(z^∗₁)≤0, then all roots with positive real parts of (19) have the same sum as those of the polynomial equation (16), for 𝜏∈[0,𝜏₀).

From Lemmas 5 and 6, we can have the following lemma.

Lemma 7. (i) The chronic infected equilibrium 𝐸^∗ of the system (5) is absolutely stable if and only if the equilibrium 𝐸^∗ 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 equilibrium 𝐸^∗ 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 𝑖𝜔₀.
Then, one turns to an investigation of local stability of the chronic infected equilibrium 𝐸^∗ in the case of 𝑅₀<1.

Theorem 8. For 𝑅₀<1 holds, there exists a sequence of values for 𝜏: 𝜏₀<𝜏₁<⋯<𝜏_𝑘<⋯, (21) such that (10) has a pair of purely imaginary roots 𝑖𝜔₀ when 𝜏=𝜏_𝑘, 𝑘=0,1,2,…. That is, the chronic infected equilibrium 𝐸^∗ 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 𝑖𝜔₀(𝜔₀>0) for characteristic equation (10).
Assume that for some 𝜏>0, 𝑖𝜔(𝜔>0) is a root of (10), which implies that −𝑖𝜔³−(2𝑑+𝐴+𝑎)𝜔²+[𝑑()]()2𝑎+𝐴+𝑑+𝑎𝐴𝑖𝜔+𝑎𝑑𝑑+𝐴=𝐵𝑒^{−𝑖𝜔𝜏}𝑒^−𝑑𝜏−𝜔²+2𝑑𝑖𝜔+𝑑². (22)
Note that 𝐴>0, 𝐵>0 because of the positivities of the parameters 𝜆, 𝜇, 𝜅, 𝛽, 𝑑 and the properties 𝑥^∗(𝑡)+(𝜆/𝑑)≤[𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡)], 𝑦^∗(𝑡)≤[𝑥^∗(𝑡)+(𝜆/𝑑)+𝑦^∗(𝑡)].
Separating the real and imaginary parts and using Euler’s formula give −(2𝑑+𝐴+𝑎)𝜔²()𝑑+𝑎𝑑𝑑+𝐴=𝐵²−𝜔²𝑒^−𝑑𝜏cos𝜔𝜏+2𝐵𝑑𝜔𝑒^−𝑑𝜏sin𝜔𝜏,−𝜔³+[𝑑(2𝑎+𝐴+𝑑)+𝑎𝐴]𝜔=2𝐵𝑑𝜔𝑒^−𝑑𝜏cos𝜔𝜏−𝐵𝑑²−𝜔²𝑒^−𝑑𝜏sin𝜔𝜏, (23) which is equivalent to 𝜔⁶+𝑝𝜔⁴+𝑞𝜔²+𝑟=0, (24) where 𝑝=(2𝑑+𝐴+𝑎)²−2[𝑑(2𝑎+𝐴+𝑑)+𝑎𝐴]=𝐴²+2𝑑²+𝑎²[𝑑()]+2𝑑𝐴>0,𝑞=2𝑎+𝐴+𝑑+𝑎𝐴²−2(2𝑑+𝐴+𝑎)𝑎𝑑(𝑑+𝐴),𝑟=−𝑒^−2𝑑𝜏𝐵²𝑑²−𝜔²²+8𝐵²𝑑²<0. (25)
In order to solve (24), we first consider the following: 𝜔⁶+𝑝𝜔⁴+𝑞𝜔²+𝑟₀=0, (26) where 𝑟₀=−𝐵²𝑑²−𝜔²²+8𝐵²𝑑²<0. (27)
By Lemma 5, there is a unique positive 𝜔₀ satisfying (26). From (26), we get the corresponding 𝜏^′_𝑘>0 such that (26) has a pair of purely imaginary roots 𝜏^′_𝑘=1𝜔₀arccos𝑎𝑑(𝑑+𝐴)−(2𝑑+𝐴+𝑎)𝜔²₀−𝜔³₀+[𝑑(2𝑎+𝐴+𝑑)+𝑎𝐴]𝜔₀×𝐵²𝑑²−𝜔²²+8𝐵²𝑑²⁻¹+2𝑘𝜋𝜔₀(𝑘=0,1,2,…). (28)
Therefore, by using Rouché’s theorem [3], there is a unique positive 𝜏_𝑘=𝜏^′_𝑘+𝑜(1/𝑘) satisfying (26), that is, the characteristic equation (10) has a pair of purely imaginary roots of the form ±𝑖𝜔₀ as 𝑘→∞. By Lemma 7, we complete the proof of Theorem 8.
Next, we turn to show that d(Re𝜆)d𝜏_{𝜏=𝜏𝑘}>0. (29)
This will signify that there exists at least one eigenvalue with positive real part 𝜏=𝜏_𝑘. We first consider the following: 𝑧³+(2𝑑+𝐴+𝑎)𝑧²+[𝑑(2𝑎+𝐴+𝑑)+𝑎𝐴]𝑧+𝑎𝑑(𝑑+𝐴)=𝐵𝑒^−𝑧𝜏𝑧²+2𝑑𝑧+𝑑². (30)
Differentiating (30) with respect to 𝜏, we have d𝜆d𝜏⁻¹=3𝑧²+2(𝐴+2𝑑+𝑎)𝑧+[𝑑(2𝑎+𝐴+𝑑)+𝑎𝐴]−𝐵𝑧(𝑑+𝑧)²𝑒^−𝑧𝜏+𝜏𝑒^−𝑧𝜏𝐵(𝑧+𝑑)²−𝐵𝑧(𝑑+𝑧)²𝑒^−𝑧𝜏−2𝑒^−𝑧𝜏𝐵(𝑧+𝑑)−𝐵𝑧(𝑑+𝑧)²𝑒^−𝑧𝜏=(𝑑+𝑧)3𝑧²()𝑧+[𝑑()×()+2𝐴+2𝑑+𝑎2𝑎+𝐴+𝑑+𝑎𝐴]−𝐵𝑧𝑑+𝑧³𝑒^{−𝑧𝜏−1}−2𝑒^−𝑧𝜏𝐵(𝑧+𝑑)²−𝐵𝑧(𝑑+𝑧)³𝑒^−𝑧𝜏−𝜏𝑧=−𝑑3𝑧²()𝑧+[𝑑()×𝑧+2𝐴+2𝑑+𝑎2𝑎+𝐴+𝑑+𝑎𝐴](𝑑+𝑧)³+(2𝑑+𝐴+𝑎)𝑧²+[𝑑()]𝑧(2𝑎+𝐴+𝑑+𝑎𝐴+𝑎𝑑𝑑+𝐴)⁻¹−𝑧²−[𝑑()]()×(𝑧2𝑎+𝐴+𝑑+𝑎𝐴𝑧−2𝑎𝑑𝑑+𝐴𝑑+𝑧)³+(2𝑑+𝐴+𝑎)𝑧²+[𝑑()]𝑧(2𝑎+𝐴+𝑑+𝑎𝐴+𝑎𝑑𝑑+𝐴)⁻¹−𝜏𝑧. (31)
For the sake of simplicity, let 𝑎₁=𝐴+2𝑑+𝑎, 𝑎₂=𝑑(2𝑎+𝐴+𝑑)+𝑎𝐴, and 𝑎₃=𝑎𝑑(𝑑+𝐴), and (28) can be written as follows: d𝜆d𝜏⁻¹𝑑=−3𝑧²+2𝑎₁𝑧+𝑎₂(𝑑+𝑧)𝑧³+𝑎₁𝑧²+𝑎₂𝑧+𝑎₃−𝑧²−𝑎₂𝑧−2𝑎₃(𝑑+𝑧)𝑧³+𝑎₁𝑧²+𝑎₂𝑧+𝑎₃−𝜏𝑧=−(3𝑑+1)𝑧²−2𝑎₁𝑑−𝑎₂𝑧−𝑎₂𝑑−2𝑎₃(𝑑+𝑧)𝑧³+𝑎₁𝑧²+𝑎₂𝑧+𝑎₃−𝜏𝑧. (32)
Therefore, d()signRe𝜆d𝜏_{𝜏=𝜏𝑘}=signRed𝜆d𝜏⁻¹_{𝜏=𝜏𝑘}−()𝑧=signRe3𝑑+1²−2𝑎₁𝑑−𝑎₂𝑧−𝑎₂𝑑−2𝑎₃(𝑑+𝑧)𝑧³+𝑎₁𝑧²+𝑎₂𝑧+𝑎₃−𝜏𝑧=signRe(3𝑑+1)𝜔²₀−𝑖2𝑎₁𝑑−𝑎₂𝜔₀−𝑎₂𝑑−2𝑎₃𝑑−𝑖𝜔₀×𝑑²+𝜔²₀×𝑎²₃−𝑎₁𝜔²₀²+𝑎₂𝜔₀−𝜔³₀²⁻¹×𝑎²₃−𝑎₁𝜔²₀−𝑖𝑎₂𝜔₀−𝜔³₀()𝜔=sign3𝑑+1⁶₀+𝑝²𝜔⁴₀+𝐵𝑑²𝜔²₀𝑑²+𝜔²₀𝑎²₃−𝑎₁𝜔²₀²+𝑎₂𝜔₀−𝜔³₀²>0. (33)
This root of characteristic equation (9) crosses the imaginary axis from the left to the right as 𝜏 continuously varies from a number less than 𝜏_𝑘 to one greater than 𝜏_𝑘 again by Rouché’s theorem [3]. Therefore, the transversality condition holds, and the conditions for Hopf bifurcation [11] are then satisfied at 𝜏=𝜏_𝑘. In conclusion, we have the following stability and bifurcation results to (5).

Theorem 9. Suppose that 𝑅₀<1 holds. Then, for each fixed 𝜏>0, there exists a sequence of values for 𝜏: 𝜏₀<𝜏₁<⋯<𝜏_𝑘<⋯, (34) such that the positive equilibrium 𝐸^∗ is asymptotically stable when 𝜏∈[0,𝜏₀), and unstable when 𝜏>𝜏₀. Furthermore, (5) undergoes a Hopf bifurcation at 𝐸^∗ when 𝜏=𝜏_𝑘=𝜏^′_𝑘+𝑜(1/𝑘),𝑘=0,1,2,…, where 𝜏^′_𝑘 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 𝑥₁=𝑥−𝑥^∗, 𝑥₂=𝑒−𝑒^∗, 𝑥₃=𝑦−𝑦^∗, 𝑥^∼_𝑖(𝑡)=𝑥_𝑖(𝜏𝑡), 𝜏=𝜏_𝑘+𝜗, and 𝜗∈𝑅=(−∞,∞). 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 𝐶([−1,0],𝑅³) as 𝑥̇(𝑡)=𝐿_𝜗𝑥_𝑡+𝑓𝜗,𝑥_𝑡, (35) where 𝐶([−1,0],𝑅³) is the Banach space of continuous functions mapping the interval [−𝜏,0] into 𝑅³, 𝑥(𝑡)=(𝑥₁(𝑡),𝑥₂(𝑡),𝑥₃(𝑡))^𝑇∈𝑅³, 𝑥_𝑡(𝑠)=𝑥(𝑡+𝑠) for 𝑠∈[−𝜏,0] and 𝐿_𝜗∶𝐶→𝑅³, 𝑓∶𝑅×𝐶→𝑅³ are read, respectively, as 𝐿_𝜗𝜙=𝜏_𝑘−()𝜙+𝜗𝑑+𝐴0−𝐵𝐴−𝑑𝐵00−𝑎¹(0)𝜙₂(0)𝜙₃(0)+𝜏_𝑘+𝜗000−𝐴𝑒^{−𝑑𝜏𝑘}0𝐵𝑒^{−𝑑𝜏𝑘}𝐴𝑒^{−𝑑𝜏𝑘}0−𝐵𝑒^{−𝑑𝜏𝑘}𝜙₁()𝜙−1₂()𝜙−1₃(),𝑓−1𝜗,𝑥_𝑡=𝜏_𝑘×−𝜙+𝜗𝛽𝑘₁(0)+(𝜆/𝑑)𝜙₃(0)𝜇𝜙₁(0)+(𝜆/𝑑)+𝜙₃(0𝜙)𝛽𝑘₁(0)+(𝜆/𝑑)𝜙₃(0)𝜇𝜙₁(0)+(𝜆/𝑑)+𝜙₃(0−)𝛽𝑘𝑒^{−𝑑𝜏𝑘}𝜙₁(−1)+(𝜆/𝑑)𝜙₃(−1)𝜇𝜙₁(−1)+(𝜆/𝑑)+𝜙₃(−1)𝛽𝑘𝑒^{−𝑑𝜏𝑘}𝜙₁(−1)+(𝜆/𝑑)𝜙₃(−1)𝜇𝜙₁(−1)+(𝜆/𝑑)+𝜙₃(.−1) (36)(37)Obviously, 𝐿_𝜗 is a continuous linear function mapping 𝐶 into 𝑅³, by the Riesz representation theorem, there exists a 3×3 matrix function 𝜂(𝑠,𝜇) of bounded variation for 𝑠∈[−1,0], such that 𝐿_𝜗𝜙=⁰₋₁d𝜂𝑠,𝜇𝜙(𝑠). (38)

In fact, we can choose 𝜂(𝑠,𝜗)=𝜏_𝑘−()𝛿(𝑠)+𝜏+𝜗𝑑+𝐴0−𝐵𝐴−𝑑𝐵00−𝑎_𝑘+𝜗000−𝐴𝑒^{−𝑑𝜏𝑘}0𝐵𝑒^{−𝑑𝜏𝑘}𝐴𝑒^{−𝑑𝜏𝑘}0−𝐵𝑒^{−𝑑𝜏𝑘}𝛿(𝑠+1), (39) where 𝛿(⋅) denotes the Dirac delta function.

If 𝜙 is any given function in 𝐶 and 𝑥(𝜙) is the unique solution of the linearized equation 𝑥̇(𝑡)=𝐿_𝜗(𝑥_𝑡) of (35) with initial function 𝜙 at zero, then the solution operator 𝑇(𝑡)∶𝐶→𝐶 is defined by 𝑇(𝑡)𝜙=𝑥_𝑡(𝜙), for all 𝑡≥0. It is obvious that 𝑇(𝑡), 𝑡≥0, is a strongly continuous semigroup of linear transformation on [0,+∞) and the infinitesimal generator 𝐴_𝜗 of 𝑇(𝑡),𝑡≥0, is 𝐴_𝜗=d𝜙(𝑠)[)∫d𝑠,𝑠∈−1,0₀₋₁d𝜂𝑠,𝜇𝜙(𝑠),𝑠=0. (40) For 𝜙∈𝐶¹([−1,0],𝑅³), the space of functions mapping the interval [−1,0] into 𝑅³ which have a continuous first derivative also defines 𝑅_𝜗𝜙=[)𝑓0,𝑠∈−1,0𝜗,𝜙,𝑠=0. (41)

Then, system (35) is equivalent to 𝑥̇(𝑡)=𝐴_𝜗𝑥_𝑡+𝑅_𝜗𝑥_𝑡. (42)

For 𝜓∈𝐶¹([−1,0],(𝑅³)^∗), the space of functions mapping interval [0,1] into the three-dimensional row vectors which have continuous first derivative defines 𝐴^∗_𝜗=−d𝜓(𝑠)(]∫d𝑠,𝑠∈0,1₀₋₁d𝜂^𝑇(𝑡,0)𝜙(−𝑡),𝑠=0, (43) and a bilinear inner product 𝜓(𝜁),𝜙(𝑠)=𝜓(0)𝜙(0)−⁰₋₁^𝑠_𝜉=0𝜓(𝜉−𝑠)d𝜂(𝑠)𝜙(𝜉)d𝜉, (44) where 𝜂(𝑠)=𝜂(𝑠,0). Then, 𝐴(0) and 𝐴^∗ are adjoint operators. By the discussion in Section 2, we know that 𝑧 are eigenvalues of 𝐴(0). Hence, they are also eigenvalues of 𝐴^∗. We first need to compute the eigenvectors of 𝐴(0) and 𝐴^∗ corresponding to 𝑖𝜔₀𝜏_𝑘 and −𝑖𝜔₀𝜏_𝑘, respectively. Suppose that 𝑞(𝑠)=(1,𝑎,𝑏)^𝑇𝑒^{𝑖𝜔0𝜏𝑘𝑠} is the eigenvectors of 𝐴(0) corresponding to 𝑖𝜔₀𝜏_𝑘, then 𝐴(0)𝑞(𝑠)=𝑖𝜔₀𝜏_𝑘𝑞(𝑠). Then, from the definition of 𝐴(0) and (36), (38), and (39), we have 𝜏_𝑘−()𝑞(0)𝑑+𝐴0−𝐵𝐴−𝑑𝐵00−𝑎+𝜏_𝑘000−𝐴𝑒^{−𝑑𝜏𝑘}0𝐵𝑒^{−𝑑𝜏𝑘}𝐴𝑒^{−𝑑𝜏𝑘}0−𝐵𝑒^{−𝑑𝜏𝑘}𝑞(−1)=𝑖𝜔₀𝜏_𝑘𝑞(0). (45)

For 𝑞(−1)=𝑞(0)𝑒^{−𝑖𝜔0𝜏𝑘}, then we obtain 1𝜌=−1+𝐵+𝐴𝑑−𝑖𝜔₀𝑎+𝑖𝜔₀𝐵d²+𝜔²₀,𝜎=−𝐴+𝑑+𝑖𝜔₀𝐵. (46)

Similarly, we can obtain the eigenvector 𝑞^∗(𝑠)=𝐷(1,𝜌^∗,𝜎^∗)𝑒^{𝑖𝜔0𝜏𝑘} of 𝐴^∗ corresponding to −𝑖𝜔₀𝜏_𝑘, where 𝜌^∗1=𝜌=−1+𝐵+𝐴𝑑−𝑖𝜔₀𝑎+𝑖𝜔₀𝐵𝑑²+𝜔²₀,𝜎^∗=𝜎=−𝐴+𝑑+𝑖𝜔₀𝐵. (47)

We need to determine the value of 𝐷 to ensure that ⟨𝑞^∗(𝜍),𝑞(𝑠)⟩=1. By (44), we have 𝑞^∗(𝜍)(𝑠=,𝑞)𝐷1,𝜌^∗,𝜎^∗1,𝜌,𝜎^𝑇−⁰₋₁^𝑠_𝜉=0𝐷1,𝜌^∗,𝜎^∗𝑒^{−𝑖𝜔0𝜏𝑘(𝜉−𝑠)}d𝜂(𝑠)1,𝜌,𝜎^𝑇×𝑒^{𝑖𝜔0𝜏𝑘𝜉}=d𝜉𝐷1+𝜌²+|𝜎|²−⁰₋₁1,𝜌^∗,𝜎^∗𝑠𝑒^{𝑖𝜔0𝜏𝑘𝑠}d𝜂(𝑠)1,𝜌,𝜎^𝑇=𝐷1+𝜌²+|𝜎|²+𝑒^{−𝑑𝜏𝑘}𝑒^{−𝑖𝜔0𝜏𝑘}(𝐴−𝐵𝜎)𝜌^∗−𝜎^∗. (48)

Therefore, we can choose 𝐷 as 1𝐷=1+𝜌²+|𝜎|²+𝑒^{−𝑑𝜏𝑘}𝑒^{−𝑖𝜔0𝜏𝑘}(𝐴−𝐵𝜎)𝜌^∗−𝜎^∗. (49)

Next, we will compute the coordinate to describe the center manifold 𝐶₀ at 𝜗=0. Let 𝑥_𝑡 be the solution of (42) when 𝜗=0. Define𝜐(𝑡)=𝑞^∗(𝜍),𝑥_𝑡,𝑊(𝑡,𝑠)=𝑥_𝑡(𝑠)−2Re𝜐(𝑡)𝑞(𝑠). (50)

On the center manifold 𝐶₀, we have 𝑊(𝑡,𝑠)=𝑊(𝜐(𝑡),𝜐(𝑡),𝑠), (51) where 𝑊(𝜐(𝑡),𝜐(𝑡),𝑠)=𝑊₂₀(𝑠)𝜐²2+𝑊₁₀(𝑠)𝜐𝜐+𝑊₀₂(𝑠)𝜐²2+⋯. (52)𝜐 and 𝜐 are local coordinates for the center manifold 𝐶₀ in the direction of 𝑞^∗ and 𝑞^∗. Note that 𝑊 is real if 𝑥_𝑡 is real. We only consider real solutions. For solution 𝑥_𝑡∈𝐶₀ of (42), since 𝜗=0, we have 𝜐̇(𝑡)=𝑖𝜔₀𝜏_𝑘𝜐(𝑡)+𝑞^∗(0)𝑓0,𝑊(𝜐,𝜐,0)+2Re𝜐𝑞(𝑠)∶=𝑖𝜔₀𝜏_𝑘𝜐(𝑡)+𝑞^∗(0)𝑓₀(𝜐,𝜐). (53)

For convenience of calculation, let 𝑞^∗(0)𝑓₀(𝜐,𝜐)=𝑓₂₀𝜐²2+𝑓₁₁(𝑠)𝜐𝜐+𝑓₀₂(𝑠)𝜐²2+𝑓₂₁𝜐²𝜐2+⋯. (54)

On the other hand, it follows from (50) and (52) that 𝑥_𝑡(𝑠)=𝑊(𝜐,𝜐,0)+2Re𝜐𝑞(𝑠)=𝑊₂₀(𝑠)𝜐²2+𝑊₁₀(𝑠)𝜐𝜐+𝑊₀₂(𝑠)𝜐²2+(1,𝑎,𝑏)^𝑇𝑒^{𝑖𝜔0𝜏𝑘𝑠}𝜐+1,𝑎,𝑏^𝑇𝑒^{−𝑖𝜔0𝜏𝑘𝑠}𝜐+⋯. (55)

Notice that 𝜙₁(0)=𝑞₁(0)𝑣+𝑞₁(0)𝑣+𝑊⁽¹⁾₂₀(0)𝑣²2+𝑊⁽¹⁾₁₁(0)𝑣𝑣+𝑊⁽¹⁾₀₂(0)𝑣²2𝜙+⋯,₂(0)=𝑞₂(0)𝑣+𝑞₂(0)𝑣+𝑊⁽²⁾₂₀(0)𝑣²2+𝑊⁽²⁾₁₁(0)𝑣𝑣+𝑊⁽²⁾₀₂(0)𝑣²2𝜙+⋯,₃(0)=𝑞₃(0)𝑣+𝑞₃(0)𝑣+𝑊⁽³⁾₂₀(0)𝑣²2+𝑊⁽³⁾₁₁(0)𝑣𝑣+𝑊⁽³⁾₀₂(0)𝑣²2𝜙+⋯,₁(−1)=𝑞₁(0)𝑣𝑒^{−𝑖𝜔∗0𝜏∗0}+𝑞₁(0)𝑣𝑒^{𝑖𝜔∗0𝜏∗0}+𝑊⁽¹⁾₂₀()𝑣−1²2+𝑊⁽¹⁾₁₁(−1)𝑣𝜙𝑣+⋯,₂(−1)=𝑞₂(0)𝑣𝑒^{−𝑖𝜔∗0𝜏∗0}+𝑞₂(0)𝑣𝑒^{𝑖𝜔∗0𝜏∗0}+𝑊⁽²⁾₂₀()𝑣−1²2+𝑊⁽²⁾₁₁(−1)𝑣𝜙𝑣+⋯,₃(−1)=𝑞₃(0)𝑣𝑒^{−𝑖𝜔∗0𝜏∗0}+𝑞₃(0)𝑣𝑒^{𝑖𝜔∗0𝜏∗0}+𝑊⁽³⁾₂₀()𝑣−1²2+𝑊⁽³⁾₁₁(−1)𝑣𝑣+⋯, (56)(57)(58)(59)(60)(61)

Taking (37) into account, we can obtain that𝑞^∗(0)𝑓₀(𝜐,𝜐)=𝑞^∗(0)𝑓0,𝑥_𝑡=𝜏_𝑘𝐷1,𝜌^∗,𝜎^∗−𝛽𝑘𝜙₁(0)+(𝜆/𝑑)𝜙₃(0)𝜇𝜙₁(0)+(𝜆/𝑑)+𝜙₃(0𝜙)𝛽𝑘₁(0)+(𝜆/𝑑)𝜙₃(0)𝜇𝜙₁(0)+(𝜆/𝑑)+𝜙₃(0−)𝛽𝑘𝑒^{−𝑑𝜏𝑘}𝜙₁(−1)+(𝜆/𝑑)𝜙₃(−1)𝜇𝜙₁(−1)+(𝜆/𝑑)+𝜙₃(−1)𝛽𝑘𝑒^{−𝑑𝜏𝑘}𝜙₁(−1)+(𝜆/𝑑)𝜙₃(−1)𝜇𝜙₁(−1)+(𝜆/𝑑)+𝜙₃(−1)=𝜏_𝑘𝐷−𝛽𝑘𝜙₁(0)+(𝜆/𝑑)𝜙₃(0)𝜇𝜙₁(0)+(𝜆/𝑑)+𝜙₃(0+)𝜌^∗𝛽𝑘𝜙₁(0)+(𝜆/𝑑)𝜙₃(0)𝜇𝜙₁(0)+(𝜆/𝑑)+𝜙₃(0−)𝜌^∗𝛽𝑘𝑒^{−𝑑𝜏𝑘}𝜙₁(−1)+(𝜆/𝑑)𝜙₃(−1)𝜇𝜙₁(−1)+(𝜆/𝑑)+𝜙₃(+−1)𝜎^∗𝛽𝑘𝑒^{−𝑑𝜏𝑘}𝜙₁(−1)+(𝜆/𝑑)𝜙₃(−1)𝜇𝜙₁(−1)+(𝜆/𝑑)+𝜙₃(=𝜌−1)_∗𝜏−1^𝑘𝐷𝛽𝑘𝜇Δ₀𝜆𝑑+𝑞₁(0)𝑣+𝑞₁(0)𝑣+𝑊⁽¹⁾₂₀(0)𝑣²2+𝑊⁽¹⁾₁₁(0)𝑣𝑣+𝑊⁽¹⁾₀₂(0)𝑣²2×𝑞+⋯₃(0)𝑣+𝑞₃(0)𝑣+𝑊⁽³⁾₂₀(0)𝑣²2+𝑊⁽³⁾₁₁(0)𝑣𝑣+𝑊⁽³⁾₀₂(0)𝑣²2×𝜆+⋯𝑑+𝑞₁(0)𝑣+𝑞₁(0)𝑣+𝑊⁽¹⁾₂₀(0)𝑣²2+𝑊⁽¹⁾₁₁(0)𝑣𝑣+𝑊⁽¹⁾₀₂(0)𝑣²2+⋯+𝑞₃(0)𝑣+𝑞₃(0)𝑣+𝑊⁽³⁾₂₀(0)𝑣²2+𝑊⁽³⁾₁₁(0)𝑣𝑣+𝑊⁽³⁾₀₂(0)𝑣²2+𝜎+⋯_∗−𝜌^∗𝜏_𝑘𝐷𝛽𝑘𝑒^{−𝑑𝜏𝑘}𝜇Δ₋₁×𝜆𝑑+𝑞₁(0)𝑣𝑒^{−𝑖𝜔∗0𝜏∗0}+𝑞₁(0)𝑣𝑒^{𝑖𝜔∗0𝜏∗0}+𝑊⁽¹⁾₂₀(0)𝑣²2+𝑊⁽¹⁾₁₁(0)𝑣×𝑞𝑣+⋯₃(0)𝑣𝑒^{−𝑖𝜔∗0𝜏∗0}+𝑞₃(0)𝑣𝑒^{𝑖𝜔∗0𝜏∗0}+𝑊⁽³⁾₂₀(0)𝑣²2+𝑊⁽³⁾₁₁(0)𝑣𝑣+𝑊⁽³⁾₀₂(0)𝑣²2×𝜆+⋯𝑑+𝑞₁(0)𝑒^{−𝑖𝜔∗0𝜏∗0}+𝑞₁(0)𝑣𝑒^{𝑖𝜔∗0𝜏∗0}+𝑊⁽¹⁾₂₀(0)𝑣²2+𝑊⁽¹⁾₁₁(0)𝑣𝑣+𝑊⁽¹⁾₀₂(0)𝑣²2+⋯+𝑞₃(0)𝑣𝑒^{−𝑖𝜔∗0𝜏∗0}+𝑞₃(0)𝑣𝑒^{𝑖𝜔∗0𝜏∗0}+𝑊⁽³⁾₂₀(0)𝑣²2+𝑊⁽³⁾₁₁(0)𝑣𝑣+𝑊⁽³⁾₀₂(0)𝑣²2+⋯, (62)

where Δ₀=𝜆𝑑+𝑞₁(0)𝑣+𝑞₁(0)𝑣+𝑊⁽¹⁾₂₀(0)𝑣²2+𝑊⁽¹⁾₁₁(0)𝑣𝑣+𝑊⁽¹⁾₀₂(0)𝑣²2+⋯+𝑞₃(0)𝑣+𝑞₃(0)𝑣+𝑊⁽³⁾₂₀(0)𝑣²2+𝑊⁽³⁾₁₁(0)𝑣𝑣+𝑊⁽³⁾₀₂(0)𝑣²2+⋯²,Δ₋₁=𝜆𝑑+𝑞₁(0)𝑣𝑒^{−𝑖𝜔∗0𝜏∗0}+𝑞₃(0)𝑣𝑒^{𝑖𝜔∗0𝜏∗0}+𝑊⁽¹⁾₂₀(0)𝑣²2+𝑊⁽¹⁾₁₁(0)𝑣𝑣+𝑊⁽¹⁾₀₂(0)𝑣²2+⋯+𝑞₃(0)𝑣+𝑞₃(0)𝑣+𝑊⁽³⁾₂₀(0)𝑣²2+𝑊⁽³⁾₁₁(0)𝑣𝑣+𝑊⁽³⁾₀₂(0)𝑣²2+⋯². (63)

Straightforward calculation leads to 𝑞^∗(0)𝑓₀(𝜐,𝜐)=𝜌^∗−1𝛽𝑘𝜇Δ₀×𝜆𝑑𝑞₁(0)+𝑞₃(0)𝑞₃(0)𝑣²+𝜆𝑑𝑞₁(0)𝑞₃(0)𝑣²+𝜆𝑑²𝑊⁽³⁾₂₀(0)𝑣²2+𝜆𝑑𝑞₁(0)+𝑞₃(0)𝑞₃(0)𝑣𝑣𝜆+2Re𝑑𝑞₁(0)𝑞₃(0)𝑣𝑣+𝜆𝑑²𝑊⁽³⁾₁₁(0)𝑣+𝜎𝑣+⋯_∗−𝜌^∗𝜏_𝑘𝐷𝛽𝑘𝑒^{−𝑑𝜏𝑘}𝜇Δ₋₁×𝜆𝑑𝑞₁(0)+𝑞₃(0)𝑞₃(0)𝑣²+𝜆𝑑𝑞₁(0)𝑞₃(0)𝑣²+𝜆𝑑²𝑊⁽³⁾₂₀()𝑣−1²2+𝜆𝑑𝑞₁(0)+𝑞₃(0)𝑞₃(0)𝑣𝑣𝜆+2Re𝑑𝑞₁(0)𝑞₃(0)𝑣𝑣+𝜆𝑑²𝑊⁽³⁾₁₁(−1)𝑣𝑣+⋯. (64)

Comparing the coefficients with (54), we have 𝑓₂₀=𝜌^∗−1𝛽𝑘𝜇Δ₀×𝜆𝑑𝑞₁(0)+𝑞₃(0)𝑞₃(0)+𝜆𝑑𝑞₁(0)𝑞₃(0)+12𝜆𝑑²𝑊⁽³⁾₂₀(0+𝜎_∗)−𝜌^∗𝜏_𝑘𝐷𝛽𝑘𝑒^{−𝑑𝜏𝑘}𝜇Δ₋₁×𝜆𝑑𝑞₁(0)+𝑞₃(0)𝑞₃(0)+𝜆𝑑𝑞₁(0)𝑞₃(0)+12𝜆𝑑²𝑊⁽³⁾₂₀(,𝑓−1)₁₁=𝜆𝑑𝑞₁(0)+𝑞₃(0)𝑞₃(0)𝜆+2Re𝑑𝑞₁(0)𝑞₃(0+𝜆)𝑑²𝑊⁽³⁾₁₁()+𝜆−1𝑑𝑞₁(0)+𝑞₃(0)𝑞₃(0)𝜆+2Re𝑑𝑞₁(0)𝑞₃(0+𝜆)𝑑²𝑊⁽³⁾₁₁(),𝑓−1₀₂=𝜌^∗−1𝛽𝑘𝜇Δ₀×𝜆𝑑𝑞₁(0)+𝑞₃(0)𝑞₃(0)+𝜆𝑑𝑞₁(0)𝑞₃(0)+𝜆𝑑²𝑊⁽³⁾₂₀(0)𝑣²2+𝜎^∗−𝜌^∗𝜏_𝑘𝐷𝛽𝑘𝑒^{−𝑑𝜏𝑘}𝜇Δ₋₁×𝜆𝑑𝑞₁(0)+𝑞₃(0)𝑞₃(0)+𝜆𝑑𝑞₁(0)𝑞₃(0)+12𝜆𝑑²𝑊⁽³⁾₂₀(,𝑓−1)₂₁=𝜌^∗−1𝛽𝑘𝜇Δ₀×𝜆𝑑𝑞₁(0)+𝑞₃(0)𝑞₃(0)+𝜆𝑑𝑞₁(0)𝑞₃(0)+𝜆𝑑²𝑊⁽³⁾₂₀(0)𝑣²2+𝜎^∗−𝜌^∗𝜏_𝑘𝐷𝛽𝑘𝑒^{−𝑑𝜏𝑘}𝜇Δ₋₁×𝜆𝑑𝑞₁(0)+𝑞₃(0)𝑞₃(0)+𝜆𝑑𝑞₁(0)𝑞₃(0)+12𝜆𝑑²𝑊⁽³⁾₂₀(−1). (65)

Since 𝑊₂₀(𝑠) and 𝑊₁₁(𝑠) are in 𝑓₂₁, we still need to compute them. From (40) and (44), we have 𝑊̇=𝑥̇_𝑡−𝜐̇𝑞−𝜐̇𝑞=𝐴𝑊−2Re𝑞^∗(0)𝑓₀𝑞(𝑠),𝑠∈[−1,0),𝐴𝑊−2Re𝑞^∗(0)𝑓₀𝑞(𝑠)+𝑓₀,𝑠=0, (66) that is, 𝑊̇=𝐴𝑊+𝐻(𝜐,𝜐,𝑠), (67) where 𝐻(𝜐,𝜐,𝑠)=𝐻₂₀(𝑠)𝜐²2+𝐻₁₀(𝑠)𝜐𝜐+𝐻₀₂(𝑠)𝜐²2+⋯. (68)

Substituting the corresponding series into (67) and comparing the coefficients, we have 𝐴−2𝑖𝜔₀𝜏_𝑘𝑊₂₀(𝑠)=−𝐻₂₀(𝑠),𝐴𝑊₁₁(𝑠)=−𝐻₁₁(𝑠),…. (69)

From (67), we know that for 𝑠∈[−1,0), 𝐻(𝜐,𝜐,𝑠)=−𝑞^∗(𝑠)𝑓₀𝑞(𝑠)−𝑞^∗(s)𝑓₀𝑞(𝑠)=−𝑓(𝜐,𝜐)𝑞(𝑠)−𝑓(𝜐,𝜐)𝑞(𝑠). (70)

Comparing the coefficients with (68) gives 𝐻₂₀(𝑠)=−𝑓₂₀𝑞(𝑠)−𝑓₀₂𝑞(𝑠),𝐻₁₁(𝑠)=−𝑓₁₁𝑞(𝑠)−𝑓₁₁𝑞(𝑠). (71)

From the definition of 𝐴 and (69) and (71), we obtain 𝑊̇₂₀(s)=2𝑖𝜔₀𝜏_𝑘𝑊₂₀(𝑠)+𝑓₂₀𝑞(𝑠)+𝑓₀₂𝑞(𝑠). (72)

For 𝑞(𝑠)=(1,𝜌,𝜎)^𝑇𝑒^{𝑖𝜔0𝜏𝑘𝑠}, hence 𝑊₂₀(𝑠)=𝑖𝑓₂₀𝜔₀𝜏_𝑘𝑞(0)𝑒^{𝑖𝜔0𝜏𝑘𝑠}+𝑖𝑓₀₂3𝜔₀𝜏_𝑘𝑞(0)𝑒^{−𝑖𝜔0𝜏𝑘𝑠}𝑞(𝑠)+𝐸₁𝑒^{2𝑖𝜔0𝜏𝑘𝑠}, (73) where 𝐸₁=(𝐸₁⁽¹⁾,𝐸₁⁽²⁾,𝐸₁⁽³⁾) is a constant vector. Similarly, we know that 𝑊₁₁(𝑠)=−𝑖𝑓₁₁𝜔₀𝜏_𝑘𝑞(0)𝑒^{𝑖𝜔0𝜏𝑘𝑠}+𝑖𝑓₁₁𝜔₀𝜏_𝑘𝑞(0)𝑒^{−𝑖𝜔0𝜏𝑘𝑠}𝑞(𝑠)+𝐸₂, (74) where 𝐸₂=(𝐸₂⁽¹⁾,𝐸₂⁽²⁾,𝐸₂⁽³⁾) is a constant vector.

In what follows, we will calculate 𝐸₁ and 𝐸₂. From the definition of 𝐴 and (69), we have ⁰₋₁d𝜂^𝑇(𝑠)𝑊₂₀(𝑠)=2𝑖𝜔₀𝜏_𝑘𝑊₂₀(0)−𝐻₂₀(0),₀₋₁d𝜂^𝑇(𝑠)𝑊₁₁(𝑠)=−𝐻₁₁(0), (75)(76) where 𝜂(𝑠)=𝜂(𝑠,0). When 𝑠=0, we obtain that 𝐻(𝜐,𝜐,0)=−2Re𝑞^∗(0)𝑓₀𝑞(0)+𝑓₀=−𝑞^∗(0)𝑓₀𝑞(0)−𝑞^∗(0)𝑓₀𝑞(0)+𝑓₀, (77) that is, 𝐻₂₀(𝑠)𝜐²2+𝐻₁₁(𝑠)𝜐𝜐+𝐻₀₂(𝑠)𝜐²2(0)𝑓+⋯=−𝑞₂₀𝜐²2+𝑓₁₁(𝑠)𝜐𝜐+𝑓₀₂(𝑠)𝜐²2+𝑓₂₁𝜐²𝜐2⋯−𝑞(0)𝑓₂₀𝜐²2+𝑓₁₁(𝑠)𝜐𝜐+𝑓₀₂(𝑠)𝜐²2+𝑓₂₁𝜐²𝜐2⋯+𝑓₀, (78) with 𝑓₀=𝜏_𝑘−𝛽𝑘𝜙₁(0)+(𝜆/𝑑)𝜙₃(0)𝜇𝜙₁(0)+(𝜆/𝑑)+𝜙₃(0𝜙)𝛽𝑘₁(0)+(𝜆/𝑑)𝜙₃(0)𝜇𝜙₁(0)+(𝜆/𝑑)+𝜙₃(0−)𝛽𝑘𝑒^{−𝑑𝜏𝑘}𝜙₁(−1)+(𝜆/𝑑)𝜙₃(−1)𝜇𝜙₁(−1)+(𝜆/𝑑)+𝜙₃(−1)𝛽𝑘𝑒^{−𝑑𝜏𝑘}𝜙₁(−1)+(𝜆/𝑑)𝜙₃(−1)𝜇𝜙₁(−1)+(𝜆/𝑑)+𝜙₃(−1). (79)

From (50), we have 𝑓₀=𝑓¹_𝑣2𝑓²_𝑣2𝑓³_𝑣2𝜐²+𝑓¹_𝑣𝑣𝑓²_𝑣𝑣𝑓³_𝑣𝑣𝑣𝑣+⋯, (80) where 𝑓¹_𝑣2=−𝛽𝑘𝜇Δ₀𝜆𝑑𝑞₁(0)+𝑞₃(0)𝑞₃(0)+𝜆𝑑𝑞₁(0)𝑞₃(0)+12𝜆𝑑²×𝑖𝑓₂₀𝜔₀𝜏_𝑘𝑞(0)𝑒^{𝑖𝜔0𝜏𝑘𝑠}+𝑖𝑓₀₂3𝜔₀𝜏_𝑘𝑞(0)𝑒^{−𝑖𝜔0𝜏𝑘𝑠}𝑞(𝑠)+𝐸₁𝑒^{2𝑖𝜔0𝜏𝑘𝑠},𝑓²_𝑣2=𝛽𝑘𝜇Δ₀𝜆𝑑𝑞₁(0)+𝑞₃(0)𝑞₃(0)+𝜆𝑑𝑞₁(0)𝑞₃(0)+12𝜆𝑑²×𝑖𝑓₂₀𝜔₀𝜏_𝑘𝑞(0)𝑒^{𝑖𝜔0𝜏𝑘𝑠}+𝑖𝑓₀₂3𝜔₀𝜏_𝑘𝑞(0)𝑒^{−𝑖𝜔0𝜏𝑘𝑠}𝑞(𝑠)+𝐸₁𝑒^{2𝑖𝜔0𝜏𝑘𝑠}+𝜆𝑑𝑞₁(0)+𝑞₃(0)𝑞₃(0)𝜆+2Re𝑑𝑞₁(0)𝑞₃(0+𝜆)𝑑²𝑊⁽³⁾₁₁()+𝜆−1𝑑𝑞₁(0)+𝑞₃(0)𝑞₃(0)𝜆+2Re𝑑𝑞₁(0)𝑞₃(0+𝜆)𝑑²𝑊⁽³⁾₁₁(),𝑓−1_3𝑣2=−𝛽𝑘𝑒^{−𝑑𝜏𝑘}𝜇Δ₋₁𝜆𝑑𝑞₁(0)+𝑞₃(0)𝑞₃(0)+𝜆𝑑𝑞₁(0)𝑞₃(0)+12𝜆𝑑²×𝑖𝑓₂₀𝜔₀𝜏_𝑘𝑞(0)𝑒^{𝑖𝜔0𝜏𝑘𝑠}+𝑖𝑓₀₂3𝜔₀𝜏_𝑘𝑞(0)𝑒^{−𝑖𝜔0𝜏𝑘𝑠}𝑞(𝑠)+𝐸₁𝑒^{2𝑖𝜔0𝜏𝑘𝑠}. (81)(82)(83)

By (70) and (71), we have 𝐻₂₀(𝑠)=−𝑓₂₀𝑞(𝑠)−𝑓₀₂𝑞(𝑠)+𝑓_𝑣2,𝐻₁₁(𝑠)=−𝑓₁₁𝑞(𝑠)−𝑓₁₁𝑞(𝑠)+𝑓_𝑣𝑣. (84)

For 𝑖𝜔₀𝜏_𝑘 is the eigenvalue of 𝐴(0) and 𝑞(0) is the corresponding eigenvector, we obtain 𝑖𝜔₀𝜏_𝑘𝐼−⁰₋₁𝑒^{𝑖𝜔0𝜏𝑘𝑠}d𝜂(𝑠)𝑞(0)=0,−𝑖𝜔₀𝜏_𝑘𝐼−⁰₋₁𝑒^{−𝑖𝜔0𝜏𝑘𝑠}d𝜂(𝑠)𝑞(0)=0. (85)

So substituting (66) and (73) into (68), we obtain 2𝑖𝜔₀𝜏_𝑘𝐼−⁰₋₁𝑒^{2𝑖𝜔0𝜏𝑘𝑠}d𝜂(s)𝐸₁=𝑓¹_𝑣2𝑓²_𝑣2𝑓³_𝑣2, (86) that is, 2𝑖𝜔₀𝜏_𝑘+()𝐴𝑒𝑑+𝐴0𝐵_{−𝑑𝜏𝑘}𝑒^{−2𝑖𝜔0𝜏𝑘}−12𝑖𝜔₀𝜏_𝑘+𝑑𝐵𝑒^{−𝑑𝜏𝑘}𝑒^{−2𝑖𝜔0𝜏𝑘}−1−𝐴𝑒^{−𝑑𝜏𝑘}𝑒^{−2𝑖𝜔0𝜏𝑘}02𝑖𝜔₀𝜏_𝑘+𝑎+𝐵𝑒^{−𝑑𝜏𝑘}𝑒^{−2𝑖𝜔0𝜏𝑘}𝐸₁=𝑓¹_𝑣2𝑓²_𝑣2𝑓³_𝑣2. (87)

It follows that𝐸₁⁽¹⁾=1𝑀₁𝑓det𝑓^1𝑣20𝐵_2𝑣22𝑖𝜔₀𝜏_𝑘+𝑑𝐵𝑒^{−𝑑𝜏𝑘}𝑒^{−2𝑖𝜔0𝜏𝑘}𝑓−1_3𝑣202𝑖𝜔₀𝜏_𝑘+𝑎+𝐵𝑒^{−𝑑𝜏𝑘}𝑒^{−2𝑖𝜔0𝜏𝑘},𝐸₁⁽²⁾=1𝑀₁det2𝑖𝜔₀𝜏_𝑘+(𝑑+𝐴)𝑓¹_𝑣2𝐵𝐴𝑒^{−𝑑𝜏𝑘}𝑒^{−2𝑖𝜔0𝜏𝑘}−1𝑓²_𝑣2𝐵𝑒^{−𝑑𝜏𝑘}𝑒^{−2𝑖𝜔0𝜏𝑘}−1−𝐴𝑒^{−𝑑𝜏𝑘}𝑒^{−2𝑖𝜔0𝜏𝑘}𝑓³_𝑣22𝑖𝜔₀𝜏_𝑘+𝑎+𝐵𝑒^{−𝑑𝜏𝑘}𝑒^{−2𝑖𝜔0𝜏𝑘},𝐸₁⁽³⁾=1𝑀₁det2𝑖𝜔₀𝜏_𝑘+(𝑑+𝐴)0𝑓¹_𝑣2𝐴𝑒^{−𝑑𝜏𝑘}𝑒^{−2𝑖𝜔0𝜏𝑘}−12𝑖𝜔₀𝜏_𝑘+𝑑𝑓²_𝑣2−𝐴𝑒^{−𝑑𝜏𝑘}𝑒^{−2𝑖𝜔0𝜏𝑘}0𝑓³_𝑣2, (88)

where𝑀₁=det2𝑖𝜔₀𝜏_𝑘+()𝐴𝑒𝑑+𝐴0𝐵_{−𝑑𝜏𝑘}𝑒^{−2𝑖𝜔0𝜏𝑘}−12𝑖𝜔₀𝜏_𝑘+𝑑𝐵𝑒^{−𝑑𝜏𝑘}𝑒^{−2𝑖𝜔0𝜏𝑘}−1−𝐴𝑒^{−𝑑𝜏𝑘}𝑒^{−2𝑖𝜔0𝜏𝑘}02𝑖𝜔₀𝜏_𝑘+𝑎+𝐵𝑒^{−𝑑𝜏𝑘}𝑒^{−2𝑖𝜔0𝜏𝑘}. (89)

Similarly, substituting (67) and (74) into (69), we get𝐸₂⁽¹⁾=1𝑀₂𝑓det𝑓^1𝑣𝑣0𝐵_2𝑣𝑣𝑑𝐵𝑒^{−𝑑𝜏𝑘}𝑒^{−𝑖𝜔0𝜏𝑘}𝑓−1_3𝑣𝑣0𝑎+𝐵𝑒^{−𝑑𝜏𝑘}𝑒^{−𝑖𝜔0𝜏𝑘},𝐸₂⁽²⁾=1𝑀₂()𝑓det𝑑+𝐴¹_𝑣𝑣𝐵𝐴𝑒^{−𝑑𝜏𝑘}𝑒^{−𝑖𝜔0𝜏𝑘}−1𝑓²_𝑣𝑣𝐵𝑒^{−𝑑𝜏𝑘}𝑒^{−𝑖𝜔0𝜏𝑘}−1−𝐴𝑒^{−𝑑𝜏𝑘}𝑒^{−𝑖𝜔0𝜏𝑘}𝑓³_𝑣𝑣𝑎+𝐵𝑒^{−𝑑𝜏𝑘}𝑒^{−𝑖𝜔0𝜏𝑘},𝐸₂⁽³⁾=1𝑀₂()det𝑑+𝐴0𝑓¹_𝑣𝑣𝐴𝑒^{−𝑑𝜏𝑘}𝑒^{−𝑖𝜔0𝜏𝑘}−1𝑑𝑓²_𝑣𝑣−𝐴𝑒^{−𝑑𝜏𝑘}𝑒^{−𝑖𝜔0𝜏𝑘}0𝑓³_𝑣𝑣, (90)

where𝑀₂()=det𝐴𝑒𝑑+𝐴0𝐵_{−𝑑𝜏𝑘}𝑒^{−𝑖𝜔0𝜏𝑘}−1𝑑𝐵𝑒^{−𝑑𝜏𝑘}𝑒^{−𝑖𝜔0𝜏𝑘}−1−𝐴𝑒^{−𝑑𝜏𝑘}𝑒^{−𝑖𝜔0𝜏𝑘}0𝑎+𝐵𝑒^{−𝑑𝜏𝑘}𝑒^{−𝑖𝜔0𝜏𝑘}. (91)

Thus, we can determine 𝑊₂₀(𝑠) and 𝑊₁₁(𝑠) from (66) and (67). Furthermore, we can compute 𝑓₂₁ by (55). Thus, we can compute the following values: 𝑐₍₁₎(0)=𝑖2𝜔₀𝜏_𝑘𝑓₂₀𝑓₁₁−2𝑓₁₁²−𝑓₀₂²3+𝑓₂₁2,𝜗₂𝑐=−Im₁(0)Re𝜆^′𝜏_𝑘,𝜎₂=2Re𝑐₍₁₎(0,𝑇₂)𝑐=−Im₁(0)+𝜗₂Re𝜆^′𝜏_𝑘𝜔₀𝜏_𝑘. (92)

Then from [8], we can give the properties of bifurcating periodic solution in the center manifold at the critical values 𝜏_𝑘. More specifically, we have the following theorem.

Theorem 10. The properties of the Hopf bifurcation are determined by the values in (75).(1)𝜗₂ determines the directions of the Hopf bifurcation: if 𝜗₂>0(𝜗₂<0), then the Hopf bifurcation is supercritical (subcritical), and the bifurcating periodic solutions exist for 𝜏>𝜏_𝑘(𝜏<𝜏_𝑘).(2)𝜗₂ determines the stability of the bifurcating periodic solutions: the bifurcating periodic solutions are stable (unstable) if 𝜗₂>0(𝜗₂<0).(3)𝑇₂ determines the period of the bifurcating periodic solutions: the period increases (decreases) if 𝑇₂>0  (𝑇₂<0).

4. Numerical Simulations

As an example, we consider the system given in [1] with 𝜆=1.5×10⁵, 𝑑=0.011, 𝛽=9.8×10⁻⁷, 𝑘=150, 𝑎=0.1, and  𝜇=0.67, then, there is a chronic infected equilibrium 𝐸^∗ which satisfies the conditions indicated in Theorem 9. When 𝜏=0, the chronic infected equilibrium 𝐸^∗=(−1.3515×10⁷,−1.0120×10⁵,1.2090×10⁵) 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.6029, 𝜏₀=8.1786 and 𝜆(𝜏₀)=0.2103−0.3942𝑖. From the formulae (92) in Section 3, it follows that 𝑐₁(0)=−0.1052−0.3061𝑖, 𝜗₂=8.2963>0,  𝜎₂=−0.2103<0, and 𝑇₂=1.6737>0. Thus, the chronic infected equilibrium 𝐸^∗ is asymptotically stable when 0≤𝜏<𝜏₀ as illustrated by the computer simulations (see Figures 2 and 3).

(a)

(b)

Figure 1 

The response of the system (5) initial data (1.2𝑒7,1.8𝑒6,6𝑒5)𝑇, when 𝜏=0.

Figure 2 

The response of the system (5) initial data (1.2𝑒7,1.8𝑒6,6𝑒5)𝑇, when 𝜏=8.1786<𝜏0.

Figure 3 

3-phase portrait of the system (5) initial data (1.2𝑒7,1.8𝑒6,6𝑒5)𝑇, when 𝜏=8.1786<𝜏0.

When 𝜏 passes through the critical value 𝜏₀, 𝐸^∗ loses its stability, and a Hopf bifurcation occurs, that is, a family of periodic solutions bifurcate from 𝐸^∗. Since 𝜗₂ and 𝜎₂=−0.2103<0, the Hopf bifurcation is supercritical, and the direction of the bifurcation is 𝜏>𝜏₀, and these bifurcating periodic solutions from 𝐸^∗ at 𝜏₀ are stable, which are depicted in Figures 2 and 3.

In Figure 4, we try to reflect the changes of stability of chronic infected equilibrium 𝐸^∗ as 𝜏 increases from 0 to 9. In Figure 4, each vertical blue strip corresponds to component of  𝑥(𝑡),𝑒(𝑡), and 𝑦(𝑡) at each 𝜏, respectively. From Figure 3, we see that if 𝜏∈(0,8.17864293), approximately, the vertical amplitudes of 𝑥(𝑡), 𝑒(𝑡), and 𝑦(𝑡) are as small as a point, suggesting that 𝐸^∗ is asymptotically stable; if 𝜏 increases, the vertical amplitudes of 𝑥(𝑡), 𝑒(𝑡), and 𝑦(𝑡) will become larger and larger, showing that 𝐸^∗ 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.2𝑒7,1.8𝑒6,6𝑒5)𝑇, when 𝜏=8.21>𝜏0.

Figure 5 

3-phase portrait of the system (5) initial data (1.2𝑒7,1.8𝑒6,6𝑒5)𝑇, 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.

Acknowledgments

This work was partially supported by the National Natural Science Foundation of China (Grant nos. 11201179 and 11026112). This work was partially supported by the Foundation of University of Jinan (Grant no. XKY1133) and the Shandong Province Young and Middle-Aged Scientists Research Awards Fund (Grant no. BS2011SF001). The authors would like to thank the reviewers for their careful reading of the paper and for their valuable comments and suggestions.