# Bifurcation Behavior Analysis in a Predator-Prey Model

Authors: Nan Wang; Min Zhao; Hengguo Yu; Chuanjun Dai; Beibei Wang; Pengfei Wang
Journal: Discrete Dynamics in Nature and Society
 (2016)

---

## Abstract

A predator-prey model is studied mathematically and numerically. The aim is to explore how some key factors influence dynamic evolutionary mechanism of steady conversion and bifurcation behavior in predator-prey model. The theoretical works have been pursuing the investigation of the existence and stability of the equilibria, as well as the occurrence of bifurcation behaviors (transcritical bifurcation, saddle-node bifurcation, and Hopf bifurcation), which can deduce a standard parameter controlled relationship and in turn provide a theoretical basis for the numerical simulation. Numerical analysis ensures reliability of the theoretical results and illustrates that three stable equilibria will arise simultaneously in the model. It testifies the existence of Bogdanov-Takens bifurcation, too. It should also be stressed that the dynamic evolutionary mechanism of steady conversion and bifurcation behavior mainly depend on a specific key parameter. In a word, all these results are expected to be of use in the study of the dynamic complexity of ecosystems.


---

## Body

## 1. Introduction

The dynamical behaviors between different populations and their complex properties have been given close attention by biologists and ecologists. Since the pioneering work of Lotka and Volterra, the research interest in predator-prey dynamics has achieved constant attention. It is well known that these models can directly reflect changes in the size of populations. Considerable improvements are that the relevant theories become more and more complete in this category in recent years [1–5].The predator-prey models have extensive applicability in the field of biological problems. The biologist can use them to study the relationship between species in different domains [6–10]. Yang and Zhao [11] have established a fish-algae consumption model to explore how to apply the complex dynamics between fish-algae populations to expound the mechanism of algae blooms; these results will be helpful in controlling algae bloom. González-Olivares and Rojas-Palma [12] have established a Gause type predator-prey model with Allee effects and considered three standard functional responses, respectively. They found that different types of functional responses will lead to a model’s dynamic behavior change. Their results perfectly explore that the expression of one interaction term has a significant effect on the stability and persistence of the population model, which is meaningful in establishing biological mathematical model. Dhar et al. [13] have proposed a mathematical model to study how instantaneous nutrient recycling affects the dynamic characteristic of aquatic ecosystem. The nutrient supply rate was found to affect the local stability of equilibrium; this work was significant for models involving flowing waters. Luo [14] has considered a mathematical model to study how the periodic environment influences the internal operating characteristics of the aquatic ecosystem. He pointed out that the temperature has a certain time periodicity, the fluctuation of temperature can lead to changes in intrinsic carrying capacity, and the growth rate of prey populations can also be influenced by time periodicity; these results in [14] are more accordant with the actual situation.It is a common phenomenon in nature that one predator lives on multiple prey species. To address this issue, some researchers have begun to consider alternative prey to depict the dynamic predator mechanism [15–18]. It is easy to find out that alternative prey can variously affect the dynamic capture feature. On one hand, alternative prey can increase the predation quantity for the focal prey, because more prey biomass may result in higher predation rates for both prey items. On the other hand, the alternative prey population can also lower predation on the focal prey because of predator preference for the alternative prey resources [19]. At the same time, many studies show that the alternative population can intensively influence the dynamic behavior of the aquatic ecosystem [20, 21]. Based on this mechanism,  Kar and Chattopadhyay [22] have developed a two-species predator-prey model which includes the effect of alternative prey. The model can be depicted as(1)dNdT=r1N1-Nk1-a1Nb1+NP,dPdT=a1e1Nb1+NP+s1P1-Nk1-d1P,where N and P represent prey and predator population densities or biomass at time T, respectively. Here, r1 stands for the intrinsic growth rate of the prey without any environment limitations, k1 is the environmental carrying capacity of the prey in the absence of the predator, a1N/b1+N is the Holling type-II functional response [23], which is used to depict the average feeding rate of the predator when the predator spends time seeking prey, where b1 is the half saturation constant for the Holling type-II, a1 is the grazing rate of the predator population, e1 and d1 are the conversional rate and mortality rate of predator, respectively, and s1P1-N/k1 indicates the portion of biomass of predator increments from the alternative prey, where s1 represents the growth rate of the predator on account of the alternative prey. From the formula, we can see that when the quantity of focal prey N approaches the environmental carrying capacity k1, the amount of alternative prey consumed by the predator will tend to be zero [24].The concept of Allee effect was firstly derived from the research of Allee and Bowen [25]. Since then, the Allee effect received the attention of many researchers [26–29]. Allee effects can be roughly classified into two types: strong and weak [30]. For these two forms, there is a critical value that is referred to as the Allee threshold, respectively. The fist form means that if the population size is below the threshold, the species will become extinct. When the growth rate gradually decreases but remains positive with a low population size, the Allee effect is described as weak. Aulisa and Jang [31] have established a continuous-time predator-prey model to study influences in dynamical behaviors when the prey population possesses Allee effect. They pointed out that both species will become extinct if the prey population size falls below a certain threshold. Pan et al. [32] have considered a reaction-diffusion phytoplankton-zooplankton model with double Allee effects on prey population. They pointed out that the Allee effects can make the dynamical behaviors of a system increasingly complex.The strong Allee effect can be depicted by the following form:(2)dNdT=r1N1-Nk1Nm1-1,where m1 is Allee effect threshold. If population density or size is below the threshold, this population is doomed to extinction. Here 0<m1<k1, the other parameters’ significance is the same as model (1).Now we will establish a predator-prey model with strong Allee effects:(3)dNdT=r1N1-Nk1Nm1-1-a1Nb1+NP,dPdT=a1e1Nb1+NP+s1P1-Nk1-d1P.The parameters are all greater than zero and have the same significance as above. For simplicity, we write the above model in dimensionless form as follows: we take the scalingN=k1n, P=r1p/a1, and T=t/r1; then, model (3) can be simplified as(4a)dndt=n1-nmn-1-nb+np=F1,(4b)dpdt=enb+np+sp1-n-dp=F2,where e=a1e1/r1,  b=b1/k1,  s=s1/r1,  d=d1/r1,  m=k1/m1, and apparently m>1.In the succedent subsections of this paper, we introduce the problem of determining the number of equilibria. Stability analysis of equilibrium point is also presented. Then, we provide a demonstration of several types of bifurcation and give a set of parameter values to prove the existence of Bogdanov-Takens (BT) bifurcation. In Section3, numerical simulations are given to illustrate the theoretical analysis results, followed by conclusions in Section 4.


## 2. Qualitative Analysis

### 2.1. Equilibria

In this subsection, we mainly concentrate on the existence of positive equilibrium of model (4a) and (4b). Model (4a) and (4b) has three boundary equilibria: E00,0, E11,0, and E21/m,0 and they are always existent. The interior equilibria are the intersection points of the vertical isocline and horizontal isocline in the interior of the first quadrant. The expressions of vertical and horizontal isoclines are as follows:(5a)1-nmn-1-pb+n=0,(5b)enb+n+s1-n-d=0.The horizontal isocline is vertical line and the number of perpendiculars is decided by (5b). We denote (5b) by H(n)=Gn/b+n=0, where G(n)=-sn2+e+s-sb-dn+b(s-d). G(n)=0 have the same solutions as (5b), and two solutions exist in the equation at most. We denoted them by n3 and n4, where n4<n3. From the geometric properties of isoclines we can know that model (4a) and (4b) has interior equilibrium points if and only if 1/m<n3<1 or 1/m<n4<1. Then, the maximum number of equilibria for model (4a) and (4b) is five. We use E3n3,p3 and E4n4,p4 to signify the interior equilibria, the existence and stability of which are shown as follows.The Jacobian matrix of model (4a) and (4b) at the equilibrium point Ei is(6)Jni,pi=Ani,piBni,piCni,piDni,pi,(7)Ani,pi=-3mni2+2m+2ni-1-bpib+ni2,Bni,pi=-nib+ni,Cni,pi=bepib+ni2-spi,Dni,pi=enib+ni+s-d-sni.Theorem 1.
(1)E00,0 is an asymptotically stable node point if s<d and a saddle point for s>d; if s=d, it is a high order singularity. The Jacobian matrix around E0 shows that(8)J0,0=-100s-d.(2) E1(1,0) is a asymptotically stable node point when e/b+1<d while it is a saddle point when e/b+1>d. The Jacobian matrix around E1 is (9)J1,0=1-m-1b+10eb+1-d.(3) E21/m,0 is always unstable. The Jacobian matrix around E2 is(10)J1/m,0=1-1m-1bm+10ebm+1+s-d-sm.One of the eigenvalues is 1-1/m>0, so E2 is unstable.Theorem 2.
Ife-d-bd>0 and -s/m2+e+s-sb-d/m+bs-d<0 hold, model (4a) and (4b) has only one interior equilibrium E3n3,p3. Furthermore, E3(n3,p3) is always unstable, especially when d=s+e+sb-2bse; it is a high-order singularity, where(11a)n3=s+e-sb-d+sb-s+d-e2+4bss-d2s,(11b)p3=1-n3mn3-1b+n3.Proof.
According to (6),(12)det⁡Jn3,p3=-Cn3,p3Bn3,p3=n3p3b+n3beb+n32-s.Let ϕ=n3-be/s+b; substituting (11a) into ϕ, we get(13)ϕ=s+e-d+sb-2bes+s+e-d-sb2+4bss-d2s=s+e-d+sb-2bes+s+e-d+sb2-4bse2s≥0.
Then, we know whenn3≥be/s-b and detJn3,p3≤0, E3 is a saddle point or a high-order singularity.Theorem 3.
Ife-d-bd<0 and -s/m2+e+s-sb-d/m+bs-d>0, model (4a) and (4b) has only one interior equilibrium E4n4,p4. Furthermore, E4(n4,p4) is locally asymptotically stable when n4∈L, where L=n∣An<0,n∈1/m,1 and(14a)n4=s+e-sb-d-sb-s+d-e2+4bss-d2s,(14b)p4=1-n4mn4-1b+n4.Proof.
According to (6), (15)trace⁡Jn4,p4=An4=-3mn42+2m+2n4-1-b1-n4mn4-1b+n4,det⁡Jn4,p4=-Cn4Bn4=n4b+n4bep4b+n42-sp4.
Imitating the proof of Theorem2, we can know(16)n4<bes-b;then, detJn4,p4=-Bn4Cn4>0, which means two eigenvalues have the same sign.
Since(17)A1m=1-1m>0,A1=1-m<0.Solving A(n)=0, we can get that there exists a set L⊂1/m,1, such that A(n)<0 for all n∈L, where L=1+m-mb+1-m+mb+m21+b+b2/3m,1. Then, we know that if n4∈L, traceJn4,p4=An4<0 and both eigenvalues of J(n4,p4) have negative real parts, and E4(n4,p4) is locally asymptotically stable.
From the above discussion, we find that model (4a) and (4b) may possess two interior equilibria E3 and E4 simultaneously if e-d-bd<0,  -s/m2+e+s-sb-s/m+bs-d<0,   sb-s+d-e2+4bss-d>0, and 2s/m<e+s-sb-d<2s. The stability condition of E4 is still n4∈L.
Considering the existence and stability conditions integrated, we know that there will appear three stable equilibria synchronously if the scope of those parameters satisfies the above conditions. Based on the above demonstration, some cases where the stability of equilibria may occur are cited in Table1, where the existence of equilibria is classified according to the magnitude of s and d.Table 1
E
0
0,0
E
1
1,0
E
2
1
/
m
,
0
E
3
n
3
,
p
3
E
4
n
4
,
p
4
s
>
d
Saddle point
Stable node point
Unstable node point
Saddle point
Nonexistent
s
=
d
Unstable (high-order singularity)
Stable node point
Unstable node point
Saddle point
Nonexistent
s
<
d
Δ
=
0
Stable node point
Stable node point
Saddle point
E
3 and E4 coincide and are a saddle-node point
Δ
>
0
Stable node point
Stable node point
Unstable node point
Saddle point
Nonexistent
Stable node point
Saddle point
Saddle point
Nonexistent
Stable node point
Stable node point
Stable node point
Saddle point
Saddle point
Stable node point


### 2.2. Local Bifurcation

From Table1 we know that the variation of parameter value will lead to the number change of interior equilibria. We take d as variable parameter and find that changing the value of d will vary equilibrium’s number; when the value of d increases across the threshold dTC1=e/b+1, the interior equilibrium E4 bifurcates from E1, and when the value of d is across dTC2=e/bm+1+s1-1/m, another interior equilibrium E3 can bifurcate from E2. Then, there exist two transcritical bifurcations, which are denoted by TC1 and TC2, respectively.Theorem 4.
(1) Model (4a) and (4b) undergoes transcritical bifurcation at E11,0 when the value of parameter d equals the transcritical bifurcation threshold dTC1=e/b+1.
(2) Model (4a) and (4b) undergoes another transcritical bifurcation at E21/m,0 when the value of parameter d equals the transcritical bifurcation threshold dTC2=e/bm+1+s1-1/m.Proof.
It can be easily seen thatE11,0 coincides with E4n4,p4 when d=e/b+1. According to the theorems in [33, 34], it can be found that one interior equilibrium point branches off from E1 when d passes the threshold dTC1=e/b+1 and also conforms with the transversality condition for transcritical bifurcation. The same notation we followed in this paper is mentioned in [33].
We can calculate the value of Jacobian matrix of model (4a) and (4b) as DetJ1,0dTC1, and then 1,0 is a nonhyperbolic equilibrium point when d=dTC1. Let the eigenvectors v=1,b+11-mT and w=0,1T indicate the eigenvectors corresponding to zero eigenvalues of J1,0,dTC1 and J1,0,dTC1T, respectively. Next the transversality conditions for the transcritical bifurcation are satisfied be verified,where F=F1,F2T:(18)wTFd1,0;dTC1=0,10,0T=0,wTDFd1,0;dTC1v=b+1m-1≠0,wTD2F1,0;dTC1v,v=eb1+b2-sb+11-m.
We notice thatm>1 and wTD2F1,0;dTC1v,v<0 if e>sb+12/b and then the transcritical bifurcation is supercritical, which means that an interior equilibrium point arises through E1 under this condition. Another aspect is that if e<sb+12/b, the transcritical bifurcation is subcritical and an interior equilibrium vanishes across E1.
In the following, we demonstrate that another interior equilibrium point bifurcates fromE21/m,0 through transcritical bifurcation at the threshold dTC2=e/bm+1+s1-1/m. Similarly, we can get DetJ1/m,0dTC2=0, v=1,b+1/mm-1T, and w=0,1Tare the eigenvectors corresponding to J1/m,0,dTC2 and J1/m,0,dTC2T, respectively. We have (19)wTFd1m,0;dTC2=0,10,0T=0,wTDFd1m,0;dTC2v=1-mb+1m≠0,wTD2F1m,0;dTC2v,v=beb+1/m2-sb+1mm-1.
Ife<sb+1/m2/b and wTD2F1/m,0;dTC2v,v<0, the transcritical bifurcation is supercritical and an interior equilibrium point appears across E2 under this condition. Another aspect is that if e>sb+1/m2/b, the transcritical bifurcation is subcritical and an interior equilibrium vanishes across E2.Theorem 5.
Model (4a) and (4b) undergoes saddle-node bifurcation when d=dSN, where dSN=s+e+bs-2bse.Proof.
It is easy to calculate that the discriminant ofGn=0 is equal to zero when d=s+e+bs-2bse, which means Gn=0 has a twofold root. We denote this root by n∗. Model (4a) and (4b) has only one interior equilibrium point E∗n∗,p∗ correspondingly and the constituents are given by (20)n∗=ebs-b,p∗=1-n∗mn∗-1b+n∗.
The Jacobian matrix evaluate atE∗ is given by (21)Jn∗,p∗,dSN=An∗,p∗,dSNBn∗,p∗,dSNCn∗,p∗,dSNDn∗,p∗,dSN=Mbse-100,where M=A(n∗,p∗,dSN). Obviously, matrix Jn∗,p∗,dSN has a double zero eigenvalues when DetJn∗,p∗,dSN=0. We can count out the eigenvectors corresponding to the zero eigenvalue, which are v=1,M/1-bs/eT and w=0,1T. Utilizing the expressions of vectors v and w, as well as n∗∈1/m,1, we can get (22)wTFdn∗,p∗;dSN=-1-n∗mn∗-1b+n∗<0,wTD2Fn∗,p∗;dSNv,v=-eb1-n∗mn∗-1b+n∗2.
Given the value range ofn∗, we know that wTD2Fn∗,p∗;dSNv,v cannot be equal to zero. Then, by Sotomayor’s theorem we can prove that the model undergoes a saddle-node bifurcation when the parameter d goes via the critical threshold d=dSN.Theorem 6.
In the case ofs<d, the interior equilibrium point E4 changes its stability through the Hopf-bifurcation at the threshold b=bH, where bH=-3mn42-2m+1n4+1/2mn4-m-1.Proof.
Because the interior equilibrium pointE3 is always a saddle, then Hopf bifurcation can only take place at E4. Parameter b can drive equilibrium E4 into an unstable state when b>bH, so b=bH is the critical value where the stability of E4 changes. Next we will prove the necessary condition for Hopf bifurcation to occur. To testify the transversality condition of the Hopf bifurcation of the model’s solution, we take λ=αb+βbi, where λ is an eigenvalue of the Jacobian matrix J(n4,p4,bH). If the Hopf bifurcation occurs at b=bH, λ is a purely imaginary number, such that TrJ(n4,p4,bH)=0, DetJ(n4,p4,bH)≠0, and dTrJ(n4,p4,bH)/db≠0. Substituting b by bH=-3mn42-2m+1n4+1/2mn4-m-1, we can get (23)Tr⁡Jn4,p4,bH=0,DetJn4,p4,bH=n41-n4mn4-1bHebH+n42-s,ddbTr⁡Jn4,p4,bH=-n41-n4mn4-1b+n42.
Using (16), we know DetJ(n4,p4,bH)>0 and dTrJ(n4,p4,bH)/db≠0. Then, the transversality condition of a Hopf bifurcation is satisfied [35].Theorem 7.
Model (4a) and (4b) undergoes a Bogdanov-Takens (BT) bifurcation of codimension two.
As we know,TrJ=0 and DetJ=0 are the necessary conditions of the occurrence of Hopf and saddle-node, respectively. If both Hopf and saddle-node conditions hold, there will be a new bifurcation called Bogdanov-Takens bifurcation. In this situation, the Jacobian matrix has a double zero eigenvalue [36]. Since the explicitly analytical expressions for thresholds of BT bifurcation are quite difficult to determine, we give a numerical example to confirm the system exhibit BTs bifurcation. We fix m=1.5, s=0.1, and d=0.4. We find that, at eBT,bBT=0.4930734388,0.2386077213, TrJE∗eBT,bBT=0 and DetJE∗eBT,bBT=0. Moreover, ∂2F/∂n2-∂F/∂n∂2F/∂n∂p/∂2F/∂n∂p+∂2F/∂n∂p=-2.800482957 and 1/2∂F/∂n∂2F/∂n2-∂F/∂n2(∂2F/∂n∂p)/∂F/∂p+1/2∂F/∂p∂2F/∂n2-∂F/∂n(∂2F/∂n∂p)=-0.000073428878; these expressions prove the transversality conditions for a BT bifurcation [37].


## 2.1. Equilibria

In this subsection, we mainly concentrate on the existence of positive equilibrium of model (4a) and (4b). Model (4a) and (4b) has three boundary equilibria: E00,0, E11,0, and E21/m,0 and they are always existent. The interior equilibria are the intersection points of the vertical isocline and horizontal isocline in the interior of the first quadrant. The expressions of vertical and horizontal isoclines are as follows:(5a)1-nmn-1-pb+n=0,(5b)enb+n+s1-n-d=0.The horizontal isocline is vertical line and the number of perpendiculars is decided by (5b). We denote (5b) by H(n)=Gn/b+n=0, where G(n)=-sn2+e+s-sb-dn+b(s-d). G(n)=0 have the same solutions as (5b), and two solutions exist in the equation at most. We denoted them by n3 and n4, where n4<n3. From the geometric properties of isoclines we can know that model (4a) and (4b) has interior equilibrium points if and only if 1/m<n3<1 or 1/m<n4<1. Then, the maximum number of equilibria for model (4a) and (4b) is five. We use E3n3,p3 and E4n4,p4 to signify the interior equilibria, the existence and stability of which are shown as follows.The Jacobian matrix of model (4a) and (4b) at the equilibrium point Ei is(6)Jni,pi=Ani,piBni,piCni,piDni,pi,(7)Ani,pi=-3mni2+2m+2ni-1-bpib+ni2,Bni,pi=-nib+ni,Cni,pi=bepib+ni2-spi,Dni,pi=enib+ni+s-d-sni.Theorem 1.
(1)E00,0 is an asymptotically stable node point if s<d and a saddle point for s>d; if s=d, it is a high order singularity. The Jacobian matrix around E0 shows that(8)J0,0=-100s-d.(2) E1(1,0) is a asymptotically stable node point when e/b+1<d while it is a saddle point when e/b+1>d. The Jacobian matrix around E1 is (9)J1,0=1-m-1b+10eb+1-d.(3) E21/m,0 is always unstable. The Jacobian matrix around E2 is(10)J1/m,0=1-1m-1bm+10ebm+1+s-d-sm.One of the eigenvalues is 1-1/m>0, so E2 is unstable.Theorem 2.
Ife-d-bd>0 and -s/m2+e+s-sb-d/m+bs-d<0 hold, model (4a) and (4b) has only one interior equilibrium E3n3,p3. Furthermore, E3(n3,p3) is always unstable, especially when d=s+e+sb-2bse; it is a high-order singularity, where(11a)n3=s+e-sb-d+sb-s+d-e2+4bss-d2s,(11b)p3=1-n3mn3-1b+n3.Proof.
According to (6),(12)det⁡Jn3,p3=-Cn3,p3Bn3,p3=n3p3b+n3beb+n32-s.Let ϕ=n3-be/s+b; substituting (11a) into ϕ, we get(13)ϕ=s+e-d+sb-2bes+s+e-d-sb2+4bss-d2s=s+e-d+sb-2bes+s+e-d+sb2-4bse2s≥0.
Then, we know whenn3≥be/s-b and detJn3,p3≤0, E3 is a saddle point or a high-order singularity.Theorem 3.
Ife-d-bd<0 and -s/m2+e+s-sb-d/m+bs-d>0, model (4a) and (4b) has only one interior equilibrium E4n4,p4. Furthermore, E4(n4,p4) is locally asymptotically stable when n4∈L, where L=n∣An<0,n∈1/m,1 and(14a)n4=s+e-sb-d-sb-s+d-e2+4bss-d2s,(14b)p4=1-n4mn4-1b+n4.Proof.
According to (6), (15)trace⁡Jn4,p4=An4=-3mn42+2m+2n4-1-b1-n4mn4-1b+n4,det⁡Jn4,p4=-Cn4Bn4=n4b+n4bep4b+n42-sp4.
Imitating the proof of Theorem2, we can know(16)n4<bes-b;then, detJn4,p4=-Bn4Cn4>0, which means two eigenvalues have the same sign.
Since(17)A1m=1-1m>0,A1=1-m<0.Solving A(n)=0, we can get that there exists a set L⊂1/m,1, such that A(n)<0 for all n∈L, where L=1+m-mb+1-m+mb+m21+b+b2/3m,1. Then, we know that if n4∈L, traceJn4,p4=An4<0 and both eigenvalues of J(n4,p4) have negative real parts, and E4(n4,p4) is locally asymptotically stable.
From the above discussion, we find that model (4a) and (4b) may possess two interior equilibria E3 and E4 simultaneously if e-d-bd<0,  -s/m2+e+s-sb-s/m+bs-d<0,   sb-s+d-e2+4bss-d>0, and 2s/m<e+s-sb-d<2s. The stability condition of E4 is still n4∈L.
Considering the existence and stability conditions integrated, we know that there will appear three stable equilibria synchronously if the scope of those parameters satisfies the above conditions. Based on the above demonstration, some cases where the stability of equilibria may occur are cited in Table1, where the existence of equilibria is classified according to the magnitude of s and d.Table 1
E
0
0,0
E
1
1,0
E
2
1
/
m
,
0
E
3
n
3
,
p
3
E
4
n
4
,
p
4
s
>
d
Saddle point
Stable node point
Unstable node point
Saddle point
Nonexistent
s
=
d
Unstable (high-order singularity)
Stable node point
Unstable node point
Saddle point
Nonexistent
s
<
d
Δ
=
0
Stable node point
Stable node point
Saddle point
E
3 and E4 coincide and are a saddle-node point
Δ
>
0
Stable node point
Stable node point
Unstable node point
Saddle point
Nonexistent
Stable node point
Saddle point
Saddle point
Nonexistent
Stable node point
Stable node point
Stable node point
Saddle point
Saddle point
Stable node point


## 2.2. Local Bifurcation

From Table1 we know that the variation of parameter value will lead to the number change of interior equilibria. We take d as variable parameter and find that changing the value of d will vary equilibrium’s number; when the value of d increases across the threshold dTC1=e/b+1, the interior equilibrium E4 bifurcates from E1, and when the value of d is across dTC2=e/bm+1+s1-1/m, another interior equilibrium E3 can bifurcate from E2. Then, there exist two transcritical bifurcations, which are denoted by TC1 and TC2, respectively.Theorem 4.
(1) Model (4a) and (4b) undergoes transcritical bifurcation at E11,0 when the value of parameter d equals the transcritical bifurcation threshold dTC1=e/b+1.
(2) Model (4a) and (4b) undergoes another transcritical bifurcation at E21/m,0 when the value of parameter d equals the transcritical bifurcation threshold dTC2=e/bm+1+s1-1/m.Proof.
It can be easily seen thatE11,0 coincides with E4n4,p4 when d=e/b+1. According to the theorems in [33, 34], it can be found that one interior equilibrium point branches off from E1 when d passes the threshold dTC1=e/b+1 and also conforms with the transversality condition for transcritical bifurcation. The same notation we followed in this paper is mentioned in [33].
We can calculate the value of Jacobian matrix of model (4a) and (4b) as DetJ1,0dTC1, and then 1,0 is a nonhyperbolic equilibrium point when d=dTC1. Let the eigenvectors v=1,b+11-mT and w=0,1T indicate the eigenvectors corresponding to zero eigenvalues of J1,0,dTC1 and J1,0,dTC1T, respectively. Next the transversality conditions for the transcritical bifurcation are satisfied be verified,where F=F1,F2T:(18)wTFd1,0;dTC1=0,10,0T=0,wTDFd1,0;dTC1v=b+1m-1≠0,wTD2F1,0;dTC1v,v=eb1+b2-sb+11-m.
We notice thatm>1 and wTD2F1,0;dTC1v,v<0 if e>sb+12/b and then the transcritical bifurcation is supercritical, which means that an interior equilibrium point arises through E1 under this condition. Another aspect is that if e<sb+12/b, the transcritical bifurcation is subcritical and an interior equilibrium vanishes across E1.
In the following, we demonstrate that another interior equilibrium point bifurcates fromE21/m,0 through transcritical bifurcation at the threshold dTC2=e/bm+1+s1-1/m. Similarly, we can get DetJ1/m,0dTC2=0, v=1,b+1/mm-1T, and w=0,1Tare the eigenvectors corresponding to J1/m,0,dTC2 and J1/m,0,dTC2T, respectively. We have (19)wTFd1m,0;dTC2=0,10,0T=0,wTDFd1m,0;dTC2v=1-mb+1m≠0,wTD2F1m,0;dTC2v,v=beb+1/m2-sb+1mm-1.
Ife<sb+1/m2/b and wTD2F1/m,0;dTC2v,v<0, the transcritical bifurcation is supercritical and an interior equilibrium point appears across E2 under this condition. Another aspect is that if e>sb+1/m2/b, the transcritical bifurcation is subcritical and an interior equilibrium vanishes across E2.Theorem 5.
Model (4a) and (4b) undergoes saddle-node bifurcation when d=dSN, where dSN=s+e+bs-2bse.Proof.
It is easy to calculate that the discriminant ofGn=0 is equal to zero when d=s+e+bs-2bse, which means Gn=0 has a twofold root. We denote this root by n∗. Model (4a) and (4b) has only one interior equilibrium point E∗n∗,p∗ correspondingly and the constituents are given by (20)n∗=ebs-b,p∗=1-n∗mn∗-1b+n∗.
The Jacobian matrix evaluate atE∗ is given by (21)Jn∗,p∗,dSN=An∗,p∗,dSNBn∗,p∗,dSNCn∗,p∗,dSNDn∗,p∗,dSN=Mbse-100,where M=A(n∗,p∗,dSN). Obviously, matrix Jn∗,p∗,dSN has a double zero eigenvalues when DetJn∗,p∗,dSN=0. We can count out the eigenvectors corresponding to the zero eigenvalue, which are v=1,M/1-bs/eT and w=0,1T. Utilizing the expressions of vectors v and w, as well as n∗∈1/m,1, we can get (22)wTFdn∗,p∗;dSN=-1-n∗mn∗-1b+n∗<0,wTD2Fn∗,p∗;dSNv,v=-eb1-n∗mn∗-1b+n∗2.
Given the value range ofn∗, we know that wTD2Fn∗,p∗;dSNv,v cannot be equal to zero. Then, by Sotomayor’s theorem we can prove that the model undergoes a saddle-node bifurcation when the parameter d goes via the critical threshold d=dSN.Theorem 6.
In the case ofs<d, the interior equilibrium point E4 changes its stability through the Hopf-bifurcation at the threshold b=bH, where bH=-3mn42-2m+1n4+1/2mn4-m-1.Proof.
Because the interior equilibrium pointE3 is always a saddle, then Hopf bifurcation can only take place at E4. Parameter b can drive equilibrium E4 into an unstable state when b>bH, so b=bH is the critical value where the stability of E4 changes. Next we will prove the necessary condition for Hopf bifurcation to occur. To testify the transversality condition of the Hopf bifurcation of the model’s solution, we take λ=αb+βbi, where λ is an eigenvalue of the Jacobian matrix J(n4,p4,bH). If the Hopf bifurcation occurs at b=bH, λ is a purely imaginary number, such that TrJ(n4,p4,bH)=0, DetJ(n4,p4,bH)≠0, and dTrJ(n4,p4,bH)/db≠0. Substituting b by bH=-3mn42-2m+1n4+1/2mn4-m-1, we can get (23)Tr⁡Jn4,p4,bH=0,DetJn4,p4,bH=n41-n4mn4-1bHebH+n42-s,ddbTr⁡Jn4,p4,bH=-n41-n4mn4-1b+n42.
Using (16), we know DetJ(n4,p4,bH)>0 and dTrJ(n4,p4,bH)/db≠0. Then, the transversality condition of a Hopf bifurcation is satisfied [35].Theorem 7.
Model (4a) and (4b) undergoes a Bogdanov-Takens (BT) bifurcation of codimension two.
As we know,TrJ=0 and DetJ=0 are the necessary conditions of the occurrence of Hopf and saddle-node, respectively. If both Hopf and saddle-node conditions hold, there will be a new bifurcation called Bogdanov-Takens bifurcation. In this situation, the Jacobian matrix has a double zero eigenvalue [36]. Since the explicitly analytical expressions for thresholds of BT bifurcation are quite difficult to determine, we give a numerical example to confirm the system exhibit BTs bifurcation. We fix m=1.5, s=0.1, and d=0.4. We find that, at eBT,bBT=0.4930734388,0.2386077213, TrJE∗eBT,bBT=0 and DetJE∗eBT,bBT=0. Moreover, ∂2F/∂n2-∂F/∂n∂2F/∂n∂p/∂2F/∂n∂p+∂2F/∂n∂p=-2.800482957 and 1/2∂F/∂n∂2F/∂n2-∂F/∂n2(∂2F/∂n∂p)/∂F/∂p+1/2∂F/∂p∂2F/∂n2-∂F/∂n(∂2F/∂n∂p)=-0.000073428878; these expressions prove the transversality conditions for a BT bifurcation [37].


## 3. Numerical Results

In order to verify the correctness and feasibility of the theoretical results, a series of numerical simulations will be depicted in detail. Many phase diagrams are given to display the dynamics properties of model (4a) and (4b), which are based on pplane8 routines [38] in Matlab 7.1. In fact, according to [39], we can know pplane8 is a powerful tool for studying planar autonomous systems of differential equations, which can rapidly and accurately draw trajectories of each phase plane, count each critical point, and correctly characterize each equilibrium point of the studied systems. It is easy to find from Table 1 that model (4a) and (4b) only has an interior equilibrium point E3 if s≥d; then, the premise of parametric ranges in Figure 1 is s>d. It should be pointed out from Figure 1(a) that two vertical lines passing through points e/b+1,0 and e/bm+1+s1-s/m,0 are the transcritical bifurcation curves, which have been named TC1 and TC2, respectively. Furthermore, if d<e/b+1 is feasible, that is to say, model (4a) and (4b) does not have any interior equilibrium point when the value of d is positioned within region I of Figure 1(a), three boundary equilibria E0(0,0), E1(1,0), and E2(2/3,0) exist.  E0(0,0) is asymptotically stable, while E1(1,0) and E2(2/3,0) are unstable, the results of which have been shown in Figure 2(a) (d=0.131). We know those three boundary equilibria E0(0,0), E1(1,0), and E2(2/3,0) always exist no matter what the value of d is. Thus, in the following cases we do not introduce the existence of boundary equilibria separately. However, if the value of d is gradually increasing across the transcritical bifurcation TC1 and finally enters into region II, model (4a) and (4b) has an interior equilibrium point E3. It is worth emphasizing that the critical value of transcritical bifurcation is dTC1=e/b+1=0.13167 and the interior equilibrium point E30.97064,0.039763 is a saddle point, boundary equilibria E0 and E1 are stable, and E2 is an unstable node point, which has been shown in Figure 2(b) (d=0.132). As the value of d gradually increases to dTC2=e/bm+1+s1-1/m=0.132083, which is a critical value of another transcritical bifurcation and finally enters into region III, model (4a) and (4b) has two interior equilibria E30.91328,0.093458 and E40.71172,0.052831, which are unstable; E0 and E1 are stable, and E2 is a saddle point (see Figure 2(c), d=0.1325).Figure 1
Bifurcation diagram of model (4a) and (4b) in n-d plane constructed for m=1.5, b=2.0, s=0.1, e=0.395 (a), and e=0.405 (b, c). Red perpendicular dashed lines represent the critical value for different bifurcation occurrence. Three horizontal lines stand for three boundary equilibria: E0 (golden), E1 (gray), and E2 (black). Two interior equilibria E3 and E4 are presented by blue and green curves, respectively. The solid lines indicate equilibrium in stable state and dotted line indicates unstable state, where, according to the stability theorem, these results are calculated and drawn based on Maple 14 platform.
(a)
(b)
(c)Figure 2
Phase portraits for different concomitant case of the model’s equilibria. The horizontal axis is prey populationn and the vertical axis is predator population p. The green curves are stable or unstable orbits; the red point is the equilibrium point. We take m=1.5, b=2.0, s=0.1, e=0.395 in (a–c), and e=0.405 in (d–h).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)To further fully explore the existence of the equilibria and the occurrence of the bifurcation behavior in model (4a) and (4b), we take another set of parameter values in Figure 1(b) and repeat the steps in Figure 1(a). The bifurcation diagram in Figure 1(b) has been depicted in detail as follows. It is obvious to find that model (4a) and (4b) has an interior equilibrium point E4 if the value of d exceeds the threshold value dTC3=e/bm+1+s1-1/m=0.134583, which is transcritical bifurcation threshold TC3. It is worthwhile to point out that the interior equilibrium point E40.69128,0.030678 is unstable on the account of n4=0.69128∉L0.83822,1. The boundary equilibrium E0(0,0) is stable, and E1(1,0) and E22/3,0 are saddle points (see Figure 2(d), d=0.1349). We take another set of parameters values as the supplement to display the case that there exists only one interior equilibrium point and it is locally asymptotically stable (see Figure 3(a)), where E1(1,0) and E2(2/3,0) are saddle points, E0(0,0) and E40.886762,0.040628 are stable, and L=0.84652,1. Furthermore, if the value of d increases beyond the threshold value of transcritical bifurcation threshold dTC4=e/b+1=0.135 and finally enters into region V, te model (4a) and (4b) has three boundary equilibria E0(0,0), E1(1,0), and E22/3,0 and two interior equilibria E30.97862,0.029801 and E40.71938,0.060341 (see Figure 2(e), d=0.1352). It is necessary to underline that E00,0 and E11,0 are stable and n4=0.71938∉L0.83822,1, so the interior equilibria E4 and E3 are all unstable.Figure 3
Phase portraits form=1.5, s=0.1, and b=0.2; (a) e=0.844, d=0.7; (b) e=0.599, d=0.5.
(a)
(b)In order to clearly explore the steady characteristic of the interior equilibrium pointE4, Figure 1(c) will be given, which is the partially enlarged view of Figure 1(b). It should be stressed that the interior equilibrium E4 will change the stable state if the value of d increased beyond the line L1 and finally enters into region VI, which suggests that a Hopf bifurcation can lead to the appearance of a limit cycle in the vicinity of the interior equilibrium E40.83816,0.11816 if d=0.135787825659 (see Figure 2(f)). From Figure 2(f), we can observe that this limit cycle around E4 is stable as it attracts two neighboring trajectories: the trajectory (bottle green curve) lying inside the limit cycle and the trajectory (blue curve) lying outside; these two trajectories move ectad and entad, respectively, and converge on the limit cycle. But at this time the interior equilibrium point E30.85397,0.11709 and boundary equilibrium point E22/3,0 are unstable, and E00,0 and E11,0 are stable. When the value of d enters into domain VI, E40.86298,0.042889 becomes stable, E30.92702,0.032122 remains unstable, and E0 and E1 are all stable in this domain. Thus, it is interesting to know that model (4a) and (4b) will show three stationary phenomenon (see Figure 2(g), d=0.13579). Due to the parameter values, the nature of E3 and E4 only can be seen clearly on the enlarged view (see Figure 2(g)). Therefore, we take another group of values to exhibit the three stable states in overall view in Figure 3(b). If the value of d gradually increases and reaches dSN=bs+e+2bse+s=0.1357900212, it is easy to find that E30.84605,0.11789 and E40.84605,0.11789 coincide at line L2 and the coincident point is called saddle-node point, which has features of both saddle and node points (see Figure 2(h)). However, it will disappear if the value of d is increased higher than dSN, which signifies that model (4a) and (4b) will undergo a saddle-node bifurcation if the value of d increases across the threshold value dSN=bs+e+2bse+s. In a word, it is worthy of our summary that model (4a) and (4b) can show a complex dynamic evolutionary process of steady conversion and bifurcation behavior with increase of key parameter d.In addition, model (4a) and (4b) has three boundary equilibria E0(0,0), E1(1,0), and E22/3,0 and an interior equilibrium point E3, which is unstable if s>d or s=d. However, it is worth stressing that the boundary equilibrium point E0(0,0) is a saddle point if s>d and is an unstable high-order singularity if s=d (see Figures 4(a) and 4(b)). At the same time, it can be known from Theorem 7 that model (4a) and (4b) can undergo a Bogdanov-Takens bifurcation if TrJE∗eBT,bBT=0 and DetJE∗eBT,bBT=0, which is shown in Figure 4(c).Figure 4
Phase portraits form=1.5, s=0.1, e=0.28, b=2.0, (a) d=0.1, and (b) d=0.095. The parameter values in image (c) are introduced in Theorem 7.
(a)
(b)
(c)Based on the above analysis, the key parameterd can impose influence on dynamic evolutionary mechanism of steady conversion and bifurcation behavior and lead to model (4a) and (4b) having three stationary phenomena and multiple bifurcation behaviors, which can in turn prove that the theoretical results are correct and the complex dynamics of model (4a) and (4b) mainly depend on some key parameters. Moreover, these results show that the method of using mathematical model to study the ecological problems is feasible.


## 4. Conclusions

On the basis of the theories and methods of ecology, a predator-prey model is studied numerically and analytically in this paper. The aim is to probe how some key factors influence dynamic evolutionary mechanism of steady conversion and bifurcation behavior in predator-prey model. The theoretical works have been promoting the investigation of the existence and stability of the equilibria, as well as some conditions of some bifurcations behaviors, such as transcritical bifurcation, saddle-node bifurcation, and Hopf bifurcation, which can deduce a standard parameter controlled relationship and in turn provide a theoretical basis for the numerical simulation.Numerical analysis indicates that the dynamic evolutionary mechanism of steady conversion and bifurcation behavior mainly depend on a specific key parameterd. Within this framework, the direct and indirect effects caused by the specific key parameter d are investigated by means of bifurcation analysis and phase diagram. It is obvious to find that the existence and stability of interior equilibria E3 and E4 mainly depend on a key parameter d. These results suggest that the key parameter d plays an important role in the prey-predator model. In addition, when the value of the key parameter d is larger than some critical value, the model can possess multiple bifurcation behaviors, such as transcritical bifurcation, saddle-node bifurcation, and Hopf bifurcation. Thus, it is worthwhile to remark that the key parameter d has a profound effect on the population bifurcation dynamical behaviors. In a word, some key parameters can alter population dynamics and features in prey-predator model. In addition, it is our hope that all these results can be applied in the study of the dynamic complexity of ecosystems.



---

*Source: 3565316-2016-03-07.xml*