The Dynamical Behavior of a Predator-Prey System with Holling Type II Functional Response and Allee Effect

Shuangte Wang¹, Hengguo Yu^1,2, Chuanjun Dai^2,3, Min Zhao^2,3

¹College of Mathematics and Physics, Wenzhou University, Wenzhou, China.
²Key Laboratory for Subtropical Oceans & Lakes Environment & Biological Resources Utilization Technology of Zhejiang, Wenzhou University, Wenzhou, China.
³School of Life and Environmental Science, Wenzhou University, Wenzhou, China.

DOI: 10.4236/am.2020.115029                

Abstract

In this paper, we mainly considered the dynamical behavior of a predator-prey system with Holling type II functional response and Allee-like effect on predator, including stability analysis of equilibria and Hopf bifurcation. Firstly, we gave some sufficient conditions to guarantee the existence, the local and global stability of equilibria as well as non-existence of limit cycles. By using the cobweb model, some cases about the existence of interior equilibrium are also illustrated with numerical outcomes. These existence and stability conclusions of interior equilibrium are also suitable in corresponding homogeneous reaction-diffusion system subject to the Neumann boundary conditions. Secondly, we theoretically deduced that our system has saddle-node bifurcation, transcritical bifurcation and Hopf bifurcation under certain conditions. Finally, for the Hopf bifurcation, we choose d as the bifurcation parameter and presented some numerical simulations to verify feasibility and effectiveness of the theoretical derivation corresponding to the existence of y_k, respectively. The Hopf bifurcations are supercritical and limit cycles generated by the critical points are stable.

Keywords

Predator-Prey System, Holling Type II Functional Response, Allee Effect, Stability, Hopf Bifurcation

1. Introduction

In this paper, we consider a predator-prey system with Holling type II functional response and Allee effect on predator, which is described by the following nonlinear ordinary differential equations (ODEs)

$\stackrel{\dot{}}{x}=r_{1}x\left(1-\frac{x}{K_1}\right)-\frac{qxy}{a+x}-m_{1}x-d{x}^2:=f\left(x,y\right),$ (1a)

$\stackrel{\dot{}}{y}=r_{2}y\left(1-\frac{y}{K_2}\right)\frac{y}{y+e}+\frac{e_{1}qxy}{a+x}-m_{2}y:=g\left(x,y\right)$ (1b)

subject to initial conditions $x\left(0\right)\mathrm{,}y\left(0\right)\ge 0$. Here, functions $x=x\left(t\right)$ and $y=y\left(t\right)$ are the prey and predator densities at time t, respectively. All above positive constants have biological considerations. Parameters $r_1$ and $r_2$ denote the intrinsic growth rate of prey and predator, respectively; $K_1$ and $K_2$ represent the carrying capacity of the environment for prey and predator, respectively; a is the half-saturation constant; q is the search efficiency of predator for prey; $m_1$ and $m_2$ are the mortality rate of prey and predator species, respectively; $e_1$ is the biomass conversion; d is the intra-specific competition coefficient; e is the Allee effect constant. For convenience, we denote $e_{1}q$ as c. The specific growth term $r_{1}x\left(1-\frac{x}{K_1}\right)$ governs the increase of prey in the lack of predator, while the specific growth term $r_{2}y\left(1-\frac{y}{K_2}\right)$ governs the increase of predator. The square term $d{x}^2$ denotes intrinsic decrease of prey. The coupled term $\frac{xy}{a+x}$, named Holling type II functional response, describes how predators transform harvested prey into the growth of itself, and also refers to the change in the density of prey attached per unit time per predator as prey density changes. It was proposed by C.S. Holling in 1965 as well as other Holling type functional responses and extended the classical Lotka-Volterra predator-prey systems in biomathematics [1] [2] [3] [4].

The term $\frac{y}{y+e}$ on predator is called the Allee effect. The Allee effect, named by Warder Clyde Alee, plays a significant role in determining dynamical behavior of predator-prey systems, even ecological and social models [5] [6] [7]. These phenomena occur in small or sparse populations and are widely accepted to be common in nature. They are related to a positive correlation between population size and fitness at very low population size and another phenomenological feature that may induce extinction of populations [8] [9] [10]. The Allee effects are called strong if they cause critical population sizes, while they are called weak if they do not result in critical sizes [7]. The strong Allee effect causes the so-called Allee threshold that populations need to surpass it in order to avoid extinction and survive permanently [9] [11]. On the contrary, a population with weak Allee effect does not have such above threshold [12]. Other names are positive density dependence (in contrast to classical negative density dependence) and depensatory dynamics (in contrast to classical compensatory dynamics) [7]. Most references concentrate on the strong Allee effect.

Like Holling type functional responses, predator-prey systems incorporating Allee effect have also received much attention from ecologists and mathematicians [13] - [18]. For Allee effect incorporating into the prey population with multiply form, in reference [13], the authors extended a predator-prey system with strong Allee effect and a functional response of the ratio of prey to predator. By means of bifurcation analysis and advance numerical techniques for the calculation of invariant manifolds of equilibria, they stated the consequences of the (dis)appearance of limit cycles, homoclinic orbits, and heterclinic connections in the global arrangement of phase plane near a Bogdanov-Takens bifurcation. The reference [14] examined global behavior of a Gause-type predator-prey system with Holling type III functional response and weak Allee effect on the prey growth. They also proved that the origin is a saddle point and obtained existence of two limit cycles surround a stable interior equilibria: just like that with strong Allee effect. In [15], Zu et al. studied a predator-prey system with Allee effect on prey and investigated local asymptotic stability of equilibria. Meanwhile, the authors found that Allee effect of prey population can bring about unstable or stable periodic fluctuations. Based on above reference, Zu and Mimura [16] considered the Holling type II functional response and local asymptotic stability of equilibria. Compared with the model without Allee effect, they found that the effect on prey increases extinction risk of predator and prey. An explicit algorithm was also obtained to determine the direction of Hopf bifurcation as well as the stability of periodic solutions. The phenomenon of periodicity in this reference is similar to that in [15].

While the reference [17] showed that, when the Leslie-Gower predator-prey model with additive Allee effect on prey has two positive equilibria, there exists a separatrix curve that separates the behavior of trajectories, i.e. the model is highly sensitive to initial conditions. In [18], the authors began with local population dynamics and then constructed a model including both local population and metapopulation dynamics. Their results indicated that the Allee-like effect in a metapopulation may emerge from the imposed Allee effect at the local population level and severe demographic stochasticity may compound the metapopulation extinction risk posed by the Allee effect.

The Allee effect was also researched in predator-prey systems with delay, impulse or diffusion. Xiao et al. [19] considered local asymptotic stability and Hopf bifurcation of a Holling type II predator-prey model incorporating time delay and Allee effect in prey. By using the delay as a bifurcation parameter, they showed that if the birth rate is small enough or the Allee effect is large enough, then both prey and predator extinct, i.e. the Allee effect can influence stability of the system. In [20], by using the Mawhins continuation theorem of coincidence degree and analysis techniques, the authors obtained some sufficient conditions of the existence of periodic solutions in a nonautonomous predator-prey system with Holling type II functional response, strong Allee effect and impulsive perturbation. They proved that their system has at least one positive $\omega$ -periodic solution. While in reference [21], Cui et al. considered a diffusive predator-prey system with strong Allee effect and a protection zone for the prey. If the Allee effect threshold is low and the protection zone is large, they showed that the over exploitation phenomenon can be avoided.

This paper mainly concentrate on the dynamical behavior analysis of a complicated and realistic predator-prey system (1) with Allee-like effect [18] [22] on the specific growth term of predator with multiply form and an intrinsic decrease term on prey, which is different from above references involving ODEs systems with the additive Allee effect on prey or Allee-like effect on prey. The rest of this paper is organized as follows. Preliminaries, such as boundedness and permanence, are given in Section 2. In section 3, we give sufficient conditions for stability analysis of equilibria by using linearization technique and non-existence of limit cycles. In Section 4, bifurcation analysis and numerical simulations of Hopf bifurcation are presented. In Section 5, we present summary and some remarks.

2. Preliminaries

In this section, we devote to give priori foundations for our system. Firstly, we denote the first quadrant as $R^{+2}$, and its closure is denoted as $R_{+}^2=\overline{R^{+2}}$. For biological consideration, the system (1) is defined on the domain $R_{+}^2$ and all the solutions are nonnegative with initial conditions $x\left(t_0\right)\mathrm{,}y\left(t_0\right)\ge 0$, i.e. $R_{+}^2$ is an invariant set. Thus, any trajectory starting from $R_{+}^2$ cannot cross the coordinate axes. Furthermore, all the solutions are bounded. Now we will prove following theorems.

Theorem 1 (Uniform boundedness) All the solutions of system (1) are uniformly bounded in $R_{+}^2$ with initial conditions $x\left(t_0\right)\mathrm{,}y\left(t_0\right)\ge 0$.

Proof of theorem 1. Here we introduce an auxiliary function $z=e_{1}x+y$ and a constant $m=\min \left\{m_1\mathrm{,}{m}_2\right\}$, it is obvious to see that

$\begin{array}{c}\stackrel{\dot{}}{z}\le e_1\left[r_{1}x\left(1-\frac{x}{K_1}\right)-m_{1}x\right]+\left[\frac{r_2{y}^2}{K_2\left(y+e\right)}\left(K_2-y\right)-m_{2}y\right]\\ \le \frac{1}{4}\left(e_1{r}_1{K}_1+r_2{K}_2\right)-mz\mathrm{.}\end{array}$

By applying the theory of differential inequality [23], we obtain

$\underset{t\to \infty }{\lim \sup }z\le \frac{1}{4m}\left(e_1{r}_1{K}_1+r_2{K}_2\right).$ (2)

All the solutions of the system (1) are confined in the region

$\left\{\left(x,y\right)|0\le z\left(x,y\right)\le \frac{1}{4m}\left(e_1{r}_1{K}_1+r_2{K}_2\right)+ϵ,\text{for}\text{any}ϵ>0\right\}.$

Thus we complete the proof. □

Theorem 2 (Permanence) If parameters satisfy $\left(c-m_2\right)\left(1-\lambda \right){\omega }_1>m_{2}a$ and

${\omega }_1=\frac{r_1-m_1-\frac{q{M}_2}{a}}{\frac{r_1}{K_1}+d}>0,$ (3)

where $\lambda \in \left(\mathrm{0,1}\right)$ and $M_2$ is a positive upper bound of y, then the system (1) is permanence in $R_{+}^2$.

Proof of theorem 2. From the theorem 1, there exist a positive upper bound ${\xi }_1$, such that

$\max \left\{\underset{t\to \infty }{\lim \sup }x\mathrm{,}\underset{t\to \infty }{\lim \sup }y\right\}\le {\xi }_1\mathrm{.}$ (4)

For the Equation (1a), we have $M_2>0$ and sufficiently large $T_0$, such that

$\stackrel{\dot{}}{x}\ge x\left[\left(r_1-m_1-\frac{q{M}_2}{a}\right)-\left(\frac{r_1}{K_1}+d\right)x\right]\mathrm{,}\forall t\ge T_0\mathrm{.}$

Similarly, by using lemmas in [24], we have $\underset{t\to \infty }{\lim \inf }x\ge {\omega }_1$. That is to say, there exist sufficiently large T, such that $x\ge \left(1-\lambda \right){\omega }_1$, $\forall t\ge T$. By substituting it into Equation (1b), we derive an inequality

$\stackrel{\dot{}}{y}\ge -\frac{r_2}{K_2}{y}^2+\frac{\left(1-\lambda \right){\omega }_1\left(c-m_2\right)-a{m}_2}{\left(1-\lambda \right){\omega }_1+a}y\mathrm{.}$

By using above lemmas again, we complete the proof.

3. Equilibria

In this section we will discuss equilibria of system (1) with their existence conditions and stability analysis. It is obvious that the system (1) has equilibria: $E_0\mathrm{<mo>\; :\; </mo>}=\left(\mathrm{0,0}\right)$, $E_1\mathrm{<mo>\; :\; </mo>}=\left(x_1\mathrm{,0}\right)$, $E_2^{\left(k\right)}\mathrm{<mo>\; :\; </mo>}=\left(\mathrm{0,}{y}_k\right)$, $k=\mathrm{1,2}$, where $x_1\mathrm{<mo>\; :\; </mo>}=\frac{r_1-m_1}{\frac{r_1}{K_1}+d}$ and

$y_{1,2}:=\frac{r_2-m_2\pm \sqrt{\Delta }}{\frac{2r_2}{K_2}},\Delta ={\left(r_2-m_2\right)}^2-\frac{4m_{2}e{r}_2}{K_2}.$ (5)

The point $E_1$ enquiries $r_1>m_1$, while the points $E_2^{\left(k\right)}$ all exist when $r_2>m_2$ and $\Delta >0$. If $r_2>m_2$ and $\Delta =0$, then $y_1=y_2$, i.e. the two equilibria $E_2^{\left(1\right)}$ and $E_2^{\left(2\right)}$ collide with each other and the system (1) has an unique boundary equilibrium when $y>0$. In this special case, we denote it as $E_2=\left(\mathrm{0,}{y}_2\right)$.

3.1. Existence of Interior Equilibria

Here we firstly denote interior equilibrium as $E_{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; :\; </mo>}=\left(x_{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}_{\mathrm{<mo>\; *\; </mo>}}\right)$ or $\left(s_1\mathrm{,}{s}_2\right)$. The point $E_{\mathrm{<mo>\; *\; </mo>}}$ must satisfies following algebraic equations $f\left(x,y\right)=g\left(x,y\right)=0$ or

$r_1\left(1-\frac{x}{K_1}\right)-\frac{qy}{a+x}-m_1-dx=\mathrm{0,}$ (6a)

$r_2\left(1-\frac{y}{K_2}\right)\frac{y}{y+e}+\frac{cx}{a+x}-m_2=0.$ (6b)

The Equation (6a) implies that $r_1>m_1$, and we have: $y\to \frac{\left(r_1-m_1\right)a}{q}$ as $x\to 0$ ; $x\to x_1$ as $y\to 0$. For the Equation (6b), we have: $y\to y_k$ ( $k=1$ or 2) as $x\to 0$ ; $x\to \frac{a{m}_2}{c-m_2}$ or $\infty$ as $y\to 0$. The derivative of implicit function $y=y\left(x\right)$ from the Equation (6b) is

$y^{\prime }\left(x\right)=\frac{K_{2}ac{\left(y+e\right)}^2}{r_2\left[y\left(y+2e\right)-K_{2}e\right]{\left(a+x\right)}^2}\mathrm{,}$ (7)

thus we denote the positive root of the quadratic equation $y\left(y+2e\right)-K_{2}e=0$ as $y_{s_1}=-e+\sqrt{e^2+K_{2}e}$. Based on above analysis, some cases will be given to illustrate the existence of the interior equilibria $E_{\mathrm{<mo>\; *\; </mo>}}$ when the curve of left hand side function from the Equation (6b) falls in $R_{+}^2$.

Case 1. If following conditions

$0<\frac{\left(r_1-m_1\right)a}{q}<y_k<y_{s_1},0<\frac{a{m}_2}{c-m_2}<x_1$ (8)

hold, then an interior equilibrium exist. For instance, let $r_1=1$, $r_2=1$, $K_1=20$, $K_2=80$, $q=0.3$, $c=15$, $a=1$, $m_1=0.3$, $m_2=0.5$, $e=9$ and $d=6.497942687$, then equilibrium $E_3\approx \left(\mathrm{0.02257205799,1.882212454}\right)$.

Case 2. If following condition

$y_{s_1}<y_k<\frac{\left(r_1-m_1\right)a}{q}$ (9)

holds, then an interior equilibrium exist. For instance, let $r_1=1$, $r_2=1$, $K_1=20$, $K_2=80$, $q=0.01$, $c=15$, $a=1$, $m_1=0.3$, $m_2=0.5$, $e=9$ and $d=6.5$, then equilibrium $E_4\approx \left(\mathrm{0.01846631478,58.97384772}\right)$.

Case 3. If following conditions

$0<y_k<\frac{\left(r_1-m_1\right)a}{q}<y_{s_1},x_1<\frac{a{m}_2}{c-m_2}$ (10)

hold, then an interior equilibrium exist. For instance, let $r_1=1$, $r_2=6$, $K_1=20$, $K_2=80$, $q=1$, $c=15$, $a=10$, $m_1=0.3$, $m_2=0.5$, $e=9$ and $d=6.5$, then equilibrium $E_5\approx \left(\mathrm{0.09827333849,0.5686300643}\right)$.

Case 4. If following conditions

$0<y_k<y_{s_1}<\frac{\left(r_1-m_1\right)a}{q},x_1<\frac{a{m}_2}{c-m_2}$ (11)

hold, then an interior equilibrium exist. For instance, let $r_1=1$, $r_2=6$, $K_1=20$, $K_2=80$, $q=0.1$, $c=15$, $a=10$, $m_1=0.3$, $m_2=0.5$, $e=9$ and $d=6.5$, then equilibrium $E_6\approx \left(\mathrm{0.1060408106,0.5490299197}\right)$.

Case 5. If following conditions

$c=m_2,0<y_k<\frac{\left(r_1-m_1\right)a}{q}<y_{s_1}$ (12)

hold, then an interior equilibrium exist. For instance, let $r_1=1$, $r_2=0.8$, $K_1=20$, $K_2=10$, $q=0.3$, $c=0.05$, $a=1$, $m_1=0.3$, $m_2=0.05$, $e=10$, $d=0.4$, then equilibrium $E_7\approx \left(\mathrm{1.484367905,0.2652844730}\right)$.

Case 6. If $y_k$ not exist and following condition

$0<\frac{a{m}_2}{c-m_2}<x_1$ (13)

holds, then an interior equilibrium exist. For instance, let $r_1=0.4$, $r_2=0.2$, $K_1=2.5$, $K_2=0.5$, $q=0.3$, $c=0.4$, $a=0.22$, $m_1=0.35$, $m_2=0.15$, $e=0.1$ and $d=0.1750056774$, then equilibrium $E_8\approx \left(\mathrm{0.09221932871,0.01988420968}\right)$.

Case 7. If $y_1\mathrm{<mo>\; <\; </mo>0}$ and the condition (13) holds, then an interior equilibrium exist. For instance, letting $r_1=0.4$, $r_2=0.2$, $K_1=10$, $K_2=20$, $q=0.3$, $c=1$, $a=0.4$, $m_1=0.3$, $m_2=0.5$, $e=0.4$ and $d=0.1220933748$, then equilibrium $E_9\approx \left(\mathrm{0.3360884070,0.1116947808}\right)$.

From above cases, we will also give some cases to illustrate the phenomena that two interior equilibria exist at the same time.

Case 8. If following conditions

$0<y_2<y_{s_1}<y_1<\frac{\left(r_1-m_1\right)a}{q},x_1<\frac{a{m}_2}{c-m_2}$ (14)

hold, then two interior equilibria coexist. For instance, let $r_1=1$, $r_2=1$, $K_1=20$, $K_2=80$, $q=0.05$, $c=1$, $a=3$, $m_1=0.3$, $m_2=0.5$, $e=4.24$ and $d=1$, then equilibrium $E_{10}^{\left(1\right)}\approx \left(\mathrm{0.08140728865,37.87187278}\right)$ and $E_{10}^{\left(2\right)}\approx \left(\mathrm{0.6389354884,2.119151346}\right)$. Here $y_{10}^{\left(1\right)}>y_1$ and $y_{10}^{\left(2\right)}<y_2$.

Case 9. If following conditions

$0<y_2<y_{s_1}<y_1<\frac{\left(r_1-m_1\right)a}{q},c=m_2$ (15)

hold, then two interior equilibria coexist. For instance, let $r_1=1$, $r_2=1$, $K_1=20$, $K_2=80$, $q=0.05$, $c=0.5$, $a=3$, $m_1=0.3$, $m_2=0.5$, $e=4.24$ and $d=1$, then equilibrium $E_{11}^{\left(1\right)}\approx \left(\mathrm{0.1008042045,36.84720274}\right)$ and $E_{11}^{\left(2\right)}\approx \left(\mathrm{0.6244397702,3.214025704}\right)$. Here $y_{11}^{\left(1\right)}>y_1$ and $y_{11}^{\left(2\right)}<y_2$.

3.2. Stability Analysis

In this subsection, we use the Routh-Hurwitz criterion and the Perron’s theorems to analyze local stability of above equilibria. Recall the system (1) again, its Jacobian matrix takes the following form $J={\left(J_{ij}\right)}_{2\times 2}$, where components are

$\begin{array}{l}{J}_{11}=r_1-m_1-2x\left(\frac{r_1}{K_1}+d\right)-\frac{qay}{{\left(a+x\right)}^2},J_{12}=-\frac{qx}{a+x},J_{21}=\frac{cay}{{\left(a+x\right)}^2},\\ J_{22}=\frac{-r_{2}y}{K_2{\left(y+e\right)}^2}\left[2y^2+\left(3e-K_2\right)y-2K_{2}e\right]+\frac{cx}{a+x}-m_2.\end{array}$

For the trivial equilibrium $E_0$ and its Jacobian matrix $J\left(E_0\right)$, we have: (a) $E_0$ is an asymptotically stable node if $r_1<m_1$ ; (b) $E_0$ is a saddle point if $r_1>m_1$.

Notice that Jacobian matrix at the axial equilibrium $E_1$ takes the form

$J\left(E_1\right)=\left[\begin{array}{cc}{m}_1-r_1& \frac{-q{x}_1}{x_1+a}\\ 0& \frac{c{x}_1}{x_1+a}-m_2\end{array}\right]\mathrm{.}$ (16)

Combining the existence condition of $E_1$, we have: (a) $E_1$ is a saddle point if $\frac{c{x}_1}{x_1+a}>m_2$ ; (b) $E_1$ is an asymptotically stable node if $\frac{c{x}_1}{x_1+a}<m_2$.

In the case that two axial equilibrium $E_2^{\left(k\right)}\left(k=\mathrm{1,2}\right)$ all exist, the Jacobian matrices are

$J\left(E_2^{\left(k\right)}\right)=\left[\begin{array}{cc}{r}_1-m_1-\frac{q{y}_k}{a}& 0\\ \frac{c{y}_k}{a}& \frac{\left(m_2-r_2\right)y_k+2m_{2}e}{y_k+e}\end{array}\right]\mathrm{.}$ (17)

Since $J_{22}\left(E_2^{\left(1\right)}\right)<0$, we have: (a) $E_2^{\left(1\right)}$ is a saddle point if $r_1-m_1>\frac{q{y}_1}{a}$ ; (b) $E_2^{\left(1\right)}$ is an asymptotically stable node if $r_1-m_1<\frac{q{y}_1}{a}$. For the equilibrium $E_2^{\left(2\right)}$, it is obvious that $J_{22}\left(E_2^{\left(2\right)}\right)>0$, we have: (a) $E_2^{\left(2\right)}$ is an unstable node if $r_1-m_1>\frac{q{y}_2}{a}$ ; (b) $E_2^{\left(2\right)}$ is a saddle point if $r_1-m_1<\frac{q{y}_2}{a}$ ; (c) $E_2^{\left(2\right)}$ is an unstable higher order singular point if $r_1-m_1=\frac{q{y}_2}{a}$. In the special case that two axial equilibria $E_2^{\left(k\right)}\left(k=\mathrm{1,2}\right)$ collide with each other, the Jacobian matrix is

$J\left(E_2\right)=\left[\begin{array}{cc}{r}_1-m_1-\frac{q{y}_2}{a}& 0\\ \frac{c{y}_2}{a}& 0\end{array}\right]\mathrm{,}$ (18)

thus $E_2$ is a higher order singular point. If $r_1-m_1>\frac{q{y}_2}{a}$, it is unstable.

For the interior equilibrium $E_{\mathrm{<mo>\; *\; </mo>}}$, its Jacobian matrix is

$J\left(E_{\mathrm{<mo>\; *\; </mo>}}\right)=\left[\begin{array}{cc}-\frac{r_1{x}_{\mathrm{<mo>\; *\; </mo>}}}{K_1}-d{x}_{\mathrm{<mo>\; *\; </mo>}}+\frac{q{x}_{\mathrm{<mo>\; *\; </mo>}}{y}_{\mathrm{<mo>\; *\; </mo>}}}{{\left(a+x_{\mathrm{<mo>\; *\; </mo>}}\right)}^2}& \frac{-q{x}_{\mathrm{<mo>\; *\; </mo>}}}{a+x_{\mathrm{<mo>\; *\; </mo>}}}\\ \frac{ca{y}_{\mathrm{<mo>\; *\; </mo>}}}{{\left(a+x_{\mathrm{<mo>\; *\; </mo>}}\right)}^2}& \frac{-r_2{y}_{\mathrm{<mo>\; *\; </mo>}}}{K_2{\left(y_{\mathrm{<mo>\; *\; </mo>}}+e\right)}^2}\left(y_{\mathrm{<mo>\; *\; </mo>}}^2+2e{y}_{\mathrm{<mo>\; *\; </mo>}}-K_{2}e\right)\end{array}\right]\mathrm{.}$ (19)

Denoting a new discriminant ${\Delta }_{\mathrm{<mo>\; *\; </mo>}}=A_1^2-4A_2$ with $A_1\mathrm{<mo>\; :\; </mo>}=tr\left[J\left(E_{\mathrm{<mo>\; *\; </mo>}}\right)\right]$ and $A_2:=\det \left[J\left(E_{*}\right)\right]$, we have:

(a) If $A_1<0$ and

(a1) $A_2>0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is an asymptotically stable node; (a2) $A_2<0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is a saddle point ;

(b) If $A_1=0$ and

(b1) $A_2>0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is a center or a focus; (b2) $A_2<0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is a saddle point;

(c) If $A_1>0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is unstable and

(c1) ${\Delta }_{*}=0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is a node; (c2) ${\Delta }_{*}<0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is a focus; (c3) ${\Delta }_{*}>0$ and $A_2>0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is a node; (c4) ${\Delta }_{*}>0$ and $A_2<0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is a saddle point.

Due to the further consideration of stability at the point $E_{\mathrm{<mo>\; *\; </mo>}}$, we will give following theorem to explain its global stability.

Theorem 3 (Global asymptotic stability) If the interior equilibrium $E_{\mathrm{<mo>\; *\; </mo>}}$ exist and parameters satisfy

$\frac{q{y}_{*}}{a\left(a+x_{*}\right)}<\frac{r_1}{K_1}+d,y_{*}>K_2,$ (20)

then $E_{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotically stable.

Proof of theorem 3. Here we take an unbounded positive definite Lyapunov function

$V=x-x_{\mathrm{<mo>\; *\; </mo>}}-x_{\mathrm{<mo>\; *\; </mo>}}\ln \frac{x}{x_{\mathrm{<mo>\; *\; </mo>}}}+A\left(y-y_{\mathrm{<mo>\; *\; </mo>}}-y_{\mathrm{<mo>\; *\; </mo>}}\ln \frac{y}{y_{\mathrm{<mo>\; *\; </mo>}}}\right)$ (21)

with $A=\frac{q\left(x_{\mathrm{<mo>\; *\; </mo>}}+a\right)}{ac}$. Introducing new variables $\hat{x}=x-x_{*}$, $\hat{y}=y-y_{*}$ and computing derivative along orbits of the system (1), we have time derivative ${\frac{\text{d}V}{\text{d}t}|}_{\left(1\right)}=-\left[\frac{r_1}{K_1}+d-\frac{q{y}_{*}}{\left(a+x\right)\left(a+x_{*}\right)}\right]{\hat{x}}^2+\frac{A{r}_2}{\left(y+e\right)\left(y_{*}+e\right)}\left[e-\frac{y{y}_{*}+e\left(y+y_{*}\right)}{K_2}\right]{\hat{y}}^2.$

It is obvious that ${\frac{\text{d}V}{\text{d}t}|}_{\left(1\right)}$ is negative definite. Consequently, the Lyapunov function V satisfies the asymptotic stability theorem in [25]. Thus we complete the proof. □

Remark. If $\frac{q{y}_{\mathrm{<mo>\; *\; </mo>}}}{a\left(a+x_{\mathrm{<mo>\; *\; </mo>}}\right)}\le \frac{r_1}{K_1}+d$ and $y_{\mathrm{<mo>\; *\; </mo>}}\ge K_2$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is asymptotically stable since $A_1<0$ and $A_2>0$.

3.3. Closed Orbits and Limit Cycles

In this subsection, we consider non-existence of closed orbits and limit cycles of system (1). Firstly, taking a diffeomorphism $\phi \mathrm{<mo>\; :\; </mo>}u=x$, $v=y$, $\text{d}\tau =\frac{\text{d}t}{\left(a+x\right)\left(y+e\right)}$ which preserving the orientation of time, the system (1) is topologically equivalent to following system [14] [26] [27]

$\stackrel{\dot{}}{x}=P\left(x,y\right):=\left[r_{1}x\left(1-\frac{x}{K_1}\right)-m_{1}x-d{x}^2\right]\left(a+x\right)\left(y+e\right)-qxy\left(y+e\right),$ (22a)

$\stackrel{\dot{}}{y}=Q\left(x,y\right):=r_2{y}^2\left(1-\frac{y}{K_2}\right)\left(a+x\right)+cxy\left(y+e\right)-m_{2}y\left(a+x\right)\left(y+e\right).$ (22b)

Notice that we still denote u, v and $\tau$ as x, y and t. Above system is a $C^{\infty }$ -qualitatively equivalent polynomial extension of the system (1) to $R_{+}^2$ and more convenient to study limit cycles throughout [13].

Theorem 4 (Non-existence of limit cycles) If $ad>c$ , $m_1>r_1$ and $m_2>r_2$ , then for system (22), there are no closed orbits and limit cycles in $R^{+2}$.

Proof of theorem 4. Here we take a Dulac function $B\left(x\mathrm{,}y\right)=\frac{1}{xy}$ and calculate following partial derivative:

$\frac{\partial \left(BP\right)}{\partial x}+\frac{\partial \left(BQ\right)}{\partial y}=\frac{1}{K_1{K}_{2}xy}{\displaystyle \underset{1\le i+j\le 3}{\sum }}{a}_{ij}{x}^i{y}^j\mathrm{,}$

where coefficients

$\begin{array}{l}{a}_{11}=-\left[\left(ad-c+m_1+m_2-r_1-r_2\right)K_1+a{r}_1\right]K_2,\\ a_{10}=-\left[\left(ad+m_1-r_1\right)K_1+a{r}_1\right]e{K}_2,\\ a_{01}=-K_1{K}_{2}a\left(m_2-r_2\right)\end{array}$

and other unlisted coefficients are all non-positive. One can also use the Bendixson criteria or Dulac function $B\left(x\mathrm{,}y\right)=\frac{1}{\left(a+x\right)\left(y+e\right)}$ to complete the proof. □

4. Bifurcations

In this section, we will consider and give sufficient conditions to show the existence of saddle-node bifurcation, transcritical bifurcation and Hopf bifurcation of the system (1). Firstly, we denote the system (1) as the following form for simplicity and convenient:

$\left(\begin{array}{c}\stackrel{\dot{}}{x}\\ \stackrel{\dot{}}{y}\end{array}\right)=\left(\begin{array}{c}f\\ g\end{array}\right):=F\left(x,y\right).$ (23)

4.1. Saddle-Node Bifurcation

As we see, when $r_2>m_2$ and $y_1=y_2$, or

$m_2=m_2^{\left[SN\right]}:=r_2+\frac{2r_{2}e}{K_2}\left(1-\sqrt{1+\frac{K_2}{e}}\right),$ (24)

the two equilibria $E_2^{\left(k\right)}\left(k=\mathrm{1,2}\right)$ collide with each other and the system (23) has an unique boundary equilibrium $E_2$ when $y>0$. Hence there is a chance of bifurcation around this unique equilibrium. Here we choose $m_2$ as bifurcation parameter and eigenvector $v=\left(\begin{array}{c}0\\ 1\end{array}\right)$ corresponding to the zero eigenvalue for the matrix (18). The eigenvector corresponding to the zero eigenvalue for the transpose of matrix (18) is $w=\left(\begin{array}{c}1\\ w_2\end{array}\right)$, where

$w_2=-\frac{r_1-m_1-\frac{q{y}_2}{a}}{\frac{c{y}_2}{a}}\ne 0.$

Suppose $r_1-m_1-\frac{q{y}_2}{a}<0$, then the following transversality conditions are hold:

$\begin{array}{l}{w}^{\text{T}}{F}_{m_2}\left(E_2,m_2^{\left[SN\right]}\right)=-w_2{y}_2\ne 0,\\ w^{\text{T}}\left[D^{2}F\left(E_2,m_2^{\left[SN\right]}\right)\left(v,v\right)\right]\frac{2w_2{r}_2}{K_2}\left(\sqrt{\frac{e}{K_2+e}}-1\right)\ne 0.\end{array}$ (25)

Thus we have following theorem by using the Sotomayor’s theorem [28] [29].

Theorem 5 (Saddle-node bifurcation) Suppose the point $E_2$ exist, if $r_1-m_1-\frac{q{y}_2}{a}<0$ , then the system (23) undergoes a saddle-node bifurcation around point $E_2$ with respect to the bifurcation parameter $m_2$.

4.2. Transcritical Bifurcation

Keeping in mind that there is an exchange of stability between the points $E_2^{\left(1\right)}$ and $E_2^{\left(2\right)}$ when $r_1$ crosses the threshold $r_1^{\left[TC\right]}\mathrm{<mo>\; :\; </mo>}=m_1+\frac{q{y}_1}{a}$. Suppose $J_{22}\left(E_2^{\left(1\right)}\right)<0$. Here we choose $r_1$ as bifurcation parameter and eigenvector $v=\left(\begin{array}{c}1\\ v_2\end{array}\right)$ corresponding to the zero eigenvalue for the Jacobian matrix $J\left(E_2^{\left(1\right)}\right)$ when $r_1=r_1^{\left[TC\right]}$, where $v_2=-\frac{c{y}_1}{a{J}_{22}\left(E_2^{\left(1\right)}\right)}>0$. The eigenvector corresponding to the zero eigenvalue for the transpose of matrix (18) is $w=\left(\begin{array}{c}1\\ 0\end{array}\right)$, then the following transversality conditions are

$\begin{array}{l}{w}^{\text{T}}{F}_{r_1}\left(E_2^{\left(1\right)}\mathrm{,}{r}_1^{\left[TC\right]}\right)=\mathrm{0,}{w}^{\text{T}}D{F}_{r_1}\left(E_2^{\left(1\right)}\mathrm{,}{r}_1^{\left[TC\right]}\right)v=\mathrm{1,}\\ w^{\text{T}}{D}^{2}F\left(E_2^{\left(1\right)}\mathrm{,}{r}_1^{\left[TC\right]}\right)\left(v\mathrm{,}v\right)=\frac{-2r_1^{\left[TC\right]}}{K_1}-2d+\frac{2q{y}_1}{a^2}+\frac{2qc{y}_1}{a^2{J}_{22}\left(E_2^{\left(1\right)}\right)}\ne 0.\end{array}$ (26)

Thus we have following theorem by using the Sotomayor’s theorem [28] [29].

Theorem 6 (Transcritical bifurcation) Suppose the two axial equilibria $E_2^{\left(k\right)}\left(k=\mathrm{1,2}\right)$ coexist and $J_{22}\left(E_2^{\left(1\right)}\right)<0$ , if $a\ge K_1$ or $J_{22}\left(E_2^{\left(1\right)}\right)\ge -c$ , then the system (23) undergoes a transcritical bifurcation around the point $E_2^{\left(1\right)}$ with respect to the bifurcation parameter $r_1$.

4.3. Hopf Bifurcation

Here the point $E_{\mathrm{<mo>\; *\; </mo>}}$ exist and we choose d as bifurcation parameter. Suppose that $\lambda \left(d\right)=\alpha \left(d\right)\pm i\omega \left(d\right)$ are a pair of conjugate eigenvalues of matrix $J\left(E_{\mathrm{<mo>\; *\; </mo>}}\right)$, where $\alpha \left(d\right)=A_1\left(d\right)/2$. The critical value $d^{\left[H\right]}$ satisfies

$A_1\left(d^{\left[H\right]}\right)=0,A_2\left(d^{\left[H\right]}\right)>0,{\alpha }^{\prime }\left(d^{\left[H\right]}\right)\ne 0.$ (27)

Then the system (23) undergoes a Hopf bifurcation around the point $E_{\mathrm{<mo>\; *\; </mo>}}$ with respect to the bifurcation parameter d.

We will calculate the first Lyapunov number $\sigma$ at the point $E_{\mathrm{<mo>\; *\; </mo>}}$, which is used to determine the stability of limit cycles. Therefore, translating the point $E_{\mathrm{<mo>\; *\; </mo>}}$ to the $E_0$ by a linear transformation in the proof of theorem 3, the system (23) in power series around the origin (drop the hats for the sake of convenience as usual) is

$\stackrel{\dot{}}{x}={\displaystyle \underset{1\le i+j\le 3}{\sum }}{a}_{ij}{x}^i{y}^j+F_1\left(x,y\right),\stackrel{\dot{}}{y}={\displaystyle \underset{1\le i+j\le 3}{\sum }}{b}_{ij}{x}^i{y}^j+F_2\left(x,y\right),$ (28)

where coefficients are

$\begin{array}{l}{a}_{10}=-\frac{r_1{s}_1}{K_1}-d{s}_1+\frac{q{s}_1{s}_2}{{\left(a+s_1\right)}^2},a_{01}=-\frac{q{s}_1}{a+s_1},a_{20}=-\frac{r_1}{K_1}+\frac{q{s}_{2}a}{{\left(a+s_1\right)}^3}-d,\\ a_{02}=0,a_{11}=-\frac{qa}{{\left(a+s_1\right)}^2},a_{30}=-\frac{q{s}_{2}a}{{\left(a+s_1\right)}^4},a_{03}=0,a_{21}=\frac{qa}{{\left(a+s_1\right)}^3},a_{12}=0;\\ b_{10}=\frac{c{s}_{2}a}{{\left(a+s_1\right)}^2},b_{01}=-\frac{r_2{s}_2}{K_2{\left(s_2+e\right)}^2}\left(s_2^2+2e{s}_2-K_{2}e\right),b_{20}=-\frac{c{s}_{2}a}{{\left(a+s_1\right)}^3},\\ b_{02}=\frac{\left[\left(K_2-3s_2\right)e^2-3s_2^{2}e-s_2^3\right]r_2}{K_2{\left(s_2+e\right)}^3},b_{11}=\frac{ca}{{\left(a+s_1\right)}^2},\\ b_{30}=\frac{c{s}_{2}a}{{\left(a+s_1\right)}^4},b_{03}=-\frac{r_2{e}^2\left(K_2+e\right)}{K_2{\left(s_2+e\right)}^4},b_{21}=-\frac{ca}{{\left(a+s_1\right)}^3},b_{12}=0,\end{array}$

and $F_1\left(x\mathrm{,}y\right)\mathrm{,}{F}_2\left(x\mathrm{,}y\right)=O\left({\left|x\mathrm{,}y\right|}^4\right)$ are smooth functions. From references [28] [30], the first Lyapunov number for a planar system is given by

$\begin{array}{c}\sigma =-\frac{3\pi }{2a_{01}{A}_2^{3/2}}\{[a_{10}{b}_{10}\left(a_{11}^2+a_{11}{b}_{02}+a_{02}{b}_{11}\right)+a_{10}{a}_{01}\left(b_{11}^2+a_{20}{b}_{11}+a_{11}{b}_{02}\right)\\ +b_{10}^2\left(a_{11}{a}_{02}+2a_{02}{b}_{02}\right)-2a_{10}{b}_{10}\left(b_{02}^2-a_{20}{a}_{02}\right)-2a_{10}{a}_{01}\left(a_{20}^2-b_{20}{b}_{02}\right)\\ -a_{01}^2\left(2a_{20}{b}_{20}+b_{11}{b}_{20}\right)+\left(a_{01}{b}_{10}-2a_{10}^2\right)\left(b_{11}{b}_{02}-a_{11}{a}_{20}\right)]\\ -\left(a_{10}^2+a_{01}{b}_{10}\right)\left[3\left(b_{10}{b}_{03}-a_{01}{a}_{30}\right)+2a_{10}\left(a_{21}+b_{12}\right)+\left(b_{10}{a}_{12}-a_{01}{b}_{21}\right)\right]\}.\end{array}$ (29)

We have to give following numerical simulations by computing the first Lyapunov number at the point $E_{\mathrm{<mo>\; *\; </mo>}}$, since above expression is much too complicated.

Theorem 7 (Hopf bifurcation) Assume that the equilibrium $E_{\mathrm{<mo>\; *\; </mo>}}$ exist and parameters satisfy conditions (27), then the system (23) undergoes a Hopf bifurcation around the equilibrium $E_{\mathrm{<mo>\; *\; </mo>}}$ as parameter d passes through the critical value $d^{\left[H\right]}$. The Hopf bifurcation is supercritical and limit cycles are stable if $\sigma <0$ ; the Hopf bifurcation is subcritical and limit cycles are unstable if $\sigma >0$.

Example 4.1 Combining the theorem 4, we consider the case 1 in subsection 3.1 again. To investigate how the control parameter d affects the dynamical behavior of our system (1), Figure 1 and Figure 2 depict the time series diagrams and phase diagrams corresponding to values $d=6.539$ and $d=6.48$ , respectively. This implies that the Hopf bifurcation occurs in the system (1). The first Lyapunov number $\sigma <0$ , thus the Hopf bifurcation is supercritical and a limit cycle generated by the critical point is stable.

Example 4.2 We consider the case 6 in subsection 3.1 again and notice the critical value $d^{\left[H\right]}$. Figure 3 and Figure 4 depict the time series diagrams and phase diagrams corresponding to values $d=0.1775$ and $d=0.175$ , respectively. The first Lyapunov number is also found to be negative.

Example 4.3 We consider the case 7 in subsection 3.1 again and notice the critical value $d^{\left[H\right]}$. Figure 5 and Figure 6 depict the time series diagrams and phase diagrams corresponding to values $d=0.125$ and $d=0.1215$ , respectively. The first Lyapunov number is also negative.

5. Summary and Remarks

In summary, we consider a predator-prey system with Holling type II functional response and Allee effect on the specific growth term of predator with multiply form. In Subsection 3.1, some cases about the existence of interior equilibria $E_{\mathrm{<mo>\; *\; </mo>}}$ are derived with the help of the cobweb model but we neglect the monotonicity

(a) (b) (c) (d)

Figure 1. Figures in the Example 4.1 with d = 6.539. (a) Time series diagram of population x; (b) Time series diagram of population y; (c) Phase diagram of populations x and y; (d) Enlarged phase diagram of populations x and y.

(a) (b) (c) (d)

Figure 2. Figures in the Example 4.1 with d = 6.48. (a) Time series diagram of population x; (b) Time series diagram of population y; (c) Phase diagram of populations x and y; (d) Enlarged phase diagram of populations x and y.

(a) (b) (c) (d)

Figure 3. Figures in the Example 4.2 with d = 0.1775. (a) Time series diagram of the prey x; (b) Time series diagram of the predator y; (c) Phase diagram of two populations x and y; (d) Enlarged phase diagram of two populations x and y.

(a) (b) (c) (d)

Figure 4. Figures in the Example 4.2 with d = 0.175. (a) Time series diagram of the prey x; (b) Time series diagram of the predator y; (c) Phase diagram of two populations x and y; (d) Enlarged phase diagram of two populations x and y.

(a) (b) (c) (d)

Figure 5. Figures in the Example 4.3 with d = 0.125. (a) Time series for population x; (b) Time series for population y; (c) shows phase diagram of two populations x and y; (d) shows enlarged phase diagram of two populations x and y.

(a) (b) (c) (d)

Figure 6. Figures in the Example 4.3 with d = 0.1215. (a) Time series for population x; (b) Time series for population y; (c) shows phase diagram of two populations x and y; (d) shows enlarged phase diagram of two populations x and y.

of functions sometimes, and we also extend the zero-point theorem to the plane $R_{+}^2$ ; for instance, when $y_k$ do not exist, a point from the Equation (6a) should belongs to the domain $\Sigma$ which is surrounded by the isocline in the Equation (6b) and two axes, while another point from the Equation (6a) should belong to ${\Sigma }^c$. It is our expectancy that this method can also be available in other complicated predator-prey systems in the form of ordinary/partial differential equations.

In Subsection 3.2, we give conclusions of stability when equilibria are hyperbolic, but omit critical cases of non-hyperbolic equilibria. These need to be researched further by using qualitative analysis and topologically equivalent systems. The saddle-node, transcritical and Hopf bifurcations in Section 4 correspond to the three critical cases: (a) $y_1=y_2$ or $\det \left[J\left(E_2\right)\right]=0$, $tr\left[J\left(E_2\right)\right]<0$ ; (b) $y_1\ne y_2$, $\det \left[J\left(E_2^{\left(1\right)}\right)\right]=0$, $tr\left[J\left(E_2^{\left(1\right)}\right)\right]<0$ ; (c) $\det \left[J\left(E_{\mathrm{<mo>\; *\; </mo>}}\right)\right]>0$, $tr\left[J\left(E_{\mathrm{<mo>\; *\; </mo>}}\right)\right]=0$, respectively. In Subsection 3.3, we obtained a theorem of non-existence conditions of limit cycles, while the existence, uniqueness and number of limit cycles need to be considered further. This theorem is suitable in the consideration of existence of limit cycles and Hopf bifurcation which generates a stable (unstable) limit cycle (see Example 4.1). For the corresponding homogeneous reaction-diffusion system subject to the Neumann boundary conditions:

$u_t=D_1{u}_{xx}+r_{1}u\left(1-\frac{u}{K_1}\right)-\frac{quv}{a+u}-m_{1}u-d{u}^2,t>0,x\in \Omega ,$ (30a)

$v_t=D_2{v}_{xx}+r_{2}v\left(1-\frac{v}{K_2}\right)\frac{v}{v+e}+\frac{cuv}{a+u}-m_{2}v,t>0,x\in \Omega ,$ (30b)

${\partial }_{\nu }u=0,{\partial }_{\nu }v=0,t\le 0,x\in \partial \left(\Omega \right),$ (30c)

$u\left(0,x\right)=u_0\left(x\right)\ge 0,v\left(0,x\right)=v_0\left(x\right)\ge 0,x\in \Omega .$ (30d)

where $\Omega =\left(0,l\pi \right)\left(l>0\right)$ is an one-dimensional domain; $D_1$ and $D_2$ are two positive diffusive constants; the interior equilibrium $E_{\mathrm{<mo>\; *\; </mo>}}=\left(u_{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}_{\mathrm{<mo>\; *\; </mo>}}\right)$ is asymptotically stable and Turing instability will not occur under the conditions (or the remark) of the theorem 3; thus we may get potential Hopf bifurcation points and transversality condition $D_n\left({\delta }_j^H\right)\ne 0$ may holds as we choose $u_{\mathrm{<mo>\; *\; </mo>}}=\delta$ as the Hopf bifurcation parameter [31] [32]. Moreover, for the corresponding homogeneous diffusion system (without terms $u_t$ and $v_t$ in above system) subject to the Neumann boundary conditions:

$-D_1{u}_{xx}=r_{1}u\left(1-\frac{u}{K_1}\right)-\frac{quv}{a+u}-m_{1}u-d{u}^2,t>0,x\in \Omega ,$ (31a)

$-D_2{v}_{xx}=r_{2}v\left(1-\frac{v}{K_2}\right)\frac{v}{v+e}+\frac{cuv}{a+u}-m_{2}v,t>0,x\in \Omega ,$ (31b)

${\partial }_{\nu }u=0,{\partial }_{\nu }v=0,t\ge 0,x\in \partial \left(\Omega \right),$ (31c)

$u\left(0,x\right)=u_0\left(x\right)\ge 0,v\left(0,x\right)=v_0\left(x\right)\ge 0,x\in \Omega .$ (31d)

Subsection 3.2 about local stability conclusions is also useful for this diffusive system (31) in local stability, steady state bifurcation and Hopf bifurcation analysis [33]. Finally, in Subsection 4.3, the bifurcation parameter d reveals Hopf bifurcation phenomenon as it crosses $d^{\left[H\right]}$ by using numerical simulations corresponding to the existence of $y_k$, respectively. The Hopf bifurcations are supercritical and limit cycles generated by the critical points are stable.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant No.31570364) and the National Key Research and Development Program of China (Grant No.2018YFE0103700). We thank the Editor and the referee(s) for their works.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Holling, C.S. (1965) The Functional Response of Predator Density and Its Role in Mimicry and Population Regulations. Memoirs of the Entomological Society of Canada, 45, 3-60. https://doi.org/10.4039/entm9745fv
[2]	Pei, Y.Z., Chen, L.S., Zhang, Q.R. and Li, C.G. (2005) Extinction and Performance of One-Prey Multi-Predators of Holling Type II Function Response System with Impulsive Biological Control. Journal of Theoretical Biology, 235, 495-503. https://doi.org/10.1016/j.jtbi.2005.02.003
[3]	Misha, P. and Raw, S.N. (2019) Dynamical Complexities in a Predator-Prey System Involving Teams of Two Prey and One Predator. Journal of Applied Mathematics and Computing, 61, 1-24. https://doi.org/10.1007/s12190-018-01236-9
[4]	Huang, J.C., Ruan, S.G. and Song, J. (2014) Bifurcation in a Predator-Prey System of Leslie Type with Generalized Holling Type III Functional Response. Journal of Differential Equations, 257, 1721-1752. https://doi.org/10.1007/s12190-018-01236-9
[5]	Allee, W.C. (1927) Animal Aggregations. The Quarterly Review of Biology, 2, 367-398. https://doi.org/10.1086/394281
[6]	Allee, W.C. (1931) Animal Aggregations: A Study in General Sociology. University of Chicago Press, Chicago. https://doi.org/10.5962/bhl.title.7313
[7]	Drake, J.M. and Kramer, A.M. (2011) Allee Effects. Nature Education Knowledge, 3, 2. https://www.nature.com/scitable/knowledge/library/allee-effects-19699394/
[8]	Dos Santos, L.S., Cabella, B.C.T. and Martinez, A.S. (2014) Generalized Allee Effect Model. Theory in Biosciences, 133, 117-124. https://doi.org/10.1007/s12064-014-0199-6
[9]	Stephens, P.A. and Sutherland, W.J. (1999) Consequences of the Allee Effect for Behaviour, Ecology and Conservation. Trends in Ecology & Evolution, 14, 401-405. https://doi.org/10.1016/S0169-5347(99)01684-5
[10]	Stephens, P.A., Sutherland, W.J. and Freckleton, R.P. (1999) What Is the Allee Effect? Oikos, 87, 185-190. https://doi.org/10.2307/3547011
[11]	Bascompte, J. (2003) Extinction Thresholds: Insights from Simple Models. Annales Zoologici Fennici, 40, 99-114.
[12]	Wang, M.H. and Kot, M. (2001) Speeds of Invasion in a Model with Strong or Weak Allee Effects. Mathematical Biosciences, 171, 83-97. https://doi.org/10.1016/S0025-5564(01)00048-7
[13]	Aguirre, P., Flores, J.D. and Gonzalez-Olivares, E. (2014) Bifurcations and Global Dynamics in a Predator-Prey Model with a Strong Allee Effect on the Prey, and a Ratio-Dependent Functional Response. Nonlinear Analysis: Real World Applications, 16, 235-249. https://doi.org/10.1016/j.nonrwa.2013.10.002
[14]	González-Olivares, E. and Rojas-Palma, A. (2012) Limit Cycles in a Gause-type Predator-Prey Model with Sigmoid Functional Response and Weak Allee Effect on Prey. Mathematical Methods in the Applied Sciences, 35, 963-975. https://doi.org/10.1002/mma.2509
[15]	Zu, J. (2013) Global Qualitative Analysis of a Predator-Prey System with Allee Effect on the Prey Species. Mathematics and Computers in Simulation, 94, 33-54. https://doi.org/10.1016/j.matcom.2013.05.009
[16]	Zu, J. and Mimura, M. (2010) The Impact of Allee Effect on a Predator-Prey System with Holling Type II Functional Response. Applied Mathematics and Computation, 217, 3542-3556. https://doi.org/10.1016/j.amc.2010.09.029
[17]	Cai, Y.L., Zhao, C.D., Wang, W.M. and Wang, J.F. (2015) Dynamics of a Leslie-Gower Predator-prey Model with Additive Allee Effect. Applied Mathematical Modelling, 39, 2092-2106. https://doi.org/10.1016/j.apm.2014.09.038
[18]	Zhou, S.R. and Wang, G. (2004) Allee-Like Effects in Metapopulation Dynamics. Mathematical Biosciences, 189, 103-113. https://doi.org/10.1016/j.mbs.2003.06.001
[19]	Xiao, Z.W., Xie, X.D. and Xue, Y.L. (2018) Stability and Bifurcation in a Holling Type II Predator-Prey Model with Allee Effect and Time Delay. Advances in Difference Equations, Article No.: 288. https://doi.org/10.1186/s13662-018-1742-4
[20]	Sun, Y.S., Wang, K.H. and Gui, Z.J. (2017) Periodic Solution of a Nonautonomous Predator-Prey System of Holling Type II with Strong Allee Effect and Impulsive Perturbation. 6th International Conference on Energy, Environment and Sustainable Development, DEStech Transactions on Environment, Energy and Earth Sciences, Phuket, Thailand, April 2017, 159-164. https://doi.org/10.12783/dteees/eesd2017/11992
[21]	Cui, R.H., Shi, J.P. and Wu, B.Y. (2014) Strong Allee Effect in a Diffusive Predator-Prey System with a Protection Zone. Journal of Differential Equations, 256, 108-129. https://doi.org/10.1016/j.jde.2013.08.015
[22]	Wang, G., Liang, X.G. and Wang, F.Z. (1999) The Competitive Dynamics of Populations Subject to an Allee Effect. Ecological Modelling, 124, 183-192. https://doi.org/10.1016/S0304-3800(99)00160-X
[23]	Birkhoff, G. and Rota, G.C. (1962) Ordinary Differential Equations Introductions to Higher Mathematics. Ginn and Company, Boston.
[24]	Chen, F. (2005) On a Nonlinear Nonautonomous Predator-Prey Model with Diffusion and Distributed Delay. Journal of Computational and Applied Mathematics, 180, 33-49. https://doi.org/10.1016/j.cam.2004.10.001
[25]	Merkin, D.R. Afagh, F.F. and Smirnov, A.L. (1997) Introduction to the Theory of Stability. Springer, New York.
[26]	Chicone, C. (2006) Ordinary Differential Equations with Applications. World Scientific, Springer-Verlag, New York.
[27]	Andronov, A. (1973) Qualitative Theory of Second-Order Dynamic Systems. Halsted Press.
[28]	Perko, L. (2001) Differential Equations and Dynamical Systems. 3rd Edition, Springer-Verlag, New York. https://doi.org/10.1007/978-1-4613-0003-8
[29]	Pirayesh, B., Pazirandeh, A. and Akbari, M. (2016) Local Bifurcation Analysis in Nuclear Reactor Dynamics by Sotomayor’s Theorem. Annals of Nuclear Energy, 94, 716-731. https://doi.org/10.1016/j.anucene.2016.04.021
[30]	Kuznetsov, Y. (1998) Elements of Applied Bifurcation Theory. 2nd Edition, Springer-Verlag, New York.
[31]	Yi, F.Q., Wei, J.J. and Shi, J.P. (2009) Bifurcation and Spatiotemporal Patterns in a Homogeneous Diffusive Predator-Prey System. Journal of Differential Equations, 246, 1944-1977. https://doi.org/10.1016/j.jde.2008.10.024
[32]	Wan, A.Y., Song, Z.Q., Zheng, L.F. (2016) Patterned Solutions of a Homogenous Diffusive Predator-Prey System of Holling Type-III. Acta Mathematicae Applicatae Sinica (English Series), 32, 1073-1086. https://doi.org/10.1007/s10255-016-0628-z
[33]	Casten, R.G. and Holland, C.J. (1977) Stability Properties of Solutions to Systems of Reaction-Diffusion Equations. SIAM Journal on Applied Mathematics, 33, 353-364. https://doi.org/10.1137/0133023