Complex dynamics of a predator–prey model with opportunistic predator and weak Allee effect in prey

ABSTRACT

In this work, we first modify a Lotka–Volterra predator–prey system to incorporate an opportunistic predator and weak Allee effect in prey. The prey will be extinct if the combined effect of hunting and other food resources of predator is large. Otherwise, the dynamic behaviour of the system is extremely rich. A series of bifurcations such as saddle-node bifurcation, Hopf bifurcation, and Bogdanov–Takens bifurcation can happen. The validity of the theoretical results are supported with numerical simulations.

KEYWORDS:

• Opportunistic predator
• Allee effect
• predator–prey
• Hopf bifurcation
• Bogdanov–Takens bifurcation

MATHEMATICS SUBJECT CLASSIFICATIONS:

• 34C23
• 34D20
• 92D40

1. Introduction

Predation is an important interaction among species in nature, where each species in a food chain must capture or be captured by other species. Mathematical models described by differential equations have provided insights on understanding the evolution and fitness of species. For example, results on persistence [1,9,19], stability of equilibria [22], and bifurcation [32] have corresponding ecological significance.

Predators can be divided into two types: specialist predator and generalist predator [14–16]. A specialist predator completely depends on the prey while a generalist predator has alternative food resources other than specialized prey. In the former, predators die out in the absence of the prey whereas, in the latter, predators can survive in the absence of the specified prey. What's interesting is that there exists opportunistic predators [6] among specialist predator. Different from traditional specialist predator, this type predator such as Jackal and Martial eagle will snatch the prey that other predators have already captured or steal the leftovers of other predators. Therefore, a simple model on predation with an opportunistic predator is described by the following Lotka-Volterra predator–prey system, (1) $\begin{array}{rl}\dot{x}& =rx(1-\frac{x}{K})-bxy,\\ \dot{y}& =cxy+my-\mathrm{d}{y}^2,\end{array}$(1) where the term my represents the intrinsic growth due to other food resources in the absence of prey. The dynamical behaviour of system (1(1) $\begin{array}{rl}\dot{x}& =rx(1-\frac{x}{K})-bxy,\\ \dot{y}& =cxy+my-\mathrm{d}{y}^2,\end{array}$(1) ) was analysed completely in Chen [8]. It admits at most a positive equilibrium, which is globally asymptotically stable when exists.

Nevertheless, a predator–prey system is affected by many factors. One currently studied is the Allee effect, which can be regarded as a negative density dependence when the population growth rate is decreased at low population size [13]. Allee effect might be caused by difficulties in finding mates, social dysfunction at small population sizes, inbreeding depression [11]. Allee effect is observed not only in animals but also in microorganisms and plants. As the density of the species increases, the influence of Allee effect will be weakened and even vanishes. Allee effect can be divided into two types, strong Allee effect [5,30] and weak Allee effect [7,29]. The strong Allee effect implies that there is a threshold for the population size. If the population size is lower than this size, the per capita growth rate is negative. But the weak Allee effect means that the per capita growth rate decreases but still remains positive.

In recent years, many scholars have extensively studied predator–prey models with both generalist predators and Allee effect [2–4,23–25,27,28]. In [24], Sen et al. studied the impact of Allee effect in a predator–prey model with a generalist predator. The model exhibits simpler dynamics than the corresponding one without Allee effect. It undergoes the same local and global bifurcations. However, the tristability [17] for the one without Allee effect will not occur. Meanwhile, Arancibia-Ibarra and Flores [3] studied a modified Leslie-Gower predator–prey model with Holling type II functional response, strong and weak Allee effects in the prey, and a generalist predator. It is shown that the model with strong Allee effect has at most two positive equilibria, one always being a saddle and the other exhibiting multi-stability phenomenon since the equilibrium can be stable or unstable. The model with weak Allee effect has at most three positive equilibria, where one is always a saddle and the other two can be stable or unstable nodes.

However, predator–prey models with opportunistic predators and Allee effect in prey have been rarely studied. Therefore we modify the following Lotka-Volterra type predator–prey system involving Allee effect studied by Merdan [20], (2) $\begin{array}{rl}\dot{x}& =rx(1-x)\frac{x}{\text{β}+x}-axy,\\ \dot{y}& =ay(x-y).\end{array}$(2) Taking into account other food resources for predators, we modify (2(2) $\begin{array}{rl}\dot{x}& =rx(1-x)\frac{x}{\text{β}+x}-axy,\\ \dot{y}& =ay(x-y).\end{array}$(2) ) as (3) $\begin{array}{rl}\dot{x}& =rx(1-\frac{x}{K})\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-\mathrm{d}{y}^2,\end{array}$(3) where x and y stand for the densities of the prey and predators at time t, respectively. Here r is the intrinsic fertility rate of the prey in the absence of predators; K is the carrying capacity of prey; the term $\frac{x}{A+x}$ describes the Allee effect in prey with A reflecting the intensity of Allee effect; b is the maximal predator per capita capture rate; c is the conversion coefficient; m stands for the intrinsic fertility rate of predators in the absence of prey; d is the intraspecific competition coefficient of predators. All r, K, A, b, c, m, and d are positive constants.

For simplicity of forthcoming calculations and to reduce the number of parameters, we introduce the dimensionless changes of variables $\overline{x}=Kx,\overline{y}=\frac{r}{d}y,\overline{t}=rt$and denote $\overline{A}=\frac{A}{K},\overline{b}=\frac{b}{dr},\overline{c}=\frac{cK}{r},\overline{m}=\frac{m}{d}.$After dropping the bars, system (3(3) $\begin{array}{rl}\dot{x}& =rx(1-\frac{x}{K})\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-\mathrm{d}{y}^2,\end{array}$(3) ) becomes (4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) From the biological point of view, the initial conditions of (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) are $\left(x\right(0),y(0\left)\right)\in {\mathbb{R}}_{+}^2$. For each such initial condition, it is easy to see that (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) has a unique global solution, which stays in ${\mathbb{R}}_{+}^2$. Moreover, it follows from the first equation of (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) that $\dot{x}\le x(1-x),$which implies that $\underset{x\to \mathrm{\infty }}{lim sup}x\left(t\right)\le 1$. Then, for any $ϵ>0$, there exists $T\ge 0$ such that $x\left(t\right)\le 1+ϵ$ for $t\ge T$. This, together with the second equation of (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ), implies that $\dot{y}\left(t\right)\le \left(c\right(1+ϵ)+m)y-y^2\text{for}t\ge T,$which yields $\underset{t\to \mathrm{\infty }}{lim sup}y\left(t\right)\le c(1+ϵ)+m$. Since ε is arbitrary, we immediately get $\underset{t\to \mathrm{\infty }}{lim sup}y\left(t\right)\le c+m$ by letting $ϵ\to 0$. Thus we have shown that every solution is bounded and, furthermore, the above differential inequalities also produce the attracting and positively invariant set $\left\{\right(x,y)\in {\mathbb{R}}_{+}^2:0\le x\le 1,0\le y\le c+m\}.$The goal of this paper is to analyse the dynamics of (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ). We first discuss the existence and stability of equilibria in Section 2. The results indicate the occurrence of various bifurcations. Thus, in Section 3, we provide details on saddle-node bifurcation, Hopf bifurcation, and Bogdanov–Takens bifurcation. The paper concludes with some discussion.

2. Existence and stability of equilibria

2.1. Existence of equilibria

We start with the existence of equilibria of (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ).

System (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) always has the three boundary equilibria $E_0(0,0)$, $E_1(1,0)$, and $E_2(0,m)$. Other equilibria are positive ones, which are determined by $\begin{array}{rl}& (1-x)\frac{x}{A+x}-by=0,\\ & y=cx+m.\end{array}$Thus x is a positive zero of $g\left(x\right)$ in $(0,1)$, where (5) $g\left(x\right)=(bc+1)x^2+(Abc+bm-1)x+Abm.$(5) Note that g is quadratic in x with $g\left(0\right)=Abm>0$. Therefore, it is necessary that 1−bm>0 and $A<A_0=\frac{1-bm}{bc}$ if $g\left(x\right)=0$ has positive roots. The discriminant of $g\left(x\right)=0$ is $\mathrm{\Delta }=b^2{c}^2{A}^2-2(b^{2}cm+bc+2bm)A+(bm-1{)}^2.$Regarded as a quadratic function of A, Δ has two zeros $A_{\pm }=\frac{cmb+c+2m\pm 2\sqrt{m(c+m)(bc+1)}}{b{c}^2}.$It is easy to check that $A_{-}<A_0<A_{+}$. Thus- when bm<1,

1. if $A>A_{-}$ then $g\left(x\right)>0$ for x>0;

2. if $0<A<A_{-}$ then $g\left(x\right)$ has two distinct positive roots, (6) $x_3=\frac{1-Abc-bm-\sqrt{\mathrm{\Delta }}}{2(bc+1)}\text{and}{x}_4=\frac{1-Abc-bm+\sqrt{\mathrm{\Delta }}}{2(bc+1)};$(6)

3. if $A=A_{-}$ then $g\left(x\right)=0$ has a unique positive real root (7) $x_{\ast }=\frac{1-Abc-bm}{2(bc+1)}.$(7)

Note that $x_3<x_{\ast }<x_4<1$.

Summarizing the above discussion, we have obtained the following result on equilibria of (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ).

Theorem 2.1

The following statements on equilibria of (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) hold.

1. System (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) always has the boundary equilibria $E_0(0,0)$, $E_1(1,0)$, and $E_2(0,m)$.

2. System (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) can have positive equilibria only if bm<1. In this case,

   1. if $A<A_{-}$ then there are two positive equilibria $E_3(x_3,y_3)$ and $E_4(x_4,y_4)$, where $x_3$ and $x_4$ are given by (6(6) $x_3=\frac{1-Abc-bm-\sqrt{\mathrm{\Delta }}}{2(bc+1)}\text{and}{x}_4=\frac{1-Abc-bm+\sqrt{\mathrm{\Delta }}}{2(bc+1)};$(6) ), $y_{3,4}=c{x}_{3,4}+m$;

   2. if $A=A_{-}$ then there is a unique positive equilibrium $E_{\ast }(x_{\ast },y_{\ast })$, where $x_{\ast }$ is given by (7(7) $x_{\ast }=\frac{1-Abc-bm}{2(bc+1)}.$(7) ) and $y_{\ast }=c{x}_{\ast }+m$;

   3. if $A>A_{-}$ then there is no positive equilibrium.

In the following, we analyse the stability of the equilibria stated in Theorem 2.1.

2.2. Stability of the boundary equilibria $E_0$, $E_1$, and $E_2$

We first consider the boundary equilibria.

Theorem 2.2

The three boundary equilibria, $E_0$, $E_1$, and $E_2$ of (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) are repelling saddle-node, saddle, and stable node, respectively.

Proof.

First, the Jacobian matrix of (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) at $E_0(0,0)$ is $J_{E_0}=\left(\begin{array}{cc}0& 0\\ 0& m\end{array}\right)$. As 0 is an eigenvalue of $J_{E_0}$, we make a Taylor expansion at $(0,0)$ and the time scale $\mathrm{d}t=m\mathrm{d}\tau$ to change (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) into $\begin{array}{rl}\dot{x}& =\frac{1}{Am}{x}^2-\frac{b}{m}xy+O(|x,y{|}^3),\\ \dot{y}& =y+\frac{c}{m}xy-\frac{1}{m}{y}^2,\end{array}$(here we still use dot to denote the derivative with respect to τ). Since the coefficient of $x^2$ is $\frac{1}{Am}>0$, by Zhang et al. [31, Theorem 7.1 in Chapter 2], $E_0$ is a saddle-node and the repelling parabolic is located at the right half plane.

Next, the Jacobian matrices at $E_1(1,0)$ and $E_2(0,m)$ are respectively $J_{E_1}=\left(\begin{array}{cc}-{\displaystyle \frac{1}{A+1}}& -b\\ 0& c+m\end{array}\right)\text{and}{J}_{E_2}=\left(\begin{array}{cc}-bm& 0\\ cm& -m\end{array}\right).$It follows immediately that $E_1$ is a saddle and $E_2$ is a stable node.

Note that the x-axis and the y-axis are positively invariant. By Theorem 2.1, if $bm\ge 1$ or (bm<1 and $A>A_{-}$) then (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) has no positive equilibrium and hence there is no closed orbit in the first quadrant. At the meantime, both $E_0$ and $E_1$ are unstable. This gives the following result.

Corollary 2.3

If $bm\ge 1$ or (bm<1 and $A>A_{-}$), then the equilibrium $E_2$ is globally asymptotically stable (see Figure 1).

Figure 1. When $A=0.7$, b=2, c=0.25, and m=0.05, the equilibrium $E_2$ of (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) is globally asymptotically stable.

2.3. Types of the positive equilibria $E_3$ and $E_4$

Next, we discuss the types of the positive equilibria $E_3$ and $E_4$. For this purpose, denote $m^{\ast }=-\frac{x_4(c{A}^2+2c{x}_{4}A+c{x_4}^2+2x_{4}A+{x_4}^2-A)}{(A+x_4{)}^2}>0.$

Theorem 2.4

Suppose bm<1 and $0<A<A_{-}$. Then the equilibrium $E_3$ of (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) is always a saddle and $E_4\text{is a}\{\begin{array}{ll}\text{source}& \text{if}m<m^{\ast },\\ \text{sink}& \text{if}m>m^{\ast },\\ \text{centre or fine focus}& \text{if}m=m^{\ast }.\end{array}$

Proof.

The Jacobian matrix at any positive equilibrium $E(x,y)$ of (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) is given by $J_E=\left(\begin{array}{cc}-{\displaystyle \frac{x(2Ax+x^2-A)}{(A+x{)}^2}}& -bx\\ c(cx+m)& -cx-m\end{array}\right).$Thus $det\left(J_E\right)=\frac{x(cx+m)\left[\right(bc+1)x^2+2A(bc+1)x+bc{A}^2-A]}{(A+x{)}^2}$and $\mathrm{t}\mathrm{r}\left(J_E\right)=\frac{x(1-2x)}{A+x}-\frac{x^2(1-x)}{(A+x{)}^2}-cx-m.$For g defined by (5(5) $g\left(x\right)=(bc+1)x^2+(Abc+bm-1)x+Abm.$(5) ), we have $g^{\prime }\left(x\right)=2(bc+1)x+Abc+bm-1.$As $g\left(x\right)=0$, we have $c=-\frac{x^2+(bm-1)x+bmA}{bx(A+x)}.$Then $det\left(J_E\right)$ and $g^{\prime }\left(x\right)$ can be simplified as $\begin{array}{rl}det\left(J_E\right)& =\frac{(bm{A}^2+2bmAx+(bm-A-1\left)x^2\right)(x-1)x}{(A+x{)}^{3}b},\\ g^{\prime }\left(x\right)& =-\frac{bm{A}^2+2bmAx+(bm-A-1)x^2}{x(A+x)},\end{array}$which give (8) $det\left(J_E\right)=\frac{(1-x)x^2}{(A+x{)}^{2}b}{g}^{\prime }\left(x\right).$(8) Recall 0<x<1. Therefore, the positive equilibrium $E(x,y)$ is an elementary equilibrium if $g^{\prime }\left(x\right)\ne 0$, a hyperbolic saddle if $g^{\prime }\left(x\right)<0$, a degenerate equilibrium if $g^{\prime }\left(x\right)=0$. As $g^{\prime }\left(x_3\right)<0$ and $g^{\prime }\left(x_4\right)>0$, we know that $E_3$ is a saddle and $E_4$ is an anti-saddle [18]. Furthermore, $\mathrm{t}\mathrm{r}\left(J_{E_4}\right)=\frac{x_4(1-2x_4)}{A+x_4}-\frac{{x_4}^2(1-x_4)}{(A+x_4{)}^2}-c{x}_4-m.$It is easy to see that $\mathrm{t}\mathrm{r}\left(J_{E_4}\right)>0$, $\mathrm{t}\mathrm{r}\left(J_{E_4}\right)=0$, and $\mathrm{t}\mathrm{r}\left(J_{E_4}\right)<0$ if $m<m^{\ast }$, $m=m^{\ast }$, and $m>m^{\ast }$ respectively. Then the type of $E_4$ follows and this completes the proof.

Remark 2.1

Suppose bm<1 and $A<A_{-}$. If $m>m^{\ast }$ then, by Theorems 2.2 and 2.4, $E_2$ is a stable node, $E_3$ is a saddle, and $E_4$ is a sink. Thus (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) admits bistable phenomenon and the stable manifold of the saddle $E_3$ is the separatrix (see Figure 2).

Figure 2. When A=0.5, b=2, c=0.25, and m=0.05, system (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) admits bistable phenomenon.

2.4. Stability of the positive equilibrium $E_{\ast }$

Finally, we consider the positive equilibrium $E_{\ast }$.

Theorem 2.5

Suppose that bm<1 and $A=A_{-}$. Then the unique positive equilibrium $E_{\ast }$ of (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) is $\{\begin{array}{ll}\text{an attracting saddle-node}& \text{if}bc{x}_{\ast }-c{x}_{\ast }-m>0;\\ \text{a repelling saddle-node}& \text{if}bc{x}_{\ast }-c{x}_{\ast }-m<0;\\ \text{a cusp of codimension two}& \text{if}c=c_{\ast },m=m_{\ast },b>1,\mathrm{a}\mathrm{n}\mathrm{d}{x}_{\ast }<\frac{b-1}{2b-1}.\end{array}$Here $c_{\ast }=\frac{(1-2b)x_{\ast }+b-1}{b^3{x}_{\ast }}$ and $m_{\ast }=\frac{(1-b)\left(\right(2b-1)x_{\ast }+1-b)}{b^3}$.

Proof.

The Jacobian matrix of system (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) evaluated at $E_{\ast }$ is given by $J_{E_{\ast }}=\left(\begin{array}{cc}-{\displaystyle \frac{x_{\ast }(2A_{-}{x}_{\ast }+x_{\ast }^2-A_{-})}{(A_{-}+x_{\ast }{)}^2}}& -b{x}_{\ast }\\ c(c{x}_{\ast }+m)& -c{x}_{\ast }-m\end{array}\right).$Note that (8(8) $det\left(J_E\right)=\frac{(1-x)x^2}{(A+x{)}^{2}b}{g}^{\prime }\left(x\right).$(8) ) still holds for $E_{\ast }$, that is, $det\left(J_{E_{\ast }}\right)=\frac{(1-x_{\ast })x_{\ast }^2}{(A_{-}+x_{\ast }{)}^{2}b}{g}^{\prime }\left(x_{\ast }\right)$. As $g^{\prime }\left(x_{\ast }\right)=0$, it follows that $\frac{x_{\ast }(2A_{-}{x}_{\ast }+x_{\ast }^2-A_{-})}{(A_{-}+x_{\ast }{)}^2}=-bc{x}_{\ast }$. We perform the transformation $(x,y)=(X+x_{\ast },Y+y_{\ast })$to transform the equilibrium $E_{\ast }$ to the origin and (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) into (9) $\begin{array}{rl}\dot{X}& =a_{10}X+a_{01}Y+a_{20}{X}^2+a_{11}XY+a_{02}{Y}^2+a_{30}{X}^3\\ & +a_{21}{X}^{2}Y+a_{12}X{Y}^2+a_{03}{Y}^3+P_3(X,Y),\\ \dot{Y}& =b_{10}X+b_{01}Y+b_{20}{X}^2+b_{11}XY+b_{02}{Y}^2+b_{30}{X}^3\\ & +b_{21}{X}^{2}Y+b_{12}X{Y}^2+b_{03}{Y}^3+Q_3(X,Y),\end{array}$(9) where $\begin{array}{rl}{a}_{10}& =bc{x}_{\ast },a_{01}=-b{x}_{\ast },a_{20}=-\frac{x_{\ast }^3+3A_{-}{x}_{\ast }^2+3A_{-}^2{x}_{\ast }-A_{-}^2}{A_{-}+x_{\ast }},\\ a_{11}& =-b,a_{30}=-\frac{A_{-}^2(A_{-}+1)}{(A_{-}+x_{\ast }{)}^4},b_{10}=c(c{x}_{\ast }+m),\\ b_{01}& =-(c{x}_{\ast }+m),b_{11}=c,b_{02}=-1,\\ b_{20}& =a_{02}=a_{21}=a_{12}=a_{03}=0,\end{array}$and $P_3(X,Y)$ and $Q_3(X,Y)$ are $C^{\mathrm{\infty }}$ functions of at least order four in terms of $(X,Y)$.

First, we assume that $a_{10}+b_{01}=bc{x}_{\ast }-(c{x}_{\ast }+m)\ne 0$. Then the eigenvalues of $J_{E_{\ast }}$ are ${\lambda }_1=0$ and ${\lambda }_2=a_{10}+b_{01}\ne 0$. Employing the transformation $(X,Y)=(u+v,cu+\frac{b_{01}}{a_{01}}v)\text{and}\mathrm{d}\tau =(a_{10}+b_{01})\mathrm{d}t,$we can rewrite (9(9) $\begin{array}{rl}\dot{X}& =a_{10}X+a_{01}Y+a_{20}{X}^2+a_{11}XY+a_{02}{Y}^2+a_{30}{X}^3\\ & +a_{21}{X}^{2}Y+a_{12}X{Y}^2+a_{03}{Y}^3+P_3(X,Y),\\ \dot{Y}& =b_{10}X+b_{01}Y+b_{20}{X}^2+b_{11}XY+b_{02}{Y}^2+b_{30}{X}^3\\ & +b_{21}{X}^{2}Y+b_{12}X{Y}^2+b_{03}{Y}^3+Q_3(X,Y),\end{array}$(9) ) as $\begin{array}{rl}\dot{u}& ={\overline{c}}_{20}{u}^2+{\overline{c}}_{11}uv+{\overline{c}}_{02}{v}^2+{\overline{P}}_4(u,v),\\ \dot{v}& =v+{\overline{d}}_{20}{u}^2+{\overline{d}}_{11}uv+{\overline{d}}_{02}{v}^2+{\overline{Q}}_4(u,v),\end{array}$where $\begin{array}{rl}{\overline{c}}_{20}& =-\frac{(c{x}_{\ast }+m)(A_{-}bc-bm-1)A_{-}}{2(A_{-}+x_{\ast }{)}^2((bc-c)x_{\ast }-m{)}^2},\\ {\overline{c}}_{11}& =\frac{1}{x_{\ast }c(A_{-}+x_{\ast }{)}^3((bc-c)x_{\ast }-m)}\left[\right(4b{A}_{-}{c}^3-b{c}^{2}m+6A_{-}{c}^2+2c^2+cm)x_{\ast }^4\\ & +(6b{A}_{-}^2{c}^3+b{A}_{-}{c}^{2}m+6A_{-}^2{c}^2+3A_{-}cm-m^2)x_{\ast }^3\\ & +(2b{A}_{-}^3{c}^3+3b{A}_{-}^2{c}^{2}m-2A_{-}^2{c}^2+3A_{-}^{2}cm-3A_{-}{m}^2)x_{\ast }^2\\ & +(b{A}_{-}^3{c}^{2}m-A_{-}^{3}cm-2A_{-}^{2}cm-3A_{-}^2{m}^2)x_{\ast }-A_{-}^3{m}^2],\\ {\overline{c}}_{02}& =\frac{c{x}_{\ast }+m}{b{x}_{\ast }{c}^2(A_{-}+x_{\ast }{)}^3((bc-c)x_{\ast }-m)}\left[\right(2b{A}_{-}{c}^2+3A_{-}c+c+m)x_{\ast }^3\\ & +(3b{A}_{-}^2{c}^2+2b{A}_{-}cm+3A_{-}^{2}c+3A_{-}m)x_{\ast }^2\\ & +(2b{A}_{-}^3{c}^2+b{A}_{-}^2{c}^2+3b{A}_{-}^{2}cm+A_{-}^{3}c+3A_{-}^{2}m)x_{\ast }+b{A}_{-}^{3}cm+A_{-}^{3}m],\\ {\overline{d}}_{20}& =-\frac{cb{x}_{\ast }}{(A_{-}+x_{\ast }{)}^3((bc-c)x_{\ast }-m)}\left[\right(bc+1)x_{\ast }^3+(3b{A}_{-}c+3A_{-})x_{\ast }^2\\ & +(3b{A}_{-}^{2}c+3A_{-}^2)x_{\ast }+b{A}_{-}^{3}c-A_{-}^2],\\ {\overline{d}}_{11}& =-\frac{1}{(A_{-}+x_{\ast }{)}^3((bc-c)x_{\ast }-m)}\left[\right(b^2{c}^2+b{c}^2+bc+c)x_{\ast }^4\\ & +(3b^2{A}_{-}{c}^2+3b{A}_{-}{c}^2+3b{A}_{-}c-bm+3A_{-}c+m)x_{\ast }^3\\ & +(3b^2{A}_{-}^2{c}^2+3b{A}_{-}^2{c}^2+3b{A}_{-}^{2}c-3bm{A}_{-}+3A_{-}^{2}c+3A_{-}m)x_{\ast }^2\\ & +(b^2{A}_{-}^3{c}^2+b{A}_{-}^3{c}^2-b{A}_{-}^{3}c-2b{A}_{-}^{2}c-3b{A}_{-}^{2}m+A_{-}^{3}c+3A_{-}^{2}m)x_{\ast }\\ & -b{A}_{-}^{3}m+A_{-}^{3}m],\\ {\overline{d}}_{02}& =\frac{1}{b{x}_{\ast }{c}^2(A_{-}+x_{\ast }{)}^3((bc-c)x_{\ast }-m)}\left[\right(b^2{c}^{2}m+2b{A}_{-}{c}^3+3A_{-}{c}^2+c^2+2cm)x_{\ast }^4\\ & +(3b^2{A}_{-}{c}^{2}m+3b{A}_{-}^2{c}^3+2b{A}_{-}{c}^{2}m+3A_{-}^2{c}^2+6A_{-}cm+m^2)x_{\ast }^3\\ & +(b^2{A}_{-}^3{c}^3+b^2{A}_{-}^2{c}^3+3b^2{A}_{-}^2{c}^{2}m+b{A}_{-}^3{c}^3+3b{A}_{-}^2{c}^{2}m+A_{-}^3{c}^2+6A_{-}^{2}cm\\ & +3A_{-}{m}^2)x_{\ast }^2+(b^2{A}_{-}^3{c}^{2}m+b{A}_{-}^3{c}^{2}m+2A_{-}^{3}cm+3A_{-}^2{m}^2)x_{\ast }+A_{-}^3{m}^2],\end{array}$and $\overline{P_4}(u,v)$, $\overline{Q_4}(u,v)$ are $C^{\mathrm{\infty }}$ functions of at least order three in terms of $(u,v)$. According to the centre manifold method [21] and noticing that ${\overline{c}}_{20}>0$, we can approximate the equation near the centre manifold by $\frac{\text{d}u}{\text{d}t}=-\frac{(c{x}_{\ast }+m)(A_{-}bc-bm-1)A_{-}}{2(A_{-}+x_{\ast }{)}^2((bc-c)x_{\ast }-m{)}^2}{u}^2+O\left(|u{|}^3\right).$Thus the degenerate equilibrium $E_{\ast }$ is a saddle-node by Theorem 7.1 in Zhang et al. [31].

Next, assume that $(a_{10}+b_{01})=0$. Recall that $x_{\ast }=\frac{1-A_{-}bc-bm}{2(bc+1)}\text{and}{A}_{-}=\frac{cmb+c+2m-2\sqrt{m(c+m)(bc+1)}}{b{c}^2}.$Solving c, m, and $A_{-}$ from these relations gives us $c=c_{\ast },m=m_{\ast },A_{-}=A_{-}^{\ast }=\frac{(b-1)x_{\ast }^2+x_{\ast }}{(1-2b)x_{\ast }+b-1}.$As all $c_{\ast }$, $m_{\ast }$, and $A_{-}^{\ast }$ are positive, we have (10) $b>1\text{and}0<x_{\ast }<\frac{b-1}{2b-1}.$(10) Since $a_{01}\ne 0$, we use the linear transformation $X=u,Y=-\frac{a_{10}}{a_{01}}u+\frac{1}{a_{01}}v$to transform (9(9) $\begin{array}{rl}\dot{X}& =a_{10}X+a_{01}Y+a_{20}{X}^2+a_{11}XY+a_{02}{Y}^2+a_{30}{X}^3\\ & +a_{21}{X}^{2}Y+a_{12}X{Y}^2+a_{03}{Y}^3+P_3(X,Y),\\ \dot{Y}& =b_{10}X+b_{01}Y+b_{20}{X}^2+b_{11}XY+b_{02}{Y}^2+b_{30}{X}^3\\ & +b_{21}{X}^{2}Y+b_{12}X{Y}^2+b_{03}{Y}^3+Q_3(X,Y),\end{array}$(9) ) into (11) $\begin{array}{rl}\dot{u}& =v+c_{20}{u}^2-c_{11}uv+O(|u,v{|}^3),\\ \dot{v}& =d_{20}{u}^2+d_{11}uv-d_{02}{v}^2,\end{array}$(11) where $\begin{array}{rl}{c}_{20}& =-\frac{(b-1)(b{x}_{\ast }-x_{\ast }+1)(2b{x}_{\ast }-b-x_{\ast }+1)}{b^3{x}_{\ast }(x_{\ast }-1)},\\ c_{11}& =x_{\ast },\\ d_{20}& =\frac{(b-1)(b{x}_{\ast }-x_{\ast }+1)(2b{x}_{\ast }-b-x_{\ast }+1{)}^2}{x_{\ast }(x_{\ast }-1)b^5},\\ d_{11}& =-\frac{(b-1)\left(\right(2b-1)x_{\ast }+1-b)}{b^3{x}_{\ast }},\\ d_{02}& =\frac{1}{b{x}_{\ast }}.\end{array}$By Remark 1 of Section 2.13 in Perko [21], we obtain an equivalent system of system (11(11) $\begin{array}{rl}\dot{u}& =v+c_{20}{u}^2-c_{11}uv+O(|u,v{|}^3),\\ \dot{v}& =d_{20}{u}^2+d_{11}uv-d_{02}{v}^2,\end{array}$(11) ) in a small neighbourhood of $(0,0)$, $\begin{array}{rl}\dot{u}& =v+O(|u,v{|}^3),\\ \dot{v}& =e_{20}{u}^2+e_{11}uv+O(|u,v{|}^3),\end{array}$where $\begin{array}{rl}{e}_{20}& =d_{20}=\frac{(b-1)(b{x}_{\ast }-x_{\ast }+1)(2b{x}_{\ast }-b-x_{\ast }+1{)}^2}{x_{\ast }(x_{\ast }-1)b^5},\\ e_{11}& =2c_{20}+d_{11}=\frac{(b-1)(2b{x}_{\ast }-x_{\ast }+1)(2b{x}_{\ast }-b-x_{\ast }+1)}{b^3{x}_{\ast }(1-x_{\ast })}.\end{array}$Obviously, $e_{20}<0$ and $e_{11}<0$. Therefore, the equilibrium $E_{\ast }$ is a cusp of codimension 2.

3. Bifurcation analysis

This section is devoted to the analysis of possible bifurcations of (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ), which include saddle-node bifurcation, Hopf bifurcation, and Bogdanov–Takens bifurcation.

3.1. Saddle-node bifurcation

Theorem 2.1 tells us that when bm<1, the number of equilibria of (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) depends on the Allee effect parameter A. Due to Theorem 2.5, the degenerate equilibrium $E_{\ast }$ is a saddle-node if $bc{x}_{\ast }-c{x}_{\ast }-m\ne 0$. Therefore, it is possible for (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) to undergo a saddle-node bifurcation at $A=A_{SN}=A_{-}$, which is confirmed in the following theorem.

Theorem 3.1

System (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) undergoes a saddle-node bifurcation at $E_{\ast }$ when $A=A_{SN}$ and $bc{x}_{\ast }-(c{x}_{\ast }+m)\ne 0$.

Proof.

We apply the Sotomayor's theorem [21] to finish the proof. We have seen before that $J_{E_{\ast }}$ has zero as a simple eigenvalue if $A=A_{SN}$ and $bc{x}_{\ast }-(c{x}_{\ast }+m)\ne 0$. Choose $V=\left(\begin{array}{c}{V}_1\\ V_2\end{array}\right)=\left(\begin{array}{c}1\\ c\end{array}\right)\text{and}W=\left(\begin{array}{c}{W}_1\\ W_2\end{array}\right)=\left(\begin{array}{c}1\\ \frac{b{x}_{\ast }}{c{x}_{\ast }+m}\end{array}\right),$which are - eigenvectors of $J_{E_{\ast }}$ and $J_{E_{\ast }}^T$ corresponding to the eigenvalue 0, respectively. Then $F_A(E_{\ast },A_{SN})=\left(\begin{array}{c}-{\displaystyle \frac{(A_{SN}c+c+m)b(A_{SN}bc+bm-1{)}^2}{(A_{SN}+1-bm{)}^2(bc+1)}}\\ 0\end{array}\right)$and $\begin{array}{rl}{D}^{2}F(E_{\ast },A_{SN})(V,V)& ={\left(\begin{array}{c}{\displaystyle \frac{{\mathrm{\partial }}^2{F}_1}{\mathrm{\partial }{x}^2}}{V}_1^2+2{\displaystyle \frac{{\mathrm{\partial }}^2{F}_1}{\mathrm{\partial }x\mathrm{\partial }y}}{V}_1{V}_2+{\displaystyle \frac{{\mathrm{\partial }}^2{F}_1}{\mathrm{\partial }{y}^2}}{V}_2^2\\ {\displaystyle \frac{{\mathrm{\partial }}^2{F}_2}{\mathrm{\partial }{x}^2}}{V}_1^2+2{\displaystyle \frac{{\mathrm{\partial }}^2{F}_2}{\mathrm{\partial }x\mathrm{\partial }y}}{V}_1{V}_2+{\displaystyle \frac{{\mathrm{\partial }}^2{F}_2}{\mathrm{\partial }{y}^2}}{V}_2^2\end{array}\right)}_{(E_{\ast };A_{SN})}\\ & =\left(\begin{array}{c}{\displaystyle \frac{A_{SN}(A_{SN}bc-bm-1)}{(A_{SN}+x_{\ast }{)}^2}}\\ 0\end{array}\right).\end{array}$Obviously, V and W satisfy the transversality conditions on the occurrence of saddle-node bifurcation at $E_{\ast }$ when $A=A_{SN}$, that is, $\begin{array}{rl}{W}^T{F}_A(E_{\ast },A_{SN})& =-\frac{(A_{SN}c+c+m)b(A_{SN}bc+bm-1{)}^2}{(A_{SN}+1-bm{)}^2(bc+1)}\ne 0,\\ W^T\left[D^{2}F\right(E_{\ast },A_{SN}\left)\right(V,V\left)\right]& =\frac{A_{SN}(A_{SN}bc-bm-1)}{(A_{SN}+x_{\ast }{)}^2}\ne 0.\end{array}$

When the parameter A varies to cross the surface from one side to the other side, the number of positive equilibria of (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) changes from zero to two and hence the saddle-node bifurcation yields two equilibria. The phase of the saddle-node bifurcation, obviously, corresponds with the ecological system having critical Allee value $A_{SN}$. When $A>A_{SN}$, the prey will face the risk of extinction; When $A<A_{SN}$, the dynamic behaviour of (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) becomes very complex since the saddle and anti-saddle come out in the phase, which means that the density of prey must be large enough for its survival.

3.2. Hopf-bifurcation

Next, we consider Hopf bifurcation. By Theorem 2.4, $E_3$ is always a saddle and $E_4$ can be stable or unstable. If $\mathrm{t}\mathrm{r}\left(J_{E_4}\right)=0$, then the eigenvalues of $E_4$ are a pair of purely imaginary complex numbers since $det\left(J_{E_4}\right)>0$. Thus $E_4$ becomes a non-hyperbolic equilibrium and a Hopf-bifurcation may occur around $E_4(x_4,y_4)$. This is what we are going to do in the following.

For simplicity, inspired by Dai et al. [12], we first do the state variable scaling, $\tilde{x}=\frac{x}{x_4},\tilde{y}=\frac{y}{y_4},\tau =y_{4}t,$and denote $\tilde{A}=\frac{A}{x_4},\tilde{p}=\frac{1}{y_4},\tilde{q}=\frac{x_4}{y_4},\tilde{b}=b,\tilde{c}=\frac{c{x}_4}{y_4}\tilde{m}=\frac{m}{y_4}.$After dropping the tildes, system (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) becomes (we still denote τ as t) (12) $\begin{array}{rl}\dot{x}& =x(p-qx){\displaystyle \frac{x}{A+x}}-bxy,\\ \dot{y}& =cxy+my-y^2,\end{array}$(12) where p, q, A, b, c, and m are all positive constants and p>q.

Since system (12(12) $\begin{array}{rl}\dot{x}& =x(p-qx){\displaystyle \frac{x}{A+x}}-bxy,\\ \dot{y}& =cxy+my-y^2,\end{array}$(12) ) has $\widetilde{E_4}(1,1)$ as an equilibrium, we have $b=\frac{p-q}{A+1}>0$ and m=1−c. Substituting them into (12(12) $\begin{array}{rl}\dot{x}& =x(p-qx){\displaystyle \frac{x}{A+x}}-bxy,\\ \dot{y}& =cxy+my-y^2,\end{array}$(12) ) and taking $\mathrm{d}t=\frac{1}{A+x}\mathrm{d}\tau$ (still denote τ as t), we get (13) $\begin{array}{rl}\dot{x}& =x^2(p-qx)-{\displaystyle \frac{p-q}{A+1}}x(A+x)y,\\ \dot{y}& =[cxy+(1-c)y-y^2](A+x),\end{array}$(13) where p>q and c<1. Since (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) and (13(13) $\begin{array}{rl}\dot{x}& =x^2(p-qx)-{\displaystyle \frac{p-q}{A+1}}x(A+x)y,\\ \dot{y}& =[cxy+(1-c)y-y^2](A+x),\end{array}$(13) ) have the same topological structures, we can investigate the Hopf bifurcation around $\widetilde{E_4}(1,1)$ in (13(13) $\begin{array}{rl}\dot{x}& =x^2(p-qx)-{\displaystyle \frac{p-q}{A+1}}x(A+x)y,\\ \dot{y}& =[cxy+(1-c)y-y^2](A+x),\end{array}$(13) ) to get the information on Hopf bifurcation around $E_4(x_4,y_4)$ in (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ).

The Jacobian matrix of (13(13) $\begin{array}{rl}\dot{x}& =x^2(p-qx)-{\displaystyle \frac{p-q}{A+1}}x(A+x)y,\\ \dot{y}& =[cxy+(1-c)y-y^2](A+x),\end{array}$(13) ) at $\widetilde{E_4}$ is $J_{\widetilde{E_4}}=\left(\begin{array}{cc}{\displaystyle \frac{pA-2qA-q}{A+1}}& -p+q\\ c(A+1)& -A-1\end{array}\right).$Then $det\left(J_{\widetilde{E_4}}\right)=(A+1)(p-q)c-pA+2qA+q$and $\mathrm{t}\mathrm{r}\left(J_{\widetilde{E_4}}\right)=-\frac{A^2-pA+2qA+2A+q+1}{A+1}.$It follows that $det\left(J_{\widetilde{E_4}}\right)>0$ for $c>max\{\frac{A(p-q)-q(A+1)}{A(p-q)+(p-q)},0\}\triangleq c^{\ast }$ and $\mathrm{t}\mathrm{r}\left(J_{\widetilde{E_4}}\right)\{\begin{array}{ll}=0& \text{if}p=p^{\ast },\\ >0& \text{if}p>p^{\ast },\\ <0& \text{if}p<p^{\ast },\end{array}$where $p^{\ast }=\frac{A^2+2qA+2A+q+1}{A}$.

Theorem 3.2

Suppose that p>q and $c^{\ast }<c<1$. Then the equilibrium $\widetilde{E_4}(1,1)$ of (13(13) $\begin{array}{rl}\dot{x}& =x^2(p-qx)-{\displaystyle \frac{p-q}{A+1}}x(A+x)y,\\ \dot{y}& =[cxy+(1-c)y-y^2](A+x),\end{array}$(13) ) is $\{\begin{array}{ll}\text{a stable hyperbolic node or focus }& \text{if}p<p^{\ast };\\ \text{an unstable hyperbolic node or focus }& \text{if}p>p^{\ast };\\ \text{a fine focus or centre}& \text{if}p=p^{\ast }.\end{array}$

Next, we study Hopf bifurcation around $\widetilde{E_4}$ in (13(13) $\begin{array}{rl}\dot{x}& =x^2(p-qx)-{\displaystyle \frac{p-q}{A+1}}x(A+x)y,\\ \dot{y}& =[cxy+(1-c)y-y^2](A+x),\end{array}$(13) ). Obviously, ${\frac{\mathrm{d}\mathrm{t}\mathrm{r}\left(J_{\widetilde{E_4}}\right)}{\mathrm{d}p}|}_{p=p^{\ast }}=\frac{A}{A+1}\ne 0,$namely, the transversality condition on the occurrence of Hopf bifurcation is satisfied. Through calculating the first Lyapunov coefficient, we can determine the stability of the limit cycles around the equilibrium $\widetilde{E_4}$.

We first translate $\widetilde{E_4}(1,1)$ to $(0,0)$ by $(\hat{x},\hat{y})=(x-1,y-1)$ and then use the linear transformation $\hat{u}=-\hat{x},\hat{v}=\frac{1}{\sqrt{D}}(\frac{p^{\ast }A-2qA-q}{A+1}\hat{x}+(q-p\left)\hat{y}\right)$and the time scaling $\mathrm{d}t\sqrt{D}=\mathrm{d}\tau$ to transform (13(13) $\begin{array}{rl}\dot{x}& =x^2(p-qx)-{\displaystyle \frac{p-q}{A+1}}x(A+x)y,\\ \dot{y}& =[cxy+(1-c)y-y^2](A+x),\end{array}$(13) ) into $\begin{array}{rl}\dot{\hat{u}}& =-\hat{v}+f(\hat{u},\hat{v}),\\ \dot{\hat{v}}& =\hat{u}+g(\hat{u},\hat{v}),\end{array}$where $\begin{array}{rl}D& =\frac{(A+1{)}^2(cA+cq-A+c)}{A}>0,\\ f(\hat{u},\hat{v})& ={\hat{c}}_{20}{\hat{u}}^2+{\hat{c}}_{11}\hat{u}\hat{v}+{\hat{c}}_{02}{\hat{v}}^2+{\hat{c}}_{30}{\hat{u}}^3+{\hat{a}}_{21}{\hat{u}}^2\hat{v}+{\hat{a}}_{12}\hat{u}{\hat{v}}^2+{\hat{a}}_{03}{\hat{v}}^3+O(|\hat{u},\hat{v}{|}^4),\\ g(\hat{u},\hat{v})& ={\hat{d}}_{20}{\hat{u}}^2+{\hat{d}}_{11}\hat{u}\hat{v}+{\hat{d}}_{02}{\hat{v}}^2+{\hat{d}}_{30}{\hat{u}}^3+{\hat{d}}_{21}{\hat{u}}^2\hat{v}+{\hat{d}}_{12}\hat{u}{\hat{v}}^2+{\hat{d}}_{03}{\hat{v}}^3+O(|\hat{u},\hat{v}{|}^4),\end{array}$and $\begin{array}{rl}{\hat{c}}_{20}& =\frac{q+1}{\sqrt{D}},{\hat{c}}_{11}=\frac{A+2}{A+1},{\hat{c}}_{30}=-\frac{q+1}{\sqrt{D}},{\hat{c}}_{21}=-\frac{1}{A+1},\\ {\hat{d}}_{20}& =-\frac{(c-1)A^3+(cq+3c+q-1)A^2+(3cq+q^2+3c+q)A+c{q}^2+2cq+c}{-(A+1)(A+q+1)(Ac+cq-A+c)},\\ {\hat{d}}_{11}& =-\frac{A\sqrt{D}(cA+cq-A+c+q+1)}{(A+1)(A+q+1)(cA+cq-A+c)},\\ {\hat{d}}_{02}& =\frac{A}{A+q+1},{\hat{d}}_{30}=\frac{A(cA+qA+cq+q^2+c+2q+1)}{(A+1)(A+q+1)(cA+cq-A+c)},\\ {\hat{d}}_{21}& =\frac{A\sqrt{D}(cA+cq-A+c+q+1)}{(A+1{)}^2(A+q+1\left)\right(cA+cq-A+c)},\\ {\hat{d}}_{12}& =-\frac{A}{(A+1)(A+q+1)},{\hat{d}}_{03}={\hat{c}}_{02}={\hat{c}}_{12}={\hat{c}}_{03}=0.\end{array}$Making use of the formal series method in [31], we get the first-order Lyapunov number, $l_1=-\frac{1}{8}\frac{A\sqrt{D}(q+1)\left(\right(A^2+2qA-q-1)c-A^2-2qA)}{(A+1{)}^3(cA+cq-A+c{)}^2}.$If $A^2+2qA-q-1\le 0$ then $l_1>0$ and if $A^2+2qA-q-1>0$ then $l_1>0$ because $\frac{A}{A+q+1}<c<1$. This shows that $l_1$ is always positive.

Theorem 3.3

System (13(13) $\begin{array}{rl}\dot{x}& =x^2(p-qx)-{\displaystyle \frac{p-q}{A+1}}x(A+x)y,\\ \dot{y}& =[cxy+(1-c)y-y^2](A+x),\end{array}$(13) ) experiences subcritical Hopf bifurcation at $\widetilde{E_4}$ with $p=p^{\ast }$ and $\widetilde{E_4}$ is surrounded by an unstable limit cycle.

Now, we present some numerical simulations to demonstrate the feasibility of the above results on Hopf bifurcation. We fix A=1, $c=\frac{2}{3}$, and q=1.

1. With p=7+0.1, $\widetilde{E_4}$ is a hyperbolic unstable focus (see Figure 3(a)).

2. With p=7, $\widetilde{E_4}$ is an unstable weak focus (see Figure 3(b)).

3. With p=7−0.1, $\widetilde{E_4}$ is a hyperbolic stable focus surrounded by an unstable limit cycle (see Figure 3(c)).

4. With p=7−0.26, system (13(13) $\begin{array}{rl}\dot{x}& =x^2(p-qx)-{\displaystyle \frac{p-q}{A+1}}x(A+x)y,\\ \dot{y}& =[cxy+(1-c)y-y^2](A+x),\end{array}$(13) ) admits an unstable homoclinic loop which consists of an unstable homoclinic orbit and a saddle (see Figure 3(d)).

5. With p=7−0.4, $\widetilde{E_4}$ is a stable focus (see Figure 3(e)).

Figure 3. Phase portraits of (13(13) $\begin{array}{rl}\dot{x}& =x^2(p-qx)-{\displaystyle \frac{p-q}{A+1}}x(A+x)y,\\ \dot{y}& =[cxy+(1-c)y-y^2](A+x),\end{array}$(13) ) when p differs with A=1, $c=\frac{2}{3}$, and q=1. (a) $\widetilde{E_4}$ is a hyperbolic unstable focus when p=7+0.1. (b) $\widetilde{E_4}$ is an unstable weak focus when p=7. (c) $\widetilde{E_4}$ is a hyperbolic stable focus surrounded by an unstable limit cycle when p=7−0.1. (d) System (13(13) $\begin{array}{rl}\dot{x}& =x^2(p-qx)-{\displaystyle \frac{p-q}{A+1}}x(A+x)y,\\ \dot{y}& =[cxy+(1-c)y-y^2](A+x),\end{array}$(13) ) admits an unstable homoclinic loop consisting of an unstable homoclinic orbit and a saddle when p=7−0.26. (e) $\widetilde{E_4}$ is a stable focus when p=7−0.4.

As the bifurcation parameter p changes, the stability of $\widetilde{E_4}$ switches from unstable to stable, accompanied by an unstable limit cycle. Therefore, system (13(13) $\begin{array}{rl}\dot{x}& =x^2(p-qx)-{\displaystyle \frac{p-q}{A+1}}x(A+x)y,\\ \dot{y}& =[cxy+(1-c)y-y^2](A+x),\end{array}$(13) ) undergoes subcritical Hopf bifurcation at p=7.

3.3. Bogdanov–Takens bifurcation

When bm<1 and $A=A_{-}^{\ast }$, by Theorem 2.5, the degenerate equilibrium $E_{\ast }$ of (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) is a cusp of codimension two if $c=c_{\ast }$, $m=m_{\ast }$, b>1, and $x_{\ast }<\frac{b-1}{2b-1}$. This indicates the possible occurrence of Bogdanov–Takens bifurcation around $E_{\ast }$.

Theorem 3.4

With A and m as the bifurcation parameters, system (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ) undergoes a Bogdanov–Takens bifurcation of codimension two around the positive equilibrium $E_{\ast }$ when $A=A_{-}^{\ast }$, $c=c_{\ast }$, $m=m_{\ast }$, and condition (10(10) $b>1\text{and}0<x_{\ast }<\frac{b-1}{2b-1}.$(10) ) holds.

Proof.

We perturb m and A with $m_{\ast }+{ϵ}_2$ and $A_{-}^{\ast }+{ϵ}_1$ respectively to get the unfolding system of (4(4) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A+x}-bxy,\\ \dot{y}& =cxy+my-y^2.\end{array}$(4) ), (14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) where $ϵ=({ϵ}_1,{ϵ}_2)$ is a parameter vector in a small neighbourhood of the origin. Obviously, when ${ϵ}_1={ϵ}_2=0$, system (14(14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) ) has a cusp of codimension two.

First, the transformation $(x,y)=(u_1+x_{\ast },u_2+y_{\ast })$ translates the equilibrium $(x_{\ast },y_{\ast })$ to $(0,0)$ and transforms (14(14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) ) into (15) $\begin{array}{rl}{\dot{u}}_1& ={\tilde{a}}_{10}\left(ϵ\right)u_1+{\tilde{a}}_{01}\left(ϵ\right)u_2+{\tilde{a}}_{00}\left(ϵ\right)+{\tilde{a}}_{20}\left(ϵ\right)u_1^2\\ & +{\tilde{a}}_{11}\left(ϵ\right)u_1{u}_2+O(|u_1,u_2{|}^3),\\ {\dot{u}}_2& ={\tilde{b}}_{10}\left(ϵ\right)u_1+{\tilde{b}}_{01}\left(ϵ\right)u_2+{\tilde{b}}_{00}\left(ϵ\right)+{\tilde{b}}_{20}\left(ϵ\right)u_1^2\\ & +{\tilde{b}}_{11}\left(ϵ\right)u_1{u}_2+{\tilde{b}}_{02}\left(ϵ\right)u_2^2+O(|u_1,u_2{|}^3),\end{array}$(15) where $\begin{array}{rl}{\tilde{a}}_{10}\left(ϵ\right)& =\frac{-3x^2+1}{(A_{-}^{\ast }+x+{ϵ}_1{)}^2}+\frac{x_{\ast }^2(1-x_{\ast })}{A_{-}^{\ast }+x+{ϵ}_1}-b{y}_{\ast },\\ {\tilde{a}}_{01}\left(ϵ\right)& =-b{x}_{\ast },{\tilde{b}}_{10}\left(ϵ\right)=c_{\ast }{y}_{\ast },{\tilde{b}}_{01}\left(ϵ\right)=-y_{\ast }+{ϵ}_2,\\ {\tilde{a}}_{00}\left(ϵ\right)& =\frac{x_{\ast }^2(1-x_{\ast })}{A_{-}^{\ast }+{ϵ}_1+x_{\ast }}-b{x}_{\ast }{y}_{\ast },\\ {\tilde{a}}_{20}\left(ϵ\right)& =\frac{-3x_{\ast }+1}{A_{-}^{\ast }+x_{\ast }+{ϵ}_1}+\frac{x_{\ast }(3A_{-}^{\ast }{x}_{\ast }+2{x_{\ast }}^2+3x_{\ast }{ϵ}_1-2A_{-}^{\ast }-x_{\ast }-2{ϵ}_1)}{(A_{-}^{\ast }+x_{\ast }+{ϵ}_1{)}^3},\\ {\tilde{a}}_{11}\left(ϵ\right)& =-b,{\tilde{b}}_{00}\left(ϵ\right)={ϵ}_2{y}_{\ast },{\tilde{b}}_{20}\left(ϵ\right)=0,\\ {\tilde{b}}_{11}\left(ϵ\right)& =c,{\tilde{b}}_{02}\left(ϵ\right)=-1,\end{array}$Next, with the $C^{\mathrm{\infty }}$ change of coordinates in a small neighbourhood of the origin, $(u_3,u_4)=(u_1,{\tilde{a}}_{10}(ϵ)u_1+{\tilde{a}}_{01}(ϵ)u_2+{\tilde{a}}_{00}(ϵ)+{\tilde{a}}_{20}(ϵ)u_1^2+{\tilde{a}}_{11}(ϵ)u_1{u}_2+O(|u_1,u_2{|}^3\left)\right),$system (15(15) $\begin{array}{rl}{\dot{u}}_1& ={\tilde{a}}_{10}\left(ϵ\right)u_1+{\tilde{a}}_{01}\left(ϵ\right)u_2+{\tilde{a}}_{00}\left(ϵ\right)+{\tilde{a}}_{20}\left(ϵ\right)u_1^2\\ & +{\tilde{a}}_{11}\left(ϵ\right)u_1{u}_2+O(|u_1,u_2{|}^3),\\ {\dot{u}}_2& ={\tilde{b}}_{10}\left(ϵ\right)u_1+{\tilde{b}}_{01}\left(ϵ\right)u_2+{\tilde{b}}_{00}\left(ϵ\right)+{\tilde{b}}_{20}\left(ϵ\right)u_1^2\\ & +{\tilde{b}}_{11}\left(ϵ\right)u_1{u}_2+{\tilde{b}}_{02}\left(ϵ\right)u_2^2+O(|u_1,u_2{|}^3),\end{array}$(15) ) becomes (16) $\begin{array}{rl}{\dot{u}}_3& =u_4,\\ {\dot{u}}_4& ={\tilde{c}}_{00}\left(ϵ\right)+{\tilde{c}}_{10}\left(ϵ\right)u_3+{\tilde{c}}_{01}\left(ϵ\right)u_4+{\tilde{c}}_{20}\left(ϵ\right)u_3^2\\ & +{\tilde{c}}_{11}\left(ϵ\right)u_3{u}_4+{\tilde{c}}_{02}\left(ϵ\right)u_4^2+O(|u_3,u_4{|}^3),\end{array}$(16) where $\begin{array}{rl}{\tilde{c}}_{00}\left(ϵ\right)& =-{\tilde{b}}_{01}\left(ϵ\right){\tilde{a}}_{00}\left(ϵ\right),\\ {\tilde{c}}_{10}\left(ϵ\right)& ={\tilde{a}}_{01}\left(ϵ\right){\tilde{b}}_{10}\left(ϵ\right)-{\tilde{a}}_{10}\left(ϵ\right){\tilde{b}}_{01}\left(ϵ\right),\\ {\tilde{c}}_{01}\left(ϵ\right)& ={\tilde{a}}_{10}\left(ϵ\right)+{\tilde{b}}_{01}\left(ϵ\right)-\frac{{\tilde{a}}_{00}\left(ϵ\right){\tilde{a}}_{11}\left(ϵ\right)}{{\tilde{a}}_{01}\left(ϵ\right)},\\ {\tilde{c}}_{20}\left(ϵ\right)& =-{\tilde{b}}_{01}\left(ϵ\right){\tilde{a}}_{20}\left(ϵ\right)+{\tilde{b}}_{10}\left(ϵ\right){\tilde{a}}_{11}\left(ϵ\right),\\ {\tilde{c}}_{11}\left(ϵ\right)& =2{\tilde{a}}_{20}\left(ϵ\right)+\frac{{\tilde{a}}_{00}\left(ϵ\right){\tilde{a}}_{11}^2\left(ϵ\right)-{\tilde{a}}_{10}\left(ϵ\right){\tilde{a}}_{11}\left(ϵ\right)}{{\tilde{a}}_{01}^2\left(ϵ\right)},\\ {\tilde{c}}_{02}\left(ϵ\right)& =\frac{{\tilde{a}}_{11}\left(ϵ\right)}{{\tilde{a}}_{01}\left(ϵ\right)}.\end{array}$Now, introducing a new time variable τ with $\mathrm{d}t=(1-{\tilde{c}}_{02}(ϵ\left)u_3\right)d\tau$ and still denoting τ as t for simplicity, we rewrite system (16(16) $\begin{array}{rl}{\dot{u}}_3& =u_4,\\ {\dot{u}}_4& ={\tilde{c}}_{00}\left(ϵ\right)+{\tilde{c}}_{10}\left(ϵ\right)u_3+{\tilde{c}}_{01}\left(ϵ\right)u_4+{\tilde{c}}_{20}\left(ϵ\right)u_3^2\\ & +{\tilde{c}}_{11}\left(ϵ\right)u_3{u}_4+{\tilde{c}}_{02}\left(ϵ\right)u_4^2+O(|u_3,u_4{|}^3),\end{array}$(16) ) as (17) $\begin{array}{rl}{\dot{u}}_3& =u_4(1-{\tilde{c}}_{02}(ϵ\left)u_3\right),\\ {\dot{u}}_4& =(1-{\tilde{c}}_{02}(ϵ\left)u_3\right)\left[{\tilde{c}}_{00}\right(ϵ)+{\tilde{c}}_{10}(ϵ)u_3+{\tilde{c}}_{01}(ϵ)u_4\\ & +{\tilde{c}}_{20}\left(ϵ\right)u_3^2+{\tilde{c}}_{11}\left(ϵ\right)u_3{u}_4+{\tilde{c}}_{02}\left(ϵ\right)u_4^2+O(|u_3,u_4{|}^3)].\end{array}$(17) Then $(u_5,u_6)=(u_3,u_4(1-{\tilde{c}}_{02}\left(ϵ\right)u_3)$ transforms (17(17) $\begin{array}{rl}{\dot{u}}_3& =u_4(1-{\tilde{c}}_{02}(ϵ\left)u_3\right),\\ {\dot{u}}_4& =(1-{\tilde{c}}_{02}(ϵ\left)u_3\right)\left[{\tilde{c}}_{00}\right(ϵ)+{\tilde{c}}_{10}(ϵ)u_3+{\tilde{c}}_{01}(ϵ)u_4\\ & +{\tilde{c}}_{20}\left(ϵ\right)u_3^2+{\tilde{c}}_{11}\left(ϵ\right)u_3{u}_4+{\tilde{c}}_{02}\left(ϵ\right)u_4^2+O(|u_3,u_4{|}^3)].\end{array}$(17) ) into (18) $\begin{array}{rl}{\dot{u}}_5& =u_6,\\ {\dot{u}}_6& ={\tilde{d}}_{00}\left(ϵ\right)+{\tilde{d}}_{10}\left(ϵ\right)u_5+{\tilde{d}}_{01}\left(ϵ\right)u_6+{\tilde{d}}_{20}\left(ϵ\right)u_5^2\\ & +{\tilde{d}}_{11}\left(ϵ\right)u_5{u}_6+O(|u_5,u_6{|}^3),\end{array}$(18) where $\begin{array}{rl}{\tilde{d}}_{00}\left(ϵ\right)& ={\tilde{c}}_{00}\left(ϵ\right),\\ {\tilde{d}}_{10}\left(ϵ\right)& ={\tilde{c}}_{10}\left(ϵ\right)-2{\tilde{c}}_{00}\left(ϵ\right){\tilde{c}}_{02}\left(ϵ\right),\\ {\tilde{d}}_{01}\left(ϵ\right)& ={\tilde{c}}_{01}\left(ϵ\right),\\ {\tilde{d}}_{20}\left(ϵ\right)& ={\tilde{c}}_{20}\left(ϵ\right)-2{\tilde{c}}_{02}\left(ϵ\right){\tilde{c}}_{10}\left(ϵ\right)+{\tilde{c}}_{00}\left(ϵ\right){\tilde{c}}_{02}(ϵ{)}^2,\\ {\tilde{d}}_{11}\left(ϵ\right)& =-{\tilde{c}}_{01}\left(ϵ\right){\tilde{c}}_{02}\left(ϵ\right)+{\tilde{c}}_{11}\left(ϵ\right).\end{array}$As ${\tilde{d}}_{20}\left(0\right)=\frac{(b-1)(b{x}_{\ast }-x_{\ast }+1)(2b{x}_{\ast }-b-x_{\ast }+1{)}^2}{b^5{x}_{\ast }(x_{\ast }-1)}<0$, we have ${\tilde{d}}_{20}\left(ϵ\right)<0$ when ${ϵ}_1$ and ${ϵ}_2$ are small enough. With the following change of coordinates and time scaling, $(u_7,u_8)=(u_5,\frac{u_6}{\sqrt{-{\tilde{d}}_{20}\left(ϵ\right)}})\text{and}\tau =\sqrt{-{\tilde{d}}_{20}\left(ϵ\right)}t,$we can transform (18(18) $\begin{array}{rl}{\dot{u}}_5& =u_6,\\ {\dot{u}}_6& ={\tilde{d}}_{00}\left(ϵ\right)+{\tilde{d}}_{10}\left(ϵ\right)u_5+{\tilde{d}}_{01}\left(ϵ\right)u_6+{\tilde{d}}_{20}\left(ϵ\right)u_5^2\\ & +{\tilde{d}}_{11}\left(ϵ\right)u_5{u}_6+O(|u_5,u_6{|}^3),\end{array}$(18) ) into (19) $\begin{array}{rl}{\dot{u}}_7& =u_8,\\ {\dot{u}}_8& ={\tilde{e}}_{00}\left(ϵ\right)+{\tilde{e}}_{10}\left(ϵ\right)u_7+{\tilde{e}}_{01}\left(ϵ\right)u_8-u_7^2+{\tilde{e}}_{11}\left(ϵ\right)u_7{u}_8+O(|u_7,u_8{|}^3),\end{array}$(19) where ${\tilde{e}}_{00}\left(ϵ\right)=-\frac{{\tilde{d}}_{00}\left(ϵ\right)}{{\tilde{d}}_{20}\left(ϵ\right)},{\tilde{e}}_{10}\left(ϵ\right)=-\frac{{\tilde{d}}_{10}\left(ϵ\right)}{{\tilde{d}}_{20}\left(ϵ\right)},{\tilde{e}}_{01}\left(ϵ\right)=\frac{{\tilde{d}}_{01}\left(ϵ\right)}{\sqrt{-{\tilde{d}}_{20}\left(ϵ\right)}},{\tilde{e}}_{11}\left(ϵ\right)=\frac{{\tilde{d}}_{11}\left(ϵ\right)}{\sqrt{-{\tilde{d}}_{20}\left(ϵ\right)}}.$Furthermore, the transformation $(u_9,u_{10})=(u_7-\frac{{\tilde{e}}_{10}\left(ϵ\right)}{2},u_8)$transforms (19(19) $\begin{array}{rl}{\dot{u}}_7& =u_8,\\ {\dot{u}}_8& ={\tilde{e}}_{00}\left(ϵ\right)+{\tilde{e}}_{10}\left(ϵ\right)u_7+{\tilde{e}}_{01}\left(ϵ\right)u_8-u_7^2+{\tilde{e}}_{11}\left(ϵ\right)u_7{u}_8+O(|u_7,u_8{|}^3),\end{array}$(19) ) into $\begin{array}{rl}{\dot{u}}_9& =u_{10},\\ {\dot{u}}_{10}& ={\tilde{f}}_{00}\left(ϵ\right)+{\tilde{f}}_{01}\left(ϵ\right)u_{10}-u_9^2+{\tilde{f}}_{11}\left(ϵ\right)u_9{u}_{10}+O(|u_9,u_{10}{|}^3),\end{array}$where ${\tilde{f}}_{00}\left(ϵ\right)={\tilde{e}}_{00}\left(ϵ\right)+\frac{1}{4}{\tilde{e}}_{10}(ϵ{)}^2,{\tilde{f}}_{01}(ϵ)={\tilde{e}}_{01}(ϵ)+\frac{1}{2}{\tilde{e}}_{10}(ϵ\left){\tilde{e}}_{11}\right(ϵ),{\tilde{f}}_{11}(ϵ)={\tilde{e}}_{11}(ϵ).$It follows from ${\tilde{f}}_{11}\left(0\right)=-\frac{(2b{x}_{\ast }-b-x_{\ast }+1)(b-1)(2b{x}_{\ast }-x_{\ast }+1)}{b^3{x}_{\ast }(x_{\ast }-1)}<0$ that ${\tilde{f}}_{11}\left(ϵ\right)\ne 0$ when ${ϵ}_1$ and ${ϵ}_2$ are sufficiently small. Employing the change of variables and time scaling, $(u_{11},u_{12})=\left({\tilde{f}}_{11}\right(ϵ{)}^2{u}_9,-{\tilde{f}}_{11}\left(ϵ{)}^3{u}_{10}\right)\text{and}\tau =-\frac{1}{{\tilde{f}}_{11}\left(ϵ\right)}t,$we obtain the versal unfolding of system (14(14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) ), $\begin{array}{rl}{\dot{u}}_{11}& =u_{12},\\ {\dot{u}}_{12}& ={\tilde{g}}_{00}\left(ϵ\right)+{\tilde{g}}_{01}\left(ϵ\right)u_{12}+u_{11}^2+u_{11}{u}_{12}+O(|u_{11},u_{12}{|}^3),\end{array}$where ${\tilde{g}}_{00}=-{\tilde{f}}_{00}\left(ϵ\right){\tilde{f}}_{11}(ϵ{)}^4,{\tilde{g}}_{01}=-{\tilde{f}}_{01}(ϵ\left){\tilde{f}}_{11}\right(ϵ).$Notice that ${\left|\frac{\mathrm{\partial }({\tilde{g}}_{00},{\tilde{g}}_{01})}{\mathrm{\partial }({ϵ}_1,{ϵ}_2)}\right|}_{{ϵ}_1={ϵ}_2=0}=-\frac{1}{2}\frac{(b-1)(2b{x}_{\ast }-x_{\ast }+1{)}^5{b}^2}{(b{x}_{\ast }-x_{\ast }+1{)}^5(x_{\ast }-1{)}^2{x}_{\ast }}\ne 0.$Therefore, the Bogdanov–Takens bifurcation theory [10,21] implies that (14(14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) ) undergoes a subcritical Bogdanov–Takens bifurcation of codimension two when the parameters ${ϵ}_1$ and ${ϵ}_2$ vary in a small neighbourhood of $(0,0)$, which includes a sequence of bifurcations with codimension one saddle-node bifurcation, subcritical Hopf bifurcation, and Homoclinic bifurcation [26].

According to the versal unfolding of system (14(14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) ) and results in [10], the Bogdanov–Takens bifurcation diagram consists of saddle-node, Hopf, and homoclinic bifurcation curves whose local expressions in the small neighbourhood of $(0,0)$ are as follows.

1. The saddle-node bifurcation curve $SN=\left\{\right({ϵ}_1,{ϵ}_2):{\tilde{g}}_{00}({ϵ}_1,{ϵ}_2)=0,{\tilde{g}}_{01}({ϵ}_1,{ϵ}_2)\ne 0\};$

2. the Hopf-bifurcation curve $H=\left\{\right({ϵ}_1,{ϵ}_2):{\tilde{g}}_{00}({ϵ}_1,{ϵ}_2)<0,{\tilde{g}}_{01}({ϵ}_1,{ϵ}_2)=\sqrt{-{\tilde{g}}_{00}({ϵ}_1,{ϵ}_2})\};$

3. the homoclinic curve $HL=\left\{\right({ϵ}_1,{ϵ}_2):{\tilde{g}}_{00}({ϵ}_1,{ϵ}_2)<0,{\tilde{g}}_{01}({ϵ}_1,{ϵ}_2)=\frac{5}{7}\sqrt{-{\tilde{g}}_{00}({ϵ}_1,{ϵ}_2)}\}.$

We give some numerical simulations on Bogdanov–Takens bifurcation to end this section. Choosing $(A,b,c,m)=(0.6,2,0.25,0.05)$, we can draw the bifurcation diagram of (14(14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) ) including three bifurcation curves H, HL, and SN, which divide the small neighbourhood of $(0,0)$ in the ${ϵ}_1{ϵ}_2$-plane into four regions I, II, III, and IV (see Figure 4(a)). Further, the values of ${ϵ}_1$ and ${ϵ}_2$ will change phase portrait of system (14(14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) ) and we have the following results:

1. when $({ϵ}_1,{ϵ}_2)=(0.1,-0.001)$ lies in region I, system (14(14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) ) has no positive equilibrium, which means that predator and prey will not coexist over time (see Figure 4(b));

2. when $({ϵ}_1,{ϵ}_2)=(0.1,-0.01135)$ lies in region II, system (14(14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) ) has an unstable focus $E_4$ and a saddle $E_3$ as shown in Figure 4(c);

3. when $({ϵ}_1,{ϵ}_2)=(0.1,-0.0116)$ lies in region III, Figure 4(d) indicates that system (14(14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) ) has an unstable limit cycle with a stable focus and a saddle, which infer that populations of predator and prey present periodical changes;

4. when $({ϵ}_1,{ϵ}_2)=(0.1,-0.01171)$ lies on the curve HL, an unstable homoclinic loop appears in the phase portrait of system (14(14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) ) (see Figure 4(e));

5. finally, let $({ϵ}_1,{ϵ}_2)=(0.1,-0.0118)$ lie in region IV. Then Figure 4(f) tells us that a stable focus and a saddle come out in the phase portrait of system (14(14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) ). From an ecological point of view, whether predator and prey coexist or not depends on their initial population values.

4. Conclusion

In this paper, we studied the dynamic behaviour of a modified Lotka-Volterra predator–prey model where predators have other food resources and the prey is influenced by weak Allee effect. As the Allee parameter, A, changes, the model can have at most two positive equilibria. When there is a unique positive equilibrium, $E_{\ast }$, the equilibrium is degenerate. It can be a saddle-node of codimension 1 or a cusp of codimension 2. When there are two positive equilibria, one is always a saddle ($E_3$) and the other is an anti-saddle ($E_4$).

Figure 4. Phase portraits of system (14(14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) ). (a) Bifurcation diagram of system (14(14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) ). (b) System (14(14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) ) admits no positive equilibrium when $({ϵ}_1,{ϵ}_2)=(0.1,-0.001)$. (c) System (14(14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) ) admits an unstable focus and a saddle when $({ϵ}_1,{ϵ}_2)=(0.1,-0.01135)$. (d) System (14(14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) ) admits an unstable limit cycle with a stable focus and a saddle when $({ϵ}_1,{ϵ}_2)=(0.1,0.0116)$. (e) System (14(14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) ) admits an unstable homoclinic orbit and a saddle when $({ϵ}_1,{ϵ}_2)=(0.1,-0.01171)$. (f) System (14(14) $\begin{array}{rl}\dot{x}& =x(1-x)\frac{x}{A_{-}^{\ast }+{ϵ}_1+x}-bxy,\\ \dot{y}& =c_{\ast }xy+(m_{\ast }+{ϵ}_2)y-y^2,\end{array}$(14) ) admits a stable focus and a saddle when $({ϵ}_1,{ϵ}_2)=(0.1,-0.0118)$.

The system is not structurally stable when parameters change values. As a result, bifurcation phenomenon is very important from the point of ecosystems. Choosing A as the bifurcation parameter, the system can have saddle-node bifurcation at $E_{\ast }$ where saddle-node appears in the phase diagram. The modification of the bifurcation parameter gives the variation in the number of equilibria. Therefore, by changing the bifurcation parameter, we should get that the nullclines intersect at two points, then these two points are the same, and finally, the nullclines do not intersect. Theorem 2.4 tells us that the stability of the anti-saddle $E_4$ will switch from unstable to stable as the intrinsic fertility rate of the predator, m, increases and crosses $m^{\ast }$. Especially, the anti-saddle $E_4$ is a weak focus when $m=m^{\ast }$ and there can be a subcritical Hopf bifurcation, where unstable limit cycles will appear. Moreover, when both A and m are regarded as bifurcation parameters, non-degenerate Bogdanvo-Takens bifurcation of codimension 2 occurs at the cusp $E_{\ast }$. Thus, as A and m vary, the system produces rich dynamic behaviour, which includes saddle-node bifurcation, Hopf bifurcation, and homoclinic bifurcation.

From the findings, we observe that if the Allee effect is strong enough then the prey will become extinct and the system has no positive equilibrium. All orbits in the first quadrant tend to the prey-free equilibrium, $E_2(0,m)$, which is globally asymptotically stable (see Figure 2). This is because predators have other food resources. Moreover, $E_2$ is always a stable node (locally asymptotically stable). For appropriate values of A and m, the system has two positive equilibria. When one of them is stable, there is a bistable phenomenon (see Figure 1).

Compared with the work of Merdan [20], the modified system has extra food resources and hence it exhibits more complicated dynamical behaviour such as the change of the number of equilibria, limit cycles bifurcating via Hopf bifurcation, and appearance of homoclinic loops.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported partially by the Natural Science Foundation of Fujian Province [2020J01499] and NSERC of Canada [RGPIN-2019-05892].