# Bifurcations of a Ratio-Dependent Holling-Tanner System with Refuge and Constant Harvesting

Authors: Xia Liu; Yepeng Xing
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

The bifurcation properties of a predator prey system with refuge and constant harvesting are investigated. The number of the equilibria and the properties of the system will change due to refuge and harvesting, which leads to the occurrence of several kinds bifurcation phenomena, for example, the saddle-node bifurcation, Bogdanov-Takens bifurcation, Hopf bifurcation, backward bifurcation, separatrix connecting a saddle-node and a saddle bifurcation and heteroclinic bifurcation, and so forth. Our main results reveal much richer dynamics of the system compared to the system with no refuge and harvesting.


---

## Body

## 1. Introduction

The Holling-Tanner predator-prey system has attracted much attentions from both theoretical and mathematical biologists, especially, in [1] the authors considered the ratio-dependent system of the form
(1)x˙=rx(1-xK)-αxyAy+x,y˙=sy(1-byx),
where x and y stand for prey and predator population (or densities) at time t, respectively. The predator growth is of logistic type with growth rate r and carrying capacity K in the absence of predation; α and A stand for the predator capturing rate and half saturation constant, respectively; s is the intrinsic growth rate of predator; however, carrying capacity x/b (b is the conversion rate of prey into predators) is the function on the population size of prey. They studied the global properties and the existence and uniqueness of limit cycle for system (1).Generally speaking, from the views of the optimal management and exploitation of bioeconomic resources, it is necessary and meaningful to consider the refuge or harvesting of populations in some bioeconomic models; one can see [2–11], and the references therein.In this paper we will analyze the system (1) with refuge and harvesting of the form
(2)x˙=rx(1-xK)-αy(x-m-)Ay+x-m--h-,y˙=sy(1-byx-m-),
where r, K, α, A, m-, h-, s, and b are positive constants. m- is a constant number of prey using refuges, and h- is the rate of prey harvesting.For simplicity, we first rescale the system (2).LetX=x-m-, Y=y; system (2) can be written as (still denote X, Y as x, y)
(3)x˙=r(x+m-)(1-x+m-K)-αxyAy+x-h-,y˙=sy(1-byx).Next, letτ=rt, X=x/K, and Y=αy/rK, then system (3) takes the form (still denote X, Y, and τ as x, y, t)
(4)x˙=(x+m)(1-x-m)-xyay+x-h=P(x,y),y˙=δy(β-yx)=Q(x,y),
where m=m-/K, a=Ar/α, δ=sh-/α, β=α/br, and h=h-/r.From the view of biology, we are only interested in the dynamics of the system (4) in the first quadrant.The organization of this paper is as follows. In Section2, we discuss the existence and properties of the equilibria of system (4). In Section 3, all possible bifurcation phenomena of the model in terms of the five parameters are presented, and the numerical simulations about every bifurcation phenomena are exhibited.


## 2. Qualitative Analysis of Equilibria

To obtain the boundary equilibria the following equation can be obtained(5)x2+(2m-1)x+h-m(1-m)=0.
Its discriminant is Δ0=1-4h.Obviously,Δ0≥0 if 0<h≤1/4 and Δ0<0 if h>1/4.Hence, (5) has two distinct positive solutions x01=(1-2m+1-4h)/2, x02=(1-2m-1-4h)/2 if 0<m<1/2,m(1-m)<h<1/4, a positive solution x01 if 0<m<1,0<h<m(1-m), a double solution x-=(1-2m)/2>0 if 0<m<1/2,h=1/4, and a solution x03=1-2m when h=m(1-m) and 0<m<1/2.One can obtain the positive equilibrium of (4) by solving the equation
(6)x2+(βaβ+1+2m-1)x+h+m(m-1)=0.We can derive that system (4) has two positive equilibria P1=(x1,y1) and P2=(x2,y2) if
(7)0<m<12(1-βaβ+1),m(1-m)<h<14(βaβ+1-1)2+mβaβ+1,
where
(8)xi=-12(βaβ+1+2m-1)+(-1)i+1Δ2,yi=βxi,i=1,2,Δ=(βaβ+1-1)2+4mβaβ+1-4h.
Moreover, we can show that system (4) just exists one positive equilibrium P1 if 0<h<m(1-m) and 0<m<1.The positive equilibriumP3=(x3,y3) (P*=(x*,y*)) of system (4) exists if 0<m<(1/2)(1-β/(aβ+1)), and h=m(1-m)(0<m<(1/2)(1-β/(aβ+1)), h = (1/4)(β/(aβ+1)-1)2 + mβ/(aβ+1)), where x3=1-2m-β/(aβ+1), y3 = βx3, x* = -(1/2)(β/(aβ+1)+2m-1), and y*=βx*.Summarizing the previous discussion, the number and location of equilibria of system (4) can be described by the following lemmas.Lemma 1.
Let1/2≤m<1.  (i)
System (4)  has no equilibria when h≥m(1-m);(ii)
System (4)  exist two equilibria E1=(x01,0) and P1 when 0<h<m(1-m).Lemma 2.
Let0<(1/2)(1-β/(aβ+1))≤m<1/2.  (i)
System (4) has no equilibria when h>1/4.(ii)
System (4) has a unique equilibrium E-=(x-,0) when h=1/4.(iii)
System (4) has two equilibria E1=(x01,0) and E2=(x02,0) when m(1-m)<h<1/4.(iv)
System (4) has an equilibrium E3=(x03,0) when h=m(1-m).(v)
System (4)  has two equilibria E1 and P1 when 0<h<m(1-m).Lemma 3.
Let0<m<(1/2)(1-β/(aβ+1)) and h~=(1/4)(β/(aβ+1)-1)2+mβ/(aβ+1).  (i)
System (4)  has no equilibria when h>1/4.(ii)
System (4)  has a unique equilibrium E-=(x-,0) when h=1/4.(iii)
System (4)  has two equilibria E1 and E2 when h~<h<1/4.(iv)
System (4)  has three equilibria E1, E2, and P* when h=h~.(v)
System (4)  has four equilibria E1, E2, P1, and P2 when m(1-m)<h<h~.(vi)
System (4)  has two equilibria E1 and P1 when 0<h<m(1-m).(vii)
System (4) has two equilibria E3 and P3 when h=m(1-m).Next we discuss the dynamics of system (4) in the neighborhood of each feasible equilibria. Firstly, the Jacobian matrix of system (4) at E1 is
(9)ME1=[-1-4h-10δβ].
It is easy to see that E1, if exists, is a hyperbolic saddle.Secondly, the Jacobian matrix of system (4) at E2 is
(10)ME2=[1-4h-10δβ].
One can see that boundary equilibrium E2, if exists, is an unstable hyperbolic node.The Jacobian matrix of system (4) at E3 is
(11)ME3=[2m-1-10δβ].
Hence, E3, if exists, is a saddle.Similarly, we assume boundary equilibriumE- exists, and the Jacobian matrix of system (4) at E- is obtained as follows:
(12)ME-=[0-10δβ].Hence,E- is a saddle node.The previous discussion can be summarized as follows.Theorem 4.
If the equilibriaE1, E2, and E3 exist, then E1 and E3 are hyperbolic saddle, and E2 is a hyperbolic unstable node. Moreover, E1 and E2 merge into a saddle node E-=(x-,0) when h=1/4.Remark 5.
Note that ifh=1/4, then x˙=-(x+m-1/2)2-xy/(x+ay)<0, if h>1/4, then x˙≤1/4-h<0. Thus, the prey species may go extinct as time increases for some initial values when h≥1/4. That is, biological over harvesting occurs.In the following, we will discuss the properties of interior equilibria of system (4).


### 2.1. The Properties of Interior Equilibria

The Jacobian matrix of system (4) at P1 is
(13)MP1=[β(aβ+1)2-Δ-1(aβ+1)2δβ2-δβ].
The characteristic equation is λ2+A1λ+A2=0, where
(14)A1=Δ-β(1+aβ)2+δβ,A2=δβΔ>0.
Denote  that (15)A12-4A2=β2δ2-2(Δ+β(1+aβ)2)βδ+(Δ-β(1+aβ)2)2:=F(δ).
The discriminant of F(δ)=0 is Δ1=16β3Δ/(1+aβ)2>0, then F(δ)=0 has two distinct solutions δ1 and δ2 denoted by
(16)δ1=(Δβ-11+aβ)2,δ2=(Δβ+11+aβ)2.
If Δ≠β/(1+aβ)2, it is easy to see that P1 is a node if 0<δ<δ1 or δ>δ2, a degenerate node if δ=δ1 or δ=δ2, and a focus or a center type nonhyperbolic if δ1<δ<δ2.IfΔ=β/(1+aβ)2,P1 is a node if δ>δ2, and a degenerate node if δ=δ2, and a focus or a center-type nonhyperbolic if 0<δ<δ2.To discuss the stability ofP1, we need to determine the sign of A1. Define δ~=1/(1+aβ)2-(1/β)Δ, then A1=β(δ-δ~).Clearly, ifβ/(1+aβ)2≤Δ then A1>0 for all δ; if β/(1+aβ)2>Δ, then A1>0 when δ>δ~, A1≤0 when 0<δ≤δ~; by simple computation, one can obtain F(δ~)=4Δ(Δ-β/(aβ+1)2)<0, hence δ1<δ~<δ2.The Jacobian matrix of system (4) at P2 is
(17)MP2=[β(aβ+1)2+Δ-1(aβ+1)2δβ2-δβ].
Its determinant is detMP2=-δβΔ<0.Through the previous discussion, about the stability ofP1 and P2, we have the following theorem.Theorem 6.
EquilibriumP2, if exists, must be a hyperbolic saddle. Equilibrium P1, if exists, may be a node or a focus when β/(1+aβ)2≤Δ. And when β/(1+aβ)2>Δ, P1 is a stable focus for δ2>δ>δ~, a stable degenerate node for δ=δ2, a stable node for δ>δ2, an unstable node for 0<δ<δ1, an unstable degenerate node for δ=δ1, an unstable focus for δ1<δ<δ~, and a weak focus or a center for δ=δ~.The Jacobian matrix of system (4) at P3 is
(18)MP3=[2m-1+β(aβ+2)(aβ+1)2-1(aβ+1)2δβ2-δβ].
then by the existence condition of P3, detMP3=-βδ(2m-1+β/(aβ+1))>0, trMP3=2m-1-δβ+β/(1+aβ)+β/(1+aβ)2. Then by taking similar methods used in estimating the properties of P1, we have the following theorem.Theorem 7.
Leth=m(1-m), δ-=(1/β)(2m-1+β/(1+aβ)+β/(1+aβ)2). Then,  (i)
assume0<m≤(1/2)(1-β/(1+aβ)-β/(1+aβ)2), then P3 is stable;(ii)
assume(1/2)(1-β/(1+aβ)-β/(1+aβ)2)<m<(1/2)(1-β/(1+aβ)), then P3 is stable if δ>δ-, is unstable if 0<δ<δ-,  and is a weak focus or a center if δ=δ-.The Jacobian matrix of system (4) at P* is
(19)MP*=[β(aβ+1)2-1(aβ+1)2δβ2-δβ].
One can see that detMP*=0, which indicates that P* is a degenerate singularity and maybe has complicated properties, see the following theorem.Theorem 8.
Let0<m<(1/2)(1-β/(aβ+1)), h=h~. Then system (4) has three equilibria, where E1 is a hyperbolic saddle, E2 is a hyperbolic unstable node, and P* is a degenerate singularity. More precisely,  1°
ifδ≠1/(aβ+1)2, then P* is a saddle node;2°
ifδ=1/(aβ+1)2, then P* is a cusp of codimension 2.Proof.
In order to discuss the properties of system (4) in the neighborhood of the equilibrium P*=(x*,y*), we first take x-=x-x*, y-=y-y*, then P* is translated to (0,0), and system (4) becomes (still denote x-, y- as x, y)
(20)x˙=β(aβ+1)2x-1(aβ+1)2y-g1x2+g2xy-g3y2+O(|x,y|3),y˙=δβ2x-δβy+g4x2-g5xy+g6y2+O(|x,y|3),
where
(21)g1=2aβ2(aβ+1)2((2m-1)(aβ+1)+β)+1,g2=4aβ(aβ+1)2((2m-1)(aβ+1)+β),g3=2a(aβ+1)2((2m-1)(aβ+1)+β),g4=2δβ2(aβ+1)(2m-1)(aβ+1)+β,g5=4δβ(aβ+1)(2m-1)(aβ+1)+β,g6=2δ(aβ+1)(2m-1)(aβ+1)+β.
Clearly, if δ≠1/(aβ+1)2, then tr(MP*)≠0. P*=(x*,y*) is a saddle-node. We finish the proof of the part 1°.
Whenδ=1/(aβ+1)2, trMP*=0, which implies that both eigenvalues of the matrix MP* are zero. We rewrite system (20) as
(22)x˙=β(aβ+1)2x-1(aβ+1)2y-q1x2+q2xy-q3y2+O(|x,y|3),y˙=β2(aβ+1)2x-β(aβ+1)2y+q4x2-q5xy+q6y2+O(|x,y|3),
where
(23)q1=g1,q2=g2,q3=g3,q4=2β2(aβ+1)2((2m-1)(aβ+1)+β),q5=4β(aβ+1)2((2m-1)(aβ+1)+β),q6=2(aβ+1)2((2m-1)(aβ+1)+β).
By introducing variableτ=(β/(aβ+1)2)t into previous system and rewriting τ as t for simplicity, then we obtain that
(24)x˙=x-1βy-ω1x2+ω2xy-ω3y2+O(|x,y|3),y˙=βx-y+ω4x2-ω5xy+ω6y2+O(|x,y|3),
where
(25)ω1=2aβ(2m-1)(aβ+1)+β+(aβ+1)2β,ω2=4a(2m-1)(aβ+1)+β,ω3=2aβ((2m-1)(aβ+1)+β),ω4=2β(2m-1)(aβ+1)+β,ω5=4(2m-1)(aβ+1)+β,ω6=2β((2m-1)(aβ+1)+β).
We take transformationX0=x, Y0=x-(1/β)y into (24), then system (24) is transformed to
(26)X˙0=Y0+η1X02+η2X0Y0-ω3β2Y02+O(|X0,Y0|3),Y˙0=η3X02+η4X0Y0-η5Y02+O(|X0,Y0|3),
where
(27)η1=-ω3β2-ω1+ω2β,η2=2ω3β2-ω2β,η3=-ω1β-ω2β2+ω3β3+ω4-ω5β+ω6β2β,η4=2ω3β2-ω2β-ω5+2ω6β,η5=β(ω3β+ω6).
In order to obtain the canonical normal forms of system (26), we will perform a series of C∞ transformations of variables for system (26) in a small neighborhood of (0,0)T.
Firstly, performing the transformation by takingX1=X0, Y1=Y0-ω3β2Y02, then (26) becomes
(28)X˙1=Y1+η1X12+η2X1Y1+O(|X1,Y1|3),Y˙1=η3X12+η4X1Y1-η5Y12+O(|X1,Y1|3).
Secondly, performing the transformation by takingX2=X1, Y2=Y1+η5X1Y1, then (28) becomes
(29)X˙2=Y2+η1X22+(η2-η5)X2Y2+O(|X2,Y2|3),Y˙2=η3X22+η4X2Y2+O(|X2,Y2|3).
We perform the final transformation of variables by
(30)X=X2-η2-η52X22,Y=Y2+η1X22+O(|X2,Y2|3).
Then, we obtain
(31)X˙=Y,Y˙=η3X2+(2η1+η4)XY+O(|X,Y|3).
Note that
(32)η3=-(aβ+1)2β≠0,2η1+η4=-2(aβ+1)2β≠0,
which indicates that the origin (0,0) of (31) is a cusp of codimension 2. We complete the proof.


## 2.1. The Properties of Interior Equilibria

The Jacobian matrix of system (4) at P1 is
(13)MP1=[β(aβ+1)2-Δ-1(aβ+1)2δβ2-δβ].
The characteristic equation is λ2+A1λ+A2=0, where
(14)A1=Δ-β(1+aβ)2+δβ,A2=δβΔ>0.
Denote  that (15)A12-4A2=β2δ2-2(Δ+β(1+aβ)2)βδ+(Δ-β(1+aβ)2)2:=F(δ).
The discriminant of F(δ)=0 is Δ1=16β3Δ/(1+aβ)2>0, then F(δ)=0 has two distinct solutions δ1 and δ2 denoted by
(16)δ1=(Δβ-11+aβ)2,δ2=(Δβ+11+aβ)2.
If Δ≠β/(1+aβ)2, it is easy to see that P1 is a node if 0<δ<δ1 or δ>δ2, a degenerate node if δ=δ1 or δ=δ2, and a focus or a center type nonhyperbolic if δ1<δ<δ2.IfΔ=β/(1+aβ)2,P1 is a node if δ>δ2, and a degenerate node if δ=δ2, and a focus or a center-type nonhyperbolic if 0<δ<δ2.To discuss the stability ofP1, we need to determine the sign of A1. Define δ~=1/(1+aβ)2-(1/β)Δ, then A1=β(δ-δ~).Clearly, ifβ/(1+aβ)2≤Δ then A1>0 for all δ; if β/(1+aβ)2>Δ, then A1>0 when δ>δ~, A1≤0 when 0<δ≤δ~; by simple computation, one can obtain F(δ~)=4Δ(Δ-β/(aβ+1)2)<0, hence δ1<δ~<δ2.The Jacobian matrix of system (4) at P2 is
(17)MP2=[β(aβ+1)2+Δ-1(aβ+1)2δβ2-δβ].
Its determinant is detMP2=-δβΔ<0.Through the previous discussion, about the stability ofP1 and P2, we have the following theorem.Theorem 6.
EquilibriumP2, if exists, must be a hyperbolic saddle. Equilibrium P1, if exists, may be a node or a focus when β/(1+aβ)2≤Δ. And when β/(1+aβ)2>Δ, P1 is a stable focus for δ2>δ>δ~, a stable degenerate node for δ=δ2, a stable node for δ>δ2, an unstable node for 0<δ<δ1, an unstable degenerate node for δ=δ1, an unstable focus for δ1<δ<δ~, and a weak focus or a center for δ=δ~.The Jacobian matrix of system (4) at P3 is
(18)MP3=[2m-1+β(aβ+2)(aβ+1)2-1(aβ+1)2δβ2-δβ].
then by the existence condition of P3, detMP3=-βδ(2m-1+β/(aβ+1))>0, trMP3=2m-1-δβ+β/(1+aβ)+β/(1+aβ)2. Then by taking similar methods used in estimating the properties of P1, we have the following theorem.Theorem 7.
Leth=m(1-m), δ-=(1/β)(2m-1+β/(1+aβ)+β/(1+aβ)2). Then,  (i)
assume0<m≤(1/2)(1-β/(1+aβ)-β/(1+aβ)2), then P3 is stable;(ii)
assume(1/2)(1-β/(1+aβ)-β/(1+aβ)2)<m<(1/2)(1-β/(1+aβ)), then P3 is stable if δ>δ-, is unstable if 0<δ<δ-,  and is a weak focus or a center if δ=δ-.The Jacobian matrix of system (4) at P* is
(19)MP*=[β(aβ+1)2-1(aβ+1)2δβ2-δβ].
One can see that detMP*=0, which indicates that P* is a degenerate singularity and maybe has complicated properties, see the following theorem.Theorem 8.
Let0<m<(1/2)(1-β/(aβ+1)), h=h~. Then system (4) has three equilibria, where E1 is a hyperbolic saddle, E2 is a hyperbolic unstable node, and P* is a degenerate singularity. More precisely,  1°
ifδ≠1/(aβ+1)2, then P* is a saddle node;2°
ifδ=1/(aβ+1)2, then P* is a cusp of codimension 2.Proof.
In order to discuss the properties of system (4) in the neighborhood of the equilibrium P*=(x*,y*), we first take x-=x-x*, y-=y-y*, then P* is translated to (0,0), and system (4) becomes (still denote x-, y- as x, y)
(20)x˙=β(aβ+1)2x-1(aβ+1)2y-g1x2+g2xy-g3y2+O(|x,y|3),y˙=δβ2x-δβy+g4x2-g5xy+g6y2+O(|x,y|3),
where
(21)g1=2aβ2(aβ+1)2((2m-1)(aβ+1)+β)+1,g2=4aβ(aβ+1)2((2m-1)(aβ+1)+β),g3=2a(aβ+1)2((2m-1)(aβ+1)+β),g4=2δβ2(aβ+1)(2m-1)(aβ+1)+β,g5=4δβ(aβ+1)(2m-1)(aβ+1)+β,g6=2δ(aβ+1)(2m-1)(aβ+1)+β.
Clearly, if δ≠1/(aβ+1)2, then tr(MP*)≠0. P*=(x*,y*) is a saddle-node. We finish the proof of the part 1°.
Whenδ=1/(aβ+1)2, trMP*=0, which implies that both eigenvalues of the matrix MP* are zero. We rewrite system (20) as
(22)x˙=β(aβ+1)2x-1(aβ+1)2y-q1x2+q2xy-q3y2+O(|x,y|3),y˙=β2(aβ+1)2x-β(aβ+1)2y+q4x2-q5xy+q6y2+O(|x,y|3),
where
(23)q1=g1,q2=g2,q3=g3,q4=2β2(aβ+1)2((2m-1)(aβ+1)+β),q5=4β(aβ+1)2((2m-1)(aβ+1)+β),q6=2(aβ+1)2((2m-1)(aβ+1)+β).
By introducing variableτ=(β/(aβ+1)2)t into previous system and rewriting τ as t for simplicity, then we obtain that
(24)x˙=x-1βy-ω1x2+ω2xy-ω3y2+O(|x,y|3),y˙=βx-y+ω4x2-ω5xy+ω6y2+O(|x,y|3),
where
(25)ω1=2aβ(2m-1)(aβ+1)+β+(aβ+1)2β,ω2=4a(2m-1)(aβ+1)+β,ω3=2aβ((2m-1)(aβ+1)+β),ω4=2β(2m-1)(aβ+1)+β,ω5=4(2m-1)(aβ+1)+β,ω6=2β((2m-1)(aβ+1)+β).
We take transformationX0=x, Y0=x-(1/β)y into (24), then system (24) is transformed to
(26)X˙0=Y0+η1X02+η2X0Y0-ω3β2Y02+O(|X0,Y0|3),Y˙0=η3X02+η4X0Y0-η5Y02+O(|X0,Y0|3),
where
(27)η1=-ω3β2-ω1+ω2β,η2=2ω3β2-ω2β,η3=-ω1β-ω2β2+ω3β3+ω4-ω5β+ω6β2β,η4=2ω3β2-ω2β-ω5+2ω6β,η5=β(ω3β+ω6).
In order to obtain the canonical normal forms of system (26), we will perform a series of C∞ transformations of variables for system (26) in a small neighborhood of (0,0)T.
Firstly, performing the transformation by takingX1=X0, Y1=Y0-ω3β2Y02, then (26) becomes
(28)X˙1=Y1+η1X12+η2X1Y1+O(|X1,Y1|3),Y˙1=η3X12+η4X1Y1-η5Y12+O(|X1,Y1|3).
Secondly, performing the transformation by takingX2=X1, Y2=Y1+η5X1Y1, then (28) becomes
(29)X˙2=Y2+η1X22+(η2-η5)X2Y2+O(|X2,Y2|3),Y˙2=η3X22+η4X2Y2+O(|X2,Y2|3).
We perform the final transformation of variables by
(30)X=X2-η2-η52X22,Y=Y2+η1X22+O(|X2,Y2|3).
Then, we obtain
(31)X˙=Y,Y˙=η3X2+(2η1+η4)XY+O(|X,Y|3).
Note that
(32)η3=-(aβ+1)2β≠0,2η1+η4=-2(aβ+1)2β≠0,
which indicates that the origin (0,0) of (31) is a cusp of codimension 2. We complete the proof.


## 3. Bifurcation Analysis

From previous analysis, we can see the equilibria of system (4) may be hyperbolic or degenerate singularities under appropriate conditions, which indicate that some bifurcations may occur for system (4). It is interesting to investigate what kinds of bifurcations system (4) can undergo with the original parameters varying.


### 3.1. Hopf Bifurcation

Theorem6 shows that P1, if exists, is a weak focus or a center when
(33)V={β(1+aβ)2(a,m,β,δ,h):0<m<m~,m(1-m)<h<h~,β(1+aβ)2>Δ,δ=δ~},
where m~=(1/2)(1-β/(aβ+1)).To determine the direction of Hopf bifurcation and stability ofP1 in this case, we need to compute the Liapunov coefficients of the equilibrium P1. Let δ=δ~ and by the variable u=x-x1, v=y-y1. Then we rewrite system (4) (still denote u, v as x, y)  as follows:
(34)x˙=δ~βx-1(aβ+1)2y+f1(x,y),y˙=δ~β2x-δ~βy+f2(x,y).
We perform the transformations
(35)X=x,Y=δ~βx-1(aβ+1)2y,u=X,v=(aβ+1)Yβδ~(1-δ~(aβ+1)2)
and rewrite u, v as x, y. Then  the previous system can be transformed to
(36)x˙=k1y+a20x2+a11xy+a02y2+a30x3+a21x2y+a12xy2+a03y3+O(|(x,y)|4),y˙=-k1x+d11xy+b20x2+b11xy+b02y2+b30x3+b21x2y+b12xy2+b03y3+O(|(x,y)|4),
where the expressions of a20,a11,a02,a30,a21,a12,  a03,b20,b11,b02,b30,b21,b12, and b03 depend on the parameters a, β, δ, h, and m, and k1 = βδ~(1-δ~(aβ+1)2)/(aβ+1)>0.Using the formula, the first Liapunov number is(37)σ=-18βk1(aβ3(2aβ-1)[δ~(aβ+1)2-1]x12(aβ+1)42(aβ+1)2δ~(aβ+1)2-1+2β(2aβ+1)(aβ+1)x1+aβ3(2aβ-1)[δ~(aβ+1)2-1]x12(aβ+1)4).
Therefore, there exists a surface Hb (Hp) in the parameter space (a,m,β,δ,h) which satisfies
(38)Hb={(a,m,β,δ,h):σ>0,(a,m,β,δ,h)∈V},(Hp={(a,m,β,δ,h):σ<0,(a,m,β,δ,h)∈V}).Hence, when the parameter(a,m,β,δ,h) is in Hb (Hp), the equilibrium P1 of system (4) is a weak focus of multiplicity 1 and is unstable (stable) (see [8]). Hb (Hp) is called the subcritical (supercritical) Hopf bifurcation surface of system (4).From Theorem6, we know that P1 is a stable focus for δ2>δ>δ~ and (a,m,β,δ,h)∈V, an unstable focus for δ1<δ<δ~ and (a,m,β,δ,h)∈V.Theorem 9.
(i)   System (4) has at least one unstable limit cycle if (a,m,β,δ,h)∈Hb, δ2>δ>δ~, |δ-δ~|≪1.
(ii)   System (4) has at least one stable limit cycle if (a,m,β,δ,h)∈Hp, δ1<δ<δ~, |δ-δ~|≪1.Remark 10.
Whenσ=0 system (4) maybe undergoes degenerate Hopf bifurcation for some parameter values; since the expression of σ is complicated, we do not discuss this case.Note that by Theorem7, if P3 is a weak focus or a center, then we can obtain  that  its first Lyapunov number is
(39)σ1=3aβ8k1(1+aβ)(2maβ-aβ+β+2m-1)<0,
therefore, P3 is a stable weak focus.By numerical calculation, we give the parameter values(a,β,m,h)=(1.6,2.0,0.13095,0.128), then δ~=0.00235, δ2=0.2220 and k1σ=0.31085>0, and  the existence condition of subcritical Hopf bifurcation is satisfied. If we keep a, β, m, h fixed and choose δ=0.00249, then a unstable limit cycle can be shown in Figure 1(a).(a) System (4) shows an unstable limit cycle when a=1.6, β=2, m=0.13095, and  h=0.128. δ2≈0.2220>δ=0.00249>δ~≈0.00235; (b) System (4) shows a stable limit cycle when a=1.2, β=2, m=0.08235, and h=0.08396. δ1≈0.0000391<δ=0.002919<δ~≈0.003639.
(a)
(b)When taking(a,β,m,h)=(1.2,2.0,0.08235,0.08396), then δ~=0.003639, δ1=0.0000391, and k1σ=-0.1741543<0 which satisfy the existence condition of supercritical Hopf bifurcation. Furthermore, we choose δ=0.002919; according to Theorem 9, there exists a stable limit cycle, which can be shown in Figure 1(b).


### 3.2. Backward Bifurcation

DefineR=m(1-m)/h,R*=m(1-m)/h~.Lemmas2–3 and Theorems 6–8 illustrate that when the parameter h varies in the range of (0,m(1-m)], system (4) just has only one positive equilibrium P1 which is stable. However, when h varies in the range (m(m-1),h~), system (4) has two distinct positive equilibria P1 and P2, where P1 is a stable node or focus and P2 is a saddle. Furthermore, when h=h~, system (4) has unique positive equilibrium P*. The previous discussion indicates the possibility of a backward bifurcation, which can be summarized as follows.Theorem 11.
Let0<m<(1/2)(1-β/(aβ+1)), δ>δ2. Then system (4) has a unique positive equilibrium P* when R=R*, has two distinct positive equilibria P1 and P2 when R*<R<1, where P1 is a stable node and P2 is a saddle, and has one positive equilibrium P1 or P3 when R≥1. Therefore, system  (4) undergoes a backward bifurcation when R=1.We give a numerical example in Figure2 which displays that system (4) has a backward bifurcation at R=1.Figure 2
The figure of preyx at equilibria versus R when a=1.6, β=2, m=0.1309,  and h=0.13, which displays a backward bifurcation at R=1.


### 3.3. Saddle-Node Bifurcations

From Lemmas2–3 and Theorem 4, we see that when 0<m<1/2, h=1/4, E1 and E2 degenerate into a saddle-node E-=(x-,0). This indicates that there is a saddle node bifurcation surface which takes the form
(40)SN1={(m,a,β,h,δ):0<m<12,h=14,δ>0,β>0,a>014}.Similarly, from Lemma3 and the part 1° of Theorem 8, we know that when 0<m<(1/2)(1-β/(aβ+1)) and h=h~, in ℝ2+, system (4) admits the double point P*=(x*,y*). And P* is a saddle node if δ≠1/(aβ+1)2.One also can see that when the parameterh varies in the range of (m(m-1),h~), system (4) has two distinct positive equilibria P1 and P2. From Theorem 6, we know that P1 may be a stable node, or a focus, and P2 is a saddle. These imply that system (4) undergoes another saddle-node bifurcation of codimension 1. That is, there is a second saddle-node bifurcation surface SN2 which is defined by
(41)SN2={1(aβ+1)2(m,a,β,h,δ):0<m<12(1-βaβ+1),h=h~,δ≠1(aβ+1)2}.


### 3.4. Bogdanov-Takens Bifurcation

From the part2° of Theorem 8, we can see that system (4) exists a cusp of codimension 2, which implies that there may exist the Bogdanov-Takens bifurcation in system (4). Clearly, there exists a parameter space
(42)BT={(βaβ+1-1)2(m,a,β,h,δ):0<m<12-β2(aβ+1),δ=1(aβ+1)2,h=14(βaβ+1-1)2+mβaβ+1},
such that system (4) has a cusp of codimension 2 when (m,a,β,h,δ)∈BT.To show  that  system (4) undergoes the Bogdanov-Takens bifurcation we choose δ and β as bifurcation, parameters. We need to find the universal unfolding of P*.Let(m,a,β,h,δ)∈BT, and consider the following unfold system
(43)x˙=(x+m)(1-x-m)-yxay+x-h,y˙=(δ+μ1)y(β+μ2-yx),
where μ1 and μ2 are small parameters and vary in the neighborhood of the origin.TranslatingP* to (0,0) by the transformation X=x-x* and Y=y-y*. Then system (43) is rewritten as
(44)X˙=α1X-α2Y-α3X2+α4XY-α5Y2+W1(X,Y),Y˙=l1+l2X+l3Y+l4X2+l5XY+l6Y2+W2(X,Y),
where W1 and W2 are smooth functions of X, Y at least of the third order. And
(45)α1=β(aβ+1)2,α2=1(aβ+1)2,α3=(aβ+1)2[(aβ+1)(2m-1)+β]+2aβ2(aβ+1)2[(aβ+1)(2m-1)+β],α4=4aβ(aβ+1)2[(aβ+1)(2m-1)+β],α5=2a(aβ+1)2[(aβ+1)(2m-1)+β],l1=-[1+μ1(aβ+1)2][(aβ+1)(2m-1)+β]βμ22(aβ+1)3,l2=(1+μ1(aβ+1)2)β2(aβ+1)2,l3=-(1+μ1(aβ+1)2)(β-μ2)(aβ+1)2,l4=2β2(1+μ1(aβ+1)2)(aβ+1)[(aβ+1)(2m-1)+β],l5=-4β(1+μ1(aβ+1)2)(aβ+1)[(aβ+1)(2m-1)+β],l6=2(1+μ1(aβ+1)2)(aβ+1)[(aβ+1)(2m-1)+β].Taking the change of variablesX-=X, Y-=α1X-α2Y-α3X2+α4XY-α5Y2+W1(X,Y) and rewriting X-, Y- as X, Y,  we  obtain
(46)X˙=Y,Y˙=n1+n2X+n3Y+n4X2+n5XY+n6Y2+W3(X,Y,μ),
where
(47)n1=-α2l1→0,n2=α4l1-α2l2-l3α1-2α5α1l1α2→0,n3=l3+α1+2α5l1α2→0,n4=-l4α2-l5α1+l3α3+α4l2+-2α5α1l2+2α5l1α3-l6α12α2-2α5l1α4α1+l3α5α12α22+2α52l1α12α23→-β(aβ+1)2,n5=l5-2α3+2α5l2+α4α1+2l6α1α2+2α5α1l3+2α5l1α4α22-4α52l1α1α23→-2,n6=-l6α2-α4α2-l3α5α22+2α52l1α23→-2(aβ+1)(aβ+1)(2m-1)+β,
with μ1→0, μ2→0.Takingu=X+n3/n5, substituting u in system (46), and rewriting u as X, we get that
(48)X˙=Y,Y˙=n1-n2n3n5+n4n32n52+(n2-2n4n3n5)X+n4X2+n5XY+n6Y2+W4(X,Y,μ),
where W4 is a smooth function of X, Y, and μ at least of order three. When μ1→0, μ2→0,
(49)n1-n2n3n5+n4n32n52→0,n2-2n4n3n5→0.Next, lets=t/(1-n6X), x=X,  and y=(1-n6X)Y into (48) and rewriting s, x,  and y as t, X, and Y yields
(50)X˙=Y,Y˙=ε1+ε2X+n5XY+ε3X2+W5(X,Y,μ),
where W5 is a smooth function of X, Y, and μ at least of order three and
(51)ε1=n1-n2n3n5+n4n32n52→0,ε2=n2-2n4n3n5-2n6ε1→0,ε3=n4-2n6(n2-2n4n3n5)+n62ε1→-1(aβ+1)2<0,
when μ1→0, μ2→0. Let x=(n52/ε3)X, y=(n53/ε32)Y, ν=(ε3/n5)t, and rewrite ν as t. Then system (50) becomes
(52)x˙=y,y˙=τ1+τ2x+xy+x2+W6(x,y,μ),
where W6 is a smooth function of x, y, and μ at least of order three and τ1=n54ε1/ε33,τ2=n52ε2/ε32.Then system (4) exists the following bifurcation curves in a small neighborhood of the origin in the (μ1,μ2) plane.Theorem 12.
Let0<m<1/2-β/2(aβ+1), δ=1/(aβ+1)2, h=(1/4)(β/(aβ+1)-1)2+mβ/(aβ+1). Then system (43) admits the following bifurcations:  (i)
there exists a saddle node bifurcation curveSN={(μ1,μ2):4ε3ε1=ε22+o(∥μ∥)2};(ii)
there is a Hopf bifurcation curveH={(μ1,μ2):ε1=0+o(∥μ∥)2,ε2<0};(iii)
there is a homoclinic bifurcation curveHL={(μ1,μ2):25ε1ε3+6ε22=0+o(∥μ∥)2}.The biological interpretation for the Bogdanov-Takens bifurcation is that if the harvesting rateh and the prey refuge value m satisfy 0<m<1/2-β/2(aβ+1), h=(1/4)(β/(aβ+1)-1)2+mβ/(aβ+1), and δ=1/(aβ+1)2, then the predator and prey coexist in the form of a positive equilibrium or a periodic orbit for different initial values, respectively. And there exist other values of parameters, such that the predator and prey coexist in the form of a positive equilibrium for all initial values lying inside the homoclinic loop, and the predator and prey coexist in the form of a periodic orbit with infinite period for all initial values on the homoclinic loop. By choosing β=2, a=1.6, m=(1/4)(1-β/(aβ+1)), h=(1/4)(β/(aβ+1)-1)2+mβ/(aβ+1), and δ=1/(aβ+1)2, the numerical simulations for the Bogdanov-Takens bifurcation in Theorem 12 can be shown in Figures 3, 4 and 5.Figure 3
Bifurcation diagram of system (4) near P* in the plane of μ1 and μ2.(a) System (4) shows a cusp of codimension 2 when a=1.6, β=2, m=(1/4)(1-β/(aβ+1)), h=(1/4)(β/(aβ+1)-1)2+mβ/(aβ+1), δ=1/(aβ+1)2, and μ1=0, μ2=0; (b) the cusp of codimension 2 break into an unstable focus E^1 and a hyperbolic saddle E^2 when μ1=-0.0132, μ2=-0.02.
(a)
(b)(a) The cusp of codimension2 breaks into a stable focus E^1 and a hyperbolic saddle E^2 when μ1=-0.0098, μ2=-0.02. The change of stability of the focus yields an unstable limit cycle. (b) The unstable limit cycle is broken when μ1≈-0.0081768, μ2=-0.02, reachs the manifold of the saddle E^2, and leads to a homoclinic loop occur.
(a)
(b)


### 3.5. Separatrix Connecting a Saddle-Node and a Saddle Bifurcation and Heteroclinic Bifurcation

From Theorem8 and Lemma 3, when 0<m<(1/2)(1-β/(aβ+1)), h=h~, δ≠1/(aβ+1)2, there may exist a separatrix connecting the saddle-node P* and the saddle E1. When 0<m<(1/2)(1-β/(aβ+1)), m(1-m)<h<h~, the saddle node P* separates into the hyperbolic node P1 and the hyperbolic saddle P2, which implies that system (4) undergoes a separatrix connecting a saddle node and a saddle bifurcation. Furthermore, the heteroclinic bifurcation may occur if there exists a heteroclinic orbit connecting the separatrix of saddle E1 and saddle P2.


## 3.1. Hopf Bifurcation

Theorem6 shows that P1, if exists, is a weak focus or a center when
(33)V={β(1+aβ)2(a,m,β,δ,h):0<m<m~,m(1-m)<h<h~,β(1+aβ)2>Δ,δ=δ~},
where m~=(1/2)(1-β/(aβ+1)).To determine the direction of Hopf bifurcation and stability ofP1 in this case, we need to compute the Liapunov coefficients of the equilibrium P1. Let δ=δ~ and by the variable u=x-x1, v=y-y1. Then we rewrite system (4) (still denote u, v as x, y)  as follows:
(34)x˙=δ~βx-1(aβ+1)2y+f1(x,y),y˙=δ~β2x-δ~βy+f2(x,y).
We perform the transformations
(35)X=x,Y=δ~βx-1(aβ+1)2y,u=X,v=(aβ+1)Yβδ~(1-δ~(aβ+1)2)
and rewrite u, v as x, y. Then  the previous system can be transformed to
(36)x˙=k1y+a20x2+a11xy+a02y2+a30x3+a21x2y+a12xy2+a03y3+O(|(x,y)|4),y˙=-k1x+d11xy+b20x2+b11xy+b02y2+b30x3+b21x2y+b12xy2+b03y3+O(|(x,y)|4),
where the expressions of a20,a11,a02,a30,a21,a12,  a03,b20,b11,b02,b30,b21,b12, and b03 depend on the parameters a, β, δ, h, and m, and k1 = βδ~(1-δ~(aβ+1)2)/(aβ+1)>0.Using the formula, the first Liapunov number is(37)σ=-18βk1(aβ3(2aβ-1)[δ~(aβ+1)2-1]x12(aβ+1)42(aβ+1)2δ~(aβ+1)2-1+2β(2aβ+1)(aβ+1)x1+aβ3(2aβ-1)[δ~(aβ+1)2-1]x12(aβ+1)4).
Therefore, there exists a surface Hb (Hp) in the parameter space (a,m,β,δ,h) which satisfies
(38)Hb={(a,m,β,δ,h):σ>0,(a,m,β,δ,h)∈V},(Hp={(a,m,β,δ,h):σ<0,(a,m,β,δ,h)∈V}).Hence, when the parameter(a,m,β,δ,h) is in Hb (Hp), the equilibrium P1 of system (4) is a weak focus of multiplicity 1 and is unstable (stable) (see [8]). Hb (Hp) is called the subcritical (supercritical) Hopf bifurcation surface of system (4).From Theorem6, we know that P1 is a stable focus for δ2>δ>δ~ and (a,m,β,δ,h)∈V, an unstable focus for δ1<δ<δ~ and (a,m,β,δ,h)∈V.Theorem 9.
(i)   System (4) has at least one unstable limit cycle if (a,m,β,δ,h)∈Hb, δ2>δ>δ~, |δ-δ~|≪1.
(ii)   System (4) has at least one stable limit cycle if (a,m,β,δ,h)∈Hp, δ1<δ<δ~, |δ-δ~|≪1.Remark 10.
Whenσ=0 system (4) maybe undergoes degenerate Hopf bifurcation for some parameter values; since the expression of σ is complicated, we do not discuss this case.Note that by Theorem7, if P3 is a weak focus or a center, then we can obtain  that  its first Lyapunov number is
(39)σ1=3aβ8k1(1+aβ)(2maβ-aβ+β+2m-1)<0,
therefore, P3 is a stable weak focus.By numerical calculation, we give the parameter values(a,β,m,h)=(1.6,2.0,0.13095,0.128), then δ~=0.00235, δ2=0.2220 and k1σ=0.31085>0, and  the existence condition of subcritical Hopf bifurcation is satisfied. If we keep a, β, m, h fixed and choose δ=0.00249, then a unstable limit cycle can be shown in Figure 1(a).(a) System (4) shows an unstable limit cycle when a=1.6, β=2, m=0.13095, and  h=0.128. δ2≈0.2220>δ=0.00249>δ~≈0.00235; (b) System (4) shows a stable limit cycle when a=1.2, β=2, m=0.08235, and h=0.08396. δ1≈0.0000391<δ=0.002919<δ~≈0.003639.
(a)
(b)When taking(a,β,m,h)=(1.2,2.0,0.08235,0.08396), then δ~=0.003639, δ1=0.0000391, and k1σ=-0.1741543<0 which satisfy the existence condition of supercritical Hopf bifurcation. Furthermore, we choose δ=0.002919; according to Theorem 9, there exists a stable limit cycle, which can be shown in Figure 1(b).


## 3.2. Backward Bifurcation

DefineR=m(1-m)/h,R*=m(1-m)/h~.Lemmas2–3 and Theorems 6–8 illustrate that when the parameter h varies in the range of (0,m(1-m)], system (4) just has only one positive equilibrium P1 which is stable. However, when h varies in the range (m(m-1),h~), system (4) has two distinct positive equilibria P1 and P2, where P1 is a stable node or focus and P2 is a saddle. Furthermore, when h=h~, system (4) has unique positive equilibrium P*. The previous discussion indicates the possibility of a backward bifurcation, which can be summarized as follows.Theorem 11.
Let0<m<(1/2)(1-β/(aβ+1)), δ>δ2. Then system (4) has a unique positive equilibrium P* when R=R*, has two distinct positive equilibria P1 and P2 when R*<R<1, where P1 is a stable node and P2 is a saddle, and has one positive equilibrium P1 or P3 when R≥1. Therefore, system  (4) undergoes a backward bifurcation when R=1.We give a numerical example in Figure2 which displays that system (4) has a backward bifurcation at R=1.Figure 2
The figure of preyx at equilibria versus R when a=1.6, β=2, m=0.1309,  and h=0.13, which displays a backward bifurcation at R=1.


## 3.3. Saddle-Node Bifurcations

From Lemmas2–3 and Theorem 4, we see that when 0<m<1/2, h=1/4, E1 and E2 degenerate into a saddle-node E-=(x-,0). This indicates that there is a saddle node bifurcation surface which takes the form
(40)SN1={(m,a,β,h,δ):0<m<12,h=14,δ>0,β>0,a>014}.Similarly, from Lemma3 and the part 1° of Theorem 8, we know that when 0<m<(1/2)(1-β/(aβ+1)) and h=h~, in ℝ2+, system (4) admits the double point P*=(x*,y*). And P* is a saddle node if δ≠1/(aβ+1)2.One also can see that when the parameterh varies in the range of (m(m-1),h~), system (4) has two distinct positive equilibria P1 and P2. From Theorem 6, we know that P1 may be a stable node, or a focus, and P2 is a saddle. These imply that system (4) undergoes another saddle-node bifurcation of codimension 1. That is, there is a second saddle-node bifurcation surface SN2 which is defined by
(41)SN2={1(aβ+1)2(m,a,β,h,δ):0<m<12(1-βaβ+1),h=h~,δ≠1(aβ+1)2}.


## 3.4. Bogdanov-Takens Bifurcation

From the part2° of Theorem 8, we can see that system (4) exists a cusp of codimension 2, which implies that there may exist the Bogdanov-Takens bifurcation in system (4). Clearly, there exists a parameter space
(42)BT={(βaβ+1-1)2(m,a,β,h,δ):0<m<12-β2(aβ+1),δ=1(aβ+1)2,h=14(βaβ+1-1)2+mβaβ+1},
such that system (4) has a cusp of codimension 2 when (m,a,β,h,δ)∈BT.To show  that  system (4) undergoes the Bogdanov-Takens bifurcation we choose δ and β as bifurcation, parameters. We need to find the universal unfolding of P*.Let(m,a,β,h,δ)∈BT, and consider the following unfold system
(43)x˙=(x+m)(1-x-m)-yxay+x-h,y˙=(δ+μ1)y(β+μ2-yx),
where μ1 and μ2 are small parameters and vary in the neighborhood of the origin.TranslatingP* to (0,0) by the transformation X=x-x* and Y=y-y*. Then system (43) is rewritten as
(44)X˙=α1X-α2Y-α3X2+α4XY-α5Y2+W1(X,Y),Y˙=l1+l2X+l3Y+l4X2+l5XY+l6Y2+W2(X,Y),
where W1 and W2 are smooth functions of X, Y at least of the third order. And
(45)α1=β(aβ+1)2,α2=1(aβ+1)2,α3=(aβ+1)2[(aβ+1)(2m-1)+β]+2aβ2(aβ+1)2[(aβ+1)(2m-1)+β],α4=4aβ(aβ+1)2[(aβ+1)(2m-1)+β],α5=2a(aβ+1)2[(aβ+1)(2m-1)+β],l1=-[1+μ1(aβ+1)2][(aβ+1)(2m-1)+β]βμ22(aβ+1)3,l2=(1+μ1(aβ+1)2)β2(aβ+1)2,l3=-(1+μ1(aβ+1)2)(β-μ2)(aβ+1)2,l4=2β2(1+μ1(aβ+1)2)(aβ+1)[(aβ+1)(2m-1)+β],l5=-4β(1+μ1(aβ+1)2)(aβ+1)[(aβ+1)(2m-1)+β],l6=2(1+μ1(aβ+1)2)(aβ+1)[(aβ+1)(2m-1)+β].Taking the change of variablesX-=X, Y-=α1X-α2Y-α3X2+α4XY-α5Y2+W1(X,Y) and rewriting X-, Y- as X, Y,  we  obtain
(46)X˙=Y,Y˙=n1+n2X+n3Y+n4X2+n5XY+n6Y2+W3(X,Y,μ),
where
(47)n1=-α2l1→0,n2=α4l1-α2l2-l3α1-2α5α1l1α2→0,n3=l3+α1+2α5l1α2→0,n4=-l4α2-l5α1+l3α3+α4l2+-2α5α1l2+2α5l1α3-l6α12α2-2α5l1α4α1+l3α5α12α22+2α52l1α12α23→-β(aβ+1)2,n5=l5-2α3+2α5l2+α4α1+2l6α1α2+2α5α1l3+2α5l1α4α22-4α52l1α1α23→-2,n6=-l6α2-α4α2-l3α5α22+2α52l1α23→-2(aβ+1)(aβ+1)(2m-1)+β,
with μ1→0, μ2→0.Takingu=X+n3/n5, substituting u in system (46), and rewriting u as X, we get that
(48)X˙=Y,Y˙=n1-n2n3n5+n4n32n52+(n2-2n4n3n5)X+n4X2+n5XY+n6Y2+W4(X,Y,μ),
where W4 is a smooth function of X, Y, and μ at least of order three. When μ1→0, μ2→0,
(49)n1-n2n3n5+n4n32n52→0,n2-2n4n3n5→0.Next, lets=t/(1-n6X), x=X,  and y=(1-n6X)Y into (48) and rewriting s, x,  and y as t, X, and Y yields
(50)X˙=Y,Y˙=ε1+ε2X+n5XY+ε3X2+W5(X,Y,μ),
where W5 is a smooth function of X, Y, and μ at least of order three and
(51)ε1=n1-n2n3n5+n4n32n52→0,ε2=n2-2n4n3n5-2n6ε1→0,ε3=n4-2n6(n2-2n4n3n5)+n62ε1→-1(aβ+1)2<0,
when μ1→0, μ2→0. Let x=(n52/ε3)X, y=(n53/ε32)Y, ν=(ε3/n5)t, and rewrite ν as t. Then system (50) becomes
(52)x˙=y,y˙=τ1+τ2x+xy+x2+W6(x,y,μ),
where W6 is a smooth function of x, y, and μ at least of order three and τ1=n54ε1/ε33,τ2=n52ε2/ε32.Then system (4) exists the following bifurcation curves in a small neighborhood of the origin in the (μ1,μ2) plane.Theorem 12.
Let0<m<1/2-β/2(aβ+1), δ=1/(aβ+1)2, h=(1/4)(β/(aβ+1)-1)2+mβ/(aβ+1). Then system (43) admits the following bifurcations:  (i)
there exists a saddle node bifurcation curveSN={(μ1,μ2):4ε3ε1=ε22+o(∥μ∥)2};(ii)
there is a Hopf bifurcation curveH={(μ1,μ2):ε1=0+o(∥μ∥)2,ε2<0};(iii)
there is a homoclinic bifurcation curveHL={(μ1,μ2):25ε1ε3+6ε22=0+o(∥μ∥)2}.The biological interpretation for the Bogdanov-Takens bifurcation is that if the harvesting rateh and the prey refuge value m satisfy 0<m<1/2-β/2(aβ+1), h=(1/4)(β/(aβ+1)-1)2+mβ/(aβ+1), and δ=1/(aβ+1)2, then the predator and prey coexist in the form of a positive equilibrium or a periodic orbit for different initial values, respectively. And there exist other values of parameters, such that the predator and prey coexist in the form of a positive equilibrium for all initial values lying inside the homoclinic loop, and the predator and prey coexist in the form of a periodic orbit with infinite period for all initial values on the homoclinic loop. By choosing β=2, a=1.6, m=(1/4)(1-β/(aβ+1)), h=(1/4)(β/(aβ+1)-1)2+mβ/(aβ+1), and δ=1/(aβ+1)2, the numerical simulations for the Bogdanov-Takens bifurcation in Theorem 12 can be shown in Figures 3, 4 and 5.Figure 3
Bifurcation diagram of system (4) near P* in the plane of μ1 and μ2.(a) System (4) shows a cusp of codimension 2 when a=1.6, β=2, m=(1/4)(1-β/(aβ+1)), h=(1/4)(β/(aβ+1)-1)2+mβ/(aβ+1), δ=1/(aβ+1)2, and μ1=0, μ2=0; (b) the cusp of codimension 2 break into an unstable focus E^1 and a hyperbolic saddle E^2 when μ1=-0.0132, μ2=-0.02.
(a)
(b)(a) The cusp of codimension2 breaks into a stable focus E^1 and a hyperbolic saddle E^2 when μ1=-0.0098, μ2=-0.02. The change of stability of the focus yields an unstable limit cycle. (b) The unstable limit cycle is broken when μ1≈-0.0081768, μ2=-0.02, reachs the manifold of the saddle E^2, and leads to a homoclinic loop occur.
(a)
(b)


## 3.5. Separatrix Connecting a Saddle-Node and a Saddle Bifurcation and Heteroclinic Bifurcation

From Theorem8 and Lemma 3, when 0<m<(1/2)(1-β/(aβ+1)), h=h~, δ≠1/(aβ+1)2, there may exist a separatrix connecting the saddle-node P* and the saddle E1. When 0<m<(1/2)(1-β/(aβ+1)), m(1-m)<h<h~, the saddle node P* separates into the hyperbolic node P1 and the hyperbolic saddle P2, which implies that system (4) undergoes a separatrix connecting a saddle node and a saddle bifurcation. Furthermore, the heteroclinic bifurcation may occur if there exists a heteroclinic orbit connecting the separatrix of saddle E1 and saddle P2.



---

*Source: 478315-2013-03-14.xml*