Bogdanov-Takens Bifurcation of a Delayed Ratio-Dependent Holling-Tanner Predator Prey System

Abstract

A delayed predator prey system with refuge and constant rate harvesting is discussed by applying the normal form theory of retarded functional differential equations introduced by Faria and Magalhães. The analysis results show that under some conditions the system has a Bogdanov-Takens singularity. A versal unfolding of the system at this singularity is obtained. Our main results illustrate that the delay has an important effect on the dynamical behaviors of the system.

1. Introduction

It is well known that the multiple bifurcations will occur when a predator prey system (ODE) with more interior positive equilibria, such as Bogdanov-Takens bifurcation, Hopf bifurcation, and backward bifurcation; see [1–5] for example. However, when the predator prey systems with delay and Bogdanov-Takens bifurcation are researched relative few (see [6–8] and the reference therein) using similar methods as [6–8], the authors of [9–11] consider the Bogdanov-Takens bifurcation of some delayed single inertial neuron or oscillator models.

Motivated by the works of [5, 6], we mainly consider the Bogdanov-Takens bifurcation of the following system: ̇𝑥𝑥=𝑟𝑥1−𝐾−𝛼𝑦(𝑥−𝑚)𝐴𝑦+𝑥−𝑚−̇(ℎ,𝑦=𝑠𝑦1−𝑏𝑦𝑡−𝜏)𝑥(𝑡−𝜏)−𝑚, (1) where 𝑥 and 𝑦 stand for prey and predator population (or densities) at time 𝑡, respectively. The predator growth is of logistic type with growth rate 𝑟 and carrying capacity 𝐾 in the absence of predation; 𝛼 and 𝐴 stand for the predator capturing rate and half saturation constant, respectively; 𝑠 is the intrinsic growth rate of predator; however, carrying capacity 𝑥/𝑏 (𝑏 is the conversion rate of prey into predators) is the function on the population size of prey. The parameters 𝛼, 𝐴, 𝑚, ℎ, 𝑠, 𝑏, and 𝜏 are all positive constants. 𝑚 is a constant number of prey using refuges; ℎ is the rate of prey harvesting.

For simplicity, we first rescale system (1). Let 𝑋=𝑥−𝑚, 𝑌=𝑦; system (1) can be written as (still denoting 𝑋, 𝑌 as 𝑥, 𝑦) 𝑥̇=𝑟(𝑥+𝑚)1−𝑥+𝑚𝐾−𝛼𝑥𝑦−𝐴𝑦+𝑥̇(ℎ,𝑦=𝑠𝑦1−𝑏𝑦𝑡−𝜏)𝑥(𝑡−𝜏). (2)

Next, let 𝑡=𝑟𝑡, 𝑋(𝑡)=𝑥(𝑡)/𝐾, and 𝑌(𝑡)=𝛼𝑦(𝑡)/𝑟𝐾; then system (2) takes the form (still denoting 𝑋, 𝑌, 𝑡 as 𝑥, 𝑦, 𝑡) 𝑥̇=(𝑥+𝑚)(1−𝑥−𝑚)−𝑥𝑦̇𝑦()𝑎𝑦+𝑥−ℎ,𝑦=𝛿𝑦𝛽−𝑡−𝜏𝑥(),𝑡−𝜏 (3) where 𝑚=𝑚/𝐾, 𝑎=𝐴𝑟/𝛼, 𝛿=𝑠ℎ/𝛼, 𝛽=𝛼/𝑏𝑟, ℎ=ℎ/𝑟, and 𝜏=𝑟𝜏.

When 𝜏=0, we have known that for some parameter values system (3) exhibits Bogdanov-Takens bifurcation (see [5]). Summarizing the methods used by [6] and the formulae in [12], the sufficient conditions which depend on delay to guarantee that system (3) has a Bogdanov-Takens singularity will be given. Therefore, the delay has effect on the occurrence of Bogdanov-Takens bifurcation.

In the next section we will compute the normal form and give the versal unfolding of system (3) at the degenerate equilibrium.

2. Bogdanov-Takens Bifurcation

System (3) can also be written as 𝑥̇=𝜏(𝑥+𝑚)(1−𝑥−𝑚)−𝑥𝑦,̇𝑦()𝑎𝑦+𝑥−ℎ𝑦=𝜏𝛿𝑦𝛽−𝑡−1𝑥().𝑡−1 (4)

Let 10<𝑚<2𝛽1−,1𝑎𝛽+1ℎ=ℎ=4𝛽𝑎𝛽+1−1²+𝑚𝛽.𝑎𝛽+1 (5)

Then 𝑃_∗=(𝑥_∗,𝑦_∗) is an interior positive equilibrium of systems (3) and (4), where 𝑥_∗=−(1/2)((𝛽/(𝑎𝛽+1))+2𝑚−1), 𝑦_∗=𝛽𝑥_∗.

In order to discuss the properties of system (4) in the neighborhood of the equilibrium 𝑃_∗=(𝑥_∗,𝑦_∗), let 𝑥=𝑥−𝑥_∗, 𝑦=𝑦−𝑦_∗; then 𝑃_∗ is translated to (0,0), and system (4) becomes (still denoting 𝑥, 𝑦 as 𝑥, 𝑦) 𝑥̇=𝜏𝛽𝑎𝛽+1²𝜏𝑥−𝑎𝛽+1²𝑦+𝑔₁𝑥²+𝑔₂𝑥𝑦+𝑔₃𝑦²̇+h.o.t,𝑦=𝜏𝛿𝛽²𝑥(𝑡−1)−𝜏𝛿𝛽𝑦(𝑡−1)+𝑔₄𝑥²(𝑡−1)+𝑔₅𝑥(𝑡−1)𝑦(𝑡−1)+𝑔₆𝑥(𝑡−1)𝑦+𝑔₇𝑦(𝑡−1)𝑦+h.o.t, (6) where h.o.t denotes the higher order terms and 𝑔₁=−2𝑎𝛿𝛽²𝜏(,𝑔𝛽+2𝑚−1)𝑎𝛽+1+𝜏₂=4𝑎𝛿𝛽𝜏(,𝑔𝛽+2𝑚−1)𝑎𝛽+1₃=−2𝑎𝛿𝜏(),𝑔𝛽+2𝑚−1𝑎𝛽+1₄=2𝜏𝛿𝛽²𝑎𝛽+1(),𝑔𝛽+2𝑚−1𝑎𝛽+1₅=−2𝜏𝛿𝛽𝑎𝛽+1(),𝑔𝛽+2𝑚−1𝑎𝛽+1₆=−2𝜏𝛿𝛽𝑎𝛽+1(),𝑔𝛽+2𝑚−1𝑎𝛽+1₇=2𝜏𝛿𝑎𝛽+1().𝛽+2𝑚−1𝑎𝛽+1 (7)

The characteristic equation of the linearized part of system (6) is 𝐹(𝜆)=𝜆²+𝜏𝛿𝛽𝑒^−𝜆−𝜏𝛽𝑎𝛽+1²𝜆. (8) Clearly, if 1𝛿=𝑎𝛽+1²,1𝜏≠𝛿𝛽,thatis,𝜏≠𝑎𝛽+1²𝛽; (9) then 𝜆=0 is double zero eigenvalue; if 𝛿=1/(𝑎𝛽+1)² and 𝜏=1/𝛿𝛽, that is, 𝜏=(𝑎𝛽+1)²/𝛽, then 𝜆=0 is triple zero eigenvalue. We will mainly discuss the first case in this paper.

According to the normal form theory developed by Faria and Magalhães [13], first, rewrite system (4) as 𝑋̇(𝑡)=𝐿(𝑋_𝑡), where 𝑋(𝑡)=(𝑥₁(𝑡),𝑥₂(𝑡)), 𝐿(𝜙)=𝐿^𝜙1_𝜙⁽⁻¹⁾₂₍₀₎, and 𝜙=(𝜙₁,𝜙₂). Take 𝐴₀ as the infinitesimal generator of system. Let Λ={0}, and denote by 𝑃 the invariant space of 𝐴₀ associated with the eigenvalue 𝜆=0; using the formal adjoint theory of RFDE in [13], the phase space 𝐶₁ can be decomposed by Λ as 𝐶₁=𝑃⊕𝑄. Define Φ and Ψ as the bases for 𝑃 and 𝑃^∗, the space associated with the eigenvalue 𝜆=0 of the adjoint equation, respectively, and to be normalized such that (Ψ,Φ)=𝐼, Φ̇=Φ𝐽, and Ψ̇=−𝐽Ψ, where Φ and Ψ are 2×2 matrices.

Next we will find the Φ(𝜃) and Ψ(𝑠) based on the techniques developed by [14].

Lemma 1 (see Xu and Huang [14]). The bases of 𝑃 and its dual space 𝑃^∗ have the following representations: 𝑃=spanΦ,Φ(𝜃)=𝜑₁(𝜃),𝜑₂(𝜃)𝑃,−1≤𝜃≤0,_∗=spanΨ,Ψ(𝑠)=col𝜓₁(𝑠),𝜓₂(𝑠),0≤𝑠≤1, (10) where 𝜑₁(𝜃)=𝜑⁰₁∈𝑅^𝑛⧵{0}, 𝜑₂(𝜃)=𝜑⁰₂+𝜑⁰₁𝜃, and 𝜑⁰₂∈𝑅^𝑛 and 𝜓₂(𝑠)=𝜓⁰₂∈𝑅^𝑛∗⧵{0}, and 𝜓₁(𝑠)=𝜓⁰₁−𝑠𝜓⁰₂, 𝜓⁰₁∈𝑅^𝑛∗, which satisfy(1)(𝐴+𝐵)𝜑⁰₁=0,(2)(𝐴+𝐵)𝜑⁰₂=(𝐵+𝐼)𝜑⁰₁,(3)𝜓⁰₂(𝐴+𝐵)=0,(4)𝜓⁰₁(𝐴+𝐵)=𝜓⁰₂(𝐵+𝐼),(5)𝜓⁰₂𝜑⁰₂−(1/2)𝜓⁰₂𝐵𝜑⁰₁+𝜓⁰₂𝐵𝜑⁰₂=1,(6)𝜓⁰₁𝜑⁰₂−(1/2)𝜓⁰₁𝐵𝜑⁰₁+𝜓⁰₁𝐵𝜑⁰₂+(1/6)𝜓⁰₂𝐵𝜑⁰₁−(1/2)𝜓⁰₂𝐵𝜑⁰₂=0.

By (6), we have 𝐴=^{𝜏𝛿𝛽−𝜏𝛿}₀₀ and 𝐵=⁰⁰_𝜏𝛿𝛽²_{−𝜏𝛿𝛽}; using Lemma 1, we obtain Φ(𝜃)=11Ψ(𝑠)=𝑚𝜏𝛿𝛽+𝜃𝛽𝛽𝜃,−1≤𝜃≤0,₁+𝑠𝛽𝑛𝑚₂−𝑠𝑛−𝛽𝑛𝑛,0≤𝑠≤1, (11) where 𝑚₁=𝜏𝛿𝛽(𝜏𝛿𝛽−2)/2(𝜏𝛿𝛽−1)²=𝜏𝛽[𝜏𝛽−2(𝑎𝛽+1)²]/2[𝜏𝛽−(𝑎𝛽+1)²]², 𝑚₂=(1+(𝜏𝛿𝛽−1)²)/2𝛽(𝜏𝛿𝛽−1)²=(1/2𝛽)(1+(𝑎𝛽+1)⁴/(𝜏𝛽−(𝑎𝛽+1))²), and 𝑛=𝜏𝛿/(𝜏𝛿𝛽−1)=𝜏/(𝜏𝛽−(𝑎𝛽+1)²).

The matrix 𝐽 satisfying Φ̇=Φ𝐽 is given by 𝐽=⁰¹₀₀. System (6) in the center manifold is equivalent to the system 𝑧̇=𝐽𝑧+Ψ(0)𝐹(Φ𝑧), where𝐹𝜙=𝑔₁𝜙²₁(0)+𝑔₂𝜙₁(0)𝜙₂(0)+𝑔₃𝜙²₂(0)𝑔+h.o.t₄𝜙²₁(−1)+𝑔₅𝜙₁(−1)𝜙₂(−1)+𝑔₆𝜙₁(−1)𝜙₂(0)+𝑔₇𝜙₂(−1)𝜙₂(0)+h.o.t. (12)It is easy to obtain Φ(𝜃)𝑧𝑧=¹+1𝑧𝜏𝛿𝛽+𝜃₂𝛽𝑧₁+𝑥₂𝜃,𝜙₁(0)=𝑧₁+1𝑧𝜏𝛿𝛽²,𝜙₁(−1)=𝑧₁+1𝑧𝜏𝛿𝛽−1₂,𝜙₂(0)=𝛽𝑧₁,𝜙₂(−1)=𝛽𝑧₁−𝑧₂. (13) Hence (6) becomes 𝑧̇₁=𝑧₂−𝜏𝑚₁𝑧²₁−2𝑚₁𝑧𝛿𝛽¹𝑧₂+𝑙₁𝑧²₂𝑧̇+h.o.t,₂=𝑛𝜏𝛽𝑧²₁+2𝑛𝛿𝑧₁𝑧₂+𝑙₂𝑧²₂+h.o.t, (14) where 𝑙₁=𝑚₁𝑔₁/𝜏²𝛿²𝛽²+2𝑚₂(𝛿+𝜇₁)(𝑎𝛽+1)(1−𝜏𝛿𝛽)/𝜏𝛿²(𝛽+(𝑎𝛽+1)(2𝑚−1)), 𝑙₂=−(𝑛𝑔₁/𝛽𝜏²𝛿²)+2𝑛(𝛿+𝜇₁)(𝑎𝛽+1)(1−𝜏𝛿𝛽)/𝜏𝛿²(𝛽+(𝑎𝛽+1)(2𝑚−1)).

After a series of transformations we obtain 𝑧̇₁=𝑧₂𝑧̇+h.o.t,₂=𝑑₁𝑧²₁+𝑑₂𝑧₁𝑧₂+h.o.t, (15) where 𝑑₁=𝛽𝑛𝜏=𝛽𝜏²/(𝜏𝛽−(𝑎𝛽+1)²), 𝑑₂=−𝜏𝑚₁+𝑛/𝛿=−(𝜏²𝛽²+2(𝑎𝛽+1)²[(𝑎𝛽+1)²−2𝜏𝛽])/(𝜏𝛽−(𝑎𝛽+1)²)²; if 𝜏≠(2±√2)/𝛿𝛽, then 𝑑₂≠0.

Hence, we have the following theorem.

Theorem 2. Let (5), (9), and 𝜏≠(2±√2)/𝛿𝛽 hold. Then the equilibrium 𝑃_∗ of system (4) is a Bogdanov-Takens singularity.

In the following, we will do versal unfolding of system (4) at 𝑃_∗: 𝑥̇=𝜏(𝑥+𝑚)(1−𝑥−𝑚)−𝑥𝑦,̇𝑎𝑦+𝑥−ℎ𝑦=𝜏𝛿+𝜇₁𝑦𝛽+𝜇₂−𝑦(𝑡−1)𝑥().𝑡−1 (16) When 𝜇₁=𝜇₂=0, system (16) has a Bogdanov-Takens singularity 𝑃_∗ and a two-dimensional center manifold exists.

The Taylor expansion of system (16) at 𝑃_∗ takes the form 𝑥̇=𝜏𝛿𝛽𝑥−𝜏𝛿𝑦+𝑔₁𝑥²+𝑔₂𝑥𝑦+𝑔₃𝑦²̇+h.o.t,𝑦=𝑤₀𝜇₂+𝜏𝛿+𝜇₁𝛽²𝑥(𝑡−1)−𝜏𝛿+𝜇₁𝛽𝑦(𝑡−1)+𝜏𝛿𝜇₂𝑦+𝑤₄𝑥²(𝑡−1)+𝑤₅𝑥(𝑡−1)𝑦(𝑡−1)+𝑤₆𝑥(𝑡−1)𝑦+𝑤₇𝑦(𝑡−1)𝑦+h.o.t, (17) where 𝑤₀()=−𝜏𝛿𝛽𝛽+2𝑚−1𝑎𝛽+12,𝑤𝑎𝛽+1₄=2𝜏𝛿+𝜇₁𝛽²𝑎𝛽+1(),𝑤𝛽+2𝑚−1𝑎𝛽+1₅=−2𝜏𝛿+𝜇₁𝛽𝑎𝛽+1(,𝑤𝛽+2𝑚−1)𝑎𝛽+1₆=−2𝜏𝛿+𝜇₁𝑎𝛽+1(),𝑤𝛽+2𝑚−1𝑎𝛽+1₇=2𝜏𝛿+𝜇₁𝑎𝛽+1().𝛽+2𝑚−1𝑎𝛽+1 (18)

We decompose the enlarged phase space BC of system as BC=𝑃⊕Ker𝜋. Then 𝑦 in system can be decomposed as 𝑦=Φ𝑧+𝑢 with 𝑧∈𝑅² and 𝑢∈𝑄^′. Hence, system is decomposed as 𝑧̇=𝐵₁+𝐵₂(0)𝐺(),̇𝑧+ΨΦ𝑧+𝑢𝑢=𝐴_𝑄′𝑢+(𝐼−𝜋)𝑋₀×𝐵₀+𝐵₂Φ₀𝑧+𝑢(0)+𝐺(Φ𝑧+𝑢), (19) where𝑋₀(𝜃)=𝐵𝐼,𝜃=00,−1≤𝜃<0,₀=0𝑤₀𝜇₂,𝐵₁=Ψ(0)𝐵₀=𝑤₀𝑚₂𝜇₂𝑤₀𝑛𝜇₂,𝐵₂=Ψ(0)𝜏𝛿𝛽𝜙₁(0)−𝜏𝛿𝜙₂(0)𝜏𝛿+𝜇₁𝛽²𝜙₁(−1)−𝜏𝛿+𝜇₁𝛽𝜙₂(−1)+𝜏𝛿𝜇₂𝜙₂(0),𝐺𝜙=𝑔₁𝜙²₁(0)+𝑔₂𝜙₁(0)𝜙₂(0)+𝑔₃𝜙²₂(0)𝑤+h.o.t₄𝜙²₁(−1)+𝑤₅𝜙₁(−1)𝜙₂(−1)+𝑤₆𝜙₁(−1)𝜙₂(0)+𝑤₇𝜙₂(−1)𝜙₂(0).+h.o.t (20)To compute the normal form of system at 𝑃_∗, consider 𝑧̇=𝐵₁+𝐵₂𝑧+Ψ(0)𝐺(Φ𝑧); together with (13) we obtain 𝑧̇₁=𝑤₀𝑚₂𝜇₂+𝑧₂+𝑚₂𝜏𝛿𝛽𝜇₂𝑧₁+𝑚₂𝛽𝜇₁𝛿𝑧₂−𝜏𝑚₁𝑧²₁−2𝑚₁𝑧𝛿𝛽¹𝑧₂+𝑙₁𝑧²₂𝑧̇+h.o.t,₂=𝑤₀𝑛𝜇₂+𝑛𝜏𝛿𝛽𝜇₂𝑧₁+𝑛𝛽𝜇₁𝛿𝑧₂+𝑛𝜏𝛽𝑧²₁+2𝑛𝛿𝑧₁𝑧₂+𝑙₂𝑧²₂+h.o.t. (21)

Following the normal form formula in Kuznetsov [12], system (21) can be reduced to 𝑧̇₁=𝑧₂𝑧̇+h.o.t,₂=𝛾₁+𝛾₂𝑧₁+𝑑₁𝑧²₁+𝑑₂𝑧₁𝑧₂+h.o.t, (22) where 𝛾₁()=−𝜏𝛿𝛽𝛽+2𝑚−1𝑎𝛽+12𝑎𝛽+1𝑛𝜇₂=𝜏²𝛽𝛽+(2𝑚−1)𝑎𝛽+12𝑎𝛽+1³𝑎𝛽+1²𝜇−𝜏𝛽²,𝛾₂=𝛽²𝑛²𝜏𝑛−𝜏𝑚₁𝛿𝜇₁+𝑛𝜏𝛿𝛽−2𝑤₀𝑛−𝑤𝛿𝛽⁰𝑛𝑙₂𝛿𝛽−2+2𝑚₂𝛽+𝑚𝛿𝛽𝑛𝜏₂𝜏𝛿²𝛽²+2𝑙₁𝑤₀𝑛𝛿𝛽+2𝑤₀𝑚₂−2𝑤₀𝑚²₂𝛽−2𝑤₀𝑚₂𝑙₂×𝛿𝛽−𝜏𝛿+𝑚₂𝜏𝛿𝛽+𝑛⁻¹𝜇₂. (23)

Hence system (4) exist the following bifurcation curves in a small neighborhood of the origin in the (𝜇₁,𝜇₂) plane.

Theorem 3. Let (5), (9), and 𝜏≠(2±√2)/𝛿𝛽 hold. Then system (4) admits the following bifurcations:(i)a saddle-node bifurcation curve 𝑆𝑁={(𝜇₁,𝜇₂)∶𝛾₁=(1/4𝑑₁)𝛾²₂}; (ii)a Hopf bifurcation curve 𝐻={(𝜇₁,𝜇₂)∶𝜇₂=0,𝛾₂<0}; (iii)a homoclinic bifurcation curve 𝐻𝐿={(𝜇₁,𝜇₂)∶𝛾₁=−(6/25𝑑₁)𝛾²₂,𝛾₂<0}.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This paper is supported by NSFC (11226142), Foundation of Henan Educational Committee (2012A110012), and Foundation of Henan Normal University (2011QK04, 2012PL03, and 1001).