Bifurcation Analysis in a Delayed Diffusive Leslie-Gower Model

Abstract

We investigate a modified delayed Leslie-Gower model under homogeneous Neumann boundary conditions. We give the stability analysis of the equilibria of the model and show the existence of Hopf bifurcation at the positive equilibrium under some conditions. Furthermore, we investigate the stability and direction of bifurcating periodic orbits by using normal form theorem and the center manifold theorem.

1. Introduction

The dynamic relationship between predators and their preys has long been and will continue to be one of dominant themes in both ecology and mathematical ecology due to its universal existence and importance. A major trend in theoretical work on prey-predator dynamics has been to derive more realistic models, trying to keep to maximum the unavoidable increase in complexity of their mathematics [1]. In this optic, recent years, the important Leslie-Gower predator-prey model [2, 3] has been extensively studied in [4–7]. A modified version of Leslie-Gower predator-prey model with Holling-type II functional response takes the form 𝑑𝐻𝑎𝑑𝑡=𝐻₁−𝑐−𝑏𝐻¹𝐻𝑃𝑘₁,+𝐻𝑑𝑃𝑎𝑑𝑡=𝑃₂−𝑐₂𝑃𝑘₂,+𝐻 (1) where 𝐻 and 𝑃 represent prey and predator population densities at time 𝑡, respectively. 𝑎₁,𝑎₂,𝑏,𝑐₁,𝑐₂,𝑘₁, and 𝑘₂ are positive constants. 𝑎₁ is the growth rate of prey 𝐻. 𝑎₂ describes the growth rate of predator 𝑃. 𝑏 measures the strength of competition among individuals of species 𝐻. 𝑐₁ is the maximum value of the per capita reduction of 𝐻 due to 𝑃, and 𝑐₂ is the maximum value of the per capita reduction of 𝑃 due to 𝐻, which is not available in abundance. 𝑘₁ measures the extent to which environment provides protection to prey 𝐻. 𝑘₂ measures the extent to which environment provides protection to the predator 𝑃.

On the other hand, time delay plays an important role in many biological dynamical systems, being particularly relevant in ecology [1]. For some predator-prey systems, the rate of the prey population depends on the predation of predator in the earlier times [8–14]. The results indicated that delay differential equations exhibit much more complicated dynamics than ordinary differential equations since a time delay could cause a stable equilibrium to become unstable and induce bifurcations.

In this paper, we will focus on the complex dynamics of the delay effect in the extended reaction-diffusion model. The reproduction of the individuals is modeled by diffusion with diffusion coefficients 𝐷₁>0 and 𝐷₂>0 for the prey and predator, respectively. This basic model is described by a system of two partial differential equations: 𝜕𝐻𝑎𝜕𝑡=𝐻₁−𝑐−𝑏𝐻¹𝐻𝑃(𝑡−𝜏)𝑘₁+𝐻+𝐷₁Δ𝐻,𝑥∈Ω,𝑡>0,𝜕𝑃𝑎𝜕𝑡=𝑃₂−𝑐₂𝑃𝑘₂+𝐻+𝐷₂Δ𝑃,𝑥∈Ω,𝑡>0,𝜕𝐻=𝜕𝐧𝜕𝑃𝜕𝐧=0,𝑥∈𝜕Ω,𝑡>0, (2) where 𝐻=𝐻(𝑡,𝑥), 𝑃=𝑃(𝑡,𝑥). Δ=𝜕²/𝜕𝑥², Ω is a bounded open domain in ℝ with boundary 𝜕Ω, 𝐧 is the outward unit normal vector on 𝜕Ω, and homogeneous Neumann boundary conditions reflect the situation where the population cannot move across the boundary of the domain. And we incorporate a single discrete delay 𝜏>0 in the negative feedback of the predator's density.

The rest of the paper is organized as follows. In Section 2, we give the stability property of the equilibria of model (1). In Section 3, we mainly analyze the distribution of the roots of the characteristic equation and show the occurrence of Hopf bifurcation at the positive equilibrium of model (2) under some conditions. In Section 4, we investigate the stability and direction of bifurcating periodic orbits by using normal form of theorem and the center manifold theorem, corresponding to theorems we also give some numerical simulations.

2. Equilibria Stability

In this section, we consider the existence and stability of the equilibria of model (1).

It is easy to verify that model (1) always has three boundary equilibria:(i)𝐸₁=(0,0) (extinction of prey and predator), which is a nodal source point;(ii)𝐸₂=(𝑎₁/𝑏,0) (extinction of the predator), which is a saddle point;(iii)𝐸₃=(0,𝑎₂𝑘₂/𝑐₂) (extinction of the prey), which is a stable node when 𝑎₂𝑘₂/𝑐₂>𝑎₁𝑘₁/𝑐₁.

For the positive equilibria, we have 𝑎₁𝑐−𝑏𝐻−¹𝑃𝑘₁+𝐻=0,𝑎₂−𝑐₂𝑃𝑘₂+𝐻=0, (3) which yields 𝑏𝑐₂𝐻²−𝑎₁𝑐₂−𝑎₂𝑐₁−𝑏𝑐₂𝑘₁𝐻+𝑎₂𝑐₁𝑘₂−𝑎₁𝑐₂𝑘₁=0. (4) For simplicity, we define 𝐴=𝑏𝑐₂,𝐵=𝑎₁𝑐₂−𝑎₂𝑐₁−𝑏𝑐₂𝑘₁,𝐶=𝑎₂𝑐₁𝑘₂−𝑎₁𝑐₂𝑘₁, (5) then (4) can be written as 𝐴𝐻²−𝐵𝐻+𝐶=0, (6) which has two roots given by ℎ₊=𝐵+𝐵²−4𝐴𝐶2𝐴,ℎ₋=𝐵−𝐵²−4𝐴𝐶.2𝐴 (7)

Case 1 (𝐶<0, i.e., 𝑘₂<𝑎₁𝑐₂𝑘₁/𝑎₂𝑐₁). Model (1) has a unique positive equilibrium 𝐸^∗=(ℎ^∗,𝑝^∗)=(ℎ₊,(𝑎₂/𝑐₂)(ℎ₊+𝑘₂)).

Case 2 (𝐶=0, i.e., 𝑘₂=𝑎₁𝑐₂𝑘₁/𝑎₂𝑐₁). (i)If 𝐵>0, that is, 𝑘₁<(𝑎₁𝑐₂−𝑎₂𝑐₁)/𝑏𝑐₂, model (1) has a unique positive equilibrium 𝐸=(ℎ,𝑝)=(ℎ₊,(𝑎₂/𝑐₂)(ℎ₊+𝑘₂));(ii)If 𝐵≤0, that is, 𝑘₁≥(𝑎₁𝑐₂−𝑎₂𝑐₁)/𝑏𝑐₂, (4) has no positive root; hence model (1) has no positive equilibrium.

Case 3 (𝐶>0, i.e., 𝑘₂>𝑎₁𝑐₂𝑘₁/𝑎₂𝑐₁). (i) Suppose that 𝐵>0, that is, 𝑘₁<(𝑎₁𝑐₂−𝑎₂𝑐₁)/𝑏𝑐₂, then(i1) if 𝐵²−4𝐴𝐶>0, that is, 𝑘₂<(𝑎₂𝑐₁−𝑎₁𝑐₂−𝑏𝑐₂𝑘₁)²/4𝑎₂𝑏𝑐₁𝑐₂+𝑘₁, model (1) has two positive equilibria 𝐸₊=(ℎ₊,𝑝₊)=(ℎ₊,(𝑎₂/𝑐₂)(ℎ₊+𝑘₂)) and 𝐸₋=(ℎ₋,𝑝₋)=(ℎ₋,(𝑎₂/𝑐₂)(ℎ₋+𝑘₂)); (i2) if 𝐵²−4𝐴𝐶=0, that is, 𝑘₂=(𝑎₂𝑐₁−𝑎₁𝑐₂−𝑏𝑐₂𝑘₁)²/4𝑎₂𝑏𝑐₁𝑐₂+𝑘₁, (4) has a unique positive root of multiplicity 2 given by ℎ_𝑒=𝐵/2𝐴=ℎ₊=ℎ₋, then model (1) has a unique positive equilibrium 𝐸_𝑒=(ℎ_𝑒,𝑝_𝑒)=(ℎ_𝑒,(𝑎₂/𝑐₂)(ℎ_𝑒+𝑘₂)); (i3) if 𝐵²−4𝐴𝐶=0, that is, 𝑘₂>(𝑎₂𝑐₁−𝑎₁𝑐₂−𝑏𝑐₂𝑘₁)²/4𝑎₂𝑏𝑐₁𝑐₂+𝑘₁, (4) has no positive root; hence model (1) has no positive equilibrium; (ii) if 𝐵<0, that is, 𝑘₁>(𝑎₁𝑐₂−𝑎₂𝑐₁)/𝑏𝑐₂, (4) has no positive root; hence model (1) has no positive equilibrium.

We show the bifurcation diagram to display the distribute of the positive roots; in Figure 1, the whole region has been divided into six parts; the number indicates the number of positive equilibria.

Figure 1 

The bifurcation diagram displays the distribute of the positive roots; the number indicates the number of positive equilibria.

Next, we analyze the stability of these positive equilibria. Let 𝐸=(ℎ,𝑝) be arbitrary positive equilibrium, and the Jacobian matrix for 𝐸=(ℎ,𝑝) is given by 𝐽(𝐸)=𝑐−𝑏ℎ+¹ℎ𝑝ℎ+𝑘₁²−𝑐₁ℎℎ+𝑘₁𝑎²₂𝑐₂−𝑎₂. (8) The corresponding characteristic equation is 𝜆²−tr(𝐽(𝐸))𝜆+det(𝐽(𝐸))=0, (9) where (𝐽(𝐸=𝑐tr))¹ℎ𝑝ℎ+𝑘₁²−𝑎₂=ℎ𝑎−𝑏ℎ₁−𝑎₂−2𝑏ℎ−𝑎₂+𝑏ℎ𝑘₁ℎ+𝑘₁,(𝐽(𝐸=𝑎det))²ℎ𝑎₂𝑐₁−𝑎₁𝑐₂+𝑏𝑐₂𝑘₁+2𝑏𝑐₂ℎ𝑐₂ℎ+𝑘₁. (10) And the sign of det(𝐽(𝐸)) is determined by 𝐹₁(ℎ)=𝑎₂𝑐₁−𝑎₁𝑐₂+𝑏𝑐₂𝑘₁+2𝑏𝑐₂ℎ=2𝐴ℎ−𝐵. (11) Thus, 𝐹₁ℎ^∗=2𝐴ℎ^∗−𝐵=𝐵²𝐹−4𝐴𝐶>0,₁ℎ=2𝐴𝐹ℎ−𝐵=𝐵>0,₁ℎ₊=2𝐴ℎ₊−𝐵=𝐵²𝐹−4𝐴𝐶>0;₁ℎ₋=2𝐴ℎ₋−𝐵=−𝐵²𝐹−4𝐴𝐶<0,₁ℎ_𝑒=2𝐴ℎ_𝑒−𝐵=0. (12) Hence, det(𝐽(𝐸^∗))>0, det(𝐽(𝐸))>0, det(𝐽(𝐸₊))>0, det(𝐽(𝐸₋))<0, and det(𝐽(𝐸_𝑒))=0. Obviously, the positive equilibrium 𝐸₋ is a saddle point.

In the following, we study the stability of other positive equilibria. The sign of tr(𝐽(𝐸)) is determined by 𝐺₁(ℎ)=ℎ𝑎₁−𝑎₂−2𝑏ℎ−𝑎₂+𝑏ℎ𝑘₁. (13) Then we can get 𝐺₁ℎ^∗=ℎ^∗𝑎₂𝑐₁−𝑐₂𝑐₂−ℎ^∗𝐵²−4𝐴𝐶𝑐₂−𝑎₂𝑘₁,𝐺₁ℎ=ℎ𝑎₂𝑐₁−𝑐₂𝑐₂−ℎ𝐵𝑐₂−𝑎₂𝑘₁,𝐺₁ℎ₊=ℎ₊𝑎₂𝑐₁−𝑐₂𝑐₂−ℎ₊𝐵²−4𝐴𝐶𝑐₂−𝑎₂𝑘₁,𝐺₁ℎ_𝑒=ℎ_𝑒𝑎₂𝑐₁−𝑐₂𝑐₂−𝑎₂𝑘₁. (14) Hence, if 𝑐₁≤𝑐₂, then tr(𝐽(𝐸^∗))<0, tr(𝐽(𝐸))<0, and tr(𝐽(𝐸₊))<0 are true. Summarizing the above, we can obtain the following theorem.

Theorem 1. For model (1),(i)if 𝑘₂<𝑎₁𝑐₂𝑘₁/𝑎₂𝑐₁ holds, the unique positive equilibrium 𝐸^∗ is locally asymptotically stable for 𝑐₁≤𝑐₂;(ii)if 𝑘₂=𝑎₁𝑐₂𝑘₁/𝑎₂𝑐₁ and 𝑘₁<(𝑎₁𝑐₂−𝑎₂𝑐₁)/𝑏𝑐₂ hold, the unique positive equilibrium 𝐸 is locally asymptotically stable for 𝑐₁≤𝑐₂;(iii)if 𝑎₁𝑐₂𝑘₁/𝑎₂𝑐₁<𝑘₂<(𝑎₂𝑐₁−𝑎₁𝑐₂−𝑏𝑐₂𝑘₁)²/4𝑎₂𝑏𝑐₁𝑐₂+𝑘₁ and 𝑘₁<(𝑎₁𝑐₂−𝑎₂𝑐₁)/𝑏𝑐₂ hold, model (1) has two positive equilibria, the positive equilibrium 𝐸₊ is locally asymptotically stable for 𝑐₁≤𝑐₂, and 𝐸₋ is a saddle point.

Figure 2 shows the dynamics of model (1). In this case, 𝐸₁ is a nodal source point; 𝐸₂ is a saddle point; 𝐸₃ is a nodal sink point, which is locally asymptotically stable; 𝐸₊ is locally asymptotically stable; 𝐸₋ is a saddle point. There exists a separatrix curve determined by the stable manifold of 𝐸₋, which divides the behavior of trajectories; that is, the stable manifold of saddle 𝐸₋ splits the feasible region into two parts such that orbits initiating inside tend to the positive equilibrium 𝐸₊, while orbits initiating outside tend to 𝐸₃ except for the stable manifolds of 𝐸₋. This means that, in this situation, the trajectories of the model can have different behavior strongly depending on the initial conditions.

Figure 2 

Phase portraits of model (1) with the parameters 𝑎1=2, 𝑎2=0.5, 𝑏=1, 𝑐1=𝑐2=2, 𝑘1=0.4, and 𝑘2=2.1. In this case, 𝐸1=(0,0) is a nodal source point; 𝐸2=(2,0) is a saddle point; 𝐸3=(0,0.525) is a nodal sink point, which is locally asymptotically stable; 𝐸+=(0.7791,0.7197) is a spiral sink, which is locally asymptotically stable; and 𝐸−=(0.3208,0.6052) is a saddle point. There exists a separatrix curve determined by the stable manifold of 𝐸−. The dashed curve is the 𝐻-nullcline 𝐻(𝑎1−𝑏𝐻)−𝑐1𝐻𝑃/(𝑘1+𝐻), and the dotted curve is the 𝑃-nullcline 𝑃(𝑎2−𝑐2𝑃/(𝑘2+𝐻)).

Theorem 2. For model (1), if the unique positive equilibrium 𝐸_𝑒=(ℎ_𝑒,(𝑎₂/𝑐₂)(ℎ_𝑒+𝑘₂)) exists, and 𝑘₁=(𝑎₁𝑐₂−𝑎₂𝑐₁)(𝑐₁−𝑐₂)/𝑏𝑐₂(𝑐₁+𝑐₂), then it is a cusp of codimension 2.

Proof. The Jacobian matrix at 𝐸_𝑒=(ℎ_𝑒,𝑝_𝑒) is 𝐽_𝐸𝑒=−𝑎₂𝑐₁𝑎₂𝑐₁−𝑎₁𝑐₂+𝑏𝑐₂𝑘₁𝑐₂𝑏𝑐₂𝑘₁−𝑎₂𝑐₁+𝑎₁𝑐₂𝑐₁𝑎₂𝑐₁−𝑎₁𝑐₂+𝑏𝑐₂𝑘₁𝑏𝑐₂𝑘₁−𝑎₂𝑐₁+𝑎₁𝑐₂𝑎²₂𝑐₂−𝑎₂, (15) we have know that det(𝐽(𝐸_𝑒))=0. Moreover, tr(𝐽(𝐸_𝑒))=0, if and only if 𝑘₁=𝑎₁𝑐₂−𝑎₂𝑐₁𝑐₁−𝑐₂𝑏𝑐₂𝑐₁+𝑐₂. (16) Then 𝐽_𝐸𝑒=𝑎₂−𝑐₂𝑎²₂𝑐₂−𝑎₂=𝑎₂𝑐1−²𝑎₂𝑎₂𝑐₂−1, (17) and the associate Jordan matrix is ̂0𝑐𝐽=¹𝑎₂𝑐₁−𝑎₁𝑐₂+𝑏𝑐₂𝑘₁𝑏𝑐₂𝑘₁−𝑎₂𝑐₁+𝑎₁𝑐₂00. (18)
Hence, following [15, 16], we know that the unique positive equilibrium 𝐸_𝑒 is a cusp of codimension 2.

3. Stability and Hopf Bifurcation Analysis in Delayed Reaction-Diffusion Model (2)

According to the previous section, for model (1), we know that 𝐸₁, 𝐸₂, and 𝐸₃ are unstable and 𝐸₋ is a saddle point, and note that a solution of the model (1) is also a solution of the model (2), so they are also unstable for model (2). In the following, we will focus on the dynamics of the positive equilibria of model (2). As an example, we only give the proof of the unique positive equilibrium 𝐸^∗ of model (2).

Introducing small perturbations 𝐻=𝐻−ℎ^∗, and 𝑃=𝑃−𝑝^∗ and dropping the hats for simplicity of notation, then we have 𝜕𝐻=𝜕𝑡𝐻+ℎ^∗𝑎₁−𝑏𝐻+ℎ^∗−𝑐₁𝐻+ℎ^∗𝑃(𝑡−𝜏)+𝑝^∗𝑘₁+𝐻+ℎ^∗+𝐷₁Δ𝐻,𝑥∈Ω,𝑡>0,𝜕𝑃=𝜕𝑡𝑃+𝑝^∗𝑎₂−𝑐₂𝑃+𝑝^∗𝑘₂+𝐻+ℎ^∗+𝐷₂Δ𝑃,𝑥∈Ω,𝑡>0,𝜕𝐻=𝜕𝐧𝜕𝑃𝜕𝐧=0,𝑥∈𝜕Ω,𝑡>0. (19)

Denote 𝑋=𝐻,𝑃∈𝑊^2,2(Ω)∶𝜕𝐻=𝜕𝐧𝜕𝑃𝜕𝐧=0,𝑥∈𝜕Ω. (20) In the abstract space 𝐶([−𝜏,0],𝑋), model (19) can be regarded as the following abstract functional differential equation. 𝜕𝑈𝑈𝜕𝑡=𝐷Δ𝑈+𝐿_𝑡+𝐹𝑈_𝑡, (21) where 𝑈=(𝐻,𝑃)^𝑇, 𝑈_𝑡=𝑈(𝑡+𝜃),𝜃∈[−𝜏,0], 𝐷=^𝐷1_𝐷2, dom(Δ)⊂𝑋, and 𝐿∶𝐶([−𝜏,0],𝑋)↦𝑋, 𝐹∶𝐶([−𝜏,0],𝑋)↦𝑋 are given by 𝐿𝜙=−𝑏ℎ^∗+𝑐₁ℎ^∗𝑝^∗ℎ^∗+𝑘₁²𝜙₁(0)+−𝑐₁ℎ^∗ℎ^∗+𝑘₁𝜙₂()𝑎−𝜏₂₂𝑐₂𝜙₁(0)−𝑎₂𝜙₂(0),𝐹𝜙=𝑐₁𝑘₁𝑝^∗𝑘₁+ℎ^∗3−𝑏𝜙²₁(0)−𝑐₁𝑘₁𝑘₁+ℎ^∗2𝜙₁(0)𝜙₂()−𝑐−𝜏₂𝑝^∗2𝑘₂+ℎ^∗3𝜙²₁(0)+2𝑐₂𝑝^∗𝑘₂+ℎ^∗2𝜙₁(0)𝜙₂(0)−𝑐₂𝑘₂+ℎ^∗𝜙²₂(0), (22) here 𝜙=(𝜙₁,𝜙₂)^𝑇=𝑈_𝑡∈𝐶([−𝜏,0],𝑋). Then the linearization of model (19) near (ℎ^∗,𝑝^∗) is 𝜕𝑈𝑈𝜕𝑡=𝐷Δ𝑈+𝐿_𝑡. (23) Following [17], we obtain that the characteristic equation for liner model (23) is 𝜆𝑦−𝐷Δ𝑦−𝐿𝑒^𝜆⋅𝑦=0,𝑦∈dom(Δ)⊂𝑋,𝑦≠0. (24) It is well known that the eigenvalue problem −Δ𝜓=𝜇𝜓,𝑥∈Ω,𝜕𝜓𝜕𝐧=0,𝑥∈𝜕Ω, (25) has eigenvalues 𝜇_𝑛={𝜇^𝑖_𝑛=−𝐷_𝑖𝑛²,𝑖=1,2,and𝑛=0,1,2,…}, with the corresponding eigenfunctions 𝜓_𝑛(𝑥)=cos𝑛𝑥   (𝑛∈ℕ={0,1,2,…}).

Substituting 𝑦=∑^∞_𝑛=0cos𝑛𝑥(𝑦_1𝑛,𝑦_2𝑛)^𝑇 into characteristic equation (24), we obtain −𝑏ℎ^∗+𝑐₁ℎ^∗𝑝^∗ℎ^∗+𝑘₁²−𝐷₁𝑛²−𝑐₁ℎ^∗ℎ^∗+𝑘₁𝑒^−𝜆𝜏𝑎²₂𝑐₂−𝑎₂−𝐷₂𝑛²𝑦_1𝑛𝑦_2𝑛𝑦=𝜆^1𝑛𝑦_2𝑛. (26) Therefore the characteristic equation (24) is equivalent to 𝜆²+𝐴_𝑛𝜆+𝐵_𝑛𝑒^−𝜆𝜏+𝐶_𝑛=0,𝑛∈ℕ, (27) where 𝐴_𝑛=𝐷₁+𝐷₂𝑛²+𝑎₂+𝑏ℎ^∗−𝑐₁ℎ^∗𝑝^∗ℎ^∗+𝑘₁²,𝐵_𝑛=𝑎²₂𝑐₁ℎ^∗𝑐₂ℎ^∗+𝑘₁,𝐶_𝑛=𝐷₁𝐷₂𝑛⁴+𝐷₁+𝐷₂+𝑎₂+𝑏ℎ^∗−𝑐₁ℎ^∗𝑝^∗ℎ^∗+𝑘₁²𝑛²+𝑎₂𝑏ℎ^∗−𝑎₂𝑐₁ℎ^∗𝑝^∗ℎ^∗+𝑘₁². (28)

The stability of the positive equilibrium 𝐸^∗ can be determined by the distribution of the roots of (27); that is, the equilibrium 𝐸^∗ is locally asymptotically stable if all the roots of (27) have negative real parts. From the result of [18], the sum of the multiplicities of the roots of (27) in the open right half plane changes only if a root appears on or crosses the imaginary axis. It can be verified that 𝜆=0 is not a root of (27) for 𝑛∈ℕ.

Theorem 3. If 𝑘₂<min{𝑎₁𝑐₂𝑘₁/𝑎₂𝑐₁,𝑏𝑐₂(ℎ^∗+𝑘₁)²/𝑎₂𝑐₁−2ℎ^∗−𝑘₁} holds, then the unique positive equilibrium 𝐸^∗ of model (2) is locally asymptotically stable.

Proof. Let ±𝑖𝜔  (𝜔>0) be a pair of roots of (27); substituting 𝑖𝜔 into (27), then we have −𝜔²+𝑖𝐴_𝑛𝜔+𝐵_𝑛𝑒^{−𝑖𝜔𝜏}+𝐶_𝑛=0. (29) Separating the real part from image part, we have −𝜔²+𝐵_𝑛cos𝜔𝜏+𝐶_𝑛𝐴=0,_𝑛𝜔−𝐵_𝑛sin𝜔𝜏=0, (30) then 𝜔⁴+𝐴²_𝑛−2𝐶_𝑛𝜔²+𝐶²_𝑛−𝐵²_𝑛=0, (31) where 𝐴²_𝑛−2𝐶_𝑛=𝑏ℎ^∗−𝑐₁ℎ^∗𝑝^∗ℎ^∗+𝑘₁²+𝐷₁𝑛²²+𝑎₂+𝐷₂𝑛²²𝐶>0,_2𝑛−𝐵²_𝑛=𝐷₁𝐷₂𝑛⁴+𝐷₁+𝐷₂+𝑎₂+𝑏ℎ^∗−𝑐₁ℎ^∗𝑝^∗ℎ^∗+𝑘₁²𝑛²−𝑎₂𝑐₁ℎ^∗𝑝^∗ℎ^∗+𝑘₁²−𝑎⁴₂𝑐²₁ℎ^∗2𝑐²₂ℎ^∗+𝑘₁²+𝑎₂𝑏ℎ^∗. (32) Obviously, if 𝑘₂<𝑏𝑐₂(ℎ^∗+𝑘₁)²/𝑎₂𝑐₁−2ℎ^∗−𝑘₁, 𝐶²_𝑛−𝐵²_𝑛>0 is true. Thus (31) has no positive roots for all 𝑛∈ℕ. Hence, all the roots of (27) have negative real part. This completes the proof.

If there exists an integer 𝑛₀∈ℕ such that for 0≤𝑛≤𝑛₀, 𝐶²_𝑛−𝐵²_𝑛<0, then (31) has a unique positive real root 𝜔^𝑛₀=−𝐴²_𝑛−2𝐶_𝑛+𝐴²_𝑛−2𝐶_𝑛²−4𝐶²_𝑛−𝐵²_𝑛2, (33) and (27) has a pair of pure imaginary roots ±𝑖𝜔^𝑛₀, and 𝜏^𝑗_𝑛=𝜏⁰_𝑛+2𝜋𝑗𝜔_𝑛,𝑗=0,1,2,…,0≤𝑛≤𝑛₀,𝑛₀∈ℕ, (34) where 𝜏⁰_𝑛=arccos((𝜔^𝑛₀²−𝐶_𝑛)/𝐵_𝑛)/𝜔^𝑛₀.

Let 𝜆(𝜏)=𝛾(𝜏)+𝑖𝜔(𝜏) be the root of (24), where 𝛾(𝜏^𝑗_𝑛)=0 and 𝜔(𝜏^𝑗_𝑛)=𝜔^𝑛₀ when 𝜏 is close to 𝜏^𝑗_𝑛. Then we have the following transversality condition.

Lemma 4. For 0≤𝑛≤𝑛₀(𝑛₀∈ℕ), if(H1)𝐷₁𝐷₂𝑛⁴+(𝐷₁+𝐷₂+𝑎₂+𝑏ℎ^∗−𝑐₁ℎ^∗𝑝^∗/(ℎ^∗+𝑘₁)²)𝑛²<𝑎₂𝑐₁ℎ^∗𝑝^∗/(ℎ^∗+𝑘₁)² +𝑎⁴₂𝑐²₁(ℎ^∗)²/𝑐²₂(ℎ^∗+𝑘₁)²−𝑎₂𝑏ℎ^∗ holds, then (𝑑𝛾/𝑑𝜏)|_𝜏=𝜏^𝑗_𝑛>0 for 𝑗=0,1,2,….

From this transversality condition, we know that when 𝜏 passes through these critical values 𝜏^𝑗_𝑛, the sum of the multiplicities of the roots of (27) in the open right half plane will increases at least two.

Summarizing the above results, we can obtain the following theorem.

Theorem 5. For 0≤𝑛≤𝑛₀(𝑛₀∈ℕ) if (H1) holds, the following statements are true:(i) if 𝜏∈[0,𝜏⁰_𝑛), then the equilibrium point 𝐸^∗ is locally asymptotically stable;(ii) if 𝜏>𝜏⁰_𝑛, then the equilibrium 𝐸^∗ is unstable;(iii)𝜏=𝜏^𝑗_𝑛(𝑗=0,1,2,…) are Hopf bifurcation values of model (2).

4. Direction and Stability of Spatial Hopf Bifurcation

In the previous section, we have obtained the conditions under which model (2) undergoes a Hopf bifurcation at the the equilibrium point 𝐸^∗ when 𝜏 crosses though the critical value 𝜏^𝑗_𝑛(0≤𝑛≤𝑛₀,𝑛₀∈ℕ,𝑗=0,1,2,…). In this section, we will study the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions by employing the center manifold theorem and normal form method [17, 19] for partial differential equations with delay. Then we compute the direction and stability of the Hopf bifurcation when 𝜏=𝜏^𝑗_𝑛 for fixed 𝑗∈{0,1,2,…}.

Define 𝑐𝜇=−Re₁(0)Re𝜆^′(𝜏),𝛽₂=2Re𝑐₁(0),𝑇₂𝑐=−Im₁(0)+𝜇₂Im𝜆^′(𝜏)𝜔^𝑛₀𝜏, (35) where 𝑐₁(0) is defined in the appendix. Then we can get the following theorem.

Theorem 6. For model (19), if (H1) holds, we have the following:(i)𝜇 determines the direction of the Hopf bifurcation: if 𝜇>0 (𝜇<0), then the bifurcating periodic solutions exist for 𝜏>𝜏 (𝜏<𝜏);(ii)𝛽₂ determines the stability of bifurcating periodic solutions: the bifurcating periodic solutions are orbitally asymptotically stable (unstable) if 𝛽₂<0 (𝛽₂>0);(iii)𝑇₂ determines the period of the bifurcating periodic solutions: the period increases (decreases) if 𝑇₂>0 (𝑇₂<0).

The proof is deferred to the appendix.

5. Conclusions and Remarks

In this paper, we have considered a modified version of Leslie-Gower model with Holling-type II functional and delayed diffusive predator-prey model under homogeneous Neumann boundary conditions. The value of this study lies in two folds. First, it presents local asymptotic stability of the equilibria of model with and without delay and the existence of Hopf bifurcation, which indicates that the dynamics induced by time delay are rich and complex. Second, it give the analysis of direction and stability of spatial Hopf bifurcation, from which one can find that small sufficiently delays cannot change the stability of the positive equilibrium and large delays cannot only destabilize the positive equilibrium but also induce oscillatory behaviors near the positive equilibrium.

In the following, we give some numerical examples to illustrate the dynamical behaviors of model (2). In Figure 3, 𝜏=2<𝜏⁰₀=3.435144529, the unique positive equilibrium 𝐸^∗=(2,1.5) remains the stability; the population of the predator and the prey will tend to a steady state. However, in Figure 4, 𝜏=4>𝜏⁰₀, the positive equilibrium 𝐸^∗ losses its stability and Hopf bifurcation occurs, which means that a family of stable periodic solutions bifurcate from 𝐸^∗ and the system goes into oscillations; it means that the predator coexists with the prey with oscillatory behaviors.

(a)

(b)

Figure 3 

The solutions of model (2) tends to aperiodically oscillatory orbit. The parameters are taken as 𝐷1=0.01, 𝐷2=0.1, 𝑎1=2, 𝑎2=1, 𝑏=0.5, 𝑐1=𝑐2=2, 𝑘1=𝑘2=1, and the initial values 𝐻(𝑥,𝑡)=2+0.1𝑡sin𝑥, 𝑃(𝑥,𝑡)=1.5+0.1𝑡cos𝑥. 𝑡∈[−2,0], 𝑥∈[0,𝜋]. In this case, 𝜏=4>𝜏00.

(a)

(b)

Figure 4 

The solutions of model (2) tends to aperiodically oscillatory orbit. The parameters are taken as 𝐷1=0.01, 𝐷2=0.1, 𝑎1=2, 𝑎2=1, 𝑏=0.5, 𝑐1=𝑐2=2, 𝑘1=𝑘2=1, and the initial values 𝐻(𝑥,𝑡)=2+0.1𝑡sin𝑥, 𝑃(𝑥,𝑡)=1.5+0.1𝑡cos𝑥. 𝑡∈[−4,0], 𝑥∈[0,𝜋]. In this case, 𝜏=4>𝜏00.

Our results show that time-delay can make a stable equilibrium to become unstable and induce Hopf bifurcation and the system goes into oscillations; that’s to say, the dynamical behaviors of the delay reaction-diffusion equations are much more complex and rich than reaction-diffusion equations.

Appendix

A. The Proof of Theorem 6

Setting 𝜏=𝜏+𝛼, then 𝛼=0 is the Hopf bifurcation of model (19). Let 𝐻(𝑡,𝑥)=𝐻(𝜏𝑡,𝑥), 𝑃(𝑡,𝑥)=𝑃(𝜏𝑡,𝑥), and drop the tilde for the sake of simplicity, then (19) can be transformed into 𝜕𝐻𝜕𝑡=𝜏𝐻+ℎ^∗𝑎₁−𝑏𝐻+ℎ^∗−𝑐₁𝐻+ℎ^∗𝑃(𝑡−1)+𝑝^∗𝑘₁+𝐻+ℎ^∗+𝐷₁Δ𝐻,𝑥∈Ω,𝑡>0,𝜕𝑃𝜕𝑡=𝜏𝑃+𝑝^∗𝑎₂−𝑐₂𝑃+𝑝^∗𝑘₂+𝐻+ℎ^∗+𝐷₂Δ𝑃,𝑥∈Ω,𝑡>0,𝜕𝐻=𝜕𝐧𝜕𝑃𝜕𝐧=0,𝑥∈𝜕Ω,𝑡>0. (A.1)

And the abstract functional differential equation can also be written in the form 𝑑𝑈=𝑈𝑑𝑡𝜏𝐷Δ𝑈+𝜏𝐿_𝑡+𝐺𝑈_𝑡,𝛼, (A.2) where 𝐿𝜙=−𝑏ℎ^∗+𝑐₁ℎ^∗𝑝^∗ℎ^∗+𝑘₁²𝜙₁(0)−𝑐₁ℎ^∗ℎ^∗+𝑘₁𝜙₂()𝑎−1₂₂𝑐₂𝜙₁(0)−𝑎₂𝜙₂(0),𝐺𝜙+()𝐹𝜙,𝐹𝜙=𝑐𝜙,𝛼=𝛼𝐷Δ𝜙+𝛼𝐿𝜏+𝛼₁𝑘₁𝑝^∗𝑘₁+ℎ^∗3−𝑏𝜙²₁(0)−𝑐₁𝑘₁𝑘₁+ℎ^∗2𝜙₁(0)𝜙₂()−𝑐−1₂𝑝^∗2𝑘₂+ℎ^∗3𝜙²₁(0)+2𝑐₂𝑝^∗𝑘₂+ℎ^∗2𝜙₁(0)𝜙₂(0)−𝑐₂𝑘₂+ℎ^∗𝜙²₂(0), (A.3) for 𝜙=(𝜙₁,𝜙₂)=𝑈_𝑡∈𝐶([−1,0],𝑋).

From Section 2, we know that ±𝑖𝜔^𝑛₀𝜏 are a pair of simple purely imaginary eigenvalues of the liner system 𝑑𝑈=𝑈𝑑𝑡𝜏𝐿_𝑡+𝜏𝐷Δ𝑈, (A.4) where 𝑈(𝑡)=(𝐻(𝑡),𝑃(𝑡))^𝑇∈𝑅², and 𝑈_𝑡(𝜃) is defined by 𝑈_𝑡(𝜃)=𝑈(𝑡+𝜃).

By using the Riesz representation theorem [20], we have a function 𝜂(𝜃,𝛼) of bounded variation for 𝜃∈[−1,0] such that 𝐿𝜙+𝐷Δ𝜙(0)=⁰₋₁()𝜙(𝜃),(𝜃)[]𝑑𝜂𝜃,𝛼for𝜙∈𝐶−1,0,𝑅². (A.5) Because in this paper, we discuss the existence of the Hopf  bifurcation when 𝜏=𝜏⁰₀, that is 𝑛=0; here we choose 𝜂(𝜃,𝛼)=−𝑏ℎ^∗+𝑐₁ℎ^∗𝑝^∗ℎ^∗+𝑘₁²0𝑎²₂𝑐₂−𝑎₂𝛿(𝜃)+0−𝑐₁ℎ^∗ℎ^∗+𝑘₁𝛿(),00𝜃+1 (A.6) where 𝛿 is the Dirac delta function. For 𝜙(𝜃)∈𝐶¹([−1,0],𝑅²), define 𝐴(𝛼) as 𝐴(𝛼)(𝜃)𝜙=𝑑𝜙[),∫𝑑𝜃,𝜃∈−1,0₀₋₁𝑑𝜂(𝜃,𝛼)𝜙(𝜃),𝜃=0, (A.7) and for 𝜓=(𝜓₁,𝜓₂)∈𝐶¹([0,1],(𝑅²)^∗), define 𝐴^∗𝜓(𝑠)=−𝑑𝜓(𝑠)[),∫𝑑𝑠,𝑠∈−1,0₀₋₁𝜓(−𝜉)𝑑𝜂(𝜃,0),𝑠=0. (A.8) Then 𝐴(0) and 𝐴^∗ are adjoint operators under the bilinear form 𝜓(𝑠)(𝜃)=,𝜙𝜓(0)𝜙(0)−⁰₋₁^𝜃₀𝜓(𝜉−𝜃)𝑑𝜂(𝜃,𝛼)𝜙(𝜉)𝑑𝜉. (A.9)

We know that ±𝑖𝜔^𝑛₀ are eigenvalues of 𝐴(0). Since 𝐴(0) and 𝐴^∗(0) are two adjoint operators, then ±𝑖𝜔^𝑛₀ are also eigenvalue of 𝐴^∗; we shall first try to obtain eigenvector of 𝐴(0) and 𝐴^∗ corresponding to the eigenvalues 𝑖𝜔^𝑛₀ and −𝑖𝜔^𝑛₀, respectively. Let 𝑞(𝜃)=(1,𝜌)^𝑇𝑒^{𝑖𝜔𝑛0𝜃} be the eigenvector of 𝐴(0) corresponding to the eigenvalue 𝑖𝜔^𝑛₀. Then by definition of eigenvector we have 𝐴(0)𝑞(𝜃)=𝑞(𝜃)𝑖𝜔^𝑛₀. Therefore, from (A.5), (A.6), and definition of 𝐴(0) we get −𝑏ℎ^∗+𝑐₁ℎ^∗𝑝^∗ℎ^∗+𝑘₁²−𝑖𝜔^𝑛₀−𝑐₁ℎ^∗ℎ^∗+𝑘₁𝑒^{−𝑖𝜔𝑛0}𝜏𝑎²₂𝑐₂−𝑎₂−𝑖𝜔^𝑛₀1𝜌=00, (A.10) we choose 𝜌=𝑎²₂/𝑐₂(𝑎₂+𝑖𝜔^𝑛₀), and then we get 𝑞(𝜃). On the other hand, 𝑞^∗(𝑠)=𝑀(1,𝛾)𝑒^{𝑖𝜔𝑛0𝑠} is the eigenvector of 𝐴^∗ corresponding to the eigenvalue −𝑖𝜔^𝑛₀. From the definition of 𝐴^∗, we have −𝑏ℎ^∗+𝑐₁ℎ^∗𝑝^∗ℎ^∗+𝑘₁²+𝑖𝜔^𝑛₀𝑎²₂𝑐₂−𝑐₁ℎ^∗ℎ^∗+𝑘₁𝑒^{−𝑖𝜔𝑛0}𝜏−𝑎₂+𝑖𝜔^𝑛₀𝑀=0𝑀𝛾0, (A.11) where 𝛾=𝑏ℎ^∗−𝑐₁ℎ^∗𝑝^∗ℎ^∗+𝑘₁²−𝑖𝜔^𝑛₀𝑐₂𝑎²₂. (A.12) We also assume that (𝑞^∗(𝑠),𝑞(𝜃))=1. To obtain the value of 𝑀, from (A.9) we have 𝑞^∗(𝑠)(𝜃)=,𝑞𝑀1,𝛾1,𝜌^𝑇−⁰₋₁^𝜃_𝜉=01,𝛾𝑒^{−𝑖(𝜉−𝜃)𝜔𝑛0}𝑑𝜂(𝜃)1,𝜌^𝑇𝑒^{𝑖𝜉𝜔𝑛0}𝑑𝜉, (A.13) then we choose 𝑀=1+𝜌𝛾+𝜌𝑐₁ℎ^∗ℎ^∗+𝑘₁⁻¹𝑒^{𝑖𝜔𝑛0} (A.14) such that (𝑞^∗(𝑠),𝑞(𝜃))=0 and (𝑞^∗(𝑠),𝑞(𝜃))=1. In other words, let Φ=(𝑞(𝜃),𝑞(𝜃)), Ψ=(𝑞^∗(𝑠),𝑞^∗(𝑠))^𝑇, then (Ψ,Φ)=𝐼, and 𝐼 is unit matrix.

Then the center subspace of model (A.4) is 𝑃=span{𝑞(𝜃),𝑞(𝜃)}, and the adjoint subspace is 𝑃^∗=span{𝑞^∗(𝑠),𝑞^∗(𝑠)}. Let 𝑣=(𝑣¹,𝑣²)^𝑇, where 𝑣¹=(1,0)^𝑇,𝑣²=(0,1)^𝑇. (A.15) Let 𝑚⋅𝑣 be defined by 𝑚⋅𝑣=𝑚₁𝑣¹+𝑚₂𝑣², (A.16) for 𝑚=(𝑚₁,𝑚₂)^𝑇∈([−1,0],𝑋). Hence the center subspace of linear system (A.4) is given by 𝑃_𝐶𝑁𝒞, where 𝑃_𝐶𝑁𝑃𝜙=ΦΨ,𝜙,𝑣⋅𝑣,𝜙∈𝒞,_𝐶𝑁𝒞=𝑞(𝜃)𝑧+𝑞(𝜃)𝑧⋅𝑣,𝑧∈ℂ, (A.17) and 𝒞=𝑃_𝐶𝑁𝒞⊕𝑃_𝑆𝒞, where 𝑃_𝑆𝒞 is the stable subspace.

From [17], we know that the infinitesimal generator 𝐴_𝑈 of linear model (A.4) satisfies 𝐴_𝑈𝜓=𝜓̇(𝜃); moreover 𝜓∈dom(𝐴_𝑈) if and only if 𝜓̇(𝜃)(0)(Δ),𝜓̇(𝜃0)∈𝒞,𝜓∈dom)(=𝜏₀Δ𝜓(0)+𝜏₀𝐿₀𝜓. (A.18)

First we define the coordinate to describe the center manifold at 𝛼=0; from center manifold we have 𝑤(𝑡,𝜃)=𝑤(𝑧(𝑡),𝑧(𝑡),𝜃)=𝑤₂₀(𝜃)𝑧²2+𝑤₁₁(𝜃)𝑧𝑧+𝑤₀₂(𝜃)𝑧²2+⋯. (A.19) The flow of model (A.2) in the center manifold can be written as follows: 𝑈_𝑡=Φ(𝑧,𝑧)^𝑇⋅𝑣+𝑤(𝑧,𝑧), (A.20) where 𝑧̇=𝑖𝜔^𝑛₀𝑧+𝑞^∗(0)𝐹Φ(𝑧,𝑧)^𝑇⋅𝑣+𝑤(𝑧,𝑧),0,𝑣=𝑖𝜔^𝑛₀𝑧+𝑔(𝑧,𝑧) (A.21) with 𝑔(𝑧,𝑧)=𝑞^∗(0)𝐹Φ(𝑧,𝑧)^𝑇⋅𝑣+𝑤(𝑧,𝑧),0,𝑣=𝑔₂₀𝑧²2+𝑔₁₁𝑧𝑧+𝑔₀₂𝑧²2+𝑔₂₁𝑧²𝑧2+⋯. (A.22)

Let 𝑓¹=𝐻+ℎ^∗𝑎₁−𝑏𝐻+ℎ^∗−𝑐₁𝐻+ℎ^∗𝑃(𝑡−1)+𝑝^∗𝑘₁+𝐻+ℎ^∗,𝑓²=𝑃+𝑝^∗𝑎₂−𝑐₂𝑃+𝑝^∗𝑘₂+𝐻+ℎ^∗. (A.23) From above equations, we have 𝑔₂₀2=𝑀𝜌𝑓¹₁₁₀+𝛾𝑓²₁₁₀+𝜌𝑒^{−𝑖𝜔𝑛0}𝑓¹₁₀₁+𝛾𝑓²₁₀₁+𝜌²𝑒^{−𝑖𝜔𝑛0}𝑓¹₀₁₁+𝛾𝑓²₀₁₁+𝑓¹₂₀₀+𝛾𝑓²₂₀₀2+𝜌²𝑓¹₀₂₀+𝛾𝑓²₀₂₀2+𝜌²𝑒^{−2𝑖𝜔𝑛0}𝑓¹₀₀₂+𝛾𝑓²₀₀₂2,𝑔₁₁(𝜃)=𝑀𝜌𝑓¹₁₁₀+𝛾𝑓²₁₁₀+𝜌𝑓¹₁₁₀+𝛾𝑓²₁₁₀+𝜌𝑒^{𝑖𝜔𝑛0}𝑓¹₁₀₁+𝛾𝑓²₁₀₁+𝜌𝑒^{−𝑖𝜔𝑛0}𝑓¹₁₀₁+𝛾𝑓²₁₀₁+𝜌𝜌𝑒^{−𝑖𝜔𝑛0}𝑓¹₀₁₁+𝛾𝑓²₀₁₁+𝜌𝜌𝑒^{𝑖𝜔𝑛0}𝑓¹₀₁₁+𝛾𝑓²₀₁₁+𝑓¹₂₀₀+𝛾𝑓²₂₀₀+𝜌𝜌𝑓¹₀₂₀+𝛾𝑓²₀₂₀+𝜌𝜌𝑓¹₀₁₁+𝑓²₀₁₁,𝑔₀₂2=𝑀𝜌𝑓¹₁₁₀+𝛾𝑓²₁₁₀+𝜌𝑒^{−𝑖𝜔𝑛0}𝑓¹₁₀₁+𝛾𝑓²₁₀₁+𝜌²𝑒^{−𝑖𝜔𝑛0}𝑓¹₀₁₁+𝛾𝑓²₀₁₁+𝑓¹₂₀₀+𝛾𝑓²₂₀₀2+𝜌²𝑓¹₀₂₀+𝛾𝑓²₀₂₀2+𝜌²𝑒^{2𝑖𝜔𝑛0}𝑓¹₀₀₂+𝛾𝑓²₀₀₂2,𝑔₂₁2=𝑀𝑓¹₁₁₀+𝛾𝑓²₁₁₀𝜌𝑤¹₁₁(0)+𝜌𝑤¹₂₀(0)2+𝑤²₂₀(0)2+𝑤²₁₁(0)+𝑓¹₁₀₁+𝛾𝑓²₁₀₁𝜌𝑒^{−𝑖𝜔𝑛0}𝑤¹₁₁(0)+𝜌𝑤¹₂₀(0)𝑒^{𝑖𝜔𝑛0}2+𝑤²₂₀(−1)2+𝑤²₁₁()+𝑓−1₁₀₁₁+𝛾𝑓²₀₁₁𝜌𝑒^{−𝑖𝜔𝑛0}𝑤²₁₁(0)+𝜌𝑤²₂₀(0)𝑒^{𝑖𝜔𝑛0}2+𝜌𝑤²₂₀(−1)2+𝑤²₁₁()𝜌+1−12𝑓¹₂₀₀+𝛾𝑓²₂₀₀2𝑤¹₁₁(0)+𝑤¹₂₀(0)+12𝑓¹₀₂₀+𝛾𝑓²₀₂₀2𝑤²₁₁(0)𝜌+𝑤²₂₀(0)𝜌+12𝑓¹₀₀₂+𝛾𝑓²₀₀₂2𝑤²₁₁(−1)𝑒^{−𝑖𝜔𝑛0}𝜌+𝑤²₂₀(−1)𝜌𝑒^{−𝑖𝜔𝑛0}. (A.24)

From [17], we can know that 𝑤(𝑧,𝑧) satisfies 𝑤̇=𝐴_𝑈𝑤+𝐻(𝑧,𝑧), (A.25) where 𝐻(𝑧,𝑧)=𝐻₂₀(𝜃)𝑧²2+𝐻₁₁(𝜃)𝑧𝑧+𝐻₀₂(𝜃)𝑧²2+⋯=𝑋₀𝐹𝑈_𝑡,0−ΦΨ,𝑋₀𝐹𝑈_𝑡,0,𝑣⋅𝑣. (A.26) Again we write near the origin on 𝐶₀𝑤̇=𝑧𝑤̇_𝑧+𝑧̇𝑤_𝑧. (A.27) By comparing (A.21) and (A.25) we get 𝐴_𝑈−2𝑖𝜔^𝑛₀𝑤₂₀(𝜃)=−𝐻₂₀(𝜃),𝐴_𝑈𝑤₁₁(𝜃)=−𝐻₁₁(𝜃). (A.28)

When −1≤𝜃<0, 𝐻(𝑧,𝑧)=−Φ(𝜃)Ψ(0)⟨𝐹(𝑈_𝑡,0),𝑣⟩⋅𝑣. Therefore, for −1≤𝜃<0, 𝐻₂₀(𝜃)=−𝑔₂₀𝑞(𝜃)+𝑔₀₂𝑞(𝜃𝐻)⋅𝑣,₂₁(𝜃)=−𝑔₂₁𝑞(𝜃)+𝑔₁₂𝑞(𝜃)⋅𝑣, (A.29) and for 𝜃=0, 𝐻(𝑧,𝑧)(0)=𝐹(𝑈_𝑡,0)−Φ(Ψ,⟨𝐹(𝑈_𝑡,0),𝑣⟩)⋅𝑣, then we obtain 𝐻₂₀(0)=−𝑔₂₀𝑞(0)+𝑔₀₂𝑞(0+𝜏𝜌𝑓)⋅𝑣₁₁₁₀+𝜌𝑒^{−𝑖𝜔𝑛0}𝑓¹₁₀₁+𝜌²𝑒^{−𝑖𝜔𝑛0}𝑓¹₀₁₁+𝑓¹₂₀₀2+𝜌²𝑓¹₀₂₀2+12𝜌²𝑒^{−2𝑖𝜔𝑛0}𝑓¹₀₀₂𝜌𝛾𝑓²₁₁₀+𝜌𝑒^{−𝑖𝜔𝑛0}𝛾𝑓²₁₀₁+𝜌²𝑒^{−𝑖𝜔𝑛0}𝛾𝑓²₀₁₁+𝛾𝑓²₂₀₀2+𝛾𝜌²𝑓²₀₂₀2+12𝜌²𝑒^{−2𝑖𝜔𝑛0}𝛾𝑓²₀₀₂,𝐻₁₁(0)=−𝑔₁₁𝑞(0)+𝑔₁₁𝑞(0+𝜏𝑓)⋅𝑣₁₁₁₀𝜌+𝜌+𝑓¹₁₀₁𝜌𝑒^{𝑖𝜔𝑛0}+𝜌𝑒^{−𝑖𝜔𝑛0}+𝑓¹₀₁₁𝜌𝜌𝑒^{𝑖𝜔𝑛0}+𝜌𝜌𝑒^{−𝑖𝜔𝑛0}+𝑓¹₂₀₀+𝑓¹₀₂₀𝜌𝜌+𝑓¹₀₀₂𝜌𝜌𝑓²₁₁₀𝜌+𝜌+𝑓²₁₀₁𝜌𝑒^{𝑖𝜔𝑛0}+𝜌𝑒^{−𝑖𝜔𝑛0}+𝑓²₀₁₁𝜌𝜌𝑒^{𝑖𝜔𝑛0}+𝜌𝜌𝑒^{−𝑖𝜔𝑛0}+𝑓²₂₀₀+𝑓²₀₂₀𝜌𝜌+𝑓²₀₀₂𝜌𝜌. (A.30)

By the definition of 𝐴_𝑈, we have from (A.28) 𝑤₂₀(𝜃)=𝑖𝑔₂₀𝜔^𝑛₀𝑞−(0)𝑒^{𝑖𝜃𝜔𝑛0}⋅𝑣−𝑖𝑔₀₂3𝜔^𝑛₀𝑞(0)𝑒^{−𝑖𝜃𝜔𝑛0}⋅𝑣+𝐸𝑒^{2𝑖𝜃𝜔𝑛0}, (A.31)where 𝐸=(𝐸¹,𝐸²)∈𝑅² is a constant vector.

Similarly, we get 𝑤₁₁(𝜃)=−𝑖𝑔₁₁𝜔^𝑛₀𝑞(0)𝑒^{𝑖𝜃𝜔𝑛0}𝑖⋅𝑣+𝑔₁₁𝜔^𝑛₀𝑞(0)𝑒^{−𝑖𝜃𝜔𝑛0}⋅𝑣+𝐸, (A.32) where 𝐸=(𝐸¹,𝐸²)∈𝑅² is a constant vector.

Combining (A.5) and (A.28), we obtain 𝐻₂₀(0)=2𝑖𝜔^𝑛₀𝜏𝑤₂₀(0)−𝜏𝐷Δ𝑤₂₀(0)−𝜏𝐿𝑤₂₀(𝜃), (A.33) therefore 𝜌𝑓¹₁₁₀+𝜌𝑒^{−𝑖𝜔𝑛0}𝑓¹₁₀₁+𝜌²𝑒^{−𝑖𝜔𝑛0}𝑓¹₀₁₁+𝑓¹₂₀₀2+𝜌²𝑓¹₀₂₀2+12𝜌²𝑒^{−2𝑖𝜔𝑛0}𝑓¹₀₀₂𝜌𝛾𝑓²₁₁₀+𝜌𝑒^{−𝑖𝜔𝑛0}𝛾𝑓²₁₀₁+𝜌²𝑒^{−𝑖𝜔𝑛0}𝛾𝑓²₀₁₁+𝛾𝑓²₂₀₀2+𝛾𝜌²𝑓²₀₂₀2+12𝜌²𝑒^{−2𝑖𝜔𝑛0}𝛾𝑓²₀₀₂=𝑔₂₀𝑞(𝜃)+𝑔₀₂𝑞(𝜃)+2𝑖𝜔^𝑛₀𝑤₂₀(𝜃)−⁰₋₁𝑑𝜂(𝜃)𝑤₂₀(𝜃)=−𝑔₂₀𝑞(0)+𝑔₀₂𝑞(0)3+𝐸2𝑖𝜔^𝑛₀−⁰₋₁𝑑𝜂(𝜃)𝑖𝑔₂₀𝜔^𝑛₀𝑞(𝜃)−𝑖𝑔₀₂3𝜔^𝑛₀𝑞(𝜃)−𝐸𝑒^{2𝑖𝜃𝜔𝑛0}=𝐸2𝑖𝜔^𝑛₀𝐼−⁰₋₁𝑑𝜂(𝜃)𝑒^{2𝑖𝜃𝜔𝑛0}=2𝑖𝐸𝜔^𝑛₀−1001−𝑏ℎ₊+𝑐₁ℎ₊𝑝₊ℎ₊+𝑘₁²−𝑐₁ℎ₊ℎ₊+𝑘₁𝑒^{2𝑖𝜔𝑛0}𝑎²₂𝑐₂−𝑎₂𝐸=2𝑖𝜔^𝑛₀+𝑏ℎ₊−𝑐₁ℎ₊𝑝₊ℎ₊+𝑘₁²𝑐₁ℎ₊ℎ₊+𝑘₁𝑒^{2𝑖𝜔𝑛0}−𝑎²₂𝑐₂2𝑖𝜔^𝑛₀+𝑎₂𝐸. (A.34)

From the above equation we can find the values of 𝐸¹ and 𝐸². From (A.28) we have ∫⁰₋₁𝑑𝜂(𝜃)𝑤₁₁(𝜃)=−𝐻₁₁(𝜃). Therefore 𝑓¹₁₁₀𝜌+𝜌+𝑓¹₁₀₁𝜌𝑒^{𝑖𝜔𝑛0}+𝜌𝑒^{−𝑖𝜔𝑛0}+𝑓¹₀₁₁𝜌𝜌𝑒^{𝑖𝜔𝑛0}+𝜌𝜌𝑒^{−𝑖𝜔𝑛0}+𝑓¹₂₀₀+𝑓¹₀₂₀𝜌𝜌+𝑓¹₀₀₂𝜌𝜌𝑓²₁₁₀𝜌+𝜌+𝑓²₁₀₁𝜌𝑒^{𝑖𝜔𝑛0}+𝜌𝑒^{−𝑖𝜔𝑛0}+𝑓²₀₁₁𝜌𝜌𝑒^{𝑖𝜔𝑛0}+𝜌𝜌𝑒^{−𝑖𝜔𝑛0}+𝑓²₂₀₀+𝑓²₀₂₀𝜌𝜌+𝑓²₀₀₂𝜌𝜌=𝑔₁₁𝑞(0)+𝑔₁₁𝑞(0)−⁰₋₁𝑑𝜂(𝜃)𝑤₁₁(𝜃)=𝑔₁₁𝑞(0)+𝑔₁₁𝑞(0)−⁰₋₁𝑑𝜂(𝜃)−𝑖𝑔₁₁𝜔^𝑛₀𝑞(0)𝑒^{𝑖𝜔𝑛0𝜃}+𝑖𝑔₁₁𝜔^𝑛₀𝑒^{−𝑖𝜔𝑛0𝜃}+𝐸=−⁰₋₁(𝜃)𝐸=𝑑𝜂𝑏ℎ₊−𝑐₁ℎ₊𝑝₊ℎ₊+𝑘₁²𝑐₁ℎ₊ℎ₊+𝑘₁−𝑎²₂𝑐₂𝑎₂𝐸. (A.35) In a similar manner we can compute the corresponding results in 𝐸¹ and 𝐸². Then 𝑔₂₁ can be determined. Based on the above analysis, we can see that each 𝑔_𝑖𝑗 can be determined by the parameters. Thus we can computer the following values which determine the direction and stability of bifurcating periodic orbits: 𝑐₁(0)=𝑖2𝜔^𝑛₀𝑔₁₁𝑔₂₀−2𝑔₁₁²−𝑔₀₂²3+𝑔₂₁2,𝑐𝜇=−Re₁(0)Re𝜆^′(𝜏),𝛽₂=2Re𝑐₁(0,𝑇₂)𝑐=−Im₁(0)+𝜇₂Im𝜆^′(𝜏)𝜔^𝑛₀𝜏. (A.36)

Acknowledgments

The authors would like to thank the anonymous referee for very helpful suggestions and comments which led to improvements of their original paper. This research was supported by Natural Science Foundation of Zhejiang Province (LY12A01014 and LQ12A01009), the National Basic Research Program of China (2012CB426510), and the Fund Project of Zhejiang Provincial Education Department (Y201223449 and Y201120383).