# Complexity Analysis of a Modified Predator-Prey System with Beddington–DeAngelis Functional Response and Allee-Like Effect on Predator

Authors: Shuangte Wang; Hengguo Yu
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

In this paper, complex dynamical behaviors of a predator-prey system with the Beddington–DeAngelis functional response and the Allee-like effect on predator were studied by qualitative analysis and numerical simulations. Theoretical derivations have given some sufficient and threshold conditions to guarantee the occurrence of transcritical, saddle-node, pitchfork, and nondegenerate Hopf bifurcations. Computer simulations have verified the feasibility and effectiveness of the theoretical results. In short, we hope that these works could provide a theoretical basis for future research of complexity in more predator-prey ecosystems.


---

## Body

## 1. Introduction

In reference [1], the authors simply considered a predator-prey system with Holling type II functional response and Allee-like effect on predator, which is described by the following nonlinear ordinary differential equations (ODEs):(1a)x˙=r1x1−xK1−qxya+x−m1x−dx2,(1b)y˙=r2y1−yK2yy+e+e1qxya+x−m2y,subject to initial conditions x0,y0≥0. Here we replace the Holling type II functional response qx/a+x with a functional response qx/a+bx+cy and denote the parameter q as q1 for later use, and thus above system (1a) and system (1b) have a modified version:(2a)x˙=r1x1−xK1−q1xya+bx+cy−m1x−dx2,(2b)y˙=r2y1−yK2yy+e+e1q1xya+bx+cy−m2y,where functions x=xt and y=yt are the densities of prey and predator at time t, respectively. In terms of biology, all above positive constants have practical considerations. Parameters r1 and r2 denote the intrinsic growth rate of the prey and predator, respectively; K1 and K2 represent the carrying capacity of the environment; a is the half-saturation constant; q1 is the search efficiency of predator for prey; m1 and m2 are the mortality rate of the prey and predator species, respectively; e1 is the biomass conversion and we denote e1q1 as q2 for convenience; d is the intraspecific competition coefficient of the prey; e is the Allee effect constant. The specific growth term r1x1−x/K1 governs the increase of prey in the lack of predator, while the specific growth term r2y1−y/K2 governs the increase of predator. The square term dx2 is an intrinsic decrease term on prey. The term y/y+e on the specific growth term of predator with multiply form is called the Allee-like effect [2, 3] and is different from predator-prey systems in [4–8] with Allee effect on prey. The coupled term q1x/a+bx+cy is named Beddington–DeAngelis (B-D) functional response (named after Beddington and DeAngelis et al.) [9, 10]. It is similar to the Holling type II functional response incorporating an extra term cy in denominator, which describes mutual interference among predators [11, 12]. This functional response has some of the same qualitative features as the ratio-dependent form but avoids some of singular behaviors of ratio-dependent models at low densities [11]. Obviously, when b=1 and c=0, system (2a) and system (2b) reduce to original system (1a) and system (1b).When parametersd=m1=e=0 (without Allee effect) and r2−m2<0, system (2a) and system (2b) reduce to System (2.1) in [13] and the author particularly conducted stability (local and global) and bifurcation (saddle-node, transcritical, Hopf–Andronov, and Bogdanov–Takens) analysis with a detailed mathematical analysis. When d=m1=r2=0, system (2a) and system (2b) become model system (3a) and model system (3b) in [14], which was also independently and originally proposed in [9–11], while in [14], the authors discussed local and global asymptotic stability behavior of various equilibria and Hopf bifurcation occurs when parameter m corresponding to reserved region crosses some critical values. To mimic the real-world scenario, they solved the inverse problem of estimation of system parameter m by using the sampled data. System (1.3) in [15] is similar to above system in [14] except the constant rate harvesting term. In this reference, the authors showed that it undergoes several kinds of bifurcations, such as the saddle-node bifurcation, the subcritical and supercritical Hopf bifurcation, and Bogdanov–Takens bifurcation by choosing the death rate of the predator and the harvesting rate of the prey as bifurcation parameters.Motivated by previous progress of predator-prey systems with B-D functional response or Allee-like effect, this paper mainly concentrates on dynamical analysis of system (2a) and system (2b). The rest of this paper is structured as follows. Preliminaries, such as boundedness, permanence, and existence of trivial equilibria, are given in Section 2. The existence of interior equilibrium is presented in Section 3 by virtues of the cobweb model and polynomial equations, respectively. In Section 4, we give stability analysis of equilibria and nonexistence of limit cycles. In Section 5, local codimension one bifurcations are analyzed, especially the Hopf bifurcation incorporating numerical simulations and Hopf bifurcation curves. In Section 6, we carry out short conclusions for our system.


## 2. Preliminaries

In this section, we devote to give some priori foundations. Before presenting the main results, we denote the first quadrant asℝ+2, and its closure is ℝ+2=ℝ+2¯. For biological consideration, system (2a) and system (2b) are defined on the domain ℝ+2 and all the solutions are non-negative with initial conditions x0,y0≥0, i.e., ℝ+2 is an invariant set and any orbits starting from it cannot cross the coordinate axes. Furthermore, all solutions are uniformly bounded. But now we only need to prove following theorems.


### 2.1. Boundedness and Permanence

Theorem 1 (uniform boundedness).
Suppose that a non-negative functionφx,y and its partial derivatives φx and φy are continuous in ℝ+2, then the system(3a)x˙=r1x1−xK1−φx,y−m1x−dx2,(3b)y˙=r2y1−yK2yy+e+e1φx,y−m2y,subject to x0,y0≥0 is uniformly bounded.This Theorem1 holds obviously after introducing an auxiliary function z=e1x+y [1, 16]. The following theorem with the help of comparison principle in ODEs is about permanence of system (2a) and system (2b).Theorem 2 (permanence).
If parameters satisfy(4)ω1=r1−m1−q1M2/a+cM2r1/K1+d>0,a0=q2ω11−λa+bω11−λ−m2>0,a1=r2e−cq2ω11−λa+bω11−λ2<0,cρ≤1,e≥ρa+bω11−λ,where M2 is a positive constant, λ∈0,1, and ρ=r2eK2+e/K2c2q2ω11−λ3, then system (2a) and system (2b) are permanent.Proof.
From above Theorem1, we have a positive upper bound ξ1 such that(5)maxlimsupt⟶∞x,limsupt⟶∞y≤ξ1.
Then, we obtainM2>0 and sufficiently large T1, such that(6)x˙≥xr1−m1−q1M2a+cM2−r1K1+dx,∀t≥T1.
By using the lemmas in [17], we admit liminft⟶∞x≥ω1. That is to say, there exist sufficiently large T2, such that x≥ω11−λ, ∀t≥T2. Then, we consider equation (2b) again. It yields y˙≥fyy; here the function fy incorporating its Taylor expression is(7)fy=r2yy+e1−yK2−m2+q2ω11−λa+bω11−λ+cy=a0+a1y+12f″ξyy2,ξy∈0,y,in which the second order derivative of fy is(8)f″y=−2r2eK2+eK2y+e3+2c2q2ω11−λa+bω11−λ+cy3.
From the conditions in this theorem, we know thatf″y≥0, when y≥0. Thus, fy≥a0+a1y, ∀y≥0. By using the lemmas in [17] again, the proof is completed.


### 2.2. Existence of Trivial Equilibria

In this section, we will discuss the existence conditions of trivial equilibria of system (2a) and system (2b). Firstly, from [1], it is obvious that system (2a) and system (2b) have trivial equilibria: E0≔0,0, E1≔x1,0, and E2k≔0,yk, k=1,2, where x1≔r1−m1/r1/K1+d and(9)y1,2≔r2−m2±Δ2r2/K2,Δ≔r2−m22−4r2m2eK2,The pointE1 exists when r1>m1, while the points E2k all exist if r2>m2 and Δ>0. If r2>m2 but Δ=0, then the two equilibria E21 and E22 collide with each other and we denote this equilibrium as E2=0,y2.


## 2.1. Boundedness and Permanence

Theorem 1 (uniform boundedness).
Suppose that a non-negative functionφx,y and its partial derivatives φx and φy are continuous in ℝ+2, then the system(3a)x˙=r1x1−xK1−φx,y−m1x−dx2,(3b)y˙=r2y1−yK2yy+e+e1φx,y−m2y,subject to x0,y0≥0 is uniformly bounded.This Theorem1 holds obviously after introducing an auxiliary function z=e1x+y [1, 16]. The following theorem with the help of comparison principle in ODEs is about permanence of system (2a) and system (2b).Theorem 2 (permanence).
If parameters satisfy(4)ω1=r1−m1−q1M2/a+cM2r1/K1+d>0,a0=q2ω11−λa+bω11−λ−m2>0,a1=r2e−cq2ω11−λa+bω11−λ2<0,cρ≤1,e≥ρa+bω11−λ,where M2 is a positive constant, λ∈0,1, and ρ=r2eK2+e/K2c2q2ω11−λ3, then system (2a) and system (2b) are permanent.Proof.
From above Theorem1, we have a positive upper bound ξ1 such that(5)maxlimsupt⟶∞x,limsupt⟶∞y≤ξ1.
Then, we obtainM2>0 and sufficiently large T1, such that(6)x˙≥xr1−m1−q1M2a+cM2−r1K1+dx,∀t≥T1.
By using the lemmas in [17], we admit liminft⟶∞x≥ω1. That is to say, there exist sufficiently large T2, such that x≥ω11−λ, ∀t≥T2. Then, we consider equation (2b) again. It yields y˙≥fyy; here the function fy incorporating its Taylor expression is(7)fy=r2yy+e1−yK2−m2+q2ω11−λa+bω11−λ+cy=a0+a1y+12f″ξyy2,ξy∈0,y,in which the second order derivative of fy is(8)f″y=−2r2eK2+eK2y+e3+2c2q2ω11−λa+bω11−λ+cy3.
From the conditions in this theorem, we know thatf″y≥0, when y≥0. Thus, fy≥a0+a1y, ∀y≥0. By using the lemmas in [17] again, the proof is completed.


## 2.2. Existence of Trivial Equilibria

In this section, we will discuss the existence conditions of trivial equilibria of system (2a) and system (2b). Firstly, from [1], it is obvious that system (2a) and system (2b) have trivial equilibria: E0≔0,0, E1≔x1,0, and E2k≔0,yk, k=1,2, where x1≔r1−m1/r1/K1+d and(9)y1,2≔r2−m2±Δ2r2/K2,Δ≔r2−m22−4r2m2eK2,The pointE1 exists when r1>m1, while the points E2k all exist if r2>m2 and Δ>0. If r2>m2 but Δ=0, then the two equilibria E21 and E22 collide with each other and we denote this equilibrium as E2=0,y2.


## 3. Existence of Interior Equilibrium

Here and below, we denote the interior equilibriumE∗ of system (2a) and system (2b) as x∗,y∗ or s1,s2. This equilibrium must satisfy the following coupled algebraic equations:(10a)r11−xK1−q1ya+bx+cy−m1−dx=0,(10b)r21−yK2yy+e+q2ya+bx+cy−m2=0.From equation (10a), we have (1) y⟶r1−m1a/q1−r1−m1c or ∞ (if q1=r1−m1c) when x⟶0; (2) x⟶x1 when y⟶0. Bearing equation (10b) in mind, we have (1) y⟶yk (if yk exists) when x⟶0; (2) x⟶am2/q2−m2b or ∞ (if q2=m2b). The implicit derivative y′x from equation (10b) is(11)y′x=q2K2a+cyy+e2r2y2+2ey−K2ea+bx+cy2+q2cK2xy+e2,and we denote the positive root of a quadratic equation y2+2ey−K2e=0 as ys1=−e+e2+K2e for later use.


### 3.1. Cobweb Model

Based on above approximate analysis and the cobweb model, some cases about the existence of the interior equilibriumE∗ will be illustrated when isoclines from equations (10a) and (10b) all fall in R+2.Case 1. If parameters satisfy(12)r1−m1aq1−cr1−m1>yk>ys1,x1>0,then an interior equilibrium exists. Here, we taker1=1, r2=1, K1=20, K2=8, q1=0.12, q2=0.06, a=4, m1=0.3, m2=0.5, e=0.3, b=1, c=0.1, and d=5; the first equilibrium is E3≈0.119085,3.688431 and the second equilibrium is E∗≈0.136769,0.323748.Case 2. If parameters satisfy(13)q1=cr1−m1,yk>ys1,x1>0,then an interior equilibrium exists. Here, we taker1=1, r2=1, K1=20, K2=8, q2=0.6, a=4, m1=0.3, m2=0.5, e=0.3, b=1, c=0.1, and d=5; the first interior equilibrium is E4≈0.126835,3.83244, and the second equilibrium is E7≈0.137621,0.298566.Case 3. If parameters satisfy(14)0<x1<am2q2−m2b,0<yk<minys1,r1−m1aq1−cr1−m1,then an interior equilibrium exists. Here, we taker1=1, r2=1.5, K1=20, K2=50, q1=5, q2=1.1, a=4, b=1, c=0.1, m1=0.1, m2=0.5, e=1, and d=3.565; an interior equilibrium is E5≈0.0916118,0.470837.Case 4. If parameters satisfy(15)0<x1,q2=m2b,0<yk<minys1,r1−m1aq1−cr1−m1,then an interior equilibrium exists. Here, we taker1=1, r2=1, K1=20, K2=8, q1=0.12, a=4, m1=0.3, m2=0.5, e=0.3, b=1, c=0.1, and d=5; an interior equilibrium is E6≈0.136885,0.303164. Notice that another interior equilibrium E∗≈0.118552,3.79777 verifies Case 1.Case 5. If parameters satisfy(16)0<x1<am2q2−m2b,0<yk<ys1,q1=cr1−m1,then an interior equilibrium exists (see the second equilibriumE7 in Case 2).Case 6. If yk does not exist and parameters satisfy(17)0<am2q2−m2b<x1,q1>cr1−m1,then an interior equilibrium exists. Here, we taker1=1, r2=1.5, K1=10, K2=2, q1=5, q2=1.1, a=4, b=1, c=0.1, m1=0.1, m2=0.5, e=1, and d=0.05389805885≈d8H; an interior equilibrium is E8≈0.399541,0.750397.Case 7. If yk does not exist and parameters satisfy(18)0<am2q2−m2b<x1,q1=cr1−m1,then an interior equilibrium exists. Here, we taker1=1, r2=1.5, K1=10, K2=2, q2=1.1, a=4, b=1, c=0.1, m1=0.1, m2=0.5, e=1, and d=0.11; an interior equilibrium is E9≈4.178700,2.09452.


### 3.2. Polynomial Equations ofx∗ and y∗

From equation (10a), an expression of y is(19)y=x1−xr1/K1+da+bxq1−cr1−m1+cr1/K1+dx.Substituting it into the equation (10b), a quintic algebraic equation of x∗ can be derived in the form of px≔∑k=05akxk=0 (see coefficients ak in Appendix A.1). Thanks to Niels Henrik Abel and Evariste Gallois’s ingenious works, quintic equations usually have no analytical form solutions. But we can make some special efforts to numerically derive positive roots for this equation px=0. For instance, if(20)px1=K2eq13K13K1d+r1am2K1d+r1+K1r1−m1bm2−q2>0,where x1>0 and a0<0, then there is a positive root x∗ such that x∗<x1.Case 1. If parameters satisfy(21)px1>0,x1>0,a0<0,q1≥cr1−m1,then an interior equilibrium exists. Here, we taker1=1, r2=1, K1=10, K2=8, m1=0.3, m2=0.5, d=0.1, q1=0.12, q2=0.16, e=0.3, a=0.2, b=0.3, and c=0.1; an interior equilibrium is E10≈3.491956,0.016748, and the following lemma is verified.Lemma 1.
Suppose thatfx=∑k=0nakxk is a polynomial with real coefficients, an≠0, n>1. If ana0<0, then the equation fx=0 has a positive root.Proof.
We only consider the special casean>0, and thus a0<0. It is clear that the polynomial fx has a decomposition(22)fx=∑k=0n−11nanxn−k+akxk.If we take a sufficiently large positive numberX such that X>max0≤k≤n−1nak/an1/n−k, then fX>0, and thus we complete the proof.Lemma 2.
Suppose thatfx=∑k=0nakxk is a real polynomial, an≠0, n>1. If there is a positive number x0>0 such that anfx0<0, then the equation fx=0 has a positive root.On the other hand, from equation (10b), we obtain an expression(23)x=a+cym2y+e−r2y1−y/K2q2−bm2y+e+br2y1−y/K2.And a quintic algebraic equation ofy∗ can be written in the form of qy≔∑k=05bkyk=0 (see coefficients bk in Appendix A.2). Suppose that the denominator and numerator of expression (23) are all positive for some y∗>0; then, x∗>0. Notice that if b5>0 and b0<0, there must exist a positive root for the quintic equation qy=0 (see Lemma 1).Case 2. If parameters satisfy b0<0, m2≥r2, q2>bm2, yu1>0, and qyu1>0, where(24)yu1=q2−bm2+br2+q2−bm2+br22+4br2e/K2q2−bm22br2/K2,is a root of the quadric equation in the denominator of (23), then an interior equilibrium exists. Here, we take r1=1, r2=0.5, K1=10, K2=8, m1=0.3, m2=0.6, d=0.1, q1=0.8, q2=0.6, e=0.3, a=0.2, b=0.2, and c=0.1, and an interior equilibrium is E11≈0.169319,0.212419.Case 3. If parameters satisfy b0<0, y2>0, qy2>0, yu1>0, and qyu1>0, then an interior equilibrium exists. Here, we take r1=1, r2=1, K1=20, K2=80, q1=3, q2=10, a=2, b=1, c=1, m1=0.1, m2=0.5, e=2, and d=2.092753894≈d12H, and an interior equilibrium is E12≈0.067460,0.695847.Remark 1.
Here we rewrite equations (10a) and (10b) as(25a)fx,y=r11−xK1−m1−dxa+bx+cy−q1y=0,(25b)gx,y=r2y1−yK2a+bx+cy+q2xy+e−m2a+bx+cyy+e=0.
Furthermore, if we sort them asfx,y=a0xy+a1x, gx,y=b0xy3+b1xy2+b2xy+b3x, the first eliminant(26)Ryf,g=a0xa1x000a0xa1x000a0xa1xb0xb1xb2xb3x,also yields above quintic equation px=0. Similarly, if we sort them as fx,y=a0yx2+a1yx+a2y, gx,y=b0yx+b1y, the second eliminant(27)Rxf,g=a0ya1ya2yb0yb1y00b0yb1y,yields the quintic equation qy=0.


## 3.1. Cobweb Model

Based on above approximate analysis and the cobweb model, some cases about the existence of the interior equilibriumE∗ will be illustrated when isoclines from equations (10a) and (10b) all fall in R+2.Case 1. If parameters satisfy(12)r1−m1aq1−cr1−m1>yk>ys1,x1>0,then an interior equilibrium exists. Here, we taker1=1, r2=1, K1=20, K2=8, q1=0.12, q2=0.06, a=4, m1=0.3, m2=0.5, e=0.3, b=1, c=0.1, and d=5; the first equilibrium is E3≈0.119085,3.688431 and the second equilibrium is E∗≈0.136769,0.323748.Case 2. If parameters satisfy(13)q1=cr1−m1,yk>ys1,x1>0,then an interior equilibrium exists. Here, we taker1=1, r2=1, K1=20, K2=8, q2=0.6, a=4, m1=0.3, m2=0.5, e=0.3, b=1, c=0.1, and d=5; the first interior equilibrium is E4≈0.126835,3.83244, and the second equilibrium is E7≈0.137621,0.298566.Case 3. If parameters satisfy(14)0<x1<am2q2−m2b,0<yk<minys1,r1−m1aq1−cr1−m1,then an interior equilibrium exists. Here, we taker1=1, r2=1.5, K1=20, K2=50, q1=5, q2=1.1, a=4, b=1, c=0.1, m1=0.1, m2=0.5, e=1, and d=3.565; an interior equilibrium is E5≈0.0916118,0.470837.Case 4. If parameters satisfy(15)0<x1,q2=m2b,0<yk<minys1,r1−m1aq1−cr1−m1,then an interior equilibrium exists. Here, we taker1=1, r2=1, K1=20, K2=8, q1=0.12, a=4, m1=0.3, m2=0.5, e=0.3, b=1, c=0.1, and d=5; an interior equilibrium is E6≈0.136885,0.303164. Notice that another interior equilibrium E∗≈0.118552,3.79777 verifies Case 1.Case 5. If parameters satisfy(16)0<x1<am2q2−m2b,0<yk<ys1,q1=cr1−m1,then an interior equilibrium exists (see the second equilibriumE7 in Case 2).Case 6. If yk does not exist and parameters satisfy(17)0<am2q2−m2b<x1,q1>cr1−m1,then an interior equilibrium exists. Here, we taker1=1, r2=1.5, K1=10, K2=2, q1=5, q2=1.1, a=4, b=1, c=0.1, m1=0.1, m2=0.5, e=1, and d=0.05389805885≈d8H; an interior equilibrium is E8≈0.399541,0.750397.Case 7. If yk does not exist and parameters satisfy(18)0<am2q2−m2b<x1,q1=cr1−m1,then an interior equilibrium exists. Here, we taker1=1, r2=1.5, K1=10, K2=2, q2=1.1, a=4, b=1, c=0.1, m1=0.1, m2=0.5, e=1, and d=0.11; an interior equilibrium is E9≈4.178700,2.09452.


## 3.2. Polynomial Equations ofx∗ and y∗

From equation (10a), an expression of y is(19)y=x1−xr1/K1+da+bxq1−cr1−m1+cr1/K1+dx.Substituting it into the equation (10b), a quintic algebraic equation of x∗ can be derived in the form of px≔∑k=05akxk=0 (see coefficients ak in Appendix A.1). Thanks to Niels Henrik Abel and Evariste Gallois’s ingenious works, quintic equations usually have no analytical form solutions. But we can make some special efforts to numerically derive positive roots for this equation px=0. For instance, if(20)px1=K2eq13K13K1d+r1am2K1d+r1+K1r1−m1bm2−q2>0,where x1>0 and a0<0, then there is a positive root x∗ such that x∗<x1.Case 1. If parameters satisfy(21)px1>0,x1>0,a0<0,q1≥cr1−m1,then an interior equilibrium exists. Here, we taker1=1, r2=1, K1=10, K2=8, m1=0.3, m2=0.5, d=0.1, q1=0.12, q2=0.16, e=0.3, a=0.2, b=0.3, and c=0.1; an interior equilibrium is E10≈3.491956,0.016748, and the following lemma is verified.Lemma 1.
Suppose thatfx=∑k=0nakxk is a polynomial with real coefficients, an≠0, n>1. If ana0<0, then the equation fx=0 has a positive root.Proof.
We only consider the special casean>0, and thus a0<0. It is clear that the polynomial fx has a decomposition(22)fx=∑k=0n−11nanxn−k+akxk.If we take a sufficiently large positive numberX such that X>max0≤k≤n−1nak/an1/n−k, then fX>0, and thus we complete the proof.Lemma 2.
Suppose thatfx=∑k=0nakxk is a real polynomial, an≠0, n>1. If there is a positive number x0>0 such that anfx0<0, then the equation fx=0 has a positive root.On the other hand, from equation (10b), we obtain an expression(23)x=a+cym2y+e−r2y1−y/K2q2−bm2y+e+br2y1−y/K2.And a quintic algebraic equation ofy∗ can be written in the form of qy≔∑k=05bkyk=0 (see coefficients bk in Appendix A.2). Suppose that the denominator and numerator of expression (23) are all positive for some y∗>0; then, x∗>0. Notice that if b5>0 and b0<0, there must exist a positive root for the quintic equation qy=0 (see Lemma 1).Case 2. If parameters satisfy b0<0, m2≥r2, q2>bm2, yu1>0, and qyu1>0, where(24)yu1=q2−bm2+br2+q2−bm2+br22+4br2e/K2q2−bm22br2/K2,is a root of the quadric equation in the denominator of (23), then an interior equilibrium exists. Here, we take r1=1, r2=0.5, K1=10, K2=8, m1=0.3, m2=0.6, d=0.1, q1=0.8, q2=0.6, e=0.3, a=0.2, b=0.2, and c=0.1, and an interior equilibrium is E11≈0.169319,0.212419.Case 3. If parameters satisfy b0<0, y2>0, qy2>0, yu1>0, and qyu1>0, then an interior equilibrium exists. Here, we take r1=1, r2=1, K1=20, K2=80, q1=3, q2=10, a=2, b=1, c=1, m1=0.1, m2=0.5, e=2, and d=2.092753894≈d12H, and an interior equilibrium is E12≈0.067460,0.695847.Remark 1.
Here we rewrite equations (10a) and (10b) as(25a)fx,y=r11−xK1−m1−dxa+bx+cy−q1y=0,(25b)gx,y=r2y1−yK2a+bx+cy+q2xy+e−m2a+bx+cyy+e=0.
Furthermore, if we sort them asfx,y=a0xy+a1x, gx,y=b0xy3+b1xy2+b2xy+b3x, the first eliminant(26)Ryf,g=a0xa1x000a0xa1x000a0xa1xb0xb1xb2xb3x,also yields above quintic equation px=0. Similarly, if we sort them as fx,y=a0yx2+a1yx+a2y, gx,y=b0yx+b1y, the second eliminant(27)Rxf,g=a0ya1ya2yb0yb1y00b0yb1y,yields the quintic equation qy=0.


## 4. Stability Analysis of System (2a) and System (2b)

In this section, we use the Routh–Hurwitz criterion and Perron’s theorems to analyze local stability of above equilibria in their existence interval, respectively. A theorem about global asymptotic stability and a theorem about nonexistence of limit cycles are also considered.


### 4.1. Local Stability Analysis

The Jacobian matrix of system (2a) and system (2b) takes the following form J=Jij2×2, where four components are(28)J11=r1−m1−2xr1K1+d−q1a+cyya+bx+cy2,J12=−q1xa+bxa+bx+cy2,J21=q2ya+cya+bx+cy2,J22=r2yK2y+e2−2y2−3ey+K2y+2K2e+q2xa+bxa+bx+cy2−m2.For the trivial equilibriumE0 with r1≠m1 and the axial equilibrium E1 with q2x1/bx1+a≠m2, we omit their stability [1]. When r1=m1, the transformation τ=−m2t and equations (2a) and (2b) yield x˙=r1/K1+d/m2x2 (we still use symbol t), and thus E0 is a saddle node and the parabolic sector is on the right half plane.In the case that two axial equilibriaE2kk=1,2 exist, the Jacobian matrices are(29)JE2k=r1−m1−q1yka+cyk0q2yka+cykm2−r2yk+2m2eyk+e.SinceJ22E21<0, we have the following: (a) E21 is a saddle point if r1−m1>q1y1/a+cy1; (b) E21 is an asymptotically stable node if r1−m1<q1y1/a+cy1. For the equilibrium E22, it is obvious that J22E22>0, and we have the following: (a) E22 is an unstable node if r1−m1>q1y2/a+cy2; (b) E22 is a saddle point if r1−m1<q1y2/a+cy2; (c) E22 is an unstable higher-order singular point if r1−m1=q1y2/a+cy2. In the special case that two axial equilibria E2kk=1,2 collide with each other, the Jacobian matrix is(30)JE2=r1−m1−q1y2a+cy20q2y2a+cy20,and thus E2 is a higher-order singular point. If r1−m1>q1y2/a+cy2, E2 is unstable.For the interior equilibriumE∗, its Jacobian matrix is(31)JE∗=J11E∗−q1x∗a+bx∗a+bx∗+cy∗2q2y∗a+cy∗a+bx∗+cy∗2J22E∗,where(32)J11E∗=−r1K1+dx∗+q1bx∗y∗a+bx∗+cy∗2,J22E∗=r2y∗K2y∗+e2−y∗2−2ey∗+K2e−q2cx∗y∗a+bx∗+cy∗2.Denoting a new discriminantΔ∗≔A12−4A2 with the trace A1≔trJE∗ and the determinant A2≔detJE∗, we have(a)
IfA1<0 and (a1) A2>0,Δ∗≥0, then E∗ is an asymptotically stable node; (a2) A2>0,Δ∗<0, then E∗ is an asymptotically stable focus; (a3) A2<0, then E∗ is a saddle point.(b)
IfA1=0 and (b1) A2>0, then E∗ is a center or a focus; (b2) A2<0, then E∗ is a saddle point.(c)
IfA1>0, then E∗ is unstable and (c1) Δ∗=0, then E∗ is a node; (c2) Δ∗<0, then E∗ is a focus; (c3) Δ∗>0 and A2>0, then E∗ is a node; (c4) Δ∗>0 and A2<0, then E∗ is a saddle point.WhenA2=0 but A1≠0, E∗ is a stable (unstable) node if A1<0 (A1>0) (see Theorem 7.1 in Zhifen Zhang’s book [18] for more details). It is probable that E∗ has a cusp case of codimension at least 2 which ensures potential Bogdanov–Takens bifurcation when A1=A2=0.


### 4.2. Global Asymptotic Stability

Combining the stability conclusions of the pointE0 in above section, the positive definite Lyapunov function V=e1x+y ensures that E0 is globally asymptotically stable if one of the following conditions holds:(33)ar1≤m1,r2≤m2;br1≤m1,Δ<0.Furthermore, conditionsr1<m1, r2+q2/b−m2≤0 and Theorem 1 deduce global asymptotical stability of E0 pronto. If r1>m1, r2+q2/b−m2≤0, and q2x1/a+bx1−m1≤0, equilibria E2k, E2, and E∗ do not exist, and E0 is unstable, then Theorem 1 ensures that E1 is globally asymptotically stable. For the further consideration of point E∗, the following theorem explains its global asymptotic stability.Theorem 3 (global asymptotic stability ofE∗).
If a unique interior equilibriumE∗ exists and parameters satisfy(34)r1K1+d>q1by∗aa+bx∗+cy∗,y∗>K2,then E∗ is globally asymptotically stable.Proof.
Here we take an unbounded positive definite Lyapunov function(35)V=Vx,y=x−x∗−x∗lnxx∗+Ay−y∗−y∗lnyy∗,with A=q1a+bx∗/q2a+cy∗. Introducing new variables x^=x−x∗ and y^=y−y∗, computing derivative along orbits of system (2a) and system (2b), we have(36)dVdt|2=−r1K1−d+q1by∗a+bx+cya+bx∗+cy∗x^2+Ar2y+ey∗+ee−yy∗+ey+y∗K2−q2cx∗a+bx+cya+bx∗+cy∗y^2.
It is obvious thatdV/dt|2 is negative definite. Consequently, the Lyapunov function V satisfies the asymptotic stability theorem in [19]. Thus, we complete the proof.
If conditionsq1by∗/aa+bx∗+cy∗≤r1/K1+d and y∗≥K2 hold, then A1<0, A2>0 and E∗ is an asymptotically stable node or focus. Additionally, we could solve a potential interior equilibrium E∗ with two control variables λ∈0,1 and μ≥1 from such conditions, where(37)x∗=−aK1d+r1K2cμ+aλ+q1bμK2K1abλK1d+r1,y∗=μK2,and parameters m1 and m2 are constrained by(38)m1=−a2K1d+r1λ2+aK1d+r1a+K2cdμ+br1K1+cK2μr1λ−q1bμK2K1abλK1,m2=1μK2K1q1b2K2μ+e·∑k=03akμk,where coefficients are(39)a3=−K1K22b2q1r2,a2=−acdλ−bq1q2−r2b2q1K1+aλcq2r1K22,a1=−aK1d+r1ce+aλ−beK1q1K2q2,a0=−a2q2λeK1d+r1.
Here we set some values of parameters asr1=1.21, r2=1.5, K1=2, K2=5, a=1, b=1.5, c=1.7, d=0.3, e=0.3, q1=1.2, and q2=0.21, and the unique globally asymptotically stable node E∗≈0.655936,5.04 with λ=0.95 and μ=1.008 is depicted in Figure 1 with characteristic direction θ=θ2≈0.101084 which is solved from the characteristic function(40)Θθ≈0.076926cosθ2−0.898114sinθcosθ+0.014025=C21θ−θ2+Oθ−θ22,after we make the polar-coordinate-transformation x=rcosθ+x∗, y=rsinθ+y∗. The residual real simple root of equation Θθ=0 in interval 0,2π is θ1≈1.555157. The power series of above characteristic function Θθ up to order two at the point θ2 admits indexes C21<0 and Rθ2<0 with same negative sign, i.e., p=1(odd number) and k=2, where function(41)Rθ≈0.898114cosθ2+0.076926sinθcosθ−1.438292.
Thus,θ=θ2 actually shows a fact that trajectories enter the stable node along this direction.Figure 1
Phase diagrams around an unique stable nodeE∗ with characteristic direction θ=θ2 (in blue).


### 4.3. Nonexistence of Limit Cycles

In this section, we consider nonexistence of closed orbits and limit cycles of system (2a) and system (2b). Firstly, taking a diffeomorphism φ:u=x,v=y,τ=t/a+bx+cyy+e which preserves the orientation of time, system (2a) and system (2b) are topologically equivalent to following system:(42a)x˙=Px,y≔r1x1−xK1−m1x−dx2a+bx+cyy+e−q1xyy+e,(42b)y˙=Qx,y≔r2y21−yK2a+bx+cy+q2xyy+e−m2ya+bx+cyy+e,since detDφu,v,τ=v+ea+bu+cv>0, where matrix(43)Dφu,v,τ=100010∗∗v+ea+bu+cv.Notice that we still denoteu, v, and τ as x, y, and t. Above system is a C∞-qualitatively equivalent polynomial form of system (2a) and system (2b), and it is more convenient to consider limit cycles [4, 5, 20, 21].Theorem 4 (nonexistence of limit cycles).
For system (42a) and system (42b), if parameters satisfy(44)m1+m2−r1−r2b+cde+ad−q2≥0,ad+bm1−br1≥0,ce+am2−ar2>0,cm2−r2K2+ar2≥0,then there are no closed orbits and limit cycles in ℝ+2.Proof.
Here we take a Dulac functionBx,y=1/xy and calculate following partial derivative:(45)∂BP∂x+∂BQ∂y=1K1K2xy∑1≤i+j≤3aijxiyj,where coefficients(46)a02=−2K1cm2−r2K2+ar2,a11=−K2−r1−r2+m1+m2b+cde+ad−q2K1+r1ce+a,a10=−ead+bm1−br1K1+ar1K2,a01=−K2ce+am2−ar2K1,and other unlisted coefficients are all nonpositive. Thus, we complete the proof.
Here we take over values of parameters from Section4.2 but m1=m2=2, and system (2a) and system (2b) do not have meaningful equilibria except a globally asymptotically stable node E0 (the origin O) with the characteristic direction θ=0 (the positive x-axis) since the characteristic equation Θθ=m1−m2−r1sinθcosθ, and conditions in Theorem 4 all hold. Hence, there are no closed orbits and limit cycles in this numerical case.
In addition, system (2a) and system (2b) merely own a saddle point E0 and a globally asymptotically stable node E1 when we use parameters from Section 4.2 but set a=2, d=3, m1=0.4, and m2=1.4. Conditions in Theorem 4 are also satisfied.


## 4.1. Local Stability Analysis

The Jacobian matrix of system (2a) and system (2b) takes the following form J=Jij2×2, where four components are(28)J11=r1−m1−2xr1K1+d−q1a+cyya+bx+cy2,J12=−q1xa+bxa+bx+cy2,J21=q2ya+cya+bx+cy2,J22=r2yK2y+e2−2y2−3ey+K2y+2K2e+q2xa+bxa+bx+cy2−m2.For the trivial equilibriumE0 with r1≠m1 and the axial equilibrium E1 with q2x1/bx1+a≠m2, we omit their stability [1]. When r1=m1, the transformation τ=−m2t and equations (2a) and (2b) yield x˙=r1/K1+d/m2x2 (we still use symbol t), and thus E0 is a saddle node and the parabolic sector is on the right half plane.In the case that two axial equilibriaE2kk=1,2 exist, the Jacobian matrices are(29)JE2k=r1−m1−q1yka+cyk0q2yka+cykm2−r2yk+2m2eyk+e.SinceJ22E21<0, we have the following: (a) E21 is a saddle point if r1−m1>q1y1/a+cy1; (b) E21 is an asymptotically stable node if r1−m1<q1y1/a+cy1. For the equilibrium E22, it is obvious that J22E22>0, and we have the following: (a) E22 is an unstable node if r1−m1>q1y2/a+cy2; (b) E22 is a saddle point if r1−m1<q1y2/a+cy2; (c) E22 is an unstable higher-order singular point if r1−m1=q1y2/a+cy2. In the special case that two axial equilibria E2kk=1,2 collide with each other, the Jacobian matrix is(30)JE2=r1−m1−q1y2a+cy20q2y2a+cy20,and thus E2 is a higher-order singular point. If r1−m1>q1y2/a+cy2, E2 is unstable.For the interior equilibriumE∗, its Jacobian matrix is(31)JE∗=J11E∗−q1x∗a+bx∗a+bx∗+cy∗2q2y∗a+cy∗a+bx∗+cy∗2J22E∗,where(32)J11E∗=−r1K1+dx∗+q1bx∗y∗a+bx∗+cy∗2,J22E∗=r2y∗K2y∗+e2−y∗2−2ey∗+K2e−q2cx∗y∗a+bx∗+cy∗2.Denoting a new discriminantΔ∗≔A12−4A2 with the trace A1≔trJE∗ and the determinant A2≔detJE∗, we have(a)
IfA1<0 and (a1) A2>0,Δ∗≥0, then E∗ is an asymptotically stable node; (a2) A2>0,Δ∗<0, then E∗ is an asymptotically stable focus; (a3) A2<0, then E∗ is a saddle point.(b)
IfA1=0 and (b1) A2>0, then E∗ is a center or a focus; (b2) A2<0, then E∗ is a saddle point.(c)
IfA1>0, then E∗ is unstable and (c1) Δ∗=0, then E∗ is a node; (c2) Δ∗<0, then E∗ is a focus; (c3) Δ∗>0 and A2>0, then E∗ is a node; (c4) Δ∗>0 and A2<0, then E∗ is a saddle point.WhenA2=0 but A1≠0, E∗ is a stable (unstable) node if A1<0 (A1>0) (see Theorem 7.1 in Zhifen Zhang’s book [18] for more details). It is probable that E∗ has a cusp case of codimension at least 2 which ensures potential Bogdanov–Takens bifurcation when A1=A2=0.


## 4.2. Global Asymptotic Stability

Combining the stability conclusions of the pointE0 in above section, the positive definite Lyapunov function V=e1x+y ensures that E0 is globally asymptotically stable if one of the following conditions holds:(33)ar1≤m1,r2≤m2;br1≤m1,Δ<0.Furthermore, conditionsr1<m1, r2+q2/b−m2≤0 and Theorem 1 deduce global asymptotical stability of E0 pronto. If r1>m1, r2+q2/b−m2≤0, and q2x1/a+bx1−m1≤0, equilibria E2k, E2, and E∗ do not exist, and E0 is unstable, then Theorem 1 ensures that E1 is globally asymptotically stable. For the further consideration of point E∗, the following theorem explains its global asymptotic stability.Theorem 3 (global asymptotic stability ofE∗).
If a unique interior equilibriumE∗ exists and parameters satisfy(34)r1K1+d>q1by∗aa+bx∗+cy∗,y∗>K2,then E∗ is globally asymptotically stable.Proof.
Here we take an unbounded positive definite Lyapunov function(35)V=Vx,y=x−x∗−x∗lnxx∗+Ay−y∗−y∗lnyy∗,with A=q1a+bx∗/q2a+cy∗. Introducing new variables x^=x−x∗ and y^=y−y∗, computing derivative along orbits of system (2a) and system (2b), we have(36)dVdt|2=−r1K1−d+q1by∗a+bx+cya+bx∗+cy∗x^2+Ar2y+ey∗+ee−yy∗+ey+y∗K2−q2cx∗a+bx+cya+bx∗+cy∗y^2.
It is obvious thatdV/dt|2 is negative definite. Consequently, the Lyapunov function V satisfies the asymptotic stability theorem in [19]. Thus, we complete the proof.
If conditionsq1by∗/aa+bx∗+cy∗≤r1/K1+d and y∗≥K2 hold, then A1<0, A2>0 and E∗ is an asymptotically stable node or focus. Additionally, we could solve a potential interior equilibrium E∗ with two control variables λ∈0,1 and μ≥1 from such conditions, where(37)x∗=−aK1d+r1K2cμ+aλ+q1bμK2K1abλK1d+r1,y∗=μK2,and parameters m1 and m2 are constrained by(38)m1=−a2K1d+r1λ2+aK1d+r1a+K2cdμ+br1K1+cK2μr1λ−q1bμK2K1abλK1,m2=1μK2K1q1b2K2μ+e·∑k=03akμk,where coefficients are(39)a3=−K1K22b2q1r2,a2=−acdλ−bq1q2−r2b2q1K1+aλcq2r1K22,a1=−aK1d+r1ce+aλ−beK1q1K2q2,a0=−a2q2λeK1d+r1.
Here we set some values of parameters asr1=1.21, r2=1.5, K1=2, K2=5, a=1, b=1.5, c=1.7, d=0.3, e=0.3, q1=1.2, and q2=0.21, and the unique globally asymptotically stable node E∗≈0.655936,5.04 with λ=0.95 and μ=1.008 is depicted in Figure 1 with characteristic direction θ=θ2≈0.101084 which is solved from the characteristic function(40)Θθ≈0.076926cosθ2−0.898114sinθcosθ+0.014025=C21θ−θ2+Oθ−θ22,after we make the polar-coordinate-transformation x=rcosθ+x∗, y=rsinθ+y∗. The residual real simple root of equation Θθ=0 in interval 0,2π is θ1≈1.555157. The power series of above characteristic function Θθ up to order two at the point θ2 admits indexes C21<0 and Rθ2<0 with same negative sign, i.e., p=1(odd number) and k=2, where function(41)Rθ≈0.898114cosθ2+0.076926sinθcosθ−1.438292.
Thus,θ=θ2 actually shows a fact that trajectories enter the stable node along this direction.Figure 1
Phase diagrams around an unique stable nodeE∗ with characteristic direction θ=θ2 (in blue).


## 4.3. Nonexistence of Limit Cycles

In this section, we consider nonexistence of closed orbits and limit cycles of system (2a) and system (2b). Firstly, taking a diffeomorphism φ:u=x,v=y,τ=t/a+bx+cyy+e which preserves the orientation of time, system (2a) and system (2b) are topologically equivalent to following system:(42a)x˙=Px,y≔r1x1−xK1−m1x−dx2a+bx+cyy+e−q1xyy+e,(42b)y˙=Qx,y≔r2y21−yK2a+bx+cy+q2xyy+e−m2ya+bx+cyy+e,since detDφu,v,τ=v+ea+bu+cv>0, where matrix(43)Dφu,v,τ=100010∗∗v+ea+bu+cv.Notice that we still denoteu, v, and τ as x, y, and t. Above system is a C∞-qualitatively equivalent polynomial form of system (2a) and system (2b), and it is more convenient to consider limit cycles [4, 5, 20, 21].Theorem 4 (nonexistence of limit cycles).
For system (42a) and system (42b), if parameters satisfy(44)m1+m2−r1−r2b+cde+ad−q2≥0,ad+bm1−br1≥0,ce+am2−ar2>0,cm2−r2K2+ar2≥0,then there are no closed orbits and limit cycles in ℝ+2.Proof.
Here we take a Dulac functionBx,y=1/xy and calculate following partial derivative:(45)∂BP∂x+∂BQ∂y=1K1K2xy∑1≤i+j≤3aijxiyj,where coefficients(46)a02=−2K1cm2−r2K2+ar2,a11=−K2−r1−r2+m1+m2b+cde+ad−q2K1+r1ce+a,a10=−ead+bm1−br1K1+ar1K2,a01=−K2ce+am2−ar2K1,and other unlisted coefficients are all nonpositive. Thus, we complete the proof.
Here we take over values of parameters from Section4.2 but m1=m2=2, and system (2a) and system (2b) do not have meaningful equilibria except a globally asymptotically stable node E0 (the origin O) with the characteristic direction θ=0 (the positive x-axis) since the characteristic equation Θθ=m1−m2−r1sinθcosθ, and conditions in Theorem 4 all hold. Hence, there are no closed orbits and limit cycles in this numerical case.
In addition, system (2a) and system (2b) merely own a saddle point E0 and a globally asymptotically stable node E1 when we use parameters from Section 4.2 but set a=2, d=3, m1=0.4, and m2=1.4. Conditions in Theorem 4 are also satisfied.


## 5. Local Bifurcations of System (2a) and System (2b)

In this section, we will give sufficient conditions to show the existence of saddle-node bifurcation, transcritical bifurcation, and Hopf bifurcation of system (2a) and system (2b). Firstly, we denote this system in the following form:(47)x˙y˙=Fx,y,for simplicity and convenience.


### 5.1. Saddle-Node Bifurcation

The two trivial equilibriaE2kk=1,2 collide with each other and system (47) has a unique boundary equilibrium E2 when y>0, if r2>m2 and y1=y2, or(48)m2=m2SN≔r2+2r2eK21−1+K2e.Then, there is a chance of bifurcation around this higher-order singular point. Here we choosem2 as a bifurcation parameter and select eigenvector v=01 corresponding to the zero eigenvalue for matrix (30). The eigenvector corresponding to the zero eigenvalue for the transpose of matrix (30) is w=1w2, where(49)w2=−r1−m1−q1y2/a+cy2q2y2/a+cy2.Supposer1−m1−q1y2/a+cy2<0, then the following transversality conditions hold:(50)wTFm2E2,m2SN=−w2y2≠0,wTD2FE2,m2SNv,v=2r2w2K2eK2+e−1≠0.Thus, we have following theorem by using Sotomayor’s theorem [22, 23].Theorem 5 (saddle-node bifurcation).
Suppose that the pointE2 exists; if r1−m1−q1y2/a+cy2<0, then system (47) undergoes a saddle-node bifurcation around point E2 with respect to the bifurcation parameter m2.


### 5.2. Transcritical and Pitchfork Bifurcation

The equilibriumE21 changes its stability when r1 crosses the threshold r1TC≔m1+q1y1/a+cy1; in other words, E21 is a higher-order equilibrium when r1=r1TC. Let J22E21<0. Thus, we choose the parameter r1 as a bifurcation parameter and an eigenvector v=1v2 corresponding to the zero eigenvalue for the Jacobian matrix JE21 when r1=r1TC, where v2=−q2y1/J22E21a+cy1>0. The eigenvector corresponding to the zero eigenvalue for the transpose of matrix JE21 is w=10; then, the transversality conditions are(51)wTFr1E21,r1TC=0,wTDFr1E21,r1TCv=1,wTD2FE21,r1TCv,v=−2r1TCK1+2q1by1a+cy12−2d+2q1q2ay1a+cy13J22E21.Thus, we have the following theorem by using Sotomayor’s theorem [22, 23].Theorem 6 (transcritical bifurcation).
Suppose that the two axial equilibria (E2kk=1,2) coexist and J22E21<0; if a+cy1≥bK1 or J22E21≥−q2a/ba+cy1, then system (47) undergoes a transcritical bifurcation around the point E21 with respect to the bifurcation parameter r1.For the special case(52)wTD2FE21,r1TCv,v=0,there is a chance the system (47) undergoes a pitchfork bifurcation. We still use the bifurcation parameter r1 and eigenvectors v and w; then, the fourth transversality condition is(53)wTD3FE21,r1PFv,v,v=−6q1bcJ22y1+aJ22b+q2cJ22b−q2y1+abJ22y1cy1+a5J222≠0,where J22=J22E21, r1PF≔m1+q1y1/a+cy1. Thus, we have another theorem by using Sotomayor’s theorem [22, 23].Theorem 7 (pitchfork bifurcation).
Suppose that the two axial equilibria (E2kk=1,2) coexist and J22E21<0; if condition (52) and J22E21≠−aq2/bcy1+a hold, then system (47) undergoes a pitchfork bifurcation around the point E21 with respect to the bifurcation parameter r1.


### 5.3. Hopf Bifurcation

In this section, we consider the Hopf bifurcation of system (2a) and system (2b). Here the point E∗ exists and we choose d as bifurcation parameter. Suppose that λd=αd±iωd are a pair of conjugate eigenvalues of matrix JE∗, where αd=1/2A1d. The critical value dH satisfies(54)αdH=0,α′dH≠0,A2dH>0.Then, system (47) undergoes a Hopf bifurcation around the point E∗ with respect to the bifurcation parameter d.We will calculate the first Lyapunov numberσ at the point E∗, which is used to determine the stability of limit cycles and Hopf bifurcation direction. The method and calculations of the first and second Lyapunov coefficients can be found in [22, 24, 25]. Therefore, translating the point E∗ to the origin O=0,0 by a linear transformation (I): X=x−x∗, Y=y−y∗, system (47) in power series around the origin is(55)X˙=f1X,Y=∑1≤i+j≤3aijXiYj+OX,Y4,Y˙=g1X,Y=∑1≤i+j≤3bijXiYj+OX,Y4,where coefficients are(56)a10=−r1s1K1−ds1+q1bs1s2bs1+cs2+a2,a01=−q1s1bs1+abs1+cs2+a2,a20=−r1K1+q1s2bcs2+abs1+cs2+a3−d,a02=q1s1cbs1+abs1+cs2+a3,a11=−a2+bs1+cs2a+2s1s2bcq1bs1+cs2+a3,a30=−q1s2b2cs2+abs1+cs2+a4,a03=−q1s1c2bs1+abs1+cs2+a4,a21=bs12cs2+ab−c2s22+a2q1bs1+cs2+a4,a12=cs22bs1+ac−b2s12+a2q1bs1+cs2+a4;b10=q2s2cs2+abs1+cs2+a2,b01=−r2s2−K2e+2es2+s22K2s2+e2−q2cs1s2bs1+cs2+a2,b20=−q2s2bcs2+abs1+cs2+a3,b02=r2K2−3s2e2−3s22e−s23K2s2+e3−q2s1cbs1+abs1+cs2+a3,b11=a2+bs1+cs2a+2s1s2bcq2bs1+cs2+a3,b30=q2s2b2cs2+abs1+cs2+a4,b03=−r2e2K2+eK2s2+e4+q2s1c2bs1+abs1+cs2+a4,b21=−bs12cs2+ab−c2s22+a2q2bs1+cs2+a4,b12=−cs22bs1+ac−b2s12+a2q2bs1+cs2+a4,and OX,Y4 stands for some smooth functions. After that, by using a transformation (II): u=Y, v=−Xb10+Ya10/β with β=A2E∗>0, the above system becomes a standard form:(57)u˙=f2u,v=−βv+∑2≤i+j≤3Aijuivj+Ou,v4,v˙=g2u,v=βu+∑2≤i+j≤3Bijuivj+Ou,v4,where(58)A20=a102b20+a10b10b11+b02b102b102,A11=−2a10b20+b10b11βb102,A02=−b20a01b10+a102b102,A30=a103b30+a102b10b21+a10b102b12+b03b103b103,A21=−3a102b30+2a10b10b21+b102b12βb103,A12=−3a10b30+b10b21a01b10+a102b103,A03=−b30β3b103;B20=−a02b103−a103b20+a102a20b10−a102b10b11+a10a11b102−a10b02b102βb102,B11=−2a102b20−2a10a20b10+a10b10b11−a11b102b102,B02=a10b20−a20b10βb102,B30=−a03b104−a104b30+a103a30b10−a103b10b21+a102a21b102−a102b102b12+a10a12b103−a10b03b103βb103,B21=−3a103b30−3a102a30b10+2a102b10b21−2a10a21b102+a10b102b12−a12b103b103,B12=β3a102b30−3a10a30b10+a10b10b21−a21b102b103,B03=a10b30−a30b10a01b10+a102b103.Letp=q=1/2i1 be two corresponding eigenvectors of a matrix A such that Aq=iβq, A⊤p=−iβp and <p,q>=1; the operation <x,y>=x†y (x,y∈ℂ) with the Hermitian transpose (upper symbol) † represents the usual inner product, and the Jacobian matrix is(59)A=0−ββ0.We should rewrite the functions in the right hand side of system (57) in the form of power series (60)F2x=f2xg2x=Ax+12!Bx,x+13!Cx,x,x+Ox4,where the components of linear functions B and C are(61)Bix,y=∑j,k=12∂2F2i0∂ξj∂ξkxjyk,Cix,y,z=∑j,k,l=12∂3F2i0∂ξj∂ξk∂ξlxjykzl,i=1,2,while x is the two-dimensional Euclidean norm of x. Define Tc be the largest subspace invariant by the matrix A and the generalized eigen subspace corresponding to the pair of purely imaginary eigenvalues ±iβ, i.e., Tc=spanq,q¯. That is to say, for any element z∈Tc, there must exist a linear expansion z=wq+w¯q¯. Now we can construct a two-dimensional parameterized center manifold(62)H=Hw,w¯=wq+w¯q¯+∑j+k≥2∞1j!·k!hjkwjw¯k,in which hjk∈ℂ2 and its complex conjugate is hjk¯=hkj. Combining equations (57) and Hw,w¯, we arrive at a complex equation without consideration:(63)∂H∂ww˙+∂H∂w¯w¯˙=F2H.The “chart”w for the central manifold H should be extracted from following differential equation:(64)w˙=iβw+12G21ww2+⋯.So, substituting it into equation (63) and comparing coefficients of these wjw¯k, we recursively derive(65)h11=−A−1Bq,q¯,h20=2iβI−A−1Bq,q,and an equality from coefficient of the cubic term ww2:(66)iβI−Ah21+G21q=Bq¯,h20+2Bq,h11+Cq,q,q¯.Here letterI represents the unit matrix with rank 2. It is quite apparent that equality (66) admits a solution:(67)G21=<p,Bq¯,h20+2Bq,h11+Cq,q,q¯>.At last, the first Lyapunov coefficient is calculated as(68)l1=12ReG21=14βB21+A12+3A30+3B03β+A11+2B02A02+A11A20−2A20B20−B11B02+B20.Thus, the first Lyapunov numberσ for the focus of planar system (57) is given by the formula σ=6π/βl1. Since above expression is much too complicated, we need to present some numerical simulations and figures around the point E∗ with computer simulations.Theorem 8 (nondegenerate Hopf bifurcation).
Assume that the equilibriumE∗ exists and parameters satisfy condition (54) and σ≠0; then, system (47) undergoes a nondegenerate Hopf bifurcation around this equilibrium as parameter δ=d passes through the critical value dH. The Hopf bifurcation is supercritical (subcritical), the interior equilibrium E∗ is a multiple stable (unstable) focus with multiplicity one, and limit cycles are stable (unstable) if σ<0 (σ>0).Remark 2.
The first Lyapunov coefficientl1 can be defined by 1/2βReG21 at times. On occasion, there are times when our system (55) exists some values of parameters such that σ=0 or the system may undergoes a degenerate Hopf bifurcation. Accompanied by a proper transformation, for planar system (55), the first Lyapunov number is given by the following formula [22]:(69)σ1=−3π2a01A23/2a10b10a112+a11b02+a02b11+a10a01b112+a20b11+a11b20+b102a11a02+2a02b02−2a10b10b022−a20a02−2a10a01a202−b20b02−a0122a20b20+b11b20+a01b10−2a102b11b02−a11a20−a102+a01b103b10b03−a01a30+2a10a21+b12+b10a12−a01b21,which is an advanced version of σ or l1.
Finally, for the Hopf bifurcation, we numerically give following examples to simulate how the parameterd controls dynamical behavior of system (2a) and system (2b).Example 1.
Consider Case 6 in Section3.1 and Theorem 4 again. Figure 2(a) reveals the existence of this interior equilibrium E8 in Case 6 by using the cobweb method. To investigate how the control parameter d affects dynamical behavior of our system (2a) and system (2b), Figures 2(c) and 3 depict phase diagrams corresponding to values d=0.0595 and d=0.053, respectively. When d=0.0595, this interior equilibrium is an asymptotically stable focus since A1≈−0.00111352489, A2≈0.07340334276, and Δ∗≈−0.2936121311. When d=0.053, it gets A1≈0.00017832135, A2≈0.07359407961, and Δ∗≈−0.2943762866, and there exists a limit cycle (by using the Poincare–Bendixson theorem) around this unstable focus. Here we notice that(70)αd1−αd2≈α′d12Hd1−d2+o10−6.
This implies that the Hopf bifurcation occurs in system (2a) and system (2b) when d=d8H. The first Lyapunov coefficient l1≈−0.122574, and thus the Hopf bifurcation is supercritical and a limit cycle generated by the critical point is stable. This interior equilibrium is a multiple stable focus with multiplicity one.
Besides, for a perturbed system with sufficiently small parameter vectorδ=δ1,δ2 in a neighbourhood of the origin O in the parameter plane(71a)x˙=r1x1−xK1−q1xya+bx+cy−m1x−d+δ1x2,(71b)y˙=r2y1−yK2yy+e+q2xya+bx+cy−m2+δ2y,and Hopf bifurcation analysis with two bifurcation parameters d and m2, we let δ≠0 and suppose that E∗=x∗,y∗ is an interior equilibrium of above perturbed system, where y∗=y8+w, w is sufficiently small, and(72)x∗=−cy∗+aδ2+m2−r2y∗+em2+δ2K2+r2y∗2δ2+m2−r2b−q2y∗+em2+δ2b−q2K2+br2y∗2.
Substituting it intoA1 and A2, the solutions δ1=δ1w and δ2=δ2w can be directly solved, which are written as the form up to third order:(73)δ1≈−1.503836w+24.893383w2+Ow3,δ2≈3.309615w−56.757836w2+Ow3,and hence the slope k=limw⟶0δ2w/δ1w of approximation straight line of the Hopf bifurcation curve Hp (see Figure 2(b)) in a small neighbourhood of the origin in parameter plane is approximately about −2.200781. The supercritical Hopf bifurcation curve Hp of system (71a) and system (71b) at the point E8 is numerically defined by solution (73), i.e., Hp=δ|δsatisfies73, while the variables δ and w both ensure the existence of E∗ and(74)A2w≈0.073568+2.087448w−47.606989w2+Ow3>0.
As a matter of fact, we perceive the phenomenon that Theorem8 describes the special case once δ lies on the horizontal line δ2=0.Figure 2
Figures in Example1. (a) Existence of interior equilibrium E8 in Case 6: red isocline (10a) and blue isocline (10b). (b) Hopf bifurcation curve Hp (in red). (c) A stable focus with d=0.0595.
(a)(b)(c)Figure 3
Figures in Example1 with d=0.053. (a) Phase diagrams around an unstable focus. (b) Enlarged phase diagrams: limit cycle.
(a)(b)Example 2.
Here we set parameters as well as Case 3 in Section3.2. Figures 4 and 5 depict phase diagrams corresponding to values d1=2.109 and d2=2.09, respectively. One should notice that(75)αd1−αd2≈α′d12Hd1−d2+o10−6.
Similar to above example, system (2a) and system (2b) undergo a Hopf bifurcation when d passes through d12H. The first Lyapunov coefficient l1≈−0.002729 is also found to be negative. Thus, the Hopf bifurcation is supercritical and a limit cycle generated by the critical point is stable. The interior equilibrium is also a multiple stable focus with multiplicity one. Figure 6 is the Hopf bifurcation curve corresponding to system (69) with slope k≈−0.094044 and(76)Hp=δ|δ1≈−40.550983w+639.229323w2+Ow3,δ2≈3.813575w−4.617377w2+Ow3.Figure 4
A stable focus in Example2 with d=2.109.Figure 5
Figures in Example2 with d=2.09. (a) Phase diagrams around an unstable focus. (b) Enlarged phase diagrams: limit cycle.
(a)(b)Figure 6
Hopf bifurcation curveHp (in red) in Example 2.


## 5.1. Saddle-Node Bifurcation

The two trivial equilibriaE2kk=1,2 collide with each other and system (47) has a unique boundary equilibrium E2 when y>0, if r2>m2 and y1=y2, or(48)m2=m2SN≔r2+2r2eK21−1+K2e.Then, there is a chance of bifurcation around this higher-order singular point. Here we choosem2 as a bifurcation parameter and select eigenvector v=01 corresponding to the zero eigenvalue for matrix (30). The eigenvector corresponding to the zero eigenvalue for the transpose of matrix (30) is w=1w2, where(49)w2=−r1−m1−q1y2/a+cy2q2y2/a+cy2.Supposer1−m1−q1y2/a+cy2<0, then the following transversality conditions hold:(50)wTFm2E2,m2SN=−w2y2≠0,wTD2FE2,m2SNv,v=2r2w2K2eK2+e−1≠0.Thus, we have following theorem by using Sotomayor’s theorem [22, 23].Theorem 5 (saddle-node bifurcation).
Suppose that the pointE2 exists; if r1−m1−q1y2/a+cy2<0, then system (47) undergoes a saddle-node bifurcation around point E2 with respect to the bifurcation parameter m2.


## 5.2. Transcritical and Pitchfork Bifurcation

The equilibriumE21 changes its stability when r1 crosses the threshold r1TC≔m1+q1y1/a+cy1; in other words, E21 is a higher-order equilibrium when r1=r1TC. Let J22E21<0. Thus, we choose the parameter r1 as a bifurcation parameter and an eigenvector v=1v2 corresponding to the zero eigenvalue for the Jacobian matrix JE21 when r1=r1TC, where v2=−q2y1/J22E21a+cy1>0. The eigenvector corresponding to the zero eigenvalue for the transpose of matrix JE21 is w=10; then, the transversality conditions are(51)wTFr1E21,r1TC=0,wTDFr1E21,r1TCv=1,wTD2FE21,r1TCv,v=−2r1TCK1+2q1by1a+cy12−2d+2q1q2ay1a+cy13J22E21.Thus, we have the following theorem by using Sotomayor’s theorem [22, 23].Theorem 6 (transcritical bifurcation).
Suppose that the two axial equilibria (E2kk=1,2) coexist and J22E21<0; if a+cy1≥bK1 or J22E21≥−q2a/ba+cy1, then system (47) undergoes a transcritical bifurcation around the point E21 with respect to the bifurcation parameter r1.For the special case(52)wTD2FE21,r1TCv,v=0,there is a chance the system (47) undergoes a pitchfork bifurcation. We still use the bifurcation parameter r1 and eigenvectors v and w; then, the fourth transversality condition is(53)wTD3FE21,r1PFv,v,v=−6q1bcJ22y1+aJ22b+q2cJ22b−q2y1+abJ22y1cy1+a5J222≠0,where J22=J22E21, r1PF≔m1+q1y1/a+cy1. Thus, we have another theorem by using Sotomayor’s theorem [22, 23].Theorem 7 (pitchfork bifurcation).
Suppose that the two axial equilibria (E2kk=1,2) coexist and J22E21<0; if condition (52) and J22E21≠−aq2/bcy1+a hold, then system (47) undergoes a pitchfork bifurcation around the point E21 with respect to the bifurcation parameter r1.


## 5.3. Hopf Bifurcation

In this section, we consider the Hopf bifurcation of system (2a) and system (2b). Here the point E∗ exists and we choose d as bifurcation parameter. Suppose that λd=αd±iωd are a pair of conjugate eigenvalues of matrix JE∗, where αd=1/2A1d. The critical value dH satisfies(54)αdH=0,α′dH≠0,A2dH>0.Then, system (47) undergoes a Hopf bifurcation around the point E∗ with respect to the bifurcation parameter d.We will calculate the first Lyapunov numberσ at the point E∗, which is used to determine the stability of limit cycles and Hopf bifurcation direction. The method and calculations of the first and second Lyapunov coefficients can be found in [22, 24, 25]. Therefore, translating the point E∗ to the origin O=0,0 by a linear transformation (I): X=x−x∗, Y=y−y∗, system (47) in power series around the origin is(55)X˙=f1X,Y=∑1≤i+j≤3aijXiYj+OX,Y4,Y˙=g1X,Y=∑1≤i+j≤3bijXiYj+OX,Y4,where coefficients are(56)a10=−r1s1K1−ds1+q1bs1s2bs1+cs2+a2,a01=−q1s1bs1+abs1+cs2+a2,a20=−r1K1+q1s2bcs2+abs1+cs2+a3−d,a02=q1s1cbs1+abs1+cs2+a3,a11=−a2+bs1+cs2a+2s1s2bcq1bs1+cs2+a3,a30=−q1s2b2cs2+abs1+cs2+a4,a03=−q1s1c2bs1+abs1+cs2+a4,a21=bs12cs2+ab−c2s22+a2q1bs1+cs2+a4,a12=cs22bs1+ac−b2s12+a2q1bs1+cs2+a4;b10=q2s2cs2+abs1+cs2+a2,b01=−r2s2−K2e+2es2+s22K2s2+e2−q2cs1s2bs1+cs2+a2,b20=−q2s2bcs2+abs1+cs2+a3,b02=r2K2−3s2e2−3s22e−s23K2s2+e3−q2s1cbs1+abs1+cs2+a3,b11=a2+bs1+cs2a+2s1s2bcq2bs1+cs2+a3,b30=q2s2b2cs2+abs1+cs2+a4,b03=−r2e2K2+eK2s2+e4+q2s1c2bs1+abs1+cs2+a4,b21=−bs12cs2+ab−c2s22+a2q2bs1+cs2+a4,b12=−cs22bs1+ac−b2s12+a2q2bs1+cs2+a4,and OX,Y4 stands for some smooth functions. After that, by using a transformation (II): u=Y, v=−Xb10+Ya10/β with β=A2E∗>0, the above system becomes a standard form:(57)u˙=f2u,v=−βv+∑2≤i+j≤3Aijuivj+Ou,v4,v˙=g2u,v=βu+∑2≤i+j≤3Bijuivj+Ou,v4,where(58)A20=a102b20+a10b10b11+b02b102b102,A11=−2a10b20+b10b11βb102,A02=−b20a01b10+a102b102,A30=a103b30+a102b10b21+a10b102b12+b03b103b103,A21=−3a102b30+2a10b10b21+b102b12βb103,A12=−3a10b30+b10b21a01b10+a102b103,A03=−b30β3b103;B20=−a02b103−a103b20+a102a20b10−a102b10b11+a10a11b102−a10b02b102βb102,B11=−2a102b20−2a10a20b10+a10b10b11−a11b102b102,B02=a10b20−a20b10βb102,B30=−a03b104−a104b30+a103a30b10−a103b10b21+a102a21b102−a102b102b12+a10a12b103−a10b03b103βb103,B21=−3a103b30−3a102a30b10+2a102b10b21−2a10a21b102+a10b102b12−a12b103b103,B12=β3a102b30−3a10a30b10+a10b10b21−a21b102b103,B03=a10b30−a30b10a01b10+a102b103.Letp=q=1/2i1 be two corresponding eigenvectors of a matrix A such that Aq=iβq, A⊤p=−iβp and <p,q>=1; the operation <x,y>=x†y (x,y∈ℂ) with the Hermitian transpose (upper symbol) † represents the usual inner product, and the Jacobian matrix is(59)A=0−ββ0.We should rewrite the functions in the right hand side of system (57) in the form of power series (60)F2x=f2xg2x=Ax+12!Bx,x+13!Cx,x,x+Ox4,where the components of linear functions B and C are(61)Bix,y=∑j,k=12∂2F2i0∂ξj∂ξkxjyk,Cix,y,z=∑j,k,l=12∂3F2i0∂ξj∂ξk∂ξlxjykzl,i=1,2,while x is the two-dimensional Euclidean norm of x. Define Tc be the largest subspace invariant by the matrix A and the generalized eigen subspace corresponding to the pair of purely imaginary eigenvalues ±iβ, i.e., Tc=spanq,q¯. That is to say, for any element z∈Tc, there must exist a linear expansion z=wq+w¯q¯. Now we can construct a two-dimensional parameterized center manifold(62)H=Hw,w¯=wq+w¯q¯+∑j+k≥2∞1j!·k!hjkwjw¯k,in which hjk∈ℂ2 and its complex conjugate is hjk¯=hkj. Combining equations (57) and Hw,w¯, we arrive at a complex equation without consideration:(63)∂H∂ww˙+∂H∂w¯w¯˙=F2H.The “chart”w for the central manifold H should be extracted from following differential equation:(64)w˙=iβw+12G21ww2+⋯.So, substituting it into equation (63) and comparing coefficients of these wjw¯k, we recursively derive(65)h11=−A−1Bq,q¯,h20=2iβI−A−1Bq,q,and an equality from coefficient of the cubic term ww2:(66)iβI−Ah21+G21q=Bq¯,h20+2Bq,h11+Cq,q,q¯.Here letterI represents the unit matrix with rank 2. It is quite apparent that equality (66) admits a solution:(67)G21=<p,Bq¯,h20+2Bq,h11+Cq,q,q¯>.At last, the first Lyapunov coefficient is calculated as(68)l1=12ReG21=14βB21+A12+3A30+3B03β+A11+2B02A02+A11A20−2A20B20−B11B02+B20.Thus, the first Lyapunov numberσ for the focus of planar system (57) is given by the formula σ=6π/βl1. Since above expression is much too complicated, we need to present some numerical simulations and figures around the point E∗ with computer simulations.Theorem 8 (nondegenerate Hopf bifurcation).
Assume that the equilibriumE∗ exists and parameters satisfy condition (54) and σ≠0; then, system (47) undergoes a nondegenerate Hopf bifurcation around this equilibrium as parameter δ=d passes through the critical value dH. The Hopf bifurcation is supercritical (subcritical), the interior equilibrium E∗ is a multiple stable (unstable) focus with multiplicity one, and limit cycles are stable (unstable) if σ<0 (σ>0).Remark 2.
The first Lyapunov coefficientl1 can be defined by 1/2βReG21 at times. On occasion, there are times when our system (55) exists some values of parameters such that σ=0 or the system may undergoes a degenerate Hopf bifurcation. Accompanied by a proper transformation, for planar system (55), the first Lyapunov number is given by the following formula [22]:(69)σ1=−3π2a01A23/2a10b10a112+a11b02+a02b11+a10a01b112+a20b11+a11b20+b102a11a02+2a02b02−2a10b10b022−a20a02−2a10a01a202−b20b02−a0122a20b20+b11b20+a01b10−2a102b11b02−a11a20−a102+a01b103b10b03−a01a30+2a10a21+b12+b10a12−a01b21,which is an advanced version of σ or l1.
Finally, for the Hopf bifurcation, we numerically give following examples to simulate how the parameterd controls dynamical behavior of system (2a) and system (2b).Example 1.
Consider Case 6 in Section3.1 and Theorem 4 again. Figure 2(a) reveals the existence of this interior equilibrium E8 in Case 6 by using the cobweb method. To investigate how the control parameter d affects dynamical behavior of our system (2a) and system (2b), Figures 2(c) and 3 depict phase diagrams corresponding to values d=0.0595 and d=0.053, respectively. When d=0.0595, this interior equilibrium is an asymptotically stable focus since A1≈−0.00111352489, A2≈0.07340334276, and Δ∗≈−0.2936121311. When d=0.053, it gets A1≈0.00017832135, A2≈0.07359407961, and Δ∗≈−0.2943762866, and there exists a limit cycle (by using the Poincare–Bendixson theorem) around this unstable focus. Here we notice that(70)αd1−αd2≈α′d12Hd1−d2+o10−6.
This implies that the Hopf bifurcation occurs in system (2a) and system (2b) when d=d8H. The first Lyapunov coefficient l1≈−0.122574, and thus the Hopf bifurcation is supercritical and a limit cycle generated by the critical point is stable. This interior equilibrium is a multiple stable focus with multiplicity one.
Besides, for a perturbed system with sufficiently small parameter vectorδ=δ1,δ2 in a neighbourhood of the origin O in the parameter plane(71a)x˙=r1x1−xK1−q1xya+bx+cy−m1x−d+δ1x2,(71b)y˙=r2y1−yK2yy+e+q2xya+bx+cy−m2+δ2y,and Hopf bifurcation analysis with two bifurcation parameters d and m2, we let δ≠0 and suppose that E∗=x∗,y∗ is an interior equilibrium of above perturbed system, where y∗=y8+w, w is sufficiently small, and(72)x∗=−cy∗+aδ2+m2−r2y∗+em2+δ2K2+r2y∗2δ2+m2−r2b−q2y∗+em2+δ2b−q2K2+br2y∗2.
Substituting it intoA1 and A2, the solutions δ1=δ1w and δ2=δ2w can be directly solved, which are written as the form up to third order:(73)δ1≈−1.503836w+24.893383w2+Ow3,δ2≈3.309615w−56.757836w2+Ow3,and hence the slope k=limw⟶0δ2w/δ1w of approximation straight line of the Hopf bifurcation curve Hp (see Figure 2(b)) in a small neighbourhood of the origin in parameter plane is approximately about −2.200781. The supercritical Hopf bifurcation curve Hp of system (71a) and system (71b) at the point E8 is numerically defined by solution (73), i.e., Hp=δ|δsatisfies73, while the variables δ and w both ensure the existence of E∗ and(74)A2w≈0.073568+2.087448w−47.606989w2+Ow3>0.
As a matter of fact, we perceive the phenomenon that Theorem8 describes the special case once δ lies on the horizontal line δ2=0.Figure 2
Figures in Example1. (a) Existence of interior equilibrium E8 in Case 6: red isocline (10a) and blue isocline (10b). (b) Hopf bifurcation curve Hp (in red). (c) A stable focus with d=0.0595.
(a)(b)(c)Figure 3
Figures in Example1 with d=0.053. (a) Phase diagrams around an unstable focus. (b) Enlarged phase diagrams: limit cycle.
(a)(b)Example 2.
Here we set parameters as well as Case 3 in Section3.2. Figures 4 and 5 depict phase diagrams corresponding to values d1=2.109 and d2=2.09, respectively. One should notice that(75)αd1−αd2≈α′d12Hd1−d2+o10−6.
Similar to above example, system (2a) and system (2b) undergo a Hopf bifurcation when d passes through d12H. The first Lyapunov coefficient l1≈−0.002729 is also found to be negative. Thus, the Hopf bifurcation is supercritical and a limit cycle generated by the critical point is stable. The interior equilibrium is also a multiple stable focus with multiplicity one. Figure 6 is the Hopf bifurcation curve corresponding to system (69) with slope k≈−0.094044 and(76)Hp=δ|δ1≈−40.550983w+639.229323w2+Ow3,δ2≈3.813575w−4.617377w2+Ow3.Figure 4
A stable focus in Example2 with d=2.109.Figure 5
Figures in Example2 with d=2.09. (a) Phase diagrams around an unstable focus. (b) Enlarged phase diagrams: limit cycle.
(a)(b)Figure 6
Hopf bifurcation curveHp (in red) in Example 2.


## 6. Conclusions

In summary, we firstly considered stability analysis and bifurcations of system (2a) and system (2b) with B-D functional response and Allee-like effect, which is a modified version of a predator-prey system in [1]. The polynomial’s method, derived from eliminants Ryf,g and Rxf,g, can be extended to more complicated polynomial systems. Some conclusions are the same as reference [1], such as the uniform boundedness (Theorem 1), the existence of equilibria E∗, the nonexistence (Theorem 4) of limit cycles, and so on. It is supposed that some methods and conclusions can be available in original system (1a) and system (1b), such as pitchfork bifurcation. Lemmas 1 and 2 are available in more complicated predator-prey systems. Some critical cases, such as r1−m1=q1y1/a+cy1 and A2=0, need to be researched further with the help of topologically equivalent systems or the “blow-up” method (horizontal and vertical blow-ups). Other parameters can also be considered as a bifurcation parameter δ in Hopf bifurcation, although it is described by Hopf bifurcation curve Hp, for instance δ=x∗. Notice that Theorem 2 (permanence) can be extended to its reaction-diffusive version [26, 27]. Under the conditions or the discussion of Theorem 3, the Turing instability of its corresponding reaction-diffusion system subject to the homogeneous Neumann boundary conditions [28]:(77a)ut−D1uxx=r1u1−uK1−q1uva+bu+cv−m1u−du2,(77b)vt−D2vxx=r2v1−vK2vv+e+q2uva+bu+cv−m2v,(77c)∂νu=∂νv=0,t≥0,x∈∂Ω,(77d)u0,x=u0x≥0≢0,v0,x=v0x≥0≢0,x∈Ω,will not occur, and thus it is available to find potential Hopf bifurcation points and consider transversality conditions if we choose u∗ as Hopf bifurcation parameter. Here D1 and D2 are two positive diffusive constants; Δ is the Laplacian operator; u=ut,x and v=vt,x are the densities of prey and predator, respectively; Ω is an one-dimensional bounded domain with smooth boundary ∂Ω; the symbol ∂ν is the outer flux, and no flux boundary condition is imposed; thus, the system is closed [29]; the admissible initial functions u0x and v0x are all continuous functions on Ω¯; to describe an environment surrounded by dispersal barriers, we take zero flux at ∂Ω [28]. The method of calculating the first Lyapunov coefficient l1 in Section 5.3 is a reference for deducing of Hopf bifurcation direction in reaction-diffusion systems [29, 30].Though the dynamical behavior of predator-prey systems in single species or multispecies has been researched by many previous literatures, we still need further study in biomathematics, especially the phytoplankton and zooplankton systems. Meanwhile, how to keep ecosystems in balance or coexistence states and avoid harmful effect is our next direction in this area.



---

*Source: 5618190-2021-02-15.xml*