Dynamics of a Prey-predator System with Foraging Facilitation in Predators

TL;DR

This paper analyzes a prey-predator model with predator foraging facilitation, using advanced algebraic methods to determine equilibrium conditions and bifurcations, supported by numerical simulations.

Contribution

It introduces a novel analytical approach combining pseudo-division and polynomial discrimination to study complex bifurcations in predator-prey dynamics.

Findings

Identified parameter conditions for equilibria

Classified bifurcation types including Hopf and Bogdanov-Takens

Validated results with numerical simulations

Abstract

The dynamics of a prey-predator system with foraging facilitation among predators are investigated. The analysis involves the computation of many semi-algebraic systems of large degrees. We apply the pseudo-division reduction, real-root isolation technique and complete discrimination system of polynomial to obtain parameter conditions for the exact number of equilibria and their qualitative properties as well as a complete investigation of bifurcations including saddle-node, transcritical, pitchfork, Hopf and Bogdanov-Takens bifurcations. Moreover, numerical simulations are presented to support our theoretical results.

Equations

$$\displaystyle\begin{array}[]{l}\frac{dN}{dt}=rN(1-\frac{N}{K})-\frac{ENP}{1+HEN},~{}~{}~{}\frac{dP}{dt}=e\frac{ENP}{1+HEN}-mP,\end{array}$$

$$\displaystyle\begin{array}[]{l}\frac{dN}{dt}=rN(1-\frac{N}{K})-\frac{E(P)NP}{1+HE(P)N},~{}~{}\frac{dP}{dt}=e\frac{E(P)NP}{1+HE(P)N}-mP\end{array}$$

$$\displaystyle\begin{array}[]{l}\frac{dx}{d\tau}=x\{\sigma(1-\frac{x}{k})-\frac{(1+\alpha y)y}{1+h(1+\alpha y)x}\},~{}~{}~{}\frac{dy}{d\tau}=y\{\frac{(1+\alpha y)x}{1+h(1+\alpha y)x}-1\},\end{array}$$

$$\displaystyle\begin{array}[]{l}\frac{dx}{dt}=x\{\sigma(k-x)(1+h(1+\alpha y)x)-ky(1+\alpha y)\},~{}~{}\frac{dy}{dt}=ky\{x(1+\alpha y)(1-h)-1\}.\end{array}$$

$$\displaystyle\left.\begin{array}[]{l}\mathcal{P}_{1}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k>k_{1},\alpha<\frac{1}{\sigma k}\},~{}\mathcal{S}_{1}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k=k_{1},\alpha<\frac{1}{\sigma k}\},\\ \mathcal{P}_{2}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k<k_{1},\alpha<\frac{1}{\sigma k}\},~{}\mathcal{S}_{2}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k>k_{1},\alpha=\frac{1}{\sigma k}\},\\ \mathcal{L}_{1}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k=k_{1},\alpha=\frac{1}{\sigma k}\},~{}\mathcal{S}_{3}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k<k_{1},\alpha=\frac{1}{\sigma k}\},\\ \mathcal{P}_{3}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k>k_{1},\alpha>\frac{1}{\sigma k}\},~{}\mathcal{S}_{4}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k=k_{1},\alpha>\frac{1}{\sigma k}\},\\ \mathcal{S}_{5}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k<k_{1},\alpha=\alpha_{1}\},~{}\mathcal{P}_{5}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k<k_{1},\alpha>\alpha_{1}\},\\ \mathcal{P}_{4}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k<k_{1},\frac{1}{\sigma k}<\alpha<\alpha_{1}\}\\ \end{array}\right.$$

$$\displaystyle\left.\begin{array}[]{l}\alpha_{1}:=\frac{-(h-1)^{2}k^{2}+18(h-1)k+27+(hk-k+9)\sqrt{(hk-k+1)(hk-k+9)}}{8\sigma(1-h)k^{2}}.\end{array}\right.$$

$$\displaystyle\begin{array}[]{l}\mathcal{P}_{11}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k\geq k_{2},\alpha<\frac{1}{\sigma k}~{}\mbox{or}~{}k_{3}\leq k<k_{2},\alpha_{2}<\alpha<\frac{1}{\sigma k}~{}\mbox{or}~{}k_{1}<k<k_{3},\sigma<\sigma_{1},\alpha_{2}<\alpha<\frac{1}{\sigma k}\},\\ \mathcal{P}_{12}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k_{3}\leq k<k_{2},\alpha<\alpha_{2}~{}\mbox{or}~{}k_{1}<k<k_{3},\sigma\geq\sigma_{1},\alpha<\frac{1}{\sigma k}~{}\mbox{or}~{}k_{1}<k<k_{3},\sigma<\sigma_{1},\\ \phantom{\mathcal{P}_{12}:=}\alpha<\alpha_{2}\},~{}\mathcal{S}_{11}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k_{3}\leq k<k_{2},\alpha=\alpha_{2}~{}\mbox{or}~{}k_{1}<k<k_{3},\sigma<\sigma_{1},\alpha=\alpha_{2}\},\\ \mathcal{S}_{21}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k\geq k_{3},\alpha=\frac{1}{\sigma k}~{}\mbox{or}~{}k_{1}<k<k_{3},\sigma<\sigma_{1},\alpha=\frac{1}{\sigma k}\},\\ \mathcal{S}_{22}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k_{1}<k<k_{3},\sigma>\sigma_{1},\alpha=\frac{1}{\sigma k}\},\mathcal{L}_{21}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k_{1}<k<k_{3},\sigma=\sigma_{1},\alpha=\frac{1}{\sigma k}\},\\ \mathcal{P}_{31}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k\geq k_{3},\alpha>\frac{1}{\sigma k}~{}\mbox{or}~{}k_{1}<k<k_{3},\sigma\geq\sigma_{1},\alpha>\alpha_{2}~{}\mbox{or}~{}k_{1}<k<k_{3},\sigma<\sigma_{1},\alpha>\frac{1}{\sigma k}\},\\ \mathcal{P}_{32}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k_{1}<k<k_{3},\sigma>\sigma_{1},\frac{1}{\sigma k}<\alpha<\alpha_{2}\},\\ \mathcal{S}_{31}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k_{1}<k<k_{3},\sigma>\sigma_{1},\alpha=\alpha_{2}\},\mathcal{S}_{41}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k=k_{1},\frac{1}{\sigma k}<\alpha<\alpha_{2}\},\\ \mathcal{L}_{41}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k=k_{1},\alpha=\alpha_{2}\},~{}\mathcal{S}_{42}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k=k_{1},\alpha>\alpha_{2}\},\\ \mathcal{P}_{51}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k<k_{1},\sigma\leq\sigma_{2},\alpha>\alpha_{1}~{}\mbox{or}~{}k<k_{1},\sigma>\sigma_{2},\alpha>\alpha_{2}\},\\ \mathcal{P}_{52}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k<k_{1},\sigma>\sigma_{2},\alpha_{1}<\alpha<\alpha_{2}\},~{}\mathcal{S}_{51}:=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k<k_{1},\sigma>\sigma_{2},\alpha=\alpha_{2}\}\end{array}$$

$$\displaystyle\left.\begin{array}[]{l}\alpha_{2}:=\frac{\{kh(h-1)+h+1\}\{k(h-1)^{2}+h\sigma+\sigma\}^{2}}{k^{2}(1-h)(h\sigma+1-h)^{2}\{k(h-1)^{2}+h+\sigma-1\}},\\ \sigma_{2}:=\frac{(1-h)(h^{2}k-hk+3h+3)(hk-k+1)+(1-h)(h^{2}k-hk+h+1)\sqrt{(hk-k+9)(hk-k+1)}}{2\{h(-h^{2}+2h-1)k-h^{2}+h+2\}}\end{array}\right.$$

$$\displaystyle\left.\begin{array}[]{l}f(\sigma):=\{hk(h-1)(h^{2}+3h+1)k+(h+1)^{3}\}\sigma^{3}+\{k(h-1)^{2}(h(3h+2)(h-1)k+3h^{2}+2h+2)\}\sigma^{2}\\ \phantom{f(\sigma):=}+\{k(h-1)^{3}(kh-k+1)(h^{2}k-hk-2h+1)\}\sigma+k(h-1)^{4}(hk-k+1)\end{array}\right.$$

$$\displaystyle\begin{array}[]{l}x\{\sigma(k-x)(1+h(1+\alpha y)x)-ky(1+\alpha y)\}=0,~{}~{}ky\{x(1+\alpha y)(1-h)-1\}=0.\end{array}$$

$$\displaystyle\left.\begin{array}[]{l}y=\frac{\sigma x(k-x)}{k},~{}~{}0<x<k.\end{array}\right.$$

$$\displaystyle\left.\begin{array}[]{l}F(x):=\alpha\sigma(h-1)x^{3}+\alpha\sigma k(1-h)x^{2}+k(1-h)x-k,\end{array}\right.$$

$$\displaystyle\left.\begin{array}[]{l}x_{*}:=\frac{k\alpha\sigma+\sqrt{k\alpha\sigma(k\alpha\sigma+3)}}{3\alpha\sigma}.\end{array}\right.$$

$$prem(F(x),F^{\prime }(x),x,‘m\text{mbox}{}^{\prime })=k{\alpha }^2{\sigma }^2(1-h{)}^2\{2(k\alpha \sigma +3)(1-h)x+k(1-h)-9\},$$

$$\displaystyle\left.\begin{array}[]{l}\tilde{F}(x):=2(k\alpha\sigma+3)(1-h)x+k(1-h)-9.\end{array}\right.$$

$$\displaystyle\left.\begin{array}[]{l}\tilde{F}(x_{*})=\frac{2(k\alpha\sigma+3)(1-h)\sqrt{k\alpha\sigma(k\alpha\sigma+3)}+\alpha\sigma\{2k^{2}\sigma\alpha(1-h)+9k(1-h)-27\}}{3\alpha\sigma},\end{array}\right.$$

$$\displaystyle\left.\begin{array}[]{l}F_{1}(\alpha):=4k^{2}\sigma^{2}(1-h)\alpha^{2}+\sigma\{(1-h)^{2}k^{2}+18(1-h)k-27\}\alpha+4k(1-h)^{2},\end{array}\right.$$

$$\displaystyle\left.\begin{array}[]{l}x_{*}=\frac{9-k(1-h)}{2(k\alpha_{1}\sigma+3)(1-h)}\end{array}\right.$$

$$J:=\left(\begin{array}[]{lr}J_{11}&x\{(hx(k-x)\sigma-2ky)\alpha-k\}\\ -ky(h-1)(\alpha y+1)&-k\{x(h-1)(2\alpha y+1)+1\}\end{array}\right),$$

$$\displaystyle\left.\begin{array}[]{l}D|_{E_{i}}=\frac{\sigma(k-x_{i})\{\alpha\sigma x_{i}(k-x_{i})+k\}x_{i}^{2}}{k}F^{\prime}(x_{i}),~{}~{}T|_{E_{i}}=\frac{(-h^{2}k+2hk-h\sigma-k-\sigma)x_{i}+k(h\sigma-h+1)}{1-h}.\end{array}\right.$$

$$\displaystyle\left.\begin{array}[]{l}x_{0}:={\frac{k\left(h\sigma-h+1\right)}{k\left(h-1\right)^{2}+\sigma\,\left(h+1\right)}}.\end{array}\right.$$

$$\displaystyle\left.\begin{array}[]{l}F(x_{0})=\frac{k\sigma\{k^{2}(1-h)(h\sigma-h+1)^{2}(k(h-1)^{2}+h+\sigma-1)\alpha+(-h^{2}k+hk-h-1)(h^{2}k-2hk+h\sigma+k+\sigma)^{2}\}}{\{k(h-1)^{2}+\sigma(h+1)\}^{3}}.\end{array}\right.$$

$$\displaystyle\left.\begin{array}[]{l}f_{1}(\sigma):=\{hk(h-1)+4h^{2}+5h+1\}\sigma^{2}+(h-1)\{(h-1)^{2}hk^{2}+(6h^{2}-5h-1)k-3h-3\}\sigma\\ \phantom{f_{1}(\sigma):=}-4k(h-1)^{3},\\ f_{2}(\sigma):=\{-h(h-1)^{2}k-h^{2}+h+2\}\sigma^{2}+(h-1)(h^{2}k-hk+3h+3)(hk-k+1)\sigma\\ \phantom{f_{1}(\sigma):=}+k(h-1)^{3}(hk-k+1).\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{l}\frac{du}{d\tau}=\alpha(h-1)^{2}u^{3}+\alpha(h-1)^{2}u^{2}v+\frac{(2\alpha\sigma+h-1)(h-1)^{2}}{\sigma}u^{2}\epsilon-\frac{(\alpha\sigma+h-1)(h-1)}{\sigma}u^{2}+\frac{(h-1)^{3}}{\sigma}uv\epsilon\\ \phantom{\frac{du}{d\tau}=}-\frac{(h-1)^{2}}{\sigma}uv+\frac{(h-1)^{3}}{\sigma}u\epsilon^{2}-\frac{(h-1)^{2}}{\sigma}u\epsilon+O(\parallel(u,v,\epsilon)\parallel^{3}),\\ \frac{dv}{d\tau}=v-(h^{2}-1)v\epsilon-(h^{2}-1)v^{2}-\frac{(\sigma-1)(h-1)^{2}}{\sigma}u\epsilon-\frac{(h-1)\{(\sigma-1)(\alpha\sigma+h)+1\}}{\sigma}u^{2}\\ \phantom{\frac{du}{d\tau}=}-\frac{\alpha h\sigma^{2}+(2h^{2}-h-1)\sigma-(h-1)^{2}}{\sigma}uv+\{((2h-1)\sigma-h+1)\alpha+h(h-1)\}(h-1)u^{3}\\ \phantom{\frac{du}{d\tau}=}+\{((4h-1)\sigma-h+1)\alpha+3h(h-1)\}(h-1)u^{2}v+\frac{(\sigma-1)(h-1)^{3}}{\sigma}u\epsilon^{2}+h(h-1)^{2}v^{3}\\ \phantom{\frac{du}{d\tau}=}+\frac{\{2\sigma(\sigma-1)\alpha+(2h-1)\sigma-h+1\}(h-1)^{2}}{\sigma}u^{2}\epsilon+h\{2\alpha\sigma+3(h-1)\}(h-1)uv^{2}+2h(h-1)^{2}v^{2}\epsilon\\ \phantom{\frac{du}{d\tau}=}+\frac{\{2\alpha h\sigma^{2}+(4h^{2}-5h+1)\sigma-(h-1)^{2}\}(h-1)}{\sigma}uv\epsilon+h(h-1)^{2}v\epsilon^{2}+O(\parallel(u,v,\epsilon)\parallel^{3}),\\ \frac{d\epsilon}{d\tau}=0.\end{array}\right.$$

$$v=h(u,ϵ)=a{u}^2+buϵ+c{ϵ}^2+O(\parallel (u,ϵ){\parallel }^3).$$

$$\displaystyle\begin{array}[]{l}\frac{du}{dt}=-{\frac{\left(-1+h\right)^{2}\epsilon\,u}{\sigma}}-{\frac{\left(\alpha\,\sigma+h-1\right)\left(-1+h\right){u}^{2}}{\sigma}}+{\frac{\left(-1+h\right)^{3}u{\epsilon}^{2}}{\sigma}}+c_{1}{u}^{2}\epsilon+c_{2}{u}^{3}+O(\parallel(u,\epsilon)\parallel^{4}),\end{array}$$

$$\displaystyle\begin{array}[]{l}T|_{E_{*}}=\frac{(hk-k+9)}{8(1-h)^{2}(\alpha_{1}k\sigma+3)}\{(h\sigma-h+1)\sqrt{(kh-k+1)(hk-k+9)}\\ \phantom{f_{1}(\sigma):=}-(kh^{2}-hk+h+4)\sigma+3(1-h)(hk-k+1)\}.\end{array}$$

$$\displaystyle\begin{array}[]{l}f_{2}(\sigma)=-\frac{1}{8}\big{\{}\{(h\sigma-h+1)\sqrt{(kh-k+1)(hk-k+9)}\}^{2}\\ \phantom{f_{2}(\sigma)=}-\{-(kh^{2}-hk+h+4)\sigma+3(1-h)(hk-k+1)\}^{2}\big{\}}.\end{array}$$

$$\displaystyle\begin{array}[]{l}f_{2}(\frac{3(1-h)(hk-k+1)}{kh^{2}-hk+h+4})=-\frac{2(h-1)^{2}(hk-k+9)(hk-k+1)(h^{2}k-hk+h+1)^{2}}{(h^{2}k-hk+h+4)^{2}}<0,\end{array}$$

$$\displaystyle\left\{\begin{array}[]{l}\frac{dx}{dt}=a_{100}x+a_{010}y+a_{001}\epsilon+a_{200}x^{2}+a_{020}y^{2}+a_{110}xy+a_{011}\epsilon y+a_{101}\epsilon x+O(\parallel(x,y,\epsilon)\parallel^{3}),\\ \frac{dy}{dt}=b_{100}x+b_{010}y+b_{001}\epsilon+b_{110}xy+b_{101}x\epsilon+b_{020}y^{2}+b_{011}y\epsilon+O(\parallel(x,y,\epsilon)\parallel^{3}),\\ \displaystyle\frac{d\epsilon}{dt}=0,\end{array}\right.$$

Full text

{history}

1 Introduction

Populations rarely exist in isolation, which results in ecological systems are characterized by the interaction between species and environment. Mathematical models play important roles in understanding population interactions (Freedman [1980]; Kot [2001]). An important type of interaction is predation, which leads to prey-predator models that have great importance in ecology. One of the classic prey-predator models, the Rosenzweig-MacArthur model(Rosenzweig & Macarthur, [1963]), is given by

[TABLE]

where $N(t)$ and $P(t)$ represent densities of the prey and predator at time $t$ respectively, $r$ stands for the intrinsic growth rate of prey, $K$ is the carrying capacity of prey, $e$ is the conversion rate, $m$ is the mortality rate of predator, $E$ is the encounter rate of predator with the prey and $H$ is the predator handling time of a prey individual. Some researchers (Hsu & Waltman, [1978]; Cheng, [1981]; Huang, [1988]; Turchin [2013]; Kot [2001]) have studied the dynamical behaviors of system (2), which has a coexistence equilibrium rose from transcritical bifurcation and a unique limit cycle induced by Hopf bifurcation. They also have shown both the prey and predator populations survive either to the coexistence equilibrium or the limit cycle. Another widespread type of interaction in ecological systems is cooperation among individuals (Dugatkin [1997]), which seems to be an important evolutionary cause of sociality and a key factor for exploring and understanding many aspects of how organisms are designed. There are a great variety of cooperative behaviors in nature such as cooperative defence against predators (Garay, [2009]), cooperative breeding (Courchamp et al. [2008]), alarm calling (Lehmann & Keller, [2006]) and cooperative hunting (Boesch, [1994]; Packer & Ruttan, [1988]). The behaviour of cooperation during prey hunting has been observed in many different species, for instance, some species of tuna hunt in a linear school and aggregate when they encounter a school of prey (Partridge et al., [1983]) and wolves can hunt animals bigger or faster than themselves by cooperative hunting (Schmidt & Mech, [1997]). Foraging facilitation or hunting cooperation embraces a number of specific mechanisms such as locating and capturing the prey in a bigger group (Cosner et al., [1999]), protecting any of members from predation (Krause [2002]) and intraspecific cooperation (Courchamp & Macdonald, [2001]). Recently, the foraging facilitation has been taken into consideration in some mathematical literatures (Berec, [2010]; Cosner et al., [1999]; Kimun et al., [2018]; Pribylova & Peniaskova, [2017]; Alves & Hilker, [2017]; Saheb et al., [2018]). Foraging facilitation can be depicted in mathematical models by functional response, which means the per capita feeding rate of predators on their prey. The independence of the Holling type II functional response in system (2) from predator density is hardly always true in reality because it reflects that any single predator affects the growth rate of prey independently of its conspecifics. Therefore, functional response might depend on predator density and is increasing with respect to predator density for the case of foraging facilitation. That is to say, when any of the foraging facilitation mechanisms operates, $E$ in Holling type II functional response no longer is a constant, but rather an increasing function of predator density.

Berec (Berec, [2010]) extended the classical Rosenzweig-MacArthur system by including foraging facilitation and proposed the following prey-predator system

[TABLE]

with the encounter-driven functional response $E(P):=e_{1}/(e_{2}+P)^{\omega}$, where $e_{1}>0$, $e_{2}\geq 0$ and $\omega\leq 0$. Clearly, the above model is exactly the Rosenzweig-MacArthur model as $\omega=0$, and it characterizes the foraging facilitation as $\omega<0$. Berec gave a brief overview on the number and stabilities of coexistence equilibria of system (4), and later Pribylova and Peniaskova (Pribylova & Peniaskova, [2017]) considered the bifurcation behaviors through qualitative analysis combined with numerical simulations. In the special case $e_{2}=0$ and $\omega=-1$, the functional response happens to be the one considered by Cosner (Cosner et al., [1999]), which actually describes the foraging facilitation in a spatially linear formation and aggregation when the predators encounter a cluster of prey. Kimun et al (Kimun et al., [2018]) analyzed system (4) with the special functional response. Furthermore, Alves and Hilker (Alves & Hilker, [2017]) investigated both of the two special cases $\omega=-1$, $e_{2}>0$ with $H=0$ and $H>0$ respectively and derived the result that the hunting cooperation in the prey-predator system induces Allee effects in predators. In the case $\omega=-1$, $e_{2}>0$ and $H=0$, they investigated the stabilities of equilibria and saddle-node, Hopf and Bogdanov-Takens bifurcations. In the case $\omega=-1$, $e_{2}>0$ and $H>0$, by dimensionless transformations $x=\frac{ee_{1}e_{2}}{m}N$, $y=\frac{e_{1}e_{2}}{m}P$, $\tau=mt$, $\sigma=\frac{r}{m}$, $k=\frac{ee_{1}e_{2}K}{m}$, $\alpha=\frac{m}{e_{1}e_{2}^{2}}$ and $h=\frac{mH}{e}$ system (4) can be written into

[TABLE]

where $\alpha$ describes the intensity of predator cooperation in hunting. System (6) is a direct extension of the Rosenzweig-MacArthur model by considering the foraging facilitation. Alves and Hilker (Alves & Hilker, [2017]) presented a two-parameter bifurcation diagram of system (6) for special parameter values $h=0.1$ and $k=0.8$. Therefore, further carrying out a detailed study of system (6) is the task of this paper.

Note that system (6) is orbitally equivalent to the following quartic system

[TABLE]

In this paper, we investigate the dynamics of the above system with positive parameters $h$, $k$, $\sigma$ and $\alpha$ in the closure of the first quadrant $\overline{\mathbb{R}_{+}^{2}}:=\{(x,y)\in\mathbb{R}^{2}:x\geq 0,y\geq 0\}$. $\overline{\mathbb{R}_{+}^{2}}$ is positively invariant under the flow generated by system (8). In fact, the origin $(0,0)$ is an equilibrium, the positive $y$-axis is an orbital and the positive $x$-axis consists of three orbitals, i.e., $0<x<k$, $x>k$ and the equilibrium $(k,0)$. Notice that the abscissas of equilibria of the above system are decided by those positive roots of a cubic polynomial with complicated coefficients. However, generically we cannot obtain the analytic expressions of those equilibria. In Section 2, we qualitatively analyse the cubic polynomial equilibrium function and investigate the relative positions of those roots for the equilibrium function and the trace of the Jacobian matrix. Consequently, we obtain the parameter conditions for the exact number of equilibria and their qualitative properties. Section 3 is devoted to equilibria with exact one zero eigenvalue. Restricting on the center manifold, we obtain parameter conditions for transcritical, pitchfork and saddle-node bifurcations. In Section 4, we apply the pseudo-division reduction (Winkler [1996]) and real-root isolation technique to determine the sign of the first quantity of focus, which is a quartic polynomial with complex coefficients. It is proved that at most one limit cycle bifurcates via Hopf bifurcation. In Section 5, we investigate the Bogdanov-Takens bifurcation and show that it is codimension 2. Furthermore, the complete discrimination system of polynomial (Yang, [1999]) is applied to verify the transversal condition. In Section 6, we verify the results by numerical simulations and end the paper with a brief biological implications.

2 Equilibria and Their Properties

In order to state our results conveniently, we consider the partition $\mathbb{R}_{+}^{3}:=\{(k,\sigma,\alpha)\in\mathbb{R}^{3}:k>0,\sigma>0,\alpha>0\}=\mathcal{P}_{1}\cup\mathcal{S}_{1}\cup\mathcal{P}_{2}\cup\mathcal{S}_{2}\cup\mathcal{L}_{1}\cup\mathcal{S}_{3}\cup\mathcal{P}_{3}\cup\mathcal{S}_{4}\cup\mathcal{P}_{4}\cup\mathcal{S}_{5}\cup\mathcal{P}_{5}$, where

[TABLE]

with $k_{1}:=1/(1-h)$ and

[TABLE]

We further consider partitions $\mathcal{P}_{1}=\mathcal{P}_{11}\cup\mathcal{P}_{12}\cup\mathcal{S}_{11}$, $\mathcal{S}_{2}=\mathcal{S}_{21}\cup\mathcal{S}_{22}\cup\mathcal{L}_{21}$, $\mathcal{P}_{3}=\mathcal{P}_{31}\cup\mathcal{P}_{32}\cup\mathcal{S}_{31}$, $\mathcal{S}_{4}=\mathcal{S}_{41}\cup\mathcal{L}_{41}\cup\mathcal{S}_{42}$ and $\mathcal{P}_{5}=\mathcal{P}_{51}\cup\mathcal{P}_{52}\cup\mathcal{S}_{51}$, where

[TABLE]

with $k_{2}:=(1+h)/h(1-h)$, $k_{3}:=(h+1)^{3}/\{h(1-h)(h^{2}+3h+1)\}$,

[TABLE]

and $\sigma_{1}$ is the unique positive root of the following function

[TABLE]

as $0<h<1$ and $k_{1}<k<k_{3}$. The following theorem is devoted to the number of equilibria of system (8) and their qualitative properties.

Theorem 2.1**.**

System (8) has at most four equilibria. The exact number and qualitative properties of equilibria are described in Table 1.

Proof 2.2**.**

Equilibria of system (8) are determined by the algebraic equations

[TABLE]

For $y=0$, we can find two equilibria $E_{0}:(0,0)$ and $E_{k}:(k,0)$. For $y>0$, from the second equation in (20) system (8) has no other equilibrium if $h\geq 1$. If $h<1$, substituting equality $1+hx(1+\alpha y)=x(1+\alpha y)$ into the first equation in (20), we conclude that all equilibria lie on the curve

[TABLE]

Substituting (22) into the second equation in (20), we obtain

[TABLE]

whose zeros in the interval $(0,k)$ determine all equilibria of system (8). The derivative of $F(x)$ is $F^{\prime}(x)=(1-h)(-3\alpha\sigma x^{2}+2\alpha\sigma kx+k)$, which has a unique positive root

[TABLE]

It is easily seen that $F(x)$ is monotonically increasing for $0<x<x_{*}$ and monotonically decreasing for $x>x_{*}$. We need to discuss the zeros of $F(x)$ in the interval $(0,k)$ for $0<h<1$ in two cases: $x_{*}\geq k$ and $x_{*}<k$. (I). For the case $x_{*}\geq k$, i.e., $\alpha\leq\frac{1}{\sigma k}$, the discussion is divided into the following two subcases. (I.1) If $F(k)>0$, i.e., $k>k_{1}$, then $F(x)=0$ has a unique root in the interval $(0,k)$ denoted by $x_{1}$ (see Fig. 1 (a)). The corresponding parameters $(k,\sigma,\alpha)$ locate in $\mathcal{P}_{1}\cup\mathcal{S}_{2}$. (I.2) If $F(k)\leq 0$, i.e., $k\leq k_{1}$, then $F(x)=0$ has no root in the interval $(0,k)$. The corresponding parameters $(k,\sigma,\alpha)$ locate in $\mathcal{S}_{1}\cup\mathcal{P}_{2}\cup\mathcal{L}_{1}\cup\mathcal{S}_{3}$. (II). For the case $x_{*}<k$, i.e., $\alpha>\frac{1}{\sigma k}$, we need to discuss in the following two subcases. (II.1) If $F(k)\geq 0$, i.e., $k\geq k_{1}$, then $F(x)=0$ has a unique root in the interval $(0,k)$ denoted by $x_{1}$ (see Fig. 1 (a)). The corresponding parameters $(k,\sigma,\alpha)$ locate in $\mathcal{P}_{3}\cup\mathcal{S}_{4}$. (II.2) If $F(k)<0$, i.e., $k<k_{1}$, it should be clear that we need only account for the sign of $F(x_{*})$ to determine the number of zeros of $F(x)$. Since $F^{\prime}(x_{*})=0$, we can use Maple command “prem” to get the pseudo-remainder of $F(x)$ divided by $F^{\prime}(x)$ at $x_{*}$, i.e.,

[TABLE]

where $m=9\alpha^{2}\sigma^{2}(h-1)^{2}$. Thus, at $x_{*}$ we have $F(x)=\frac{k}{9}\tilde{F}(x)$ with

[TABLE]

Substituting $x_{*}$ given by (26) into $\tilde{F}(x)$ leads to

[TABLE]

in which the sign of $2(k\alpha\sigma+3)(1-h)\sqrt{k\alpha\sigma(k\alpha\sigma+3)}$ is positive but that of $2k^{2}\sigma\alpha(1-h)+9k(1-h)-27$ is indeterminate. If $2k^{2}\sigma\alpha(1-h)+9k(1-h)-27\geq 0$, i.e., $\alpha\geq\frac{27-9k(1-h)}{2k^{2}\sigma(1-h)}$, then $\tilde{F}(x_{*})$ is positive. If $2k^{2}\sigma\alpha(1-h)+9k(1-h)-27<0$, i.e., $\frac{1}{\sigma k}<\alpha<\frac{27-9k(1-h)}{2k^{2}\sigma(1-h)}$, then the sign of $\tilde{F}(x_{*})$ is same as that of

[TABLE]

which is deduced from that $2(\alpha k\sigma+3)(1-h)\sqrt{k\alpha\sigma(\alpha k\sigma+3)}$ square minus $\alpha\sigma\{2k^{2}\sigma\alpha(1-h)+9k(1-h)-27\}$ square. Since the leading coefficient of $F_{1}(\alpha)$ and $F_{1}(\frac{27-9k(1-h)}{2k^{2}\sigma(1-h)})=\frac{(hk-k+9)^{3}}{2k^{2}(1-h)}$ are positive and $F_{1}(\frac{1}{\sigma k})=\frac{(hk-k+1)(5hk-5k-27)}{k}$ is negative under the conditions $0<h<1$ and $0<k<k_{1}$, $F_{1}(\alpha)$ has one root $\alpha_{1}$ given in (11) in the interval $(\frac{1}{\sigma k},\frac{27-9k(1-h)}{2k^{2}\sigma(1-h)})$. Hence, we can immediately obtain that $\tilde{F}(x_{*})<0$ if $\frac{1}{\sigma k}<\alpha<\alpha_{1}$, $\tilde{F}(x_{*})>0$ if $\alpha>\alpha_{1}$ and $\tilde{F}(x_{*})=0$ if $\alpha=\alpha_{1}$. Accordingly, the distribution of roots of $F(x)$ in the interval $(0,k)$ is displayed as follows. $F(x)$ has no root in the interval $(0,k)$ if $(k,\sigma,\alpha)\in\mathcal{P}_{4}$; $F(x)$ has two roots in the interval $(0,k)$ denoted by $x_{1}$, $x_{2}$ and $x_{1}<x_{2}$ if $(k,\sigma,\alpha)\in\mathcal{P}_{5}$ (see Fig. 1 (a)); $F(x)$ has one multiple root $x_{*}$ in the interval $(0,k)$ if $(k,\sigma,\alpha)\in\mathcal{S}_{5}$ (see Fig. 1 (b)). Furthermore, $x_{*}$ also can be expressed as

[TABLE]

if $\alpha=\alpha_{1}$ because $\tilde{F}(x_{*})=0$ in (28). Corresponding to the roots of $F(x)$ in the interval $(0,k)$, the positive equilibria of system (8) are $E_{1}:(x_{1},y_{1})$, $E_{2}:(x_{2},y_{2})$ or $E_{*}:(x_{*},y_{*})$, where $y_{i}=\frac{\sigma x_{i}(k-x_{i})}{k}$, $i=1,2,*$. From the above discussion we obtain the number of equilibria of system (8) as shown in Table 1.

In what follows, we study the dynamical behaviors of equilibria. Compute the Jacobian matrix of vector field (8)

[TABLE]

where $J_{11}:=\sigma(k-2x)\{1+h(\alpha y+1)x\}+\{h\sigma x(k-x)-yk\}(\alpha y+1)$ and let $T$ and $D$ denote its trace and determinant respectively. $E_{0}$ is a saddle because of $D|_{E_{0}}=-\sigma k^{2}<0$. At $E_{k}$, $D|_{E_{k}}=k^{2}\sigma(hk+1)(hk-k+1)$, $T|_{E_{k}}=-k\{(\sigma+1)(hk+1)-k\}$, $T|_{E_{k}}-4D|_{E_{k}}=k^{2}\{(\sigma-1)(hk+1)+k\}^{2}$. When $h\geq 1$, $D|_{E_{k}}>0$, $T|_{E_{k}}<0$ and $T|_{E_{k}}-4D|_{E_{k}}>0$, implying that $E_{k}$ is a stable node. When $h<1$, the qualitative properties of $E_{k}$ are displayed as follows. $D|_{E_{k}}<0$ if $k>k_{1}$, implying that $E_{k}$ is a saddle, $D|_{E_{k}}>0$, $T|_{E_{k}}<0$ and $T|_{E_{k}}-4D|_{E_{k}}>0$ if $k<k_{1}$, implying that $E_{k}$ is a stable node, $T|_{E_{k}}<0$ and $D|_{E_{k}}=0$ if $k=k_{1}$, implying that $E_{k}$ is degenerate. At $E_{i}$, $i=1,2,*$, we obtain determinant $D|_{E_{i}}$ and trace $T|_{E_{i}}$ of Jacobian matrix $J$ as follows

[TABLE]

To obtain the afore-given expressions $D|_{E_{i}}$ and $T|_{E_{i}}$, we have used the branch $1+hx_{i}(1+\alpha y_{i})=x_{i}(1+\alpha y_{i})$ and the expression of $y_{i}$. Furthermore, $T|_{E_{i}}$ is the pseudo-remainder of trace $\tilde{T}|_{E_{i}}$ of Jacobian matrix $J$ at $E_{i}$ divided by $F(x_{i})$, where $\tilde{T}|_{E_{i}}:=\sigma x_{i}\{\alpha\sigma(h+1)x_{i}^{3}-\alpha k(2h\sigma-h+\sigma+1)x_{i}^{2}+k(\alpha hk\sigma-\alpha hk+\alpha k-h-1)x_{i}+hk^{2}\}/k$. Using the MAPLE command “prem”, we can simplify trace $\tilde{T}|_{E_{i}}$ as $T|_{E_{i}}$ since $F(x_{i})=0$. The above discussion of the existence of equilibria shows that $F^{\prime}(x_{1})>0$, $F^{\prime}(x_{*})=0$ and $F^{\prime}(x_{2})<0$. Thus, we obtain $D|_{E_{1}}>0$, $D|_{E_{*}}=0$ and $D|_{E_{2}}<0$, which imply that $E_{*}$ is degenerate, $E_{2}$ is a saddle and $E_{1}$ can be neither a saddle nor a degenerate equilibrium. We only need to discuss the sign of $T|_{E_{1}}$ in the following lemma.

Lemma 2.3**.**

For $0<h<1$, $T|_{E_{1}}>0$ if $(k,\sigma,\alpha)\in\mathcal{P}_{11}\cup\mathcal{S}_{21}\cup\mathcal{P}_{31}\cup\mathcal{S}_{42}\cup\mathcal{P}_{51}$, $T|_{E_{1}}<0$ if $(k,\sigma,\alpha)\in\mathcal{P}_{12}\cup\mathcal{S}_{22}\cup\mathcal{P}_{32}\cup\mathcal{S}_{41}\cup\mathcal{P}_{52}$ and $T|_{E_{1}}=0$ if $(k,\sigma,\alpha)\in\mathcal{S}_{11}\cup\mathcal{L}_{21}\cup\mathcal{S}_{31}\cup\mathcal{L}_{41}\cup\mathcal{S}_{51}$.

*Proof of Lemma 2.3: *** The above discussion shows that the equilibrium $E_{1}$ exists for $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{P}_{1}\cup\mathcal{S}_{2}\cup\mathcal{P}_{3}\cup\mathcal{S}_{4}\cup\mathcal{P}_{5}$. Determining the sign of $T|_{E_{1}}$ is a difficulty because the explicit solution $x_{1}$ can not be obtained from equilibrium equation (24), which is a cubic equation. In order to overcome it, we need to discuss the sign of $T|_{E_{1}}$ indirectly via the relative position of the roots of equilibrium equation (24) and $T|_{E_{1}}$ together with the monotonicity of $T|_{E_{1}}$. Function $T|_{E_{1}}$ is monotonically decreasing and has one positive root. Let the root of $T|_{E_{1}}=0$ be

[TABLE]

Substituting $x_{0}$ into $F(x)$, we get

[TABLE]

The concrete strategy is described as follows. For $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{P}_{1}\cup\mathcal{S}_{2}\cup\mathcal{P}_{3}\cup\mathcal{S}_{4}$, $F(x)$ has a unique root $x_{1}$ in the interval $(0,k)$ as well as $F(x)<0$ in the interval $(0,x_{1})$ and $F(x)>0$ in the interval $(x_{1},k)$. The relative position of $x_{1}$ and $x_{0}$ is determined by the sign of $F(x_{0})$ together with relationship $x_{0}<k$. Concretely, $x_{1}>x_{0}$ if $F(x_{0})<0$, $x_{1}<x_{0}$ if $F(x_{0})>0$ and $x_{1}=x_{0}$ if $F(x_{0})=0$ (see Fig. 1 (a)). In addition that $T|_{E_{1}}$ is monotonically decreasing, then $T|_{E_{1}}<0$ if $F(x_{0})<0$, $T|_{E_{1}}>0$ if $F(x_{0})>0$ and $T|_{E_{1}}=0$ if $F(x_{0})=0$. For $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{P}_{5}$, $F(x)$ has two roots $x_{1}$ and $x_{2}$ ($x_{1}<x_{*}<x_{2}$) as well as $F(x)<0$ in the intervals $(0,x_{1})\cup(x_{2},k)$ and $F(x)>0$ in the interval $(x_{1},x_{2})$. The relative position of $x_{1}$ and $x_{0}$ is determined by the sign of $F(x_{0})$ together with relative position of $x_{0}$ and $x_{*}$. Concretely, $x_{1}>x_{0}$ if $F(x_{0})<0$ and $x_{0}<x_{*}$; $x_{1}<x_{0}$ if $F(x_{0})>0$ or $F(x_{0})\leq 0$ and $x_{0}>x_{*}$; $x_{1}=x_{0}$ if $F(x_{0})=0$ and $x_{0}<x_{*}$ (see Fig. 1 (a)). In addition that $T|_{E_{1}}$ is monotonically decreasing, then $T|_{E_{1}}<0$ if $F(x_{0})<0$ and $x_{0}<x_{*}$; $T|_{E_{1}}>0$ if $F(x_{0})>0$ or $F(x_{0})\leq 0$ and $x_{0}>x_{*}$; $T|_{E_{1}}=0$ if $F(x_{0})=0$ and $x_{0}<x_{*}$. Thus, to obtain the parameter condition for each case is the subsequent task.

Because of space cause, we just give the proof in detail for $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{4}$, but omit the verbose proof of the rest cases. By analyzing $F(x_{0})$ for $0<h<1$ and $(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}$ we have the sign of $F(x_{0})$ as follows. $F(x_{0})>0$ if $(k,\sigma,\alpha)\in\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k\geq k_{2}~{}\mbox{or}~{}k_{1}\leq k<k_{2},\alpha>\alpha_{2}~{}\mbox{or}~{}k<k_{1},\sigma>1-h-k(1-h)^{2},\alpha>\alpha_{2}\}$; $F(x_{0})<0$ if $(k,\sigma,\alpha)\in\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k_{1}\leq k<k_{2},\alpha<\alpha_{2}~{}\mbox{or}~{}k<k_{1},\sigma\leq 1-h-k(1-h)^{2}~{}\mbox{or}~{}k<k_{1},\sigma>1-h-k(1-h)^{2},\alpha<\alpha_{2}\}$; $F(x_{0})=0$ if $(k,\sigma,\alpha)\in\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k_{1}\leq k<k_{2},\alpha=\alpha_{2}~{}\mbox{or}~{}k<k_{1},\sigma>1-h-k(1-h)^{2},\alpha=\alpha_{2}\}$, where $\alpha_{2}$ is given by (15). Furthermore, $x_{0}<k$ for $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{P}_{1}\cup\mathcal{S}_{2}\cup\mathcal{P}_{3}\cup\mathcal{S}_{4}$ since $x_{0}-k=\frac{k\{(1-h)(hk-k+1)-\sigma\}}{k(1-h)^{2}+\sigma(1+h)}<0$. Now we need to find the intersections of set $\mathcal{S}_{4}$ and sets of $F(x_{0})>0$, $F(x_{0})=0$ and $F(x_{0})<0$ respectively. In order to compare the endpoints $\alpha_{2}$ with $\frac{1}{\sigma k}$, we denote $\alpha_{2}-\frac{1}{\sigma k}$ by $f(\sigma)$ given in (18), where $f(\sigma)=\frac{(1-h)\{(2h+1)\sigma-2h+2\}}{(h\sigma+1-h)^{2}}>0$ for $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{4}$, implying $\alpha_{2}>\frac{1}{\sigma k}$. For $0<h<1$, we can get $F(x_{0})>0$ as $(k,\sigma,\alpha)\in\mathcal{S}_{42}$, $F(x_{0})<0$ as $(k,\sigma,\alpha)\in\mathcal{S}_{41}$ and $F(x_{0})=0$ as $(k,\sigma,\alpha)\in\mathcal{L}_{41}$. Thus, we obtain the corresponding sign of $T|_{E_{1}}$ for this case.

Although the proof for the case $(k,\sigma,\alpha)\in\mathcal{P}_{1}\cup\mathcal{S}_{2}\cup\mathcal{P}_{3}\cup\mathcal{P}_{5}$ is omitted, we should account for the two quantities $\sigma_{1}$ and $\sigma_{2}$. In the case $(k,\sigma,\alpha)\in\mathcal{P}_{1}\cup\mathcal{S}_{2}\cup\mathcal{P}_{3}$, we still need to compare the endpoints $\alpha_{2}$ and $\frac{1}{\sigma k}$ so that function $f(\sigma)$ given in (18) need to be discussed for $(h,k,\sigma)\in\{(h,k,\sigma)\in\mathbb{R}_{+}^{3}:h<1,k_{1}<k<k_{2}\}$, the properties of which are displayed as follows. $f(\sigma)>0$ if $(h,k,\sigma)\in\{(h,k,\sigma)\in\mathbb{R}_{+}^{3}:h<1,k_{1}<k<k_{3},\sigma>\sigma_{1}\}$; $f(\sigma)<0$ if $(h,k,\sigma)\in\{(h,k,\sigma)\in\mathbb{R}_{+}^{3}:h<1,k_{1}<k<k_{3},\sigma<\sigma_{1}~{}\mbox{or}~{}h<1,k_{3}\leq k<k_{2}\}$; $f(\sigma)=0$ if $(h,k,\sigma)\in\{(h,k,\sigma)\in\mathbb{R}_{+}^{3}:h<1,k_{1}<k<k_{3},\sigma=\sigma_{1}\}$, where $\sigma_{1}$ is the unique positive root of $f(\sigma)$ for $0<h<1$ and $k_{1}<k<k_{3}$. In the case $(k,\sigma,\alpha)\in\mathcal{P}_{5}$, we need to compare the endpoints $\alpha_{1}$ and $\alpha_{2}$ for $(h,k,\sigma)\in\{(h,k,\sigma)\in\mathbb{R}_{+}^{3}:h<1,k<k_{1},\sigma>1-h-k(1-h)^{2}\}$. Substituting $\alpha_{2}$ into $F_{1}(\alpha)$, we get $F_{1}(\alpha_{2})=f_{1}(\sigma)f_{2}^{2}(\sigma)/\{k^{2}(1-h)(h\sigma-h+1)^{4}\{k(1-h)^{2}+h+\sigma-1\}^{2}\}$, where

[TABLE]

$f_{1}(\sigma)>0$* for $\sigma>0$ because all the coefficients are positive for $0<h<1$ and $0<k<k_{1}$. Furthermore, the leading coefficient of $f_{2}(\sigma)$ is positive and the constant term of which is negative, implying that $f_{2}(\sigma)$ is monotonically increasing for $\sigma>0$ and has a unique positive root $\sigma_{2}$ given in (15). Since $f_{2}(1-h-k(1-h)^{2})=-(1-h)^{2}(kh-k+1)(h^{2}k-hk+h+1)^{2}<0$ for $0<h<1$ and $0<k<k_{1}$, we have $\sigma_{2}>1-h-k(1-h)^{2}$. Hence, we obtain $\alpha_{2}=\alpha_{1}$ for $0<h<1$, $0<k<k_{1}$ and $\sigma=\sigma_{2}$ as well as $\alpha_{2}>\alpha_{1}$ for $0<h<1$, $0<k<k_{1}$, $\sigma>1-h-k(1-h)^{2}$ and $\sigma\neq\sigma_{2}$. The proof of Lemma 2.3 is completed.*

The determinant of $E_{1}$ is positive and the sign of the trace of $E_{1}$ is shown in Lemma 2.3, the qualitative properties of $E_{1}$ can be derived, namely, $E_{1}$ is an unstable node or focus if $T|_{E_{1}}>0$, $E_{1}$ is a stable node or focus if $T|_{E_{1}}<0$ and $E_{1}$ is center type if $T|_{E_{1}}=0$. The stability and topological classification for the equilibria are presented in Table 1. The proof of Theorem 2.1 is completed.

2.1 Bifurcations at $E_{k}$ and $E_{*}$

In this section, we show that both transcritical and pitchfork bifurcations may occur at $E_{k}$ and a saddle-node bifurcation may occur at $E_{*}$. Table 1 of Theorem 2.1 indicates that system (8) has a degenerate equilibrium $E_{k}$ with $T|_{E_{k}}<0$ and $D|_{E_{k}}=0$ if $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{1}\cup\mathcal{L}_{1}\cup\mathcal{S}_{41}\cup\mathcal{L}_{41}\cup\mathcal{S}_{42}$, i.e., $k=k_{1}$. The following theorem displays the bifurcations at $E_{k}$.

Theorem 2.4**.**

*For $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{1}\cup\mathcal{L}_{1}\cup\mathcal{S}_{41}\cup\mathcal{L}_{41}\cup\mathcal{S}_{42}$, $E_{k}$ is a saddle-node of system (8). Moreover,

(i) as $(k,\sigma,\alpha)$ crosses $\mathcal{S}_{1}$, i.e., $(k,\sigma,\alpha)$ varies from $\mathcal{P}_{11}\cup\mathcal{P}_{12}$ to $\mathcal{P}_{2}$, a transcritical bifurcation happens at $E_{k}$ such that a stable (resp., unstable) node $E_{1}$ and two saddles $E_{0}$ and $E_{k}$ change into a stable node $E_{k}$ and a saddle $E_{0}$ for $(k,\sigma,\alpha)\in\mathcal{P}_{12}$ (resp., $(k,\sigma,\alpha)\in\mathcal{P}_{11}$).

(ii) as $(k,\sigma,\alpha)$ crosses $\mathcal{S}_{41}$, i.e., $(k,\sigma,\alpha)$ varies from $\mathcal{P}_{32}$ to $\mathcal{P}_{52}$, a transcritical bifurcation happens at $E_{k}$ such that a stable node $E_{1}$ and two saddles $E_{0}$ and $E_{k}$ change into two stable nodes $E_{1}$ and $E_{k}$ and two saddles $E_{0}$ and $E_{2}$.

(iii) as $(k,\sigma,\alpha)$ crosses $\mathcal{L}_{41}$, i.e., $(k,\sigma,\alpha)$ varies from $\mathcal{S}_{31}$ to $\mathcal{S}_{51}$, a transcritical bifurcation happens at $E_{k}$ such that a center type equilibrium $E_{1}$ and two saddles $E_{0}$ and $E_{k}$ change into a center type equilibrium $E_{1}$, a stable node $E_{k}$ and two saddles $E_{0}$ and $E_{2}$.

(iv) as $(k,\sigma,\alpha)$ crosses $\mathcal{S}_{42}$, i.e., $(k,\sigma,\alpha)$ varies from $\mathcal{P}_{31}$ to $\mathcal{P}_{51}$, a transcritical bifurcation happens at $E_{k}$ such that two saddles $E_{0}$ and $E_{k}$ and an unstable node $E_{1}$ change into an unstable node $E_{1}$, a stable node $E_{k}$ and two saddles $E_{0}$ and $E_{2}$.

(v) as $(k,\sigma,\alpha)$ crosses $\mathcal{L}_{1}$, i.e., $(k,\sigma,\alpha)$ varies from $\mathcal{S}_{21}\cup\mathcal{S}_{22}$ to $\mathcal{S}_{3}$, a pitchfork bifurcation happens at $E_{k}$ such that two saddles $E_{0}$ and $E_{k}$ and a stable (resp., unstable) node $E_{1}$ change into a saddle $E_{0}$ and a stable node $E_{k}$ for $(k,\sigma,\alpha)\in\mathcal{S}_{22}$ (resp., $(k,\sigma,\alpha)\in\mathcal{S}_{21}$).*

Proof 2.5**.**

Let $\epsilon=k-k_{1}$. For sufficiently small $|\epsilon|$, consider system (8) suspected by the parameter $\epsilon$. Using the linear transformation $x=u+v+k$, $y=-\sigma u$ and time-rescaling $\tau:=\frac{-\sigma}{(h-1)^{2}}t$ to translate $E_{k}$ to the origin $(0,0)$ and diagonalize the linear part of the suspected system, we can change the system into the follows

[TABLE]

By Theorem 1 of Carr [1981], system (49) has a two-dimensional center manifold $W^{c}:v=h(u,\epsilon)$ near the origin, which is $C^{\infty}$ and tangent to the plane $v=0$ at the origin in the $(u,v,\epsilon)$-space. Let

[TABLE]

Since it is invariant to solutions $(u(t),v(t),\epsilon(t))$ of system (49), we can differentiate both sides of (50), which leads to the equality $\dot{v}=h_{u}\dot{u}+h_{\epsilon}\dot{\epsilon}$. Substituting equations of (49) into the equality and comparing the coefficients of $u^{2}$, $\epsilon^{2}$ and $u\epsilon$, we get $a={\left(\alpha\,{\sigma}^{2}-\alpha\,\sigma+h\sigma-h+1\right)\left(-1+h\right)/\sigma}$, $b=0$ and $c={\left(\sigma-1\right)\left(-1+h\right)^{2}/\sigma}$. Thus, system (49) restricted to center manifold (50) can be written as

[TABLE]

where $c_{1}:=(h-1)^{2}\{2\alpha\sigma^{2}-(h-1)(h-2)\sigma+(h-1)^{2}\}/\sigma^{2}$ and $c_{2}:=-(h-1)^{2}\{\alpha(h-2)\sigma^{2}-(h-1)(\alpha-h)\sigma-(h-1)^{2}\}/\sigma^{2}$.

When $\alpha\neq\frac{1-h}{\sigma}$, it shows that $\frac{(\alpha\sigma+h-1)(1-h)}{\sigma}\neq 0$ in (52) and the origin is the unique equilibrium as $\epsilon=0$ and another equilibrium arises from the origin as $\epsilon\neq 0$. Moreover, the stabilities of the equilibria exchange as $\epsilon$ varies from negative to positive. Thus, $E_{k}$ is a saddle-node as $\epsilon=0$ and system (8) undergoes a transcritical bifurcation at $E_{k}$ for $(k,\sigma,\alpha)\in\mathcal{S}_{1}\cup\mathcal{S}_{41}\cup\mathcal{L}_{41}\cup\mathcal{S}_{42}$ ([Guckenheimer & Holmes, 1983, p.149]). When $\alpha=\frac{1-h}{\sigma}$, it shows that $\frac{(\alpha\sigma+h-1)(1-h)}{\sigma}=0$ and $\frac{2(1-h)^{3}}{\sigma}\neq 0$ in (52) and the origin is the unique equilibrium as $\epsilon=0$ and the other two equilibria arise from the origin as $\epsilon>0$. Thus, $E_{k}$ is a saddle-node as $\epsilon=0$ and system (8) undergoes a pitchfork bifurcation at $E_{k}$ for $(k,\sigma,\alpha)\in\mathcal{L}_{1}$ ([Guckenheimer & Holmes, 1983, p.149]). The proof is completed.

As indicated in Theorem 2.1, system (8) has a degenerate equilibrium $E_{*}$ for $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{5}$, i.e., $D|_{E_{*}}=0$. To consider what bifurcation system (8) undergoes for this degenerate case, let us first discuss the sign of the trace $T|_{E_{*}}$ given in (36). Substituting $\alpha=\alpha_{1}$ and $x_{*}$ (given in (11) and (34) respectively) into $T|_{E_{*}}$, we obtain

[TABLE]

For $0<h<1$ and $0<k<k_{1}$, the sign of $(h\sigma-h+1)\sqrt{(kh-k+1)(hk-k+9)}$ is always positive, but that of $-(kh^{2}-hk+h+4)\sigma+3(1-h)(hk-k+1)$ is indeterminate. If $-(kh^{2}-hk+h+4)\sigma+3(1-h)(hk-k+1)\geq 0$ in (55), i.e., $0<\sigma\leq\frac{3(1-h)(hk-k+1)}{kh^{2}-hk+h+4}$, it is evident that $T|_{E_{*}}>0$. If $-(kh^{2}-hk+h+4)\sigma+3(1-h)(hk-k+1)<0$, i.e., $\sigma>\frac{3(1-h)(hk-k+1)}{kh^{2}-hk+h+4}$, we can derive the following relationship

[TABLE]

Based on the fact that $f_{2}(\sigma)>0$ if $\sigma>\sigma_{2}$, $f_{2}(\sigma)<0$ if $0<\sigma<\sigma_{2}$ and $f_{2}(\sigma)=0$ if $\sigma=\sigma_{2}$ as well as the following inequality

[TABLE]

we conclude that $T|_{E_{*}}>0$ if $\frac{3(1-h)(hk-k+1)}{kh^{2}-hk+h+4}<\sigma<\sigma_{2}$, $T|_{E_{*}}<0$ if $\sigma>\sigma_{2}$ and $T|_{E_{*}}=0$ if $\sigma=\sigma_{2}$. From the discussion, the sign of $T|_{E_{*}}$ is obtained for $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{5}$, namely, $T|_{E_{*}}>0$ if $0<\sigma<\sigma_{2}$, $T|_{E_{*}}<0$ if $\sigma>\sigma_{2}$ and $T|_{E_{*}}=0$ if $\sigma=\sigma_{2}$.

For $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{5}$ with $\sigma\neq\sigma_{2}$, the following theorem displays that system (8) undergoes a saddle-node bifurcation at $E_{*}$.

Theorem 2.6**.**

For $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{5}$ with $\sigma\neq\sigma_{2}$, $E_{*}$ is a saddle-node of system (8) and a saddle-node bifurcation happens at $E_{*}$ as $(k,\sigma,\alpha)$ crosses $\mathcal{S}_{5}$. Moreover, as $(k,\sigma,\alpha)$ changes from $\mathcal{P}_{4}$ to $\mathcal{P}_{51}\cup\mathcal{P}_{52}$, an unstable (resp., stable) node $E_{1}$ and a saddle $E_{2}$ arise for $(k,\sigma,\alpha)\in\mathcal{P}_{51}$ (resp., $(k,\sigma,\alpha)\in\mathcal{P}_{52}$).

Proof 2.7**.**

Let $\epsilon=\alpha-\alpha_{1}$. For sufficiently small $|\epsilon|$, consider system (8) suspected by the parameter $\epsilon$. By translating $E_{*}$ to the origin $(0,0)$ we can expand the suspected system as follows

[TABLE]

where the coefficients $a_{ijk}$ and $b_{ijk}$ are given in the Appendix with $\alpha_{1}$ and $x_{*}$ given in (11) and (34) respectively. Using the linear transformation $x=u+\frac{a_{100}}{b_{100}}v+\frac{a_{001}(a_{010}-a_{100})}{a_{100}(a_{100}+b_{010})}\epsilon$ and $y=-\frac{b_{100}}{b_{010}}u+v$ to diagonalize the linear part of the suspected system, we obtain the following form

[TABLE]

where $p_{ijk}$ and $q_{ijk}$ are displayed in the Appendix. By Theorem 1 of Carr [1981], system (65) has a two-dimensional center manifold $W^{c}:v=h_{1}(u,\epsilon)$ near the origin, which is $C^{\infty}$ and tangent to the plane $v=0$ at the origin in the $(u,v,\epsilon)$-space. Let

[TABLE]

Differentiating both sides of (66) leads to the equality $\dot{v}=h_{1_{u}}\dot{u}+h_{1_{\epsilon}}\dot{\epsilon}$. Substituting equations of (65) into the equality and comparing the coefficients of $u^{2}$, $\epsilon^{2}$ and $u\epsilon$, we obtain $c_{20}$, $c_{11}$ and $c_{02}$ given in the Appendix respectively. System (65) restricted to center manifold (66) can be written as

[TABLE]

where

[TABLE]

In the following, we prove $d_{2}(0)\neq 0$ for $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{5}$ with $\sigma\neq\sigma_{2}$, where

[TABLE]

Denote the numerator of $d_{2}(0)$ by $\tilde{d_{2}}(0)$, it follows immediately that $\tilde{d_{2}}(0)=-x_{*}^{4}\alpha_{1}^{2}\sigma^{3}(k-x_{*})^{2}(h-1)^{2}f_{3}(x_{*})/k$ with $f_{3}(x_{*}):=-\alpha_{1}\sigma(h+2)x_{*}^{3}+\alpha_{1}k\sigma(h+3)x_{*}^{2}-k(\alpha_{1}k\sigma-h-2)x_{*}-k(k-1)$, which is a cubic polynomial. Since $F^{\prime}(x)$ is a quartic polynomial and $F^{\prime}(x_{*})=0$, we use the Maple command “prem” to get the pseudo-remainder of $f_{3}(x_{*})$ divided by $F^{\prime}(x_{*})$.

[TABLE]

where $m=9\alpha_{1}^{2}\sigma^{2}(-1+h)^{2}$. Substituting $x_{*}$ given by (34) into the pseudo-remainder leads to $f_{3}(x_{*})=\frac{k}{18(k\alpha_{1}\sigma+3)(h-1)}\{-9k\sigma(hk-k+3)\alpha_{1}-36(h-1)k-162\}>0$ for $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{5}$. Thus, $\tilde{d_{2}}(0)<0$ for $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{5}$. In the denominator of $d_{2}(0)$, $b_{010}>0$ is obvious and $a_{100}+b_{010}$, i.e., the trace $T|_{E_{*}}$, has been discussed before this theorem. Hence, for $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{5}$ with $\sigma\neq\sigma_{2}$, we have $d_{2}(0)>0$ if $\sigma>\sigma_{2}$ and $d_{2}(0)<0$ if $0<\sigma<\sigma_{2}$.

Using the translation $u=w-\frac{d_{1}(\epsilon)}{2d_{2}(\epsilon)}$ and time-rescaling $\tau:=d_{2}(\epsilon)t$ to system (68), we get

[TABLE]

where $\zeta(\epsilon):=\{4d_{0}(\epsilon)d_{2}(\epsilon)-d_{1}^{2}(\epsilon)\}/4d_{2}^{2}(\epsilon)$. The computation yields $\zeta(0)=0$ and $\zeta^{\prime}(0)=d_{0}^{\prime}(0)/d_{2}(0)$, where $d_{0}^{\prime}(0)=-(a_{100}+b_{100})b_{001}b_{010}/b_{100}(a_{100}+b_{010})$. It is obvious that both $b_{100}$ and $b_{001}$ are positive. In addition, $a_{100}+b_{100}=\sigma x_{*}\{\alpha_{1}\sigma(k-x_{*})x_{*}+k\}(k-2x_{*})/k<0$ because $k-2x_{*}=-\frac{-k\alpha_{1}\sigma+2\sqrt{k\alpha_{1}\sigma(\alpha_{1}k\sigma+3)}}{\alpha_{1}\sigma}<0$. Thus, $\zeta^{\prime}(0)<0$ for $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{5}$ with $\sigma\neq\sigma_{2}$.

Hence, the origin is the unique equilibrium of (72) as $\epsilon=0$ and two equilibria arise from the origin as $\epsilon$ varies from [math] to positive when $\sigma\neq\sigma_{2}$. Therefore, for $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{5}$ with $\sigma\neq\sigma_{2}$, a saddle-node bifurcation occurs at $E_{*}$ as $\alpha$ changes from $\alpha=\alpha_{1}$ to $\alpha>\alpha_{1}$ such that an unstable (stable) node $E_{1}$ and a saddle $E_{2}$ emerge from $E_{*}$ if $\sigma<\sigma_{2}$ (resp. $\sigma>\sigma_{2}$). The proof is completed.

3 Hopf Bifurcation at $E_{1}$

As indicated in Theorem 2.1, $E_{1}(x_{1},y_{1})$ is of center type for $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{11}\cup\mathcal{L}_{21}\cup\mathcal{S}_{31}\cup\mathcal{L}_{41}\cup\mathcal{S}_{51}$, i.e., $T|_{E_{1}}=0$ and $D|_{E_{1}}>0$, where $x_{1}:=x_{0}$ given in (38) and $y_{1}:=\sigma k\{k(h-1)^{2}+h+\sigma-1\}(h\sigma-h+1)/\{k(h-1)^{2}+\sigma(h+1)\}^{2}$. In this section, we show that $E_{1}$ is a weak focus of multiplicity at most 1 and the Hopf bifurcation occurs at $E_{1}$. For convenience, let $\mathcal{D}:=\mathcal{S}_{11}\cup\mathcal{L}_{21}\cup\mathcal{S}_{31}\cup\mathcal{L}_{41}\cup\mathcal{S}_{51}=\{(k,\sigma,\alpha)\in\mathbb{R}_{+}^{3}:k_{1}\leq k<k_{2},\alpha=\alpha_{2}~{}\mbox{or}~{}k<k_{1},\sigma>\sigma_{2},\alpha=\alpha_{2}\}$.

Theorem 3.1**.**

For $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{D}$, equilibrium $E_{1}$ of system (8) is a stable weak focus of multiplicity 1 and one stable limit cycle arises near $E_{1}$ as $\alpha$ varies from $\alpha=\alpha_{2}$ to $\alpha>\alpha_{2}$.

Proof 3.2**.**

Translating equilibrium $E_{1}$ to the origin, system (8) becomes the following system

[TABLE]

where the coefficients are given in the Appendix. For $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{D}$, system (75) has a pair of purely imaginary eigenvalues $\pm\beta$, where $\beta:=k\sqrt{\sigma f_{2}(\sigma)}/\{(k(h-1)^{2}+\sigma(h+1))\sqrt{1-h}\}$. The transversal condition of Hopf bifurcation holds because

[TABLE]

In the following we compute the quantity of focus. Using the linear transformation $x=\frac{1}{b_{10}}u+\frac{a_{10}}{b_{10}\beta}v$, $y=\frac{1}{\beta}v$ and time-rescaling $\tau:=\beta t$ to normalize the linear part, we can change system (75) into the form

[TABLE]

where the coefficients are given in the Appendix. The following is devoted to the center-focus determination by the successive function method (Zhang et al. [1992]). We can obtain the first order focal value

[TABLE]

where $G(\sigma):=L_{4}\sigma^{4}+L_{3}\sigma^{3}+L_{2}\sigma^{2}+L_{1}\sigma+L_{0}$ with

[TABLE]

The sign of $g$ is determined by that of $G(\sigma)$. We first show $G(\sigma)<0$ for $(h,k,\sigma)\in\{(h,k,\sigma)\in\mathbb{R}_{+}^{3}:h<1,k_{1}\leq k<k_{2}\}$ by proving that all coefficients $L_{i}$ $(i=0,1,2,3,4)$ of $G(\sigma)$ are nonpositive. It is easy to check that $L_{4}<0$ and $L_{0}\leq 0$. In fact, the third factor is negative and the other two are positive in $L_{4}$ and the third factor is nonpositive and the other three are positive in $L_{0}$. To prove $L_{1}<0$, let the last factor of $L_{1}$ be $L_{11}(k)$. Since $L_{11}(k)$ is negative at the both endpoints of the interval $[k_{1},k_{2})$, Lemma 3.1 of Yang, [1999] indicates that the number of the roots for $L_{11}(k)$ in the interval $(k_{1},k_{2})$ is equal to that of positive roots for

[TABLE]

It is easily seen that $\Phi(z)$ has no positive root. We thus infer that $L_{11}(k)<0$, implying $L_{1}=k(h-1)^{4}L_{11}(k)<0$. To see $L_{2}<0$, let the second factor of $L_{2}$ be $L_{21}(k)$ and the derivative of which be

[TABLE]

The facts that the leading coefficient of $L^{\prime}_{21}(k)$ is negative and $L^{\prime}_{21}(k)$ is positive at the endpoints of the interval $[k_{1},k_{2})$ imply that $L^{\prime}_{21}(k)>0$ for $k_{1}\leq k<k_{2}$. Furthermore, since $L_{21}(k_{1})=\frac{6h+5}{1-h}>0$, we deduce $L_{21}(k)>0$. Hence, $L_{2}=(h-1)^{3}L_{21}(k)<0$. To show $L_{3}<0$, let the second factor of $L_{3}$ be $L_{31}(k)$. Analysis similar to that in the proof of $L_{21}(k)$ shows that $L_{31}(k)>0$. Hence, $L_{3}=(h-1)L_{31}(k)<0$. Consequently, it follows that $G(\sigma)<0$ for $(h,k,\sigma)\in\{(h,k,\sigma)\in\mathbb{R}_{+}^{3}:h<1,k_{1}\leq k<k_{2}\}$.

We proceed to show $G(\sigma)<0$ for $(h,k,\sigma)\in\{(h,k,\sigma)\in\mathbb{R}_{+}^{3}:h<1,k<k_{1},\sigma>\sigma_{2}\}$. In order to avoid discussing the monotonicity of $G(\sigma)$, we make the transformation $\sigma=\rho+\sigma_{2}$ to transform the problem of determining the sign of $G(\sigma)$ in the interval $(\sigma_{2},+\infty)$ to the issue of determining the sign of $\tilde{G}(\rho)$ in the interval $(0,+\infty)$, where $\tilde{G}(\rho):=\tilde{L}_{4}\rho^{4}+\tilde{L}_{3}\rho^{3}+\tilde{L}_{2}\rho^{2}+\tilde{L}_{1}\rho+\tilde{L}_{0}$ with

[TABLE]

in which $\tilde{L}_{2}$, $\tilde{L}_{1}$ and $\tilde{L}_{0}$ were reduced by the pseudo-division since $f_{2}(\sigma_{2})=0$. Likewise, we show $\tilde{G}(\rho)<0$ by proving that all coefficients $\tilde{L}_{i}$ $(i=0,1,2,3,4)$ are negative for $0<h<1$ and $0<k<k_{1}$. It follows immediately that $\tilde{L}_{4}<0$ because the last factor of $\tilde{L}_{4}$ is negative and the others are positive. To prove $\tilde{L}_{0}<0$, let the last factor of $\tilde{L}_{0}$ be $\tilde{L}_{01}(\sigma_{2})$, the constant term of which is positive. Let $LC_{0}(k)$ be the leading coefficient of $\tilde{L}_{01}(\sigma_{2})$. The fact that

[TABLE]

has no positive root shows, by Lemma 3.1 of Yang, [1999], that $LC_{0}(k)$ has no root in the interval $(0,k_{1})$. Since $LC_{0}(0)>0$, we immediately obtain $LC_{0}(k)>0$. Therefore, we have $\tilde{L}_{01}(\sigma_{2})>0$, which implies $\tilde{L}_{0}=(h-1)^{3}(h^{2}k-hk+h+1)^{4}\tilde{L}_{01}(\sigma_{2})<0$. In the following, we omit the details of the proof about $\tilde{L}_{i}<0$ ($i=1,2,3$). We claim that $\tilde{L}_{1}<0$. In fact, let the last factor of $\tilde{L}_{1}$ be $\tilde{L}_{11}(\sigma_{2})$ and obtain which is positive by analyzing the monotonicity. Since $-(h-1)^{2}(h^{2}k-hk+h+1)^{3}<0$, we conclude $\tilde{L}_{1}=-(h-1)^{2}(h^{2}k-hk+h+1)^{3}\tilde{L}_{11}(\sigma_{2})<0$. We claim that $\tilde{L}_{2}<0$. In fact, let the last factor of $\tilde{L}_{2}$ be $\tilde{L}_{21}(\sigma_{2})$ and obtain $\tilde{L}_{21}(\sigma_{2})>0$, which can derive $\tilde{L}_{2}=(h-1)(h^{2}k-hk+h+1)^{2}\tilde{L}_{21}(\sigma_{2})<0$. In the same manner, we can see that $\tilde{L}_{3}<0$. Consequently, we can assert that $\tilde{G}(\rho)<0$ for $\rho>0$, namely that $G(\sigma)<0$ for $(h,k,\sigma)\in\{(h,k,\sigma)\in\mathbb{R}_{+}^{3}:h<1,k<k_{1},\sigma>\sigma_{2}\}$.

We obtain the desired conclusion that the first order focal value $g$ is negative for $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{D}$. Therefore, the equilibrium $E_{1}$ of system (8) is a stable weak focus of multiplicity 1 and at most one stable limit cycle arises near $E_{1}$ from Hopf bifurcation as $\alpha$ varies from $\alpha=\alpha_{2}$ to $\alpha>\alpha_{2}$.

4 Bogdanov-Takens bifurcation at $E_{*}$

As presented before Theorem 2.6, $E_{*}(x_{*},y_{*})$ is degenerate with $D|_{E_{*}}=0$ and $T|_{E_{*}}=0$ for $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{5}$ with $\sigma=\sigma_{2}$, where $x_{*}=\frac{hk\sigma_{2}+k(1-h)}{k(1-h)^{2}+\sigma_{2}(h+1)}$ and $y_{*}=\frac{\sigma_{2}k(h\sigma_{2}-h+1)\{(h-1)(hk-k+1)+\sigma_{2}\}}{\{k(1-h)^{2}+\sigma_{2}(1+h)\}^{2}}$. Since $\alpha_{1}=\alpha_{2}$ if $\sigma=\sigma_{2}$, we let $\alpha_{*}:=\alpha_{1}=\alpha_{2}$. In the section, we display that $E_{*}$ is a cusp and the Bogdanov-Takens bifurcation may occur at $E_{*}$.

Lemma 4.1**.**

For $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{5}$ with $\sigma=\sigma_{2}$, the equilibrium $E_{*}$ of system (8) is a cusp.

Proof 4.2**.**

For $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{5}$ with $\sigma=\sigma_{2}$, system (8) can be transformed into the following form by translating $E_{*}$ to the origin

[TABLE]

where the coefficients are given in the Appendix with $\sigma_{2}$ given in (15). Using the linear transformation $x=-\frac{B_{01}}{B_{10}}u+v$ and $y=u$ combined with the time-rescaling $\tau:=B_{10}t$ to change system (89) into the canonical form

[TABLE]

where

[TABLE]

By the near-identity transformation $u_{1}:=u$ and $v_{1}:=v+\mathcal{A}_{20}u^{2}+\mathcal{A}_{11}uv\cdots$, system (92) can be written as the Kukles form

[TABLE]

Using a further transformation $u_{2}:=u_{1}$ and $v_{2}:=v_{1}-(\mathcal{A}_{11}+\mathcal{B}_{02})u_{1}v_{1}$ and the time-rescaling $t:=\{1+(\mathcal{A}_{11}+\mathcal{B}_{02})u_{2}\}\tau$ to eliminate the term of $v_{1}^{2}$ in (95), the system can be changed into

[TABLE]

We can assert that the coefficients $\mathcal{B}_{20}$ and $2\mathcal{A}_{20}+\mathcal{B}_{11}$ are nonzero for $0<h<1$ and $0<k<k_{1}$. In fact, we can obtain

[TABLE]

with $g_{1}(\sigma_{2}):=\{h(h-1)(h^{2}+2h-1)k+(h+1)(h^{2}+2h-2)\}\sigma_{2}^{2}+(h-1)\{h(h-1)^{3}k^{2}-(4h+3)(h-1)k-(h+4)(h+1)\}\sigma_{2}-k(h-1)^{3}(h^{2}k+h-k+2)$. It is easy to check that $\mathcal{B}_{20}<0$. Since the pseudo remainder of $g_{1}(\sigma_{2})$ divided by $f_{2}(\sigma_{2})$ is $2(1-h)\{h(hk-k+1)+1\}^{2}\{(h^{2}k(h-1)+(h+1)^{2})\sigma_{2}+k(h-1)^{2}\}$, which is positive, we obtain $2\mathcal{A}_{20}+\mathcal{B}_{11}>0$. Then, by the rescaling $u_{3}:=(2\mathcal{A}_{20}+\mathcal{B}_{11})^{2}u_{1}/\mathcal{B}_{20}$, $v_{3}:=-(2\mathcal{A}_{20}+\mathcal{B}_{11})^{3}v_{2}/\mathcal{B}_{20}^{2}$ and $\tau:=-\mathcal{B}_{20}t/(2\mathcal{A}_{20}+\mathcal{B}_{11})$ system (97) becomes

[TABLE]

It follows by Theorem 8.4 of Kuznetsov [1995] that $E_{*}$ is a cusp of system (8) for $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{5}$ with $\sigma=\sigma_{2}$. The proof of this lemma is completed.

We proceed to display that the Bogdanov-Takens bifurcation may occur at $E_{*}$ in the following theorem. We choose $\sigma$ and $\alpha$ as the bifurcation parameters and unfold the Bogdanov-Takens normal form of codimension 2 when the parameters $(\sigma,\alpha)$ are perturbed near the point $(\sigma_{2},\alpha_{*})$.

Theorem 4.3**.**

For $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{S}_{5}$ with $\sigma=\sigma_{2}$, there is a neighborhood $U$ of the point $(\sigma_{2},\alpha_{*})$ in the $(\sigma,\alpha)$-space and four curves

[TABLE]

such that system (8) undergoes a saddle-node bifurcation near $E_{*}$ as $(\sigma,\alpha)$ crossing $\mathcal{SN}^{+}\cup\mathcal{SN}^{-}$, a Hopf bifurcation near $E_{*}$ as $(\sigma,\alpha)$ crossing $\mathcal{H}$ and a homoclinic bifurcation near $E_{*}$ as $(\sigma,\alpha)$ crossing $\mathcal{HL}$, where $\alpha_{3}:=\alpha_{*}-\frac{\mu_{101}}{\mu_{110}}(\sigma-\sigma_{2})-\{\frac{A(0,0)(\mu_{101}^{2}\mu_{120}-\mu_{101}\mu_{110}\mu_{111}+\mu_{102}\mu_{110}^{2})}{A(0,0)\mu^{3}_{110}}+\frac{6(\mu_{101}\mu_{210}-\mu_{110}\mu_{201})^{2}}{25A(0,0)\mu^{3}_{110}}\}(\sigma-\sigma_{2})^{2}+O(|\sigma-\sigma_{2}|^{3})$ with $A(0,0)=\mathcal{B}_{20}$ and $\mu_{lij}$ displayed in the Appendix.

Proof 4.4**.**

Let $\epsilon_{1}:=\alpha-\alpha_{*}$ and $\epsilon_{2}:=\sigma-\sigma_{2}$. For sufficiently small $|\epsilon_{1}|$ and $|\epsilon_{2}|$, we can transform system (8) into the following form by translating $E_{*}$ to the origin and using the same translation as (92)

[TABLE]

where the coefficients are given in the Appendix. With the change of variables $(x,y)\rightarrow(u_{1},v_{1})$, where $u_{1}:=x$ and $v_{1}$ denotes the right side of the first equation in (104), system (104) can be written as the Kukles form, whose second order truncation is the following form

[TABLE]

where the coefficients are given in the Appendix. Since $\mathcal{F}_{11}(0,0)=2\mathcal{A}_{20}+\mathcal{B}_{11}>0$, we can use a parameter-dependent shift $u_{2}:=u_{1}+\frac{\mathcal{F}_{01}(\epsilon_{1},\epsilon_{2})}{\mathcal{F}_{11}(\epsilon_{1},\epsilon_{2})}$ and $v_{2}:=v_{1}$ to vanish the term proportional to $v_{1}$ in the second equation of system (107), which leads to the following system

[TABLE]

Using the near-identity transformation $u_{3}:=u_{2}$, $v_{3}:=v_{2}-\mathcal{F}_{02}(\epsilon_{1},\epsilon_{2})u_{2}v_{2}$ and time-rescaling $\tau:=(1+\mathcal{F}_{02}(\epsilon_{1},\epsilon_{2})u_{3})t$, system (111) can be changed into

[TABLE]

where

[TABLE]

We can check that $A(0,0)=\mathcal{B}_{20}<0$ and $B(0,0)=2\mathcal{A}_{20}+\mathcal{B}_{11}>0$. Thus, by the rescaling $u_{4}:=B^{2}(\epsilon_{1},\epsilon_{2})u_{3}/A(\epsilon_{1},\epsilon_{2})$, $v_{4}:=-B^{3}(\epsilon_{1},\epsilon_{2})v_{3}/A^{2}(\epsilon_{1},\epsilon_{2})$ and $t:=-A(\epsilon_{1},\epsilon_{2})\tau/B(\epsilon_{1},\epsilon_{2})$ system (113) can be changed into

[TABLE]

where

[TABLE]

Because the coefficients $E_{00}$, $E_{10}$, $F_{00}$, $F_{10}$ and $F_{01}$ in system (104) are equal to zero if $\epsilon_{1}=0$ and $\epsilon_{2}=0$, we can check $\mu_{1}(0,0)=0$ and $\mu_{2}(0,0)=0$. Consequently, we conclude that $\beta_{1}(0,0)=0$ and $\beta_{2}(0,0)=0$. Moreover, the Jacobian determinant of (118) at $(0,0)$ is given by

[TABLE]

where

[TABLE]

We utilize the theory of complete discrimination system for parametric polynomials in Yang, [1999] to determine the number of real roots of $g_{2}(k)$ and $g_{3}(k)$ in the interval $k\in(0,k_{1})$ with $0<h<1$. Let $k=\frac{1}{(1-h)(1+x^{2})}$, the number of real roots for $g_{2}(k)$ in the interval $(0,k_{1})$ is equal to the half number of that for $\tilde{g}_{2}(x)$ on the total real axis, where

[TABLE]

The discriminant sequence of $\tilde{g}_{2}(x)$ is $D:=\{D_{1},D_{2},D_{3},D_{4},D_{5},D_{6},D_{7},D_{8}\}$, where

[TABLE]

with $D_{31}$, $D_{41}$, $D_{51}$, $D_{61}$ and $D_{71}$ listed in the Appendix. It is obvious that $\mbox{sign}(D_{1})=1$ and $\mbox{sign}(D_{2})=-1$ for $0<h<1$. To discuss the sign of $D_{i}$ ($i=3\cdots 8$), we begin by considering the zeros of the single-variable function $D_{j}$ ($j=31,41,51,61,71$). Using the Maple command $``realroot(D_{j},0.000001)"$ to isolate the real roots of $D_{j}$ in the interval $(0,1)$. We see that $D_{31}$ has exactly one real root $h_{1}$ covered by $I_{1}:=[5686495/16777216,177703/524288]$. By computing $D_{31}$ at the endpoints of $I_{1}$, i.e., $D_{31}(5686495/16777216)<0$ and $D_{31}(177703/524288)>0$, we obtain the sign of $D_{31}$ as follows. $D_{31}<0$ if $h\in(0,h_{1})$, $D_{31}=0$ if $h=h_{1}$ and $D_{31}>0$ if $h\in(h_{1},1)$. $D_{41}$ has exactly one real root $h_{2}$ covered by $I_{2}:=[5448295/16777216,681037/2097152]$. Since $D_{41}(5448295/16777216)>0$ and $D_{41}(681037/2097152)<0$, the sign of $D_{41}$ is that $D_{41}>0$ if $h\in(0,h_{2})$, $D_{41}=0$ if $h=h_{2}$ and $D_{41}<0$ if $h\in(h_{2},1)$. $D_{51}$ has three real roots $h_{3}$, $h_{4}$ and $h_{5}$ covered by $I_{3}:=[3141071/8388608,6282143/16777216]$, $I_{4}:=[4311003/8388608,1077751/2097152]$ and $I_{5}:=[6104019/8388608,1526005/2097152]$ respectively. Similarly, by computing the sign of $D_{51}$ at the endpoints of the intervals $I_{3}$, $I_{4}$ and $I_{5}$ we obtain $D_{51}<0$ if $h\in(0,h_{3})\cup(h_{4},h_{5})$, $D_{51}=0$ if $h=h_{3}\cup h_{4}\cup h_{5}$ and $D_{51}>0$ if $h\in(h_{3},h_{4})\cup(h_{5},1)$. $D_{61}$ also has exactly three real roots $h_{6}$, $h_{7}$ and $h_{8}$ covered by $I_{6}:=[6273181/16777216,3136591/8388608]$, $I_{7}:=[4310379/8388608,1077595/2097152]$ and $I_{8}:=[6129165/8388608,3064583/4194304]$ respectively. We conclude similarly that $D_{61}>0$ if $h\in(0,h_{6})\cup(h_{7},h_{8})$, $D_{61}=0$ if $h=h_{6}\cup h_{7}\cup h_{8}$ and $D_{61}<0$ if $h\in(h_{6},h_{7})\cup(h_{8},1)$. $D_{71}$ has exactly three real roots $h_{9}$, $h_{10}$ and $h_{11}$ covered by $I_{9}:=[3160567/8388608,6321135/16777216]$, $I_{10}:=[4283093/8388608,2141547/4194304]$ and $I_{11}:=[6960901/8388608,3480451/4194304]$ respectively. The computation yields that $D_{71}<0$ if $h\in(0,h_{9})\cup(h_{10},h_{11})$, $D_{71}=0$ if $h=h_{9}\cup h_{10}\cup h_{11}$ and $D_{71}>0$ if $h\in(h_{9},h_{10})\cup(h_{11},1)$. Furthermore, we divide interval $(0,1)$ into 13 open subintervals and 12 single points, which are arranged in order as follows by comparing the endpoints of $I_{i}$ $(i=1,2\cdots 11)$

[TABLE]

Consequently, the sign of $D_{i}$ ($i=3\cdots 8$) is displayed as follows

[TABLE]

The sign lists of the discriminant sequence $D$ are given as follows

[TABLE]

For $h\in h_{2}\cup h_{1}\cup h_{6}\cup h_{3}\cup h_{7}\cup h_{8}\cup h_{4}\cup h_{5}$, the revised sign lists are

[TABLE]

Thus, the change number of the sign lists of the discriminant sequence $D$ is 4 and the number of the non-vanishing numbers of these lists is 8 for $h\in(0,1)/(h_{9}\cup\frac{1}{2}\cup h_{10}\cup h_{11})$ and the change number is 3 and the number of the non-vanishing numbers is 6 for $h\in h_{9}\cup\frac{1}{2}\cup h_{10}\cup h_{11}$. By Theorem 2.1 in Yang, [1999], $\tilde{g}_{2}(x)$ has no root on the total real axis. Hence, $g_{2}(k)$ has no real root in the interval $(0,k_{1})$. In addition, $g_{2}(k)$ is an even function and $g_{2}(0)<0$, then we obtain $g_{2}(k)<0$ for $0<k<k_{1}$ and $0<h<1$.

Similarly, by using the complete discrimination system of polynomial we can also obtain $g_{3}(k)<0$ for $0<k<k_{1}$ and $0<h<1$. Hence, we conclude that the Jacobian determinant is nonzero, i.e., (118) is locally invertible. System (116) therefore is locally equivalent to the universal unfolding system

[TABLE]

As indicated in Section 8.4 of Kuznetsov [1995], system (125) undergoes a saddle-node bifurcation as $(\xi_{1},\xi_{2})$ crossing $\mathcal{SN}^{+}\cup\mathcal{SN}^{-}$, where $\mathcal{SN}^{+}:=\{(\xi_{1},\xi_{2})\in U:\xi_{1}=\xi_{2}^{2}/4,\xi_{2}>0\}$ and $\mathcal{SN}^{-}:=\{(\xi_{1},\xi_{2})\in U:\xi_{1}=\xi_{2}^{2}/4,\xi_{2}<0\}$, a Hopf bifurcation as $(\xi_{1},\xi_{2})$ crossing $\mathcal{H}$, where $\mathcal{H}:=\{(\xi_{1},\xi_{2})\in U:\xi_{1}=0,\xi_{2}<0\}$ and a homoclinic bifurcation as $(\xi_{1},\xi_{2})$ crossing $\mathcal{HL}$, where $\mathcal{HL}:=\{(\xi_{1},\xi_{2})\in U:\xi_{1}=-6\xi_{2}^{2}/25+O(|\xi_{2}|^{3}),\xi_{2}<0\}$.

In what follows, we only need to present the bifurcation curve $\mathcal{HL}$ in terms of $\epsilon_{1}$ and $\epsilon_{2}$ because the bifurcation curves $\mathcal{SN}$ and $\mathcal{H}$ have already been shown in Theorem 2.6 and Theorem 3.1. For convenience, we denote $\mu_{lij}:=\partial^{i+j}\mu_{l}(0,0)/\partial^{i}\epsilon_{1}\partial^{j}\epsilon_{2}$, $l=1,2$ and $i,j=0,1,2$ given in the Appendix. We can solve $\epsilon_{1}$ and $\epsilon_{2}$ from (118) as follows

[TABLE]

Before expressing the bifurcation curve we need to prove $\mu_{110}\neq 0$ for $0<h<1$ and $0<k<k_{1}$, where $\mu_{110}$ is given in the Appendix and has the same sign as $\tilde{\mu}_{110}:=\zeta_{1}(k)\sigma_{2}+\zeta_{2}(k)$ with

[TABLE]

It is easy to obtain $\zeta_{2}(k)>0$ for $0<h<1$ and $0<k<k_{1}$. We next prove $\zeta_{1}(k)>0$ for $0<h<1$ and $0<k<k_{1}$. Let $k=\frac{1}{(1-h)(1+x^{2})}$, the number of real roots for $\zeta_{1}(k)$ in the interval $(0,k_{1})$ is equal to the half number of that for $\tilde{\zeta}_{1}(k)$ on the total real axis (Yang, [1999]), where $\tilde{\zeta}_{1}(k):=3(h+1)^{4}x^{6}+(h+9)(h+1)^{3}x^{4}+2(6h^{2}+8h+5)x^{2}+4$. Obviously, $\tilde{\zeta}_{1}(k)$ has no real root, implying that $\zeta_{1}(k)$ has no root in the interval $(0,k_{1})$. In addition, $\zeta_{1}(k)$ is an odd function and $\zeta_{1}(0)>0$, we obtain $\zeta_{1}(k)>0$ for $0<h<1$ and $0<k<k_{1}$.

For the bifurcation curve $\mathcal{HL}$, we consider $\Xi:=\beta_{1}+\frac{6}{25}\beta_{2}^{2}+O(|\beta_{2}|^{3})$. Since $\frac{\partial\Xi}{\partial\epsilon_{1}}=\frac{B^{4}(0,0)}{A^{3}(0,0)}\mu_{110}\neq 0$, by the implicit function theorem, there exists a unique function $\epsilon_{1}(\epsilon_{2})$ such that $\epsilon_{1}(0)=0$ and $\Xi(\epsilon_{1}(\epsilon_{2}),\epsilon_{2})=0$, which can be obtained as an expansion

[TABLE]

Further, on the curve $\Xi=0$, we have $\epsilon_{2}=-\frac{A^{2}(0,0)\mu_{110}}{B^{2}(0,0)(\mu_{101}\mu_{210}-\mu_{110}\mu_{201})}\beta_{2}+O(|\beta_{2}|^{2})$, in which the coefficient of $\beta_{2}$ is negative, implying that $\epsilon_{2}>0$ if $\beta_{2}<0$ and $\epsilon_{2}<0$ if $\beta_{2}>0$. Therefore, we obtain

[TABLE]

With the transformation $\epsilon_{1}=\alpha-\alpha_{1}$ and $\epsilon_{2}=\sigma-\sigma_{2}$, we can rewrite the above bifurcation curve $\mathcal{HL}$ as in Theorem 4.3. The proof of this theorem is completed.

5 Numerical simulation and discussion

In this paper we qualitatively investigate prey-predator system (6) with foraging facilitation among predators including the number and properties of the equilibria (Theorem 2.1) as well as the bifurcations of equilibria such as transcritical and pitchfork bifurcations (Theorem 2.4), saddle-node bifurcation (Theorem 2.6), Hopf bifurcation (Theorem 3.1) and Bogdanov-Takens bifurcation (Theorem 4.3). In spite that both the saddle-node and Hopf bifurcations are discussed above, they are also exhibited in the Bogdanov-Takens bifurcation. As indicated in Theorem 4.3, the neighborhood $U$ of point $(\sigma_{2},\alpha_{*})$ is divided into four regions, i.e., $U=\mathcal{SN}^{+}\cup\mathcal{SN}^{-}\cup\mathcal{H}\cup\mathcal{HL}\cup\mathcal{R}_{1}\cup\mathcal{R}_{2}\cup\mathcal{R}_{3}\cup\mathcal{R}_{4}$, where

[TABLE]

Accordingly, the dynamical behaviors of system (8) near the cusp $E_{*}$ for parameters in neighborhood $U$ of point $(\sigma_{2},\alpha_{*})$ are listed in Table 2.

We next offer some examples to demonstrate the dynamical behaviors of system (8). Let $h=0.5$ and $k=1$, we have $\sigma_{2}=0.55$, $\alpha_{*}=15.94$ and $E_{*}=(0.72,0.11)$. The bifurcation diagram of the Bogdanov-Takens bifurcation is displayed in Fig. 2.

When $(\sigma,\alpha)=(0.62,14.2)\in\mathcal{R}_{1}$, system (8) has no equilibrium except the saddle $E_{0}$ and the stable node $E_{k}$ (Fig. 3 (a)). For $(\sigma,\alpha)=(0.62,14.3)\in\mathcal{R}_{2}$, system (8) has four equilibria, i.e., the stable node $E_{k}$, the stable focus $E_{1}$ and the saddles $E_{0}$ and $E_{2}$ as shown in Fig. 3 (b). The rise of equilibria $E_{1}$ and $E_{2}$ is due to the saddle-node bifurcation. When $(\sigma,\alpha)=(0.62,14.42)\in\mathcal{R}_{3}$, system (8) has a stable limit cycle and four equilibria, i.e., the stable node $E_{k}$, the unstable focus $E_{1}$ and the saddles $E_{0}$ and $E_{2}$ as shown in Fig. 3 (c). The rise of the limit cycle is induced by the Hopf bifurcation. For $(\sigma,\alpha)=(0.62,14.55)\in\mathcal{R}_{4}$, system (8) only has four equilibria, i.e., the stable node $E_{k}$, the unstable focus $E_{1}$ and the saddles $E_{0}$ and $E_{2}$ as shown in Fig. 3 (d). The disappearance of the limit cycle is induced by the homoclinic bifurcation. Theorem 3.1 describes that one stable limit cycle arises near $E_{1}$ induced by the Hopf bifurcation as $0<h<1$ and $(k,\sigma,\alpha)$ varies from $\mathcal{S}_{11}\cup\mathcal{L}_{21}\cup\mathcal{S}_{31}$ to $\mathcal{P}_{11}\cup\mathcal{S}_{21}\cup\mathcal{P}_{31}$. Table 1 shows that system (8) has two saddles $E_{0}$ and $E_{k}$ and an unstable focus or node $E_{1}$ when $0<h<1$ and $(k,\sigma,\alpha)\in\mathcal{P}_{11}\cup\mathcal{S}_{21}\cup\mathcal{P}_{31}$. Thus, the stable limit cycle is the $\omega$-limit set of the positive solutions. For example, let $h=0.5$ and $(k,\sigma,\alpha)=(5.5,1,0.1)\in\mathcal{P}_{11}$, system (8) has two saddles $E_{0}$ and $E_{k}$ and an unstable focus or node $E_{1}$ surrounded by a stable limit cycle as shown in Fig. 4.

From an ecological point of view, the foraging facilitation among predators is an interesting phenomenon to understand the dynamics of the prey-predator interactions in ecosystems, and it is more realistic and reasonable to take into account this factor in the prey-predator system. The qualitative results of system (6) indicate that prey-predator system (6) with foraging facilitation has richer dynamic behaviors than the Rosenzweig-MacArthur system because system (2) only undergoes the transcritical and Hopf bifurcation. The analysis of system (6) reveals that population can be stabilized at the predator free equilibrium or the coexistence equilibrium with increasing the foraging facilitation $\alpha$ as the environmental capacity of prey is relatively low. How the population evolves in time depends on the initial conditions. The foraging facilitation is then beneficial for population persistence and promotes ecosystem diversity. Cooperative predators can survive in a less favorable and less productive environment, in which sufficient preys are available and the survival is more robust for higher levels of cooperation. The bistability of the system implies that the predator population goes extinct for low initial predator densities, which actually is the phenomenon of Allee effect in the predators (Courchamp et al. [2008]). Therefore, the foraging facilitation is a mechanism for inducing Allee effects in predators. For low environmental capacity of prey and weak foraging facilitation, the prey population is too small to sustain the predator population even though the foraging facilitation of predators exists. Nevertheless, the foraging facilitation can have not only positive but also negative effects for predators. For very strong foraging facilitation, the population goes to extinction due to the excessive hunting of prey population by predator population. The destabilization of the system appears due to the Hopf bifurcation even the homoclinic bifurcation that causes splitting of the stable cycle, thus ending the oscillation and consequently causing the extinction of the predators. The overexploitation can therefore backfire and result in the extinction of predators because of the increased predation pressure. It is well known that the Rosenzweig-MacArthur system demonstrates the paradox of enrichment caused by the Hopf bifurcation (Rosenzweig, [1971]), which means that a stable oscillation bifurcates from a stable equilibrium once the environmental carrying capacity of the prey exceeds a critical value. The qualitative results of system (6) reveal that this typical phenomenon of system (2) is inherited even if the foraging facilitation is introduced. We also can observe that the predators will go to extinct if the handing time of the predators $h$ is too long such as $h\geq 1$. By means of bifurcation analysis of prey-predator system (6), we have proved that hunting cooperation is not always beneficial for the predator population. Such studies of bifurcations may give insights into the important changes of dynamical behaviors of the system caused by small perturbation of parameters. The results of bifurcations provide some thresholds to control the qualitative properties of the prey-predator system.

\nonumsection

Acknowledgments The author is grateful to the associate editor and reviewers for their valuable comments and suggestions that have significantly improved the presentation and quality of the manuscript.