Stability and Local Hopf Bifurcation for a Predator-Prey Model with Delay

Abstract

A predator-prey system with disease in the predator is investigated, where the discrete delay 𝜏 is regarded as a parameter. Its dynamics are studied in terms of local analysis and Hopf bifurcation analysis. By analyzing the associated characteristic equation, it is found that Hopf bifurcation occurs when 𝜏 crosses some critical values. Using the normal form theory and center manifold argument, the explicit formulae which determine the stability, direction, and other properties of bifurcating periodic solutions are derived.

1. Introduction

Many models in ecology can be formulated as system of differential equations with time delays. The effect of the past history on the stability of system is also an important problem in population biology. Recently, the properties of periodic solutions arising from the Hopf bifurcation have been considered by many authors [1–4].

May [5] first proposed and discussed the delayed predator-prey system 𝑑𝑥(𝑡𝑟𝑑𝑡=𝑥)₁−𝑎₁₁𝑥(𝑡−𝜏)−𝑎₁₂𝑦(𝑡),𝑑𝑦(𝑡𝑑𝑡=𝑦)−𝑟₂+𝑎₂₁𝑦(𝑡)−𝑎₂₂𝑦(𝑡), (1.1) where 𝑥(𝑡) and 𝑦(𝑡) can be interpreted as the population densities of prey and predator at time 𝑡, respectively; 𝜏≥0 is the feedback time delay of the prey to the growth of the species itself; 𝑟₁>0 denotes the intrinsic growth rate of the prey, and 𝑟₂>0 denotes the death rate of the predator; the parameter 𝑎_𝑖𝑗(𝑖,𝑗=1,2) are all positive constants. System (1,1) shows that, in the absence of predator species, the prey species are governed by the well-known delayed logistic equation 𝑑𝑥/𝑑𝑡=𝑥(𝑡)[𝑟₁−𝑎₁₁𝑥(𝑡−𝜏)] and the predator species will decrease in the absence of the prey species.There has been an extensive literature dealing with system (1,1) or the system similar to (1.1), regarding boundedness of solutions, persistence, local and global stabilities of equilibria, and existence of nonconstant periodic solutions [6–9].

Recently, Faria [7] investigated the stability and Hopf bifurcation of the following system with instantaneous feedback control and two different discrete delays: 𝑑𝑥(𝑡𝑟𝑑𝑡=𝑥)₁−𝑎₁₁𝑥(𝑡)−𝑎₁₂𝑦𝑡−𝜏₁𝑑𝑦(𝑡𝑑𝑡=𝑦)−𝑟₂+𝑎₂₁𝑥𝑡−𝜏₂−𝑎₂₂𝑦(𝑡), (1.2) where 𝜏₁>0 and 𝜏₂>0. But, as pointed out by Kuang [8], in view of the fact that in real situations, instantaneous responses are rare, and thus, more realistic models should consist of delay differential equations without instantaneous feedbacks. Based on this idea, in the present paper, we combine the model (1.1) and (1.2) and consider the following delayed prey-predator system with a single delay: 𝑑𝑋(𝑡𝑟𝑑𝑡=𝑋)₁−𝑟₂𝑋(𝑡−𝜏)−𝑝𝑆(𝑡−𝜏)𝑑𝑆(𝑡𝑑𝑡=𝑆)−𝑐₁+𝑘𝑝𝑋(𝑡−𝜏)−𝜎𝐼(𝑡)+𝛾𝐼(𝑡)𝑑𝐼(𝑡(𝑡)𝑑𝑡=𝐼)𝜎𝑆−𝑐₂−𝛾, (1.3) where 𝑋,𝑆,𝐼 denote, respectively, the population of prey species, susceptible predator species and infected predator species. In addition, the coefficients 𝑟1,𝑟2,𝑝,𝑘,𝜎,𝑐₁,𝑐₂ in model (1.3) are all positive constants and their ecological meaning are interpreted as follows: 𝑟₁ denotes the intrinsic growth rate of prey and 𝑟₁/𝑟₂ denotes the carrying capacity of prey; 𝑝,𝑘,𝑐₁ and 𝑐₂ represent the predating coefficient of predator to prey, absorbing rate of predator to prey, and the death rate of susceptible and infected predator, respectively.

The main purpose of this paper is to investigate the effects of the delay on the dynamics of model (1.3) with the following initial conditions: 𝑋(𝑡)=𝜙₁(𝑡)>0,𝑆(𝑡)=𝜙₂(𝑡)>0,𝐼(𝑡)=𝜙₃(𝑡)[]𝜙>0,𝑡∈−𝜏,0₁(𝑡),𝜙₂(𝑡),𝜙₃(𝑡)∈𝐶[−𝜏,0],ℜ³_0,+, (1.4) where ℜ³_0,+={(𝑥,𝑦,𝑧)∣𝑥≥0,𝑦≥0,𝑧≥0}. We will take the delay 𝜏 as the bifurcation parameter and show that when 𝜏 passes through a certain critical value, the positive equilibrium loses its stability and a Hopf bifurcation will take place. Furthermore, when 𝜏 takes a sequence of critical values containing the above critical value, the positive equilibrium of system (1.3) will undergo a Hopf bifurcation. In particular, by using the normal form theory and the center manifold, the formulae determining the direction of Hopf bifurcation and the stability of bifurcating periodic solutions are also obtained.

The organization of this paper is as follows. In Section 2, we discuss the stability of the positive solutions and the existence of the Hopf bifurcations. In Section 3, the direction of the Hopf bifurcation and the stability of bifurcated periodic solutions are obtained by using the normal form theory and the center manifold theorem. In Section 4, we do some numerical simulations to validate our theoretical results.

2. Stability of Positive Equilibrium and Hopf Bifurcation

System (1.3) has a unique positive equilibrium 𝐸⁺(𝑋^∗,𝑆^∗,𝐼^∗) provided that the condition(H₁)𝜎𝑟₁>(𝑐₂+𝛾)𝑝,𝑘𝑝[𝜎𝑟₁−𝑝(𝑐₂+𝛾)]>𝜎𝑐₁𝑟₂is satisfied, where 𝑋^∗=𝜎𝑟₁−𝑝𝑐₂+𝛾𝜎𝑟₂,𝑆^∗=𝑐₂+𝛾𝜎,𝐼^∗=𝑐₂+𝛾𝑘𝑝𝜎𝑟₁−𝑝𝑐₂+𝛾−𝜎𝑐₁𝑟₂𝜎²𝑟₂𝑐₂+2𝛾. (2.1) Linearizing system (1.2) at 𝐸⁺ gives the following linear system: 𝑑𝑈𝑑𝑡=−𝑟₂𝑋^∗𝑈(𝑡−𝜏)−𝑝𝑋^∗𝑉(𝑡−𝜏)𝑑𝑉𝑑𝑡=𝑘𝑝𝑆^∗𝑈(𝑡−𝜏)−𝑐₁+𝜎𝐼^∗−𝑘𝑝𝑋^∗𝑉(𝑡)−𝜎𝑆^∗−𝛾𝑊𝑑𝑊𝑑𝑡=𝜎𝐼^∗𝑉. (2.2)

The characteristic matrix of this system (2.2) is 𝜆+𝑟₂𝑋^∗𝑒^−𝜆𝜏𝑝𝑋^∗𝑒^−𝜆𝜏0−𝑘𝑝𝑆^∗𝑒^−𝜆𝜏𝜆+𝑐₁+𝜎𝐼^∗−𝑘𝑝𝑋^∗𝜎𝑆^∗−𝛾0−𝜎𝐼^∗𝜆. (2.3) Thus, the characteristic equation of system (2.2) is given by 𝜆³+𝑐₁+𝜎𝐼^∗−𝑘𝑝𝑋^∗𝜆²+𝜎²𝑆^∗𝐼^∗−𝜎𝛾𝐼^∗𝜆+𝑟₂𝑋^∗𝑒^−𝜆𝜏𝜆²+𝑐₁+𝜎𝐼^∗−𝑘𝑝𝑋^∗𝜆𝑟₂𝑋^∗𝑒^−𝜆𝜏+𝜎²𝑆^∗𝐼^∗−𝜎𝛾𝐼^∗𝑟₂𝑋^∗𝑒^−𝜆𝜏+𝑘𝑝²𝑆^∗𝑋^∗𝜆𝑒^−2𝜆𝜏=0. (2.4) Let 𝑑₁=𝑐₁+𝜎𝐼^∗−𝑘𝑝𝑋^∗,𝑑₂=𝜎²𝑆^∗𝐼^∗−𝜎𝛾𝐼^∗,𝑑₃=𝑟₂𝑋^∗,𝑑₄=𝑐₁+𝜎𝐼^∗−𝑘𝑝𝑋^∗𝑟₂𝑋^∗,𝑑₅=𝜎²𝑆^∗𝐼^∗−𝜎𝛾𝐼^∗𝑟₂𝑋^∗,𝑑₆=𝑘𝑝²𝑆^∗𝑋^∗. (2.5) Then we rewrite (2.4) as: 𝜆³+𝑑₁𝜆²+𝑑₂𝜆𝑒^𝜆𝜏+𝑑₃𝜆²+𝑑₄𝜆+𝑑₅+𝑑₆𝜆𝑒^−𝜆𝜏=0. (2.6)

Obviously, 𝑖𝜔(𝜔>0) is a root of (2.6) if and only if 𝜔 satisfies −𝑖𝜔³−𝑑₁𝜔²+𝑑₂𝑖𝜔(cos𝜔𝜏+𝑖sin𝜔𝜏)−𝑑₃𝜔²+𝑑₄𝑖𝜔+𝑑₅+𝑑₆𝑖𝜔(cos𝜔𝜏−𝑖sin𝜔𝜏)=0. (2.7) Separating the real and imaginary parts, we have −𝑑₁𝜔²cos𝜔𝜏+𝜔³−𝑑₂𝜔+𝑑₆𝜔sin𝜔𝜏=𝑑₃𝜔²−𝑑₅−𝜔³+𝑑₂𝜔+𝑑₆𝜔cos𝜔𝜏−𝑑₁𝜔²sin𝜔𝜏=−𝑑₄𝜔. (2.8) By calculating, we have obtained 𝑑sin𝜔𝜏=³𝜔⁵+𝑑₁𝑑₄−𝑑₂𝑑₃−𝑑₃𝑑₆−𝑑₅𝜔³+𝑑₂𝑑₅+𝑑₅𝑑₆𝜔𝜔⁶+𝑑²₁−2𝑑₂𝜔⁴+𝑑²₂−𝑑²₆𝜔²𝑑cos𝜔𝜏=⁴−𝑑₁𝑑₃𝜔⁴+𝑑₁𝑑₅+𝑑₄𝑑₆−𝑑₂𝑑₄𝜔²𝜔⁶+𝑑²₁−2𝑑₂𝜔⁴+𝑑²₂−𝑑²₆𝜔². (2.9) Let 𝑒₁=𝑑²₁−2𝑑₂, 𝑒₂=𝑑²₂−𝑑²₆, 𝑒₃=𝑑₃, 𝑒₄=𝑑₁𝑑₄−𝑑₂𝑑₃−𝑑₃𝑑₆−𝑑₅, 𝑒₅=𝑑₂𝑑₅+𝑑₅𝑑₆, 𝑒₆=𝑑₄−𝑑₁𝑑₃, 𝑒₇=𝑑₁𝑑₅+𝑑₄𝑑₆−𝑑₂𝑑₄. Then sin𝜔𝜏,cos𝜔𝜏 can be written as 𝜔𝑒sin𝜔𝜏=³𝜔⁴+𝑒₄𝜔²+𝑒₅𝜔⁶+𝑒₁𝜔⁴+e₂𝜔²𝑒cos𝜔𝜏=⁶𝜔⁴+𝑒₇𝜔²𝜔⁶+𝑒₁𝜔⁴+𝑒₂𝜔². (2.10)(2.11) As sin²𝜔𝜏+cos²𝜔𝜏=1, so we have 𝜔¹⁰+𝑓₁𝜔⁸+𝑓₂𝜔⁶+𝑓₃𝜔⁴+𝑓₄𝜔²+𝑓₅=0, (2.12) where 𝑓₁=2𝑒₁−𝑒²₃, 𝑓₂=𝑒²₁+2𝑒₂−2𝑒₃𝑒₄−𝑒²₆, 𝑓₃=2𝑒₁𝑒₂−2𝑒₃𝑒₅−2𝑒₆𝑒₇−𝑒²₄, 𝑓₄=𝑒²₂−2𝑒₄𝑒₅−𝑒²₇, 𝑓₅=−𝑒²₅.

Denote 𝑧=𝜔², then (2.12) becomes 𝑧⁵+𝑓₁𝑧⁴+𝑓₂𝑧³+𝑓₃𝑧²+𝑓₄𝑧+𝑓₅=0. (2.13) Let 𝐺(𝑧)=𝑧⁵+𝑓₁𝑧⁴+𝑓₂𝑧³+𝑓₃𝑧²+𝑓₄𝑧+𝑓₅. (2.14) Since lim_𝑧→∞𝐺(𝑧)=+∞ and 𝑓₅<0, then we can get the following conclusion.(H₂) Equation (2.13) has at least one positive real root.

Without loss of generality, we assume that it has five positive roots, defined by 𝑧₁,𝑧₂,𝑧₃,𝑧₄,𝑧₅, respectively. Then (2.13) has five positive roots 𝜔₁=√𝑧₁,𝜔₂=√𝑧₂,𝜔₃=√𝑧₃,𝜔₄=√𝑧₄,𝜔₅=√𝑧₅. (2.15) By (2.11), we have cos𝜔_𝑘𝑒𝜏=⁶𝜔²_𝑘+𝑒₇𝜔⁴_𝑘+𝑒₁𝜔²_𝑘+𝑒₂. (2.16) Thus, if we denote 𝜏_𝑘^(𝑗)=1𝜔_𝑘𝑒arccos⁶𝜔²_𝑘+𝑒₇𝜔⁴_𝑘+𝑒₁𝜔²_𝑘+𝑒₂+2𝑗𝜋, (2.17) where 𝑘=1,…,5;𝑗=0,1,…, then ±𝑖𝜔_𝑘 is a pair of purely imaginary roots of (2.6) with 𝜏_𝑘^(𝑗). Define 𝜏₀=𝜏_𝑘⁽⁰⁾₀=min_{𝑘∈{1,…,5}}𝜏_𝑘⁽⁰⁾,𝜔₀=𝜔_𝑘0. (2.18) Note that when 𝜏=0, (2.6) becomes 𝜆³+𝑏𝜆²+(𝑎+𝑑)𝜆+𝑐=0. (2.19) By Routh-Hurwitz criterion, we know that all the roots of (2.19) have negative real parts, that is, the positive equilibrium 𝐸⁺ is locally asymptotically stable for 𝜏=0.

In order to give the main results, it is necessary to make the following assumption:(H₃)[𝑑(Re𝜆)/𝑑𝜏]_𝜏=𝜏0≠0.

Differentiating two sides of (2.6) in respect to 𝜏, we get 𝑑𝜆𝑑𝜏⁻¹=3𝜆²+2𝑑₁𝜆+𝑑₂𝜆𝑒^𝜆𝜏+𝜏𝜆³+𝑑₁𝜆²+𝑑₂𝜆𝑒^𝜆𝜏+2𝑑₃𝜆+𝑑₄+𝑑₆𝑒^−𝜆𝜏−𝑑₆𝜆𝜏𝑒^−𝜆𝜏−𝜆𝜆³+𝑑₁𝜆²+𝑑₂𝜆𝑒^𝜆𝜏+𝑑₆𝜆²𝑒^−𝜆𝜏=3𝜆²+2𝑑₁𝜆+𝑑₂𝜆𝑒^𝜆𝜏+2𝑑₃𝜆+𝑑₄+𝑑₆𝑒^−𝜆𝜏−𝜆𝜆³+𝑑₁𝜆²+𝑑₂𝜆𝑒^𝜆𝜏+𝑑₆𝜆²𝑒^−𝜆𝜏−𝜏𝜆=3𝜆²+2𝑑₁𝜆+𝑑₂𝜆𝑒^𝜆𝜏+2𝑑₃𝜆+𝑑₄+𝑑₆𝑒^−𝜆𝜏𝑑₃𝜆³+𝑑₄𝜆²+𝑑₅𝜆+2𝑑₆𝜆²𝑒^−𝜆𝜏−𝜏𝜆=−3𝜔²+𝑑₂+𝑑₆cos𝜔𝜏−2𝑑₁𝜔sin𝜔𝜏+𝑑₄−𝑑₄𝜔²−2𝑑₆𝜔²cos𝜔𝜏+𝑖−𝑑₃𝜔³+𝑑₅𝜔+2𝑑₆𝜔²+sin𝜔𝜏−3𝜔²+𝑑₂−𝑑₆sin𝜔𝜏+2𝑑₁𝜔cos𝜔𝜏+2𝑑₃𝜔−𝑑₄𝜔²−2𝑑₆𝜔²cos𝜔𝜏+𝑖−𝑑₃𝜔³+𝑑₅𝜔+2𝑑₆𝜔².sin𝜔𝜏 (2.20) Let 𝑄=−𝑑₄𝜔²−2𝑑₆𝜔²cos𝜔𝜏²+−𝑑₃𝜔³+𝑑₅𝜔+2𝑑₆𝜔²sin𝜔𝜏²>0𝑄Re𝑑𝜆𝑑𝜏⁻¹=−3𝜔²+𝑑₂+𝑑₆cos𝜔𝜏−2𝑑₁𝜔sin𝜔𝜏+𝑑₄−𝑑₄𝜔²−2𝑑₆𝜔²+cos𝜔𝜏−3𝜔²+𝑑₂−𝑑₆sin𝜔𝜏+2𝑑₁𝜔cos𝜔𝜏+2𝑑₃𝜔−𝑑₃𝜔³+𝑑₅𝜔+2𝑑₆𝜔²sin𝜔𝜏. (2.21)(2.22) Noting that sgn𝑑(Re𝜆)𝑑𝜏_𝜏=𝜏0=sgnRe𝑑𝜆𝑑𝜏⁻¹_𝜏=𝜏0. (2.23) Now, we can employ a result from Ruan and Wei [10] to analyze (2.6), which is, for the convenience of the reader, stated as follows.

Lemma 2.1 (see [2]). Consider the exponential polynomial 𝑃𝜆,𝑒^−𝜆𝜏1,…,𝑒^{−𝜆𝜏𝑚}=𝜆^𝑛+𝑝₁⁽⁰⁾𝜆^𝑛−1+⋯+𝑝⁽⁰⁾_𝑛−1𝜆+𝑝_𝑛⁽⁰⁾+𝑝₁⁽¹⁾𝜆^𝑛−1+⋯+𝑝⁽¹⁾_𝑛−1𝜆+𝑝_𝑛⁽¹⁾𝑒^−𝜆𝜏1+⋯+𝑝₁^(𝑚)𝜆^𝑛−1+⋯+𝑝^(𝑚)_𝑛−1𝜆+𝑝_𝑛^(𝑚)𝑒^{−𝜆𝜏𝑚}, (2.24) where 𝜏_𝑖≥0(𝑖=1,2,…,𝑚) and 𝑝_𝑗^(𝑖)(𝑖=1,2,…,𝑚;𝑗=1,2,…,𝑚) are constants. As (𝜏₁,𝜏₂,…,𝜏_𝑚) vary, the sum of the order of the zeros of 𝑃(𝜆,𝑒^−𝜆𝜏1,…,𝑒^{−𝜆𝜏𝑚}) on the open right half plane can change only if a zero appears on or crosses the imaginary axis.

Form Lemma 2.1, it is easy to obtain the following theorem.

Theorem 2.2. Suppose the condition, (𝐻₁), (𝐻₂), and (𝐻₃) are satisfied, then one has the following results:(i) if 𝜏∈[0,𝜏₀), then the positive equilibrium 𝐸⁺ of (1.2) is locally asymptotically stable and unstable when 𝜏>𝜏₀,(ii) system (1.2) undergoes a Hopf bifurcation at the positive equilibrium 𝐸⁺ when 𝜏=𝜏_𝑘(𝑘=0,1,2,…), where 𝜏_𝑘 is defined by (2.17).

3. Stability and Direction of Hopf Bifurcation

In this section, we will derive the explicit formulae determining the properties of the Hopf bifurcation at the critical value using the normal form theory and center manifold theorem introduced by Hassard et al. [11].

Without loss of generality, let 𝜏=𝜏_𝑘+𝜇, where 𝜏_𝑘 is defined by (2.17), 𝜇∈𝑅, then system (1.3) can be rewritten as 𝑢^′(𝑡)=𝐿_𝜇𝑢_𝑡+𝑓𝜇,𝑢_𝑡, (3.1) where 𝑢(𝑡)=(𝑢₁(𝑡),𝑢₂(𝑡),𝑢₃(𝑡))^𝑇=(𝑈(𝜏𝑡),𝑉(𝜏𝑡),𝑊(𝜏𝑡))^𝑇∈𝑅³, 𝑢_𝑡(𝜃)=𝑢(𝑡+𝜃) and 𝐿_𝜇∶𝐶→𝑅³, 𝑓∶𝑅×𝐶→𝑅³ are given, respectively, by 𝐿_𝜇𝑢_𝑡=𝜏_𝑘+𝜇0000𝑘𝑝𝑋^∗−𝑐₁−𝜎𝐼^∗−𝜎𝑆^∗+𝛾0𝜎𝐼^∗0𝑢_1𝑡(0)𝑢_2𝑡(0)𝑢_3𝑡(0)+𝜏_𝑘+𝜇−𝑟₂𝑋^∗−𝑝𝑋^∗0𝑘𝑝𝑆^∗𝑢00000^1𝑡()𝑢−1_2𝑡()𝑢−1_3𝑡(),𝑓−1𝜇,𝑢_𝑡=𝜏_𝑘+𝜇−𝑟₂𝑢_1𝑡(0)𝑢_1𝑡(−1)−𝑝𝑢_1𝑡(0)𝑢_2𝑡(−1)𝑘𝑝𝑢_2𝑡(0)𝑢_1𝑡(−1)−𝜎𝑢_2𝑡(0)𝑢_3𝑡(0)𝜎𝑢_2𝑡(0)𝑢_3𝑡(0). (3.2)(3.3) By the Riesz representation theorem, there exists a function 𝜂(𝜃,𝜇) of bounded variation for 𝜃∈[−1,0], such that 𝐿_𝜇𝜙=⁰₋₁𝑑𝜂𝜃,𝜇𝜙(𝜃)for𝜙∈[−1,0],𝑅³. (3.4) In fact, we can choose 𝜂𝜃,𝜇=𝜏_𝑘+𝜇0000𝑘𝑝𝑋^∗−𝑐₁−𝜎𝐼^∗−𝜎𝑆^∗+𝛾0𝜎𝐼^∗0𝛿(𝜃)−𝜏_𝑘+𝜇−𝑟₂𝑋^∗−𝑝𝑋^∗0𝑘𝑝𝑆^∗00𝛿(),000𝜃+1 (3.5) where 𝛿 denote the Dirac delta function. For 𝜙∈𝐶([−1,0],𝑅³), define 𝐴𝜇(𝜃)𝜙=𝑑𝜙[),𝑑𝜃,𝜃∈−1,0₀₋₁𝜙(𝑠)𝑅𝜇𝜙=[),𝑓𝑑𝜂𝑠,𝜇,𝜃=0,0,𝜃∈−1,0𝜇,𝜙,𝜃=0. (3.6) Then system (3.1) is equivalent to 𝑢^′(𝑡)=𝐴𝜇𝑢_𝑡+𝑅𝜇𝑢_𝑡. (3.7) For 𝜓∈𝐶([0,1],(𝑅³)^∗), define 𝐴^∗=−𝑑𝜓(𝑠)(],𝑑s,𝑠∈0,1₀₋₁𝜓(−𝑡)𝑑𝜂(𝑡,0),𝑠=0 (3.8) and a bilinear inner product 𝜓(𝑠),𝜙(𝜃)=𝜓(0)𝜙(0)−⁰₋₁^𝜃_𝜉=0𝜓(𝜉−𝜃)𝑑𝜂(𝜃)𝜙(𝜉)𝑑𝜉, (3.9) 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 𝑞(𝜃)=(1,𝑎₁,𝑎₂)^𝑇𝑒^{𝑖𝜔0𝜏𝑘𝜃} is the eigenvector of 𝐴(0) corresponding to 𝑖𝜔₀𝜏_𝑘, then 𝐴(0)𝑞(𝜃)=𝑖𝜔₀𝜏_𝑘𝑞(𝜃). It follows from the definition of 𝐴(0) and (3.2), (3.4) and (3.5), we have 𝜏₀𝑖𝜔₀+𝑟₂𝑋^∗𝑒^{−𝑖𝜔0𝜏0}𝑝𝑋^∗𝑒^{−𝑖𝜔0𝜏0}0−𝑘𝑝𝑋^∗𝑒^{−𝑖𝜔0𝜏0}𝑖𝜔₀−𝑘𝑝𝑋^∗+𝑐₁+𝜎𝐼^∗𝜎𝑆^∗−𝛾0−𝜎𝐼^∗𝑖𝜔₀𝑞(0)=000. (3.10) For 𝑞(−1)=𝑞(0)𝑒^{−𝑖𝜔0𝜏𝑘}, then we obtain 𝑎₁=𝑖𝜔₀+𝑟₂𝑋^∗𝑒^{−𝑖𝜔0𝜏0}𝑝𝑋^∗𝑒^{−𝑖𝜔0𝜏0},𝑎₂=𝜎𝐼^∗𝑖𝜔₀+𝑟₂𝑋^∗𝑒^{−𝑖𝜔0𝜏0}𝑖𝜔₀𝑝𝑋^∗𝑒^{−𝑖𝜔0𝜏0}. (3.11) On the other hand, suppose that 𝑞^∗(𝑠)=𝐽(1,𝑎^∗₁,𝑎^∗₂)𝑒^{𝑖𝜔0𝜏𝑘𝑆} is the eigenvector of 𝐴^∗ corresponding to −𝑖𝜔₀𝜏_𝑘, by the similar method, we have 𝑎^∗₁=−𝑖𝜔₀+𝑟₂𝑋^∗𝑒^{𝑖𝜔0𝜏0}𝑘𝑝𝑋^∗𝑒^{𝑖𝜔0𝜏0},𝑎^∗₂=−𝜎𝐼^∗−𝑖𝜔₀+𝑟₂𝑋^∗𝑒^{𝑖𝜔0𝜏0}𝑖𝜔₀𝑘𝑝𝑋^∗𝑒^{𝑖𝜔0𝜏0}. (3.12) In order to assure ⟨𝑞^∗(𝑠),𝑞(𝜃)⟩=1, we need to determine the value of 𝐽. From (3.9), we have 𝑞^∗(𝑠),𝑞(𝜃)=𝐽1,𝑎^∗₁,𝑎^∗₂1,𝑎₁,𝑎₂^𝑇−⁰₋₁^𝜃_𝜉=0𝐽1,𝑎^∗₁,𝑎^∗₂𝑒^{−𝑖𝜔0𝜏𝑘}(𝜉−𝜃)𝑑𝜂(𝜃)1,𝑎₁,𝑎₂^𝑇𝑒^{𝑖𝜔0𝜏𝑘𝜉}=𝑑𝜉𝐽1+𝑎₁𝑎^∗₁+𝑎₂𝑎^∗₂−⁰₋₁1,𝑎^∗₁,𝑎^∗₂𝜃𝑒^{𝑖𝜔0𝜏𝑘𝜃}𝑑𝜂(𝜃)1,𝑎₁,𝑎₂^𝑇=𝐽1+𝑎₁𝑎^∗₁+𝑎₂𝑎^∗₂+𝜏_𝑘−𝑟₂𝑋^∗−𝑎₁𝑝𝑋^∗+𝑘𝑝𝑆^∗𝑎^∗₁𝑒^{−𝑖𝜔0𝜏𝑘}. (3.13) Therefore, we can choose 𝐽 as 1𝐽=1+𝑎₁𝑎^∗₁+𝑎₂𝑎^∗₂+𝜏_𝑘−𝑟₂𝑋^∗−𝑎₁𝑝𝑋^∗+𝑘𝑝𝑆^∗𝑎^∗₁𝑒^{−𝑖𝜔0𝜏𝑘}. (3.14) Next we will compute the coordinate to describe the center manifold 𝐶₀ at 𝜇=0. Let 𝑢_𝑡 be the solution of (3.1) with 𝜇=0. Define 𝑧(𝑡)=𝑞^∗,𝑢_𝑡,𝑊(𝑡,𝜃)=𝑢_𝑡(𝜃)−𝑧(𝑡)𝑞(𝜃)−𝑧(𝑡)𝑞^∗(𝜃)=𝑢_𝑡(𝜃)−2Re𝑧(𝑡)𝑞(𝜃). (3.15) On the center manifold 𝐶₀, we have 𝑊(𝑡,𝜃)=𝑊(𝑧(𝑡),𝑧(𝑡),𝜃), (3.16) where 𝑊(𝑧(𝑡),𝑧(𝑡),𝜃)=𝑊₂₀(𝜃)𝑧²2+𝑊₁₁(𝜃)𝑧𝑧+𝑊₀₂(𝜃)𝑧²2+⋯. (3.17)𝑧 and 𝑧 are local coordinates of center manifold 𝐶₀ in the direction of 𝑞^∗ and 𝑞^∗. Note that 𝑊 is real if 𝑢_𝑡 is real. We consider only real solutions. For solution 𝑢_𝑡∈𝐶₀ of (3.7), since 𝜇=0, we have 𝑧^′(𝑡)=𝑞^∗,𝑢^′(𝑡)=𝑞^∗,𝐴𝑢_𝑡+𝑅𝑢_𝑡=𝑞^∗,𝐴𝑢_𝑡+𝑞^∗,𝑅𝑢_𝑡=𝐴^∗𝑞^∗,𝑢_𝑡+𝑞^∗,𝑅𝑢_𝑡=𝑖𝜔₀𝜏_𝑘𝑧+𝑞^∗(0)𝑓(0),𝑊(𝑧,𝑧,0)+2Re𝑧𝑞(0)≜𝑖𝜔₀𝜏_𝑘𝑧+𝑞^∗(0)𝑓₀(𝑧,𝑧). (3.18) We rewrite this equation as 𝑧^′(𝑡)=𝑖𝜔₀𝜏_𝑘𝑧(𝑡)+𝑔(𝑧,𝑧), (3.19) where 𝑔(𝑧,𝑧)=𝑞^∗(0)𝑓₀(𝑧,𝑧)=𝑔₂₀𝑧²2+𝑔₁₁𝑧𝑧+𝑔₀₂𝑧²2+𝑔₂₁𝑧²𝑧2+⋯. (3.20) Noticing 𝑢_𝑡(𝜃)=𝑊(𝑡,𝜃)+𝑧𝑞(𝜃)+𝑧𝑞(𝜃) and 𝑞(𝜃)=(1,𝑎₁,𝑎₂)^𝑇𝑒^{𝑖𝜔0𝜏𝑘𝜃}, we have 𝑢_1𝑡(0)=𝑧+𝑧+𝑊⁽¹⁾₂₀(0)𝑧²2+𝑊⁽¹⁾₁₁(0)𝑧𝑧+𝑊⁽¹⁾₀₂(0)𝑧²2𝑢+⋯,_2𝑡(0)=𝑎₁𝑧+𝑎₁𝑧+𝑊⁽²⁾₂₀(0)𝑧²2+𝑊⁽²⁾₁₁(0)𝑧𝑧+𝑊⁽²⁾₀₂(0)𝑧²2𝑢+⋯,_3𝑡(0)=𝑎₂𝑧+𝑎₂𝑧+𝑊⁽³⁾₂₀(0)𝑧²2+𝑊⁽³⁾₁₁(0)𝑧𝑧+𝑊⁽³⁾₀₂(0)𝑧²2𝑢+⋯,_1𝑡(−1)=𝑒^{−𝑖𝜔0𝜏𝑘}𝑧+𝑒^{𝑖𝜔0𝜏𝑘}𝑧+𝑊⁽¹⁾₂₀()𝑧−1²2+𝑊⁽¹⁾₁₁(−1)𝑧𝑧+𝑊⁽¹⁾₀₂(−1)𝑧²2𝑢+⋯,_2𝑡(−1)=𝑒^{−𝑖𝜔0𝜏𝑘}𝑎₁𝑧+𝑒^{𝑖𝜔0𝜏𝑘}𝑎₁𝑧+𝑊⁽²⁾₂₀()𝑧−1²2+𝑊⁽²⁾₁₁(−1)𝑧𝑧+𝑊⁽²⁾₀₂(−1)𝑧²2+⋯. (3.21) It follows together with (3.3), that 𝑔(𝑧,𝑧)=𝑞^∗(0)𝑓₀(𝑧,𝑧)=𝑞^∗(0)𝑓₀0,𝑢_𝑡.=𝜏_𝑘𝐽1,𝑎^∗₁,𝑎^∗₂−𝑟₂𝑢_1𝑡(0)𝑢_1𝑡(−1)−𝑝𝑢_1𝑡(0)𝑢_2𝑡(−1)𝑘𝑝𝑢_2𝑡(0)𝑢_1𝑡(−1)−𝜎𝑢_2𝑡(0)𝑢_3𝑡(0)𝜎𝑢_2𝑡(0)𝑢_3𝑡(0)=−𝜏_𝑘𝐽𝑟₂𝑧+𝑧+𝑊⁽¹⁾₂₀(0)𝑧²2+𝑊⁽¹⁾₁₁(0)𝑧𝑧+𝑊⁽¹⁾₀₂(0)𝑧²2∗𝑒+⋯_{−𝑖𝜔0𝜏𝑘}𝑧+𝑒^{𝑖𝜔0𝜏𝑘}𝑧+𝑊⁽¹⁾₂₀()𝑧−1²2+𝑊⁽¹⁾₁₁(−1)𝑧𝑧+𝑊⁽¹⁾₀₂(−1)𝑧²2+⋯+𝑝𝑧+𝑧+𝑊⁽¹⁾₂₀(0)𝑧²2+𝑊⁽¹⁾₁₁(0)𝑧𝑧+𝑊⁽¹⁾₀₂(0)𝑧²2∗𝑒+⋯_{−𝑖𝜔0𝜏𝑘}𝑎₁𝑧+𝑒^{−𝑖𝜔0𝜏𝑘}𝑎₁𝑧+𝑊⁽²⁾₂₀()𝑧−1²2+𝑊⁽²⁾₁₁(−1)𝑧𝑧+𝑊⁽²⁾₀₂(−1)𝑧²2+⋯+𝜏_𝑘𝐽𝑎^∗₁𝑘𝑝𝑎₁𝑧+𝑎₁𝑧+W⁽²⁾₂₀(0)𝑧²2+𝑊⁽²⁾₁₁(0)𝑧𝑧+𝑊⁽²⁾₀₂(0)𝑧²2∗𝑒+⋯_{−𝑖𝜔0𝜏𝑘}𝑧+𝑒^{𝑖𝜔0𝜏𝑘}𝑧+𝑊⁽¹⁾₂₀()𝑧−1²2+𝑊⁽¹⁾₁₁(−1)𝑧𝑧+𝑊⁽¹⁾₀₂(−1)𝑧²2𝑎+⋯−𝜎₁𝑧+𝑎₁𝑧+𝑊⁽²⁾₂₀(0)𝑧²2+𝑊⁽²⁾₁₁(0)𝑧𝑧+𝑊⁽²⁾₀₂(0)𝑧²2∗𝑎+⋯₂𝑧+𝑎₂𝑧+𝑊⁽³⁾₂₀(0)𝑧²2+𝑊⁽³⁾₁₁(0)𝑧𝑧+𝑊⁽³⁾₀₂(0)𝑧²2+⋯+𝜏_𝑘𝐽𝑎^∗₂𝜎𝑎₁𝑧+𝑎₁𝑧+𝑊⁽²⁾₂₀(0)𝑧²2+𝑊⁽²⁾₁₁(0)𝑧𝑧+𝑊⁽²⁾₀₂(0)𝑧²2∗𝑎+⋯₂𝑧+𝑎₂𝑧+𝑊⁽³⁾₂₀(0)𝑧²2+𝑊⁽³⁾₁₁(0)𝑧𝑧+𝑊⁽³⁾₀₂(0)𝑧²2+⋯. (3.22) h (3.20), we have 𝑔₂₀=2𝜏_𝑘𝐽−𝑟₂𝑒^{−𝑖𝜔0𝜏𝑘}−𝑝𝑎₁𝑒^{−𝑖𝜔0𝜏𝑘}+𝑎^∗₁𝑘𝑝𝑎₁𝑒^{−𝑖𝜔0𝜏𝑘}−𝑎^∗₁𝜎𝑎₁𝑎₂+𝑎^∗₂𝜎𝑎₁𝑎₂,𝑔₁₁=𝜏_𝑘𝐽𝑒^{−𝑖𝜔0𝜏𝑘}−𝑟₂−𝑝𝑎₁+𝑎^∗₁𝑘𝑝𝑎₁+𝑒^{𝑖𝜔0𝜏𝑘}−𝑟₂−𝑝𝑎₁+𝑎^∗₁𝑘𝑝𝑎₁−𝑎^∗₁𝜎𝑎₁𝑎₂−𝑎^∗₁𝜎𝑎₁𝑎₂+𝑎^∗₂𝜎𝑎₁𝑎₂+𝑎^∗₂𝜎𝑎₁𝑎₂,𝑔₀₂=2𝜏_𝑘𝐽−𝑟₂𝑒^{𝑖𝜔0𝜏𝑘}−𝑝𝑎₁𝑒^{𝑖𝜔0𝜏𝑘}+𝑎^∗₁𝑘𝑝𝑎₁𝑒^{𝑖𝜔0𝜏𝑘}−𝑎^∗₁𝜎𝑎₁𝑎₂+𝑎^∗₂𝜎𝑎₁𝑎₂,𝑔₂₁=𝜏_𝑘𝐽−𝑟₂2𝑊⁽¹⁾₁₁(−1)+𝑊⁽¹⁾₂₀(−1)+𝑊⁽¹⁾₂₀(0)𝑒^{𝑖𝜔0𝜏𝑘}+2𝑊⁽¹⁾₁₁(0)𝑒^{−𝑖𝜔0𝜏𝑘}−𝑝2𝑊⁽²⁾₁₁(−1)+𝑊⁽²⁾₂₀(−1)+𝑎₁𝑊⁽¹⁾₂₀(0)𝑒^{𝑖𝜔0𝜏𝑘}+2𝑎₁𝑊⁽¹⁾₁₁(0)𝑒^{−𝑖𝜔0𝜏𝑘}+𝑎^∗₁𝑘𝑝2𝑎₁𝑊⁽¹⁾₁₁(−1)+𝑎₁𝑊⁽¹⁾₂₀(−1)+𝑊⁽²⁾₂₀(0)𝑒^{𝑖𝜔0𝜏𝑘}+2𝑊⁽²⁾₁₁(0)𝑒^{−𝑖𝜔0𝜏𝑘}−𝑎^∗₁𝜎2𝑎₁𝑊⁽³⁾₁₁(0)+𝑎₁𝑊⁽³⁾₂₀(0)+𝑎₂𝑊⁽²⁾₂₀(0)+2𝑎₂𝑊⁽²⁾₁₁(0)+𝑎^∗₂𝜎2𝑎₁𝑊⁽³⁾₁₁(0)+𝑎₁𝑊⁽³⁾₂₀(0)+𝑎₂𝑊⁽²⁾₂₀(0)+2𝑎₂𝑊⁽²⁾₁₁(0). (3.24) Since there are 𝑊₂₀(𝜃) and 𝑊₁₁(𝜃) in 𝑔₂₁, we need to determine them.

From (3.7) and (3.15), we have 𝑊^′=𝑢^′_𝑡−𝑧^′𝑞−𝑧^′𝑞=𝐴(0)𝑊−2Re𝑞^∗(0)𝑓₀𝑞(𝜃[)𝐴(0)),𝜃∈−1,0𝑊−2Re𝑞^∗(0)𝑓₀𝑞(0)+𝑓₀(0)(,𝜃=0≜𝐴𝑊+𝐺𝑧,𝑧,𝜃), (3.25) where 𝐺(𝑧,𝑧,𝜃)=𝐺₂₀(𝜃)𝑧²2+𝐺₁₁(𝜃)𝑧𝑧+𝐺₀₂(𝜃)𝑧²2+⋯. (3.26) Substituting the corresponding series into (3.25) and comparing the coefficients, we have 𝐴(0)−2𝑖𝜔₀𝜏_𝑘𝐼𝑊₂₀(𝜃)=−𝐺₂₀(𝜃),𝐴(0)𝑊₁₁(𝜃)=−𝐺₁₁(𝜃). (3.27) From (3.20) and (3.25), we have, for 𝜃∈[−1,0)𝐺(𝑧,𝑧,𝜃)=−𝑞^∗(0)𝑓₀𝑞(𝜃)−𝑞^∗(0)𝑓₀𝑞(𝜃)=−𝑔(𝑧,𝑧)𝑞(𝜃)−𝑔(𝑧,𝑧)𝑞(𝜃). (3.28) Comparing the coefficients with (3.26), we have 𝐺₂₀(𝜃)=−𝑔₂₀𝑞(𝜃)−𝑔₀₂𝑞(𝜃),𝐺₁₁(𝜃)=−𝑔₁₁𝑞(𝜃)−𝑔₁₁𝑞(𝜃). (3.29)(3.30) From (3.27), (3.29) and the definition of 𝐴(0), we have 𝑊^′₂₀(𝜃)=2𝑖𝜔₀𝜏_𝑘𝑊₂₀(𝜃)+𝑔₂₀𝑞(𝜃)+𝑔₀₂𝑞(𝜃). (3.31) Notice that 𝑞(𝜃)=(1,𝑎₁,𝑎₂)^𝑇𝑒^{𝑖𝜔0𝜏𝑘𝜃}, hence 𝑊₂₀(𝜃)=𝑖𝑔₂₀𝜔₀𝜏_𝑘𝑞(0)𝑒^{𝑖𝜔0𝜏𝑘𝜃}+𝑖𝑔₀₂3𝜔₀𝜏_𝑘𝑞(0)𝑒^{−𝑖𝜔0𝜏𝑘𝜃}+𝐸₁𝑒^{2𝑖𝜔0𝜏𝑘𝜃}, (3.32) where 𝐸₁=(𝐸₁⁽¹⁾,𝐸₁⁽²⁾,𝐸₁⁽³⁾)^𝑇∈𝑅³ is a constant vector. Similarly, we obtain 𝑊₁₁(𝜃)=−𝑖𝑔₁₁𝜔₀𝜏_𝑘𝑞(0)𝑒^{𝑖𝜔0𝜏𝑘𝜃}+𝑖𝑔₁₁𝜔₀𝜏_𝑘𝑞(0)𝑒^{−𝑖𝜔0𝜏𝑘𝜃}+𝐸₂, (3.33) where 𝐸₂=(𝐸₂⁽¹⁾,𝐸₂⁽²⁾,𝐸₂⁽³⁾)^𝑇∈𝑅³ is a constant vector. In the following, we will seek the values of 𝐸₁ and 𝐸₂. From the definition of 𝐴(0) and (3.27), we have 𝑑𝜂(⁰₋₁𝜃)𝑊₂₀(𝜃)=2𝑖𝜔₀𝜏_𝑘𝑊₂₀(0)−𝐺₂₀(0),₀₋₁𝑑𝜂(𝜃)𝑊₁₁(𝜃)=−𝐺₁₁(0), (3.34)(3.35) where 𝜂(𝜃)=𝜂(0,𝜃).

By (3.25), when 𝜃=0, we have 𝐺(𝑧,𝑧,0)=2Re𝑞^∗(0)𝑓₀𝑞(0)+𝑓₀=−𝑞^∗(0)𝑓₀𝑞(0)−𝑞^∗(0)𝑓₀𝑞(0)+𝑓₀=−𝑔(𝑧,𝑧)𝑞(0)−𝑔(𝑧,𝑧)𝑞(0)+𝑓₀. (3.36) So we obtain G₂₀(𝜃)𝑧²2+𝐺₁₁(𝜃)𝑧𝑧+𝐺₀₂(𝜃)𝑧²2+⋯=−𝑞(0)𝑔₂₀𝑧²2+𝑔₁₁𝑧𝑧+𝑔₀₂𝑧²2−+⋯𝑞(0)𝑔₂₀𝑧²2+𝑔₁₁𝑧𝑧+𝑔₀₂𝑧²2+⋯+𝑓₀. (3.37) By (3.37), we have 𝐺₂₀(0)=−𝑔₂₀𝑞(0)−𝑔₀₂𝑞(0)+2𝜏_𝑘−𝑟₂𝑒^{−𝑖𝜔0𝜏𝑘}−𝑝𝑎₁𝑒^{−𝑖𝜔0𝜏𝑘}𝑘𝑝𝑎₁𝑒^{−𝑖𝜔0𝜏𝑘}−𝜎𝑎₁𝑎₂𝜎𝑎₁𝑎₂𝐺₁₁,(0)=−𝑔₁₁𝑞(0)−𝑔₁₁𝑞(0)+2𝜏_𝑘−𝑟₂−𝑝Re𝑎₁𝑘𝑝Re𝑎₁−𝜎Re𝑎₁𝑎₂𝜎Re𝑎₁𝑎₂. (3.38)(3.39) Noticing that 𝑖𝜔₀𝜏_𝑘𝐼−⁰₋₁𝑒^{𝑖𝜔0𝜏𝑘𝜃}𝑑𝜂(𝜃)𝑞(0)=0,−𝑖𝜔₀𝜏_𝑘𝐼−⁰₋₁𝑒^{−𝑖𝜔0𝜏𝑘𝜃}𝑑𝜂(𝜃)𝑞(0)=0. (3.40) So, substituting (3.32) and (3.38) into (3.34), we have 2𝑖𝜔₀𝜏_𝑘𝐼−⁰₋₁𝑒^{2𝑖𝜔0𝜏𝑘𝜃}𝑑𝜂(𝜃)𝐸₁=2𝜏_𝑘−𝑟₂𝑒^{−𝑖𝜔0𝜏𝑘}−𝑝𝑎₁𝑒^{−𝑖𝜔0𝜏𝑘}𝑘𝑝𝑎₁𝑒^{−𝑖𝜔0𝜏𝑘}−𝜎𝑎₁𝑎₂𝜎𝑎₁𝑎₂. (3.41) That is 2𝑖𝜔₀+𝑟₂𝑋^∗𝑒^{−𝑖𝜔0𝜏𝑘}𝑝𝑋^∗𝑒^{−𝑖𝜔0𝜏𝑘}0−𝑘𝑝𝑆^∗𝑒^{−𝑖𝜔0𝜏𝑘}2𝑖𝜔₀−𝑘𝑝𝑋^∗+𝑐₁+𝜎𝐼^∗𝜎𝑆^∗−𝛾0−𝜎𝐼^∗2𝑖𝜔₀𝐸₁=2−𝑟₂𝑒^{−𝑖𝜔0𝜏𝑘}−𝑝𝑎₁𝑒^{−𝑖𝜔0𝜏𝑘}𝑘𝑝𝑎₁𝑒^{−𝑖𝜔0𝜏𝑘}−𝜎𝑎₁𝑎₂𝜎𝑎₁𝑎₂. (3.42) Let 𝐿₁=2𝑖𝜔₀+𝑟₂𝑋^∗𝑒^{−𝑖𝜔0𝜏𝑘}𝑝𝑋^∗𝑒^{−𝑖𝜔0𝜏𝑘}0−𝑘𝑝𝑆^∗𝑒^{−𝑖𝜔0𝜏𝑘}2𝑖𝜔₀−𝑘𝑝𝑋^∗+𝑐₁+𝜎𝐼^∗𝜎𝑆^∗−𝛾0−𝜎𝐼^∗2𝑖𝜔₀. (3.43) It follows that 𝐸₁⁽¹⁾=2𝐿₁−𝑟₂𝑒^{−𝑖𝜔0𝜏𝑘}−𝑝𝑎₁𝑒^{−𝑖𝜔0𝜏𝑘}𝑝𝑋^∗𝑒^{−𝑖𝜔0𝜏𝑘}0𝑘𝑝𝑎₁𝑒^{−𝑖𝜔0𝜏𝑘}−𝜎𝑎₁𝑎₂2𝑖𝜔₀−𝑘𝑝𝑋^∗+𝑐₁+𝜎𝐼^∗𝜎𝑆^∗−𝛾𝜎𝑎₁𝑎₂−𝜎𝐼^∗2𝑖𝜔₀,𝐸₁⁽²⁾=2𝐿₁2𝑖𝜔₀+𝑟₂𝑋^∗𝑒^{−𝑖𝜔0𝜏𝑘}−𝑟₂𝑒^{−𝑖𝜔0𝜏𝑘}−𝑝𝑎₁𝑒^{−𝑖𝜔0𝜏𝑘}0−𝑘𝑝𝑆^∗𝑒^{−𝑖𝜔0𝜏𝑘}𝑘𝑝𝑎₁𝑒^{−𝑖𝜔0𝜏𝑘}−𝜎𝑎₁𝑎₂𝜎𝑆^∗−𝛾0𝜎𝑎₁𝑎₂2𝑖𝜔₀,𝐸₁⁽³⁾=2𝐿₁2𝑖𝜔₀+𝑟₂𝑋^∗𝑒^{−𝑖𝜔0𝜏𝑘}𝑝𝑋^∗𝑒^{−𝑖𝜔0𝜏𝑘}−𝑟₂𝑒^{−𝑖𝜔0𝜏𝑘}−𝑝𝑎₁𝑒^{−𝑖𝜔0𝜏𝑘}−𝑘𝑝𝑆^∗𝑒^{−𝑖𝜔0𝜏𝑘}2𝑖𝜔₀−𝑘𝑝𝑋^∗+𝑐₁+𝜎𝐼^∗𝑘𝑝𝑎₁𝑒^{−𝑖𝜔0𝜏𝑘}−𝜎𝑎₁𝑎₂0−𝜎𝐼^∗𝜎𝑎₁𝑎₂. (3.44) Similarly, we can get 𝑟₂𝑋^∗𝑒^{−𝑖𝜔0𝜏𝑘}𝑝𝑋^∗𝑒^{−𝑖𝜔0𝜏𝑘}0−𝑘𝑝𝑆^∗𝑒^{−𝑖𝜔0𝜏𝑘}−𝑘𝑝𝑋^∗+𝑐₁+𝜎𝐼^∗𝜎𝑆^∗−𝛾0−𝜎𝐼^∗0𝐸₂=−𝑟₂−𝑝Re𝑎₁𝑘𝑝Re𝑎₁−𝜎Re𝑎₁𝑎₂𝜎Re𝑎₁𝑎₂, (3.45) and hence 𝐸₂⁽¹⁾=1𝐿₂−𝑟₂−𝑝Re𝑎₁𝑝𝑋^∗𝑒^{−𝑖𝜔0𝜏𝑘}0𝑘𝑝Re𝑎₁−𝜎Re𝑎₁𝑎₂−𝑘𝑝𝑋^∗+𝑐₁+𝜎𝐼^∗𝜎𝑆^∗𝑎−𝛾𝜎Re₁𝑎₂−𝜎𝐼^∗0,𝐸₂⁽²⁾=1𝐿₂𝑟₂𝑋^∗𝑒^{−𝑖𝜔0𝜏𝑘}−𝑟₂−𝑝Re𝑎₁0−𝑘𝑝𝑆^∗𝑒^{−𝑖𝜔0𝜏𝑘}𝑘𝑝Re𝑎₁−𝜎Re𝑎₁𝑎₂𝜎𝑆^∗𝑎−𝛾0𝜎Re₁𝑎₂0,𝐸₂⁽³⁾=1𝐿₂𝑟₂𝑋^∗𝑒^{−𝑖𝜔0𝜏𝑘}𝑝𝑋^∗𝑒^{−𝑖𝜔0𝜏𝑘}−𝑟₂−𝑝Re𝑎₁−𝑘𝑝𝑆^∗𝑒^{−𝑖𝜔0𝜏𝑘}0𝑘𝑝Re𝑎₁−𝜎Re𝑎₁𝑎₂0−𝜎𝐼^∗𝜎Re𝑎₁𝑎₂, (3.46) where 𝐿₂=𝑟₂𝑋^∗𝑒^{−𝑖𝜔0𝜏𝑘}𝑝𝑋^∗𝑒^{−𝑖𝜔0𝜏𝑘}0−𝑘𝑝𝑆^∗𝑒^{−𝑖𝜔0𝜏𝑘}−𝑘𝑝𝑋^∗+𝑐₁+𝜎𝐼^∗𝜎𝑆^∗−𝛾0−𝜎𝐼^∗0=𝑟2𝑋^∗𝜎²𝑆^∗𝐼^∗−𝜎𝐼^∗𝛾𝑒^{−𝑖𝜔0𝜏𝑘}. (3.47) Therefore, we can determine 𝑊₂₀(0) and 𝑊₁₁(0), hence we can obtain 𝑔₂₁. Thus, we can compute these values 𝑐₁(0)=𝑖2𝜔₀𝜏_𝑘𝑔₂₀𝑔₁₁−2𝑔₁₁²−𝑔₀₂²3+𝑔₂₁2,𝜇₂𝑐=−Re₁(0)Re𝜆^′𝜏_𝑘,𝛽₂=2Re𝑐₁(0,𝑇₂)𝑐=−𝐼𝑚₁(0)+𝜇₂𝐼𝑚𝜆^′𝜏_𝑘𝜔₀𝜏_𝑘, (3.48) which determine the qualities of bifurcating periodic solution in the center manifold at the critical values 𝜏_𝑘, so we have the following results.

Theorem 3.1. (i) 𝜇₂ determines the directions of the Hopf bifurcation: if 𝜇₂>0(𝜇₂<0), then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for 𝜏>𝜏_𝑘(𝜏<𝜏_𝑘).
(ii) 𝛽₂ determines the stability of the bifurcating periodic solutions: the bifurcating periodic solutions are stable (unstable) if 𝛽₂<0(𝛽₂>0).
(iii) 𝑇₂ determines the period of the bifurcating periodic solutions: the period increases (decreases) if 𝑇₂>0(𝑇₂<0).

4. Discussion and Numerical Example

In this section,we present some numerical results of system (1.3) at different values of 𝜏. Form Section 3, we may determine the direction of a Hopf bifurcation and the stability of the bifurcating periodic solutions. We consider the following system: 𝑑𝑋(𝑡(𝑡)(𝑡𝑑𝑡=𝑋)[0.9−0.2𝑋−0.4𝑆)]𝑑𝑆(𝑡(𝑡)(𝑡(𝑡)𝑑𝑡=𝑆)[−0.2+0.3𝑋−0.6𝐼)]+0.2𝐼𝑑𝐼(𝑡(𝑡)],𝑑𝑡=𝐼)[0.6𝑆−0.3−0.2 (4.1) which has a positive equilibrium 𝐸⁺=(2.8333,0.8333,1.8056). Form (2.13) and (2.14), we are easy to get at least one positive real root 0.5977. In addition, it is easy to show that [𝑑(Re𝜆)/𝑑𝜏]_{𝜏=𝜏𝑘}=2.4377, the hypothesis of 𝐻₃ holds. Hence, 𝐸⁺ satisfies the condition of Theorem (2.2). When 𝜏=0, the positive equilibrium 𝐸⁺=(2.8333,0.8333,1.8056) is asymptotically stable. According to (2.18), we obtain 𝜏₀=1.124,𝜔₀=0.7731,𝜆^′(𝜏₀)=0.3504−0.1448𝑖. Form the formulae (3.48) in Section 3, it follows that 𝑐₁(0)=−15.6822+2.8655𝑖,𝜇₂=44.7551>0,𝛽₂=−31.3764<0 and 𝑇₂=4.1602>0. Thus, 𝐸⁺ is stable when 𝜏<𝜏₀ as is illustrated by the computer simulations (see Figures 1(a) and 1(b)).

(a)

(b)

Figure 1 

When 𝜏=1.11<𝜏0, the positive equilibrium 𝐸+=(2.8333,0.8333,1.8056) is asymptotically stable. (a) shows the trajectories graphs of the system (4.1) with initial data 𝑋(𝑡)=2,𝑆1(𝑡)=2,𝐼2(𝑡)=2. (b) shows the phase portrait of system (4.1).

When 𝜏 passes through the critical value 𝜏₀. 𝐸⁺ loses its stability and a Hopf bifurcation occurs, that is, a family of periodic solutions bifurcate from 𝐸⁺. Since 𝜇₂>0 and 𝛽₂<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(a) and 2(b).

(a)

(b)

Figure 2 

When 𝜏=1.13>𝜏0, bifurcation periodic solutions form 𝐸+. (a) shows the trajectory graphs of system (4.1) with initial data 𝑥(𝑡)=2,𝑦1(𝑡)=2,𝑦2(𝑡)=2. (b) shows the phase portrait of system (4.1).

Acknowledgment

This work is supported by the National Science Foundation of China (10471040), the National Science Foundation of Shanxi (2009011005-1).