Bifurcation analysis in a Holling-Tanner predator-prey model with strong Allee effect

1.   Introduction

Analyzing the dynamics of predator-prey models is of great importance in mathematical ecology for studying, interpreting and predicting species evolution, growth and interactions. It is well known that different types of biological models, exhibit great different dynamic behavior. Many authors discuss various types of predator-prey systems with hunting cooperation ^[1], Holling-Tanner ^[2,3], stochastic ^[4,5], and study their stability ^[6,7] and bifurcation ^[8,9,10,11]. One of the classical Leslie-Gower predator-prey model is

which was proposed by ^[12] and ^[13] to describe the relationship between predator and prey, where $x$ and $y$ represent population densities of prey and predator respectively, $r$ and $K$ are the intrinsic growth rate and the environmental capacity of the prey respectively, $s$ is intrinsic growth rate of the predator at time $t$, $n$ is a measure of the quality of the prey as food for the predator, $H(x)$ represents the functional response which is a key element to describe the number of predator consuming prey per unit time. The Leslie-Gower predator-prey model with Allee effect ^[14,15] and generalist predator ^[16,17], or modified Leslie-Gower predator-prey model ^[18,19,20] have been widely adopted in the biological model domains.

Functional responses play an important role in predator-prey systems and describe the transformation of organisms from lower to higher trophic levels in the biological chain. The differences in the dynamical behavior of predator-prey systems are partly attributable to the functional responses chosen. Three main types of functional responses were proposed in ^[21]: Holling Type I is a linear increasing function corresponding to lower animals, Holling Type II is hyperbolic in form corresponding to invertebrate and Holling Type III is sigmoid corresponding to vertebrate. S$\mathrm{\acute{a}}$ez and Gonz$\mathrm{\acute{a}}$lez-Olivares in ^[22] proposed the following Leslie-Gower predator-prey model with Holling Type II functional response

where $q$ represents the maximum capture rate per capita and $e$ is half of the saturated response level. They investigated the stability and bifurcation of system (1.2), and showed that local asymptotic stability of a positive equilibrium point does not imply global stability.

The Allee effect is an ecological phenomenon that describes population size and fitness ^[23,24,25]. In ^[26], the authors researched deeply a predator-prey system with strong Allee effect and Holling Type I functional response, and provided insights that system has a weak focus of order at least two. Furthermore, the authors explained the collapse of the positive equilibria is a cusp of codimension 2. In ^[27], they proposed a ratio-dependent Leslie-Gower predator-prey model with the Allee effect and fear effect on prey, and investigated the saddle-node bifurcation, degenerate Hopf bifurcation, and Bogdanov-Takens bifurcation. Martinez-Jeraldo and Aguirre analyzed multiplicative Allee effect and Holling Type I functional response, gave a concrete proof of a subcritical Hopf bifurcation or supercritical Hopf bifurcation, and performed a numerical simulation to prove the existence of a Bogdanov-Takens bifurcation (see ^[28]). In ^[29], they studied the stability and bifurcation of a Leslie-Gower predation model with Allee effect and an alternative food source.

Inspired by ^[22], the authors ^[30] studied a Holling-Tanner predator-prey model with strong Allee effect $(0 < m < K)$ on prey as follows:

For simplicity, making some substitutions

and still denoting $\tau$ by $t$, system (1.3) is topological equivalent to system

where $M, Q, S$ are positive constants.

The authors ^[30] proved the boundedness of system (1.5), and discussed the existence and stability of the equilibrium. They demonstrated the existence of separation lines in the phase plane separating basins of attraction associated with species extinction and coexistence. The author discussed numerous potential bifurcations such as saddle-node, Hopf, and Bogdanov-Takens bifurcations. However, they did not give a detailed proof of Hopf bifurcation and Bogdanov-Takens bifurcation.

Hence, in this paper, we also want to consider system (1.5), where the existence of positive equilibria are determined by third order polynomial, which makes the research more difficult. Using the different method with ^[30], we want to investigate the existence and stability of the system, and give the rigorous proof of Hopf bifurcation and Bogdanov-Takens bifurcation. We show that there exists bistable phenomena and two limit cycles in system (1.5).

The paper is organized as follows. In Section 2, the existence of equilibria are discussed. In Section 3, we focus on the stability and bifurcations of system (1.5), such as Hopf bifurcation and Bogdanov-Takens bifurcation. In Section 4, we give the proof of the Theorems. In Section 5, we give some numerical simulations to show the feasibility of main results. We give a brief summary in the last section.

2.   Existence of equilibria

We denote the domain of system (1.5) in phase plane

From Theorem 2 ^[30], the solutions of system (1.5) are ultimately upper bounded and eventually end up in $\Omega$ with initial values $u(0)\geq0$ and $v(0)\geq0$.

Apparently, system (1.5) has three boundary equilibria $(0, 0)$, $(1, 0)$ and $(M, 0)$. From ^[30], the origin $(0, 0)$ of system (1.5) is a non-hyperbolic attractor, $(1, 0)$ is a hyperbolic saddle and $(M, 0)$ is a hyperbolic repeller.

The positive equilibria satisfies:

Substituting $v = u$ into the first equation of (2.1), we get

and the derivative of $f(u)$ is

Denote

When $Q < Q_{*}$, the equation $f{'}(u) = 0$ has two real roots,

Clearly, when $Q\geq Q_{*}$, the equation $f(u) = 0$ has no positive root, that is system (1.5) has no positive equilibrium.

Note that $\frac{M+1-A}{3} < 1$ and $f'(1) = (1-M)(1+A)+Q > 0$, then $\overline{u} < 1$ if $Q < Q_{*}$. By computation, $f(1) = Q > 0$. Therefore, the positive solution of equation $f(u) = 0$ is in the interval $(0, 1)$.

Our next discussion is only under $Q < Q_{*}$. If $\overline{u}\leq 0$ or $\overline{u} > 0$ and $f(\overline{u}) > 0$, it is obvious that system (1.5) has no positive equilibrium. If $\overline{u} > 0$ and $f(\overline{u}) = 0$, equation $f(u) = 0$ has a unique positive real root $u_{\ast}$, that is system (1.5) has a unique equilibrium $E_{\ast} = (u_{\ast}, u_{\ast})$. If $\overline{u} > 0$ and $f(\overline{u}) < 0$, equation $f(u) = 0$ has two positive distinct real roots $u_{1, 2}$ (letting $u_{1} < u_{2}$), which implies that system (1.5) has two positive equilibria $E_{1, 2} = (u_{1, 2}, u_{1, 2})$.

Hence, we obtain the following theorem.

Theorem 2.1. The existence of positive equilibria of system (1.5) is classified as follows.

(1) System (1.5) has a unique equilibrium $E_{\ast} = (u_{\ast}, u_{\ast})$ if $Q < Q_{*}$, $\overline{u} > 0$ and $f(\overline{u}) = 0$.

(2) System (1.5) has two positive equilibria $E_{1} = (u_{1}, u_{1})$ and $E_{2} = (u_{2}, u_{2})$ if $Q < Q_{*}$, $\overline{u} > 0$ and $f(\overline{u}) < 0$.

(3) System (1.5) has no positive equilibrium if one of the following conditions holds:

(i) $Q\geq Q_{*}$; (ii) $Q < Q_{*}, \; \overline{u}\leq 0$; (iii) $Q < Q_{*}, \; \overline{u} > 0$ and $f(\overline{u}) > 0$.

In order to investigate the stability of the positive equilibria, we get the Jacobian matrix of system (1.5) at any positive equilibria and obtain

The determinant and trace of $J_{E}$ are respectively given by

and

The property of eigenvalues of $J_{E}$ plays a crucial role in determining the dynamics of each equilibria. In addition, equilibrium $E(u, u)$ of system (1.5) is an elementary equilibrium (a degenerate equilibrium, respectively) when $Det\big(J_{E}\big)\neq0$ ($= 0$, respectively). From the derivative property of $f(u)$, we have $f'(u^*) = 0, f'(u_1) < 0$ and $f'(u_2) > 0$. Hence, the positive equilibria $E_{\ast}$ is a degenerate equilibrium since $Det\big(J_{E_{*}}\big) = 0$. The positive equilibria $E_{1}$ and $E_{2}$ respectively are a hyperbolic saddle and an elementary equilibrium on account of $Det\big(J_{E_{1}}\big) < 0$ and $Det\big(J_{E_{2}}\big) > 0$.

According to $f'(u^*) = 0$, the Jacobian matrix at $E_{*}$ can be simplified to

and the trace of $J_{E_{*}}$ is

From $Det\big(J_{E_{\ast}}\big) = 0, \; Tr\big(J_{E_{\ast}}\big) = 0$ and $f(u_{\ast}) = 0$, we can express $S, \; A$ and $Q$ by $u_{\ast}$ and $M$ as follows:

Since $S > 0, \; A > 0$ and $Q > 0$, the degenerate equilibrium $E_{\ast}$ of system (1.5) satisfies the condition

Lemma 2.1. ^[31] The system

is equivalent to the system

by some nonsingular transformations in the neighborhood of $(0, 0)$.

Lemma 2.2. ^[31] The system

is equivalent to the system

by some nonsingular transformations in the neighborhood of $(0, 0)$, where $G = a_{31}-a_{30}a_{21}$.

3.   Main results

3.1.   Stability of equilibria

If (2.2) and $\sqrt{M} < u_{\ast} < \frac{M+1}{2}$ hold, it's easy to verify that the conditions of Theorem 2.1 (1) hold, that is $E_{\ast}$ exists. Then the following theorem shows that the degenerate equilibrium $E_{\ast}$ is a cusp of codimension 2 or 3.

Define

Theorem 3.1. If (2.2) and $\sqrt{M} < u_{\ast} < \frac{M+1}{2}$ are satisfied, the degenerate equilibrium $E_{\ast}$ is

(1) a cusp of codimension 2 if $M\in\big(0, 7-4\sqrt{3}\big)$ and $u_{\ast}\neq u_{2\ast}$, or $M\in\big[7-4\sqrt{3}, 1\big)$;

(2) a cusp of codimension 3 if $M\in\big(0, 7-4\sqrt{3}\big)$ and $u_{\ast} = u_{2\ast}$.

Now, we give the stability of the positive equilibria $E_{1}$ and $E_{2}$.

Theorem 3.2. If condition (2) in Theorem 2.1 holds, system (1.5) has two positive equilibria $E_{1}$ and $E_{2}$, where $E_{1}$ is always a hyperbolic saddle, and $E_{2}$ is

(1) unstable if $0 < S < S^{*}$;

(2) stable if $S > S^{*} > 0$ or $S^{*}\leq0$;

(3) may be a center or fine focus if $S = S^{*} > 0$,

where $S^{*} = \frac{u_{2}\left[-3u_{2}^2-2(A-M-1)u_{2}+AM+A-M\right]}{u_{2}+A}.$

3.2.   Hopf bifurcation

In this subsection, we will investigate the Hopf bifurcation of system (1.5).

When condition $(2)$ of Theorem 2.1 holds, system (1.5) has two different positive equilibria. From Theorem 3.2, $E_{1}$ is always a hyperbolic saddle, and $E_{2}$ is a repeller or an attractor which depends on the sign of $Tr\big(J_{E_{2}}\big)$. Changing the sign of $Tr\big(J_{E_{2}}\big)$, the stability of $E_{2}$ will change as well. In this section, we consider the condition $Tr\big(J_{E_{2}}\big) = 0$, which implies that system (1.5) may exist a Hopf bifurcation around $E_{2}$. In the process of calculating of the first-order Lyapunov number, we use the following transformation (see ^[33] and ^[34])

Dropping the bar, system (1.5) has the following form

where $A, \; Q, \; S$ are all positive constants and $M < 1 < a$. Because $\overline{E}_{2}(1, 1)$ is an equilibrium of system (3.2), we have $Q = (A+1)(a-1)(1-M)$. Hence system (3.2) becomes

Define

The Jacobian matrix of system (3.3) at $\overline{E}_{2}(1, 1)$ is

The determinant and trace of $J_{\overline{E}_{2}}$ are respectively

and

If $M < \frac{A+2-a}{Aa+1}\triangleq M_{1}$ and $1 < a < A+2$, then $Det\big(J_{\overline{E}_{2}}\big) > 0$. Then the stability of ${E}_{2}$ is determined by the sign of $Tr\big(J_{\overline{E}_{2}}\big)$. If $S_{\ast}\leq0$, i.e., $0 < M\leq \frac{2A+3-(A+2)a}{A+2-a}\triangleq M_{2}$ and $1 < a < \frac{2A+3}{A+2}$, we obtain $Tr\big(J_{\overline{E}_{2}}\big) < 0$. If $S_{\ast} > 0$, i.e., $M_{2} < M < M_{1}$ and $1 < a < \frac{2A+3}{A+2}$ or $0 < M < M_{1}$ and $\frac{2A+3}{A+2}\leq a < A+2$, then $Tr\big(J_{\overline{E}_{2}}\big) > 0\; ( = 0, \; < 0, \; \mathrm{respectively})$ when $S < S_{\ast}\; ( = S_{\ast}, \; > S_{\ast}, \; \mathrm{respectively})$. To summarize the above discussion, we have the following theorem.

Theorem 3.3. The stability of the equilibrium $\overline{E}_{2}$ of system (3.3) is classified as follows.

(1) $\overline{E}_{2}$ is a stable hyperbolic focus or a node if one of the following conditions holds:

(i) $0 < M\leq M_{2}$ and $1 < a < \frac{2A+3}{A+2}$;

(ii) $M_{2} < M < M_{1}$, $1 < a < \frac{2A+3}{A+2}$ and $S > S_{\ast}$;

(iii) $0 < M < M_{1}$, $\frac{2A+3}{A+2}\leq a < A+2$ and $S > S_{\ast}$.

(2) $\overline{E}_{2}$ is an unstable hyperbolic focus or a node if one of the following conditions holds:

(i) $M_{2} < M < M_{1}$, $1 < a < \frac{2A+3}{A+2}$ and $S < S_{\ast}$;

(ii) $0 < M < M_{1}$, $\frac{2A+3}{A+2}\leq a < A+2$ and $S < S_{\ast}$.

(3) $\overline{E}_{2}$ is maybe a fine focus or center if one of the following conditions holds:

(i) $M_{2} < M < M_{1}$, $1 < a < \frac{2A+3}{A+2}$ and $S = S_{\ast}$;

(ii) $0 < M < M_{1}$, $\frac{2A+3}{A+2}\leq a < A+2$ and $S = S_{\ast}$.

Remark 3.1. If $0 < M < M_{2}$ and $1 < a < \frac{2A+3}{A+2}$, $\overline{E}_{2}$ is a stable hyperbolic focus or a node. Using the transformations (1.4) and (3.1), the coefficients $M, \, \, A$ and $a$ of system (3.2) do not include the intrinsic growth rate of predator $s$ of the original system (1.3). Ecologically, when $M$ and $a$ are sufficiently small, we take certain initial value that two species can coexist in the form of steady state independent of the intrinsic growth rate of predator $s$.

When the case $(3)$ of Theorem 3.3 are satisfied, system (3.3) will go through a Hopf bifurcation. Let

with $k = M+a-1$. We obtain the following theorem about the Hopf bifurcation.

Theorem 3.4. Assume that the case $(3)$ of Theorem 3.3 are satisfied.

(1) System (3.3) undergoes a subcritical Hopf bifurcation and an unstable limit cycle appears around $\overline{E}_{2}$ when $R_{2} < 0$.

(2) System (3.3) undergoes a supercritical Hopf bifurcation and a stable limit cycle appears around $\overline{E}_{2}$ when $R_{2} > 0$.

(3) System (3.3) undergoes a degenerate Hopf bifurcation and multiple limit cycles may appear around $\overline{E}_{2}$ when $R_{2} = 0$.

3.3.   Bogdanov-Takens Bifurcation

From Theorem 3.1, the degenerate equilibrium $E_{\ast}$ of system (1.5) is a cusp of codimension 2 or 3. Hence, system (1.5) presumably exists a Bogdanov-Takens bifurcation of codimension 2 or 3. Then we will choose appropriate bifurcation parameters to verify the existence of the Bogdanov-Takens bifurcation.

When the equilibrium $E_{\ast}$ is a cusp of codimension 2, we select $Q$ and $S$ as bifurcation parameters, and have the following theorem.

Theorem 3.5. Assume that condition (2.2) and $\sqrt{M} < u_{\ast} < \frac{M+1}{2}$ are satisfied. If $M\in\big(0, 7-4\sqrt{3}\big)$ and $u_{\ast}\neq u_{2\ast}$, or $M\in\big[7-4\sqrt{3}, 1\big)$, the degenerate equilibrium $E_{\ast}$ of system (1.5) undergoes a Bogdanov-Takens bifurcation of codimension 2.

Now, we give the saddle-node bifurcation curve, Hopf bifurcation curve and homoclinic bifurcation curve of system (1.5) around $E_{\ast}$.

Theorem 3.6. Assume that condition (2.2) and $\sqrt{M} < u_{\ast} < \frac{M+1}{2}$ are satisfied. If $M\in\big(0, 7-4\sqrt{3}\big)$ and $u_{\ast}\neq u_{2\ast}$, or $M\in\big[7-4\sqrt{3}, 1\big)$, system (1.5) undergoes a Bogdanov-Takens bifurcation of codimension 2 around the degenerate equilibrium $E_{\ast}$ when $(\lambda_{1}, \lambda_{2})$ in a small neighborhood of $(0, 0)$. Moreover, there are three bifurcation curves.

(1) When $M\in\big(0, 7-4\sqrt{3}\big)$ and $u_{\ast}\in\left(u_{2\ast}, \frac{M+1}{2}\right)$, or $M\in\big[7-4\sqrt{3}, 1\big)$,

(2) When $M\in\big(0, 7-4\sqrt{3}\big)$ and $u_{\ast}\in\big(\sqrt{M}, u_{2\ast}\big)$,

$SN$, $H$ and $HL$ are respectively a saddle-node bifurcation curve, a Hopf bifurcation curve and a homoclinic bifurcation curve of system (1.5) around $E_{\ast}$.

Now we assume that the equilibrium $E_{\ast}$ is a cusp of codimension 3, and select $M, \; Q$ and $S$ as bifurcation parameters. In the next theorem, we prove system (1.5) undergoes a Bogdanov-Takens bifurcation of codimension 3 in a small parameter disturbance.

Theorem 3.7. Assume that condition (2.2) and $\sqrt{M} < u_{\ast} < \frac{M+1}{2}$ are satisfied. If $M\in\big(0, 7-4\sqrt{3}\big)$ and $u_{\ast} = u_{2\ast}$, the degenerate equilibrium $E_{\ast}$ of system (1.5) undergoes a Bogdanov-Takens bifurcation of codimension 3.

4.   Proof of Theorems

Proof of Theorem 3.1. When $S = (1-u_{\ast})(u_{\ast}-M), \; A = \frac{u_{\ast}^2(M+1-2u_{\ast})}{u_{\ast}^2-M}$ and $Q = \frac{(u_{\ast}-1)^2(M-u_{\ast})^2}{u_{\ast}^2-M}$, we move $E_{\ast}$ to the origin via the transformation $x = u-u_{\ast}$ and $y = v-u_{\ast}$. Then we have a Taylor expansion at the origin and get a new system

where the coefficients are given in Appendix A.

Next, we make a substitution $x = a_{01}x_{1}$ and $y = -a_{10}x_{1}+y_{1}$, then system (4.1) becomes

where the coefficients are given in Appendix A.

By Lemma 2.1, we get the equivalent system of (4.2) as follows:

where

and

When $\sqrt{M} < u_{\ast} < \frac{M+1}{2}$, $\frac{u_{\ast}^6(u_{\ast}-1)^4(u_{\ast}-M)^4}{(u_{\ast}^2-M)^3}\neq0$. Hence, the sign of $e_{20}$ is depends on $u_{\ast}^3-3Mu_{\ast}+M^2+M$.

Denote

Clearly, $h'(u_{*}) > 0$ if $\sqrt{M} < u_{\ast} < \frac{M+1}{2}$. Note that $h\big(\sqrt{M}\big) = M(\sqrt{M}-1)^2 > 0$, then $h(u_{*})\neq0$, which implies $e_{20}\neq0$.

Note that $\frac{u_{\ast}^4(u_{\ast}-1)^2(u_{\ast}-M)^2}{(u_{\ast}^2-M)^2}\neq0$ if $\sqrt{M} < u_{\ast} < \frac{M+1}{2}$. To determining the sign of $e_{11}$, we let

The discriminant of the equation $g(u_{\ast}) = 0$ is

Since $M\in(0, 1)$, a solution of $\Delta(M) = 0$ is given by $M_{1} = 7-4\sqrt3.$

If $M\in\big(0, 7-4\sqrt3\big)$, then $\Delta(M) > 0$, which implies that the equation $g(u_{\ast}) = 0$ has two different positive real roots

By calculation, we have $g(\sqrt{M}) = -\sqrt{M}(\sqrt{M}-1)^2(M-4\sqrt{M}+1) < 0$ and $g(\frac{M+1}{2}) = \frac{(M+1)(M-1)^2}{4} > 0$. Then, only $u_{2\ast}\in\big(\sqrt{M}, \frac{M+1}{2}\big)$. When $u_{\ast}\neq u_{2\ast}$, we have $g(u_{\ast})\neq0$, i.e., $e_{11}\neq0$.

If $M = 7-4\sqrt3,$ the equation $g(u_{\ast}) = 0$ has one positive real root $u_{\ast\ast} = 2-\sqrt{3} = \sqrt{M}$, i.e., $e_{11}\neq0$.

If $M\in\big(7-4\sqrt3, 1\big),$ we have $\Delta(M) < 0$. Clearly, $g(u_{\ast})\neq0$, i.e., $e_{11}\neq0$.

Hence, by the result in ^[32], $E_{\ast}$ is a cusp of codimension 2 if $M\in\big(0, 7-4\sqrt{3}\big)$ and $u_{\ast}\neq u_{2\ast}$, or $M\in\big[7-4\sqrt{3}, 1\big)$ (see Figures 1(a) and 1(b)).

On the other hand, if $M\in\big(0, 7-4\sqrt{3}\big)$ and $u_{\ast} = u_{2\ast}$, we obtain $e_{11} = 0$. Then system (4.1) becomes

where the coefficients are given in Appendix A.

Employing the transformation

system (4.5) can be written as

where the coefficients are given in Appendix A.

Next, introducing a new time variable $dt = (1-c_{02}\overline{x}_{1})d\tau$, and rewriting $\tau$ as $t$, we have

Making transformation $\overline{x}_{2} = \overline{x}_{1}$ and $\overline{y}_{2} = \overline{y}_{1}(1-\overline{c}_{02}\overline{x}_{1})$ once more, we obtain

where

Similar to the analysis of $e_{20}$, when $\sqrt{M} < u_{\ast} < \frac{M+1}{2}$, we have

Making the change of variable

and rewriting $\tau$ as $t$, system (4.8) turns to

where

By Lemma 2.2, system (4.9) is equivalent to

where $G = \overline{e}_{31}-\overline{e}_{30}\overline{e}_{21}$. After a brief calculation, we have

with

Note that

when $M\in\big(0, 7-4\sqrt{3}\big)$ and $\sqrt{M} < u_{2\ast} < \frac{M+1}{2}$. Therefore the sign of $G$ is related to $R_{1}$. As we observed in Figure 2, we have $R_{1} < 0$ if $M\in\big(0, 7-4\sqrt{3}\big)$, that is $G\neq0$. Our results demonstrate that $E_{\ast}$ is a cusp of codimension 3 if $M\in\big(0, 7-4\sqrt{3}\big)$ and $u_{\ast} = u_{2\ast}$ (see Figures 1(c) and 1(d)). The proof is completed.

Proof of Theorem 3.2. The equilibrium $E_{1}$ is a hyperbolic saddle is driven by $Det\big(J_{E_{1}}\big) < 0$. Notice that $Det\big(J_{E_{2}}\big) > 0$ and

When $S^{*}\leq0$, it is clear $Tr\big(J_{E_{2}}\big) < 0$, which implies $E_{2}$ is a hyperbolic stable node. When $S^{*} > 0$, obviously, $E_{2}$ is unstable (stable, a center or fine focus, respectively) if $S < S^{*}$ ($> S^{*}$, $= S^{*}$, respectively). The proof is completed.

Proof of Theorem 3.4. Obviously,

which satisfies the condition for the occurrence of Hopf bifurcation. Next we will calculate the first-order Lyapunov number which can determine the stability of limit cycle around $\overline{E}_{2}$.

First, we make the change to convert $\overline{E}_{2}$ to the origin by $\widetilde{x} = u-1$ and $\widetilde{y} = v-1.$ We have a Taylor expansion at the origin and system (3.3) becomes

where the coefficients are given in Appendix B.

On account of

we make a change of variables

System (4.11) becomes the following system

where the coefficients are given in Appendix B.

Subsequently, the first-order Lyapunov number in ^[32] at $\overline{E}_{2}$ is

The sign of $l_{1}$ is depends on $R_{2}$ since $8(1+A)^2(A+2-a-M(1+Aa)) > 0$. For specific conclusions, we can refer to Theorem 3.4. The proof is completed.

Proof of Theorem 3.5. Denote $(S_{0}, A_{0}, Q_{0}) = \Big((1-u_{\ast})(u_{\ast}-M), \frac{u_{\ast}^2(M+1-2u_{\ast})}{u_{\ast}^2-M}, \frac{(u_{\ast}-1)^2(u_{\ast}-M)^2}{u_{\ast}^2-M}\Big)$. Replacing $Q_{0}$ and $S_{0}$ with $Q_{0}+\lambda_{1}$ and $S_{0}+\lambda_{2}$, and substituting them into (1.5), we can obtain a new system

where $(\lambda_{1}, \lambda_{2})$ is a parameter vector in a small neighborhood of $(0, 0)$.

The first step, in order to move $E_{\ast}$ to the origin, we make a transformation $\widehat{x} = u-u_{\ast}$ and $\widehat{y} = v-u_{\ast}$, and obtain

where the coefficients are given in Appendix C and $P_{1}(\widehat{x}, \widehat{y}, \lambda), \; P_{2}(\widehat{x}, \widehat{y}, \lambda)$ are $C^{\infty}$ functions at least of third with respect to $(\widehat{x}, \widehat{y})$.

The second step, letting

system (4.14) can be written as

where the coefficients are given in Appendix C and $P_{3}(\widehat{x}_1, \widehat{y}_1, \lambda)$ is a $C^{\infty}$ function at least of third with respect to $(\widehat{x}_{1}, \widehat{y}_{1})$.

The third step, taking $dt = (1-\widehat{c}_{02}\widehat{x}_{1})d\tau$, $\widehat{x}_{2} = \widehat{x}_{1}$ and $\widehat{y}_{2} = (1-\widehat{c}_{02}\widehat{x}_{1})\widehat{y}_{1}$, system (4.15) is equivalent to

where

and $P_{4}(\widehat{x}_2, \widehat{y}_2, \lambda)$ is a $C^{\infty}$ function at least of third with respect to $(\widehat{x}_{2}, \widehat{y}_{2})$.

If $\sqrt{M} < u_{\ast} < \frac{M+1}{2}$, we have

where $\lambda_{1}$ and $\lambda_{2}$ are sufficiently small.

The fourth step, applying the following transformation

and still regarding $\tau$ as $t$, system (4.16) has the following form

where

and $P_{5}(\widehat{x}_3, \widehat{y}_3, \lambda)$ is a $C^{\infty}$ function at least of third with respect to $(\widehat{x}_{3}, \widehat{y}_{3})$.

The fifth step, making a transformations $\widehat{x}_{4} = \widehat{x}_{3}-\frac{\widehat{e}_{11}}{2}$ and $\widehat{y}_{4} = \widehat{y}_{3}$, system (4.17) turns to

where

and $P_{6}(\widehat{x}_4, \widehat{y}_4, \lambda)$ is a $C^{\infty}$ function at least of third with respect to $(\widehat{x}_{4}, \widehat{y}_{4})$.

Noting that $\lambda_{1}$ and $\lambda_{2}$ are sufficiently small, we obtain

when $\sqrt{M} < u_{\ast} < \frac{M+1}{2}$ and $u_{\ast}\neq u_{2\ast}$.

The last step, replacing

and rewriting $\tau$ as $t$, we get the universal unfolding of system (4.13)

where

and $P_{7}(\widehat{x}_5, \widehat{y}_5, \lambda)$ is a $C^{\infty}$ function at least of third with respect to $(\widehat{x}_{5}, \widehat{y}_{5})$.

We express $\mu_{1}$ and $\mu_{2}$ in terms of $\lambda_{1}$ and $\lambda_{2}$ as follows:

where the coefficients are given in Appendix C.

Note that

when $\sqrt{M} < u_{\ast} < \frac{M+1}{2}$, for $M\in\big(0, 7-4\sqrt{3}\big)$ and $u_{\ast}\neq u_{2\ast}$, or $M\in\big[7-4\sqrt{3}, 1\big)$.

Refer to ^[32], system (4.13) goes through a Bogdanov-Takens bifurcation of codimension 2 when $(\lambda_{1}, \lambda_{2})$ in a small neighborhood of $(0, 0)$. Trajectory topological classifications of Bogdanov-Takens bifurcation of codimension $2$ of system (4.13), see Figure 5. The proof is completed.

Proof of Theorem 3.6. The local bifurcation curve as in ^[32] is given by the following expression.

(i) The saddle-node bifurcation curve:

(ii) The Hopf bifurcation curve:

(iii) The homoclinic bifurcation curve:

Using the implicit function theorem, we can write $\lambda_{1}$ and $\lambda_{2}$ from $\mu_{1} = \mu_{1}(\lambda_{1}, \lambda_{2}, u_{\ast}, M)$ and $\mu_{2} = \mu_{2}(\lambda_{1}, \lambda_{2}, u_{\ast}, M)$ in (4.19) as follows:

where

where $Q_9$ can be found in Theorem 3.5.

Now, we consider the case $(1)$ of Theorem 3.6. When $M\in\big(0, 7-4\sqrt{3}\big)$ and $u_{\ast}\in\left(u_{2\ast}, \frac{M+1}{2}\right)$, or $M\in\big[7-4\sqrt{3}, 1\big)$, we obtain $s_{4} > 0$.

We consider the saddle-node bifurcation curve $\Gamma_{1} \triangleq \mu_{1}(\lambda_{1}, \lambda_{2}) = 0.$ By the implicit function theorem, we obtain a unique function $\lambda_{1}(\lambda_{2}) = 0$ which satisfies $\lambda_{1}(0) = 0$ and $\Gamma_{1}(\lambda_{1}(\lambda_{2}), \lambda_{2}) = 0$ since

In addition, on the curve $\Gamma_{1} = 0$, it folllows from (4.20) that $\lambda_{2} = s_{4}\mu_{2}+o(|\mu_{2}|)$. Then we obtain $\lambda_{2} > 0\; (<0)$ if $\mu_{2} > 0 \; (<0)$. Hence, the saddle-node bifurcation curve can be expressed as

The Hopf bifurcation curve is $\Gamma_{2} \triangleq \mu_{1}(\lambda_{1}, \lambda_{2})+\mu_{2}^2(\lambda_{1}, \lambda_{2}) = 0.$ Note that

Thus, from the implicit function theorem, we can obtain a unique function

which satisfies $\lambda_{1}(0) = 0$ and $\Gamma_{2}(\lambda_{1}(\lambda_{2}), \lambda_{2}) = 0$. On the curve $\Gamma_{2} = 0$, by (4.20), we obtain $\lambda_{2} = s_{4}\mu_{2}+o(|\mu_{2}|)$ and $\lambda_{2} > 0$ if $\mu_{2} > 0$. Thus, it is easy to obtain a Hopf bifurcation curve

Denote $\Gamma_{3} \triangleq \frac{25}{49}\mu_{1}(\lambda_{1}, \lambda_{2})+\mu_{2}^2(\lambda_{1}, \lambda_{2}) = 0$ as the homoclinic bifurcation curve. Notice that

Using the implicit function theorem, there exists a unique function

satisfying $\lambda_{1}(0) = 0$ and $\Gamma_{3}(\lambda_{1}(\lambda_{2}), \lambda_{2}) = 0$. Similarity, on the curve $\Gamma_{3} = 0$, we have $\lambda_{2} = s_{4}\mu_{2}+o(|\mu_{2}|)$ and $\lambda_{2} > 0$ if $\mu_{2} > 0$. The homoclinic bifurcation curve can be expressed as

The proof of Theorem 3.6 (2) is similar to the above proof, so we omit the detail proof here. The proof is completed.

Proof of Theorem 3.7. Denote $(S_{1}, A_{1}, Q_{1}) = \Big(-(u_{2\ast}-1)(u_{2\ast}-M), \frac{u_{2\ast}^2(M+1-2u_{2\ast})}{u_{2\ast}^2-M}, \frac{(u_{2\ast}-1)^2(u_{2\ast}-M)^2}{u_{2\ast}^2-M}\Big)$. Replacing $A_{1}, \; Q_{1}$ and $S_{1}$ with $A_{1}+\xi_{1}, \; Q_{1}+\xi_{2}$ and $S_{1}+\xi_{3}$ respectively, and substituting them into (1.5), we have the following system

where $(\xi_{1}, \xi_{2}, \xi_{3})$ is a parameter vector in a small neighborhood of $(0, 0, 0)$.

Next, inspired by ^[35] and ^[36], we intend to discuss the universal unfolding of a Bogdanov-Takens bifurcation of codimension 3.

The first step, moving the equilibrium $E_{\ast}$ of system (1.5) to the origin via the translation $X = u-u_{2\ast}, \; Y = v-u_{2\ast}$, system (4.21) is equivalent to

where the coefficients are given in Appendix D and $O_{1}(X, Y, \xi)$ is a $C^{\infty}$ function at least of fifth with respect to $(X, Y)$.

The second step, letting

system (4.22) can be written as

where the coefficients are given in Appendix D and $O_{2}(X, Y, \xi)$ is a $C^{\infty}$ function at least of fifth with respect to $(X_{1}, Y_{1})$.

The third step, we let $X_{1} = X_{2}+\frac{\gamma_{02}}{2}$ and $Y_{1} = Y_{2}+\gamma_{02}X_{2}Y_{2}$ to remove $Y_{1}^2$ from $\dot{Y_{1}}$. Then system (4.23) can be converted to

where the coefficients are given in Appendix D and

which has no effect on the bifurcation phenomenon (see ^[35] and ^[36]).

The fourth step, letting $X_{2} = X_{3}+\frac{\delta_{12}}{6}X_{3}^3$ and $Y_{2} = Y_{3}+\frac{\delta_{12}}{2}X_{3}^2Y_{3}$, system (4.24) turns into

where the coefficients are given in Appendix D and $R_{2}(X_{3}, Y_{3}, \xi)$ has the same properties as (4.25). Here we find that $X_{2}Y_{2}^2$ disappears from $\dot{Y_{3}}$ of system (4.25).

The fifth step, since

when $\xi_{1}, \; \xi_{2}$ and $\xi_{3}$ are sufficiently small, we make the following transformation

Let $\tau$ be represented by $t$. System (4.26) has a new expression

where the coefficients are given in Appendix D and $R_{3}(X_{4}, Y_{4}, \xi)$ has the same properties as (4.25).

The sixth step, since

when $\xi_{1}, \; \xi_{2}$ and $\xi_{3}$ are sufficiently small, we introduce a new transformation as follows

Therefore, system (4.27) has the following form ($\tau$ is rewritten by $t$),

where the coefficients are given in Appendix D and $R_{4}(X_{5}, Y_{5}, \xi)$ has the same properties as (4.25). Compared with (4.27), $X_{4}^3, \; X_{4}^4$ and $X_{4}^2Y_{4}$ disappear from $\dot{Y_{4}}$.

Clearly,

for $\sqrt{M} < u_{\ast} < \frac{M+1}{2}$. Figure 3(a) shows $R_{3} < 0$, where $R_{3}$ is listed in Appendix D.

By calculating, we have

and

for $\xi_{1}, \; \xi_{2}$ and $\xi_{3}$ are enough small.

The seventh step, we intend to convert $\eta_{20}$ and $\eta_{31}$ to $-1$ and $-1$ respectively. We denote

System (4.28) is rewritten as the following system

where the coefficients are given in Appendix D and $R_{5}(X_{6}, Y_{6}, \xi)$ has the same properties as (4.25).

The last step, we want to obtain the universal unfolding of the Bogdanov-Takens bifurcation of codimension 3. Our purpose is to remove $X_{6}$ from $\dot{Y_{6}}$ by letting $X_{6} = X_{7}-\frac{\theta_{10}}{2}$ and $Y_{6} = Y_{7}$. Hence, system (4.29) is equivalent to the following system

where the coefficients are given in Appendix D and $R_{6}(X_{7}, Y_{7}, \xi)$ has the same properties as (4.25).

After tedious calculations, we obtain

Actually,

for $\sqrt{M} < u_{\ast} < \frac{M+1}{2}$, and Figure 3(b) shows that $R_{4} < 0$ whose expression can be found in Appendix D. Therefore, the results indicates that system (4.13) undergoes a Bogdanov-Takens bifurcation of codimension 3 (^[37,38,39]). The proof is completed.

5.   Numerical simulations

Now, we give some numerical simulations to show the feasibility of Theorem 3.4.

Firstly, letting $a = 2, M = 0.1, A = 4, S = 0.28$, we obtain $Tr\big(J_{\overline{E}_{2}}\big) = 0$ and $l_{1} < 0$. Now we reduce $S = 0.28$ to $S = 0.27$ but do not change the other variables. Hence, system (3.3) undergoes a supercritical Hopf bifurcation and a stable limit cycle around $\overline{E}_{2}$ (see Figure 4(a)). From Figure 4(a), assuming that the other positive equilibrium of system (3.3) is $\overline{E}_{1}$, which is a saddle, we find the two stable manifolds of $\overline{E}_{1}$ can be treated as a separation curve between the basins of attraction of the origin and stable limit cycle around $\overline{E}_{2}$. In the biological sense, if the initial values above the separation curve, all trajectories converge to the origin. Thus the prey species and the predator species are all expected to be extinct. Additionally, if the initial values below the separation curve, all trajectories converge to the stable limit cycle. Hence the prey species and the predator species will oscillate and coexist.

Secondly, letting $a = 2, M = 0.3, A = 4, S = 0.44$, we have $Tr\big(J_{\overline{E}_{2}}\big) = 0$ and $l_{1} > 0$. Then, keeping the other variables without changing, we increase $S = 0.44$ to $S = 0.45$. Therefore, system (3.3) undergoes a subcritical Hopf bifurcation and an unstable limit cycle around $\overline{E}_{2}$ (see Figure 4(b)). As illustrated in Figure 4(b), we find that all trajectories inside the unstable limit cycle converge to $\overline{E}_{2}$, all trajectories outside the unstable limit cycle converge to the origin. Biologically, the two species will be extinct if the initial values lie outside the limit cycle, and the two species coexist if the initial values lie inside the limit cycle.

Thirdly, we give an example to show system (3.3) undergoes a degenerate Hopf bifurcation. When $a = 2, M = 0.2, A = 2.845, S = 0.408$, then we have $Tr\big(J_{\overline{E}_{2}}\big) = 0$ and $l_{1} = 0$. Keeping $a, M$ without changing, we perturb $A, S$ such that $A = 4, S = 0.358$. Hence, system (3.3) undergoes a degenerate Hopf bifurcation and two limit cycles (the inner one is stable and the outer one is unstable) around $\overline{E}_{2}$ (see Figure 4(c)). Figure 4(c) shows that the outside unstable limit cycle can be treated as a separation curve between the basins of attraction of the origin and stable limit cycle. Thus the prey and predator will be extinct if the initial values lie outside the unstable limit cycle. However, the prey and predator will oscillate and coexist if the initial values lie inside the unstable limit cycle.

Using (3.1), we are able to obtain the coefficients of original system (1.5) corresponding to the coefficients of system (3.3). For example, according to the coefficients of system (3.3) in Figure 4(a), i.e. $a = 2, M = 0.3, A = 4, S = 0.44$, we obtain the coefficients of the original system (1.5) such as $Q = 0.875, M = 0.15, A = 2, S = 0.11$ and two positive equilibria $E_1(0.352, 0.352)$ and $E_1(0.5, 0.5)$, where $E_1$ is a saddle and a stable limit cycle around $E_2$. The phase diagram of the original system (1.5) is similar to Figure 4(a), so we omit it.

Next, we give some numerical simulations to show the feasibility of Theorems 3.5 and 3.6. When $A = 24.2, \; Q = 5.0625, \; S = 0.1125,$ and $M = 0.3$, by Theorems 3.5 and 3.6, we obtain $E_{\ast}$ is a cusp of codimension two. and system (4.13) undergoes a Bogdanov-Takens bifurcation of codimension 2 (see Figure 5).

6.   Conclusions

We analyse the bifurcation of a Holling-Tanner predator-prey model with strong Allee effect on pery. ^[30] showed that system (1.5) undergoes the Hopf bifurcation and Bogdanov-Takens bifurcation by numerical simulation. Hence, using the different method with ^[30], we further investigate the rigorous proof of stability and bifurcations of system (1.5).

Now, using numerical simulations, we show the influence of Allee effect on the dynamical behavior of system (1.5). Let $A = 0.2, \, Q = 0.2, \, S = 0.02$. If $M = 0.5$, system (1.5) has no positive equilibrium. Then the origin is globally asymptotically stable (see Figure 6(a)), which implies that the prey and predator will be extinct. If $M = 0.225$, system (1.5) has two positive equilibria $E_{1} = (0.5591, 0.5591)$ and $E_{2} = (0.6, 0.6)$, where $E_{1}$ is a saddle and $E_{2}$ is unstable. From Figure 6(b), all solutions of system (1.5) tend to origin, which implies that the prey and predator is still extinct. If $M = 0.2$, the equilibrium $E_{2}$ changes from unstable to stable, and an unstable limit cycle occur. Form Figure 6(c), the unstable limit cycle becomes a separatrix curve. When the initial values lie outside of the limit cycle, the solutions will tend to origin. That is the prey and predator is still extinct. However if the initial values lie inside of the limit cycle, the solutions will tend to $E_{2}$, which implies that the prey and predator can coexist. If $M = 0.1$, system (1.5) has two positive equilibria $E_{1} = (0.2466, 0.2466)$ and $E_{2} = (0.7601, 0.7601)$, where $E_{1}$ is a saddle and $E_{2}$ is stable. From Figure 6(d), the unstable limit cycle disappears and the two stable manifolds of saddle of $E_{1}$ becomes a separatrix curve. That is, if the initial values lie right of the stable manifolds of saddle $E_{1}$, the solutions will tend to $E_{2}$, which implies that the prey and predator will stabilize to the value of $E_{2}$. If the initial values lie left of the stable manifolds of saddle $E_{1}$, the solutions will tend to origin, which implies that the prey and predator will be extinct. Hence, the small Allee effect leads to bistable phenomena. Finally, if $M = 0$, that is we consider the system (1.5) with out Allee effect, system has a unique positive equilibrium $(0.8385, 0.8385)$ which is globally asymptotically stable (see Figure 6(e)). Therefore, Allee effect is not conducive to the stability of the system (1.5). From Figure 6(e), the prey and predator will coexist if system (1.5) without the Allee effect. However, if the Allee effect is large, the system will be extinct (see Figure 6(a)). From Figures 6(b)-6(d), small Allee effect will lead to bistable phenomena. With the decrease of Allee effect, the region where prey and predator coexist becomes larger. Hence, compared to the larger Allee effect, the smaller Allee effect is beneficial to the survival of the prey and predator.

It follows from Theorem 2.1 that system (1.5) has no positive equilibria, a unique equilibrium $E_{\ast}$, or two positive equilibria $E_1$ and $E_2$ under some corresponding conditions. We confirm that the degenerate equilibrium $E_{\ast}$ of system can be a cusp of codimension 2 or 3. Further, Under some parameter perturbations, we also prove that system undergoes a Bogdanov-Takens bifurcation of codimension 2 or 3, and give the expressions of the saddle-node bifurcation curve, homoclinic bifurcation curve and the Hopf bifurcation curve. By translating $E_{2}$ to $\overline{E}_{2}(1, 1)$ (see ^[33] and ^[34]), we prove that system (1.5) undergoes a supercritical and subcritical Hopf bifurcation. Specially, system (1.5) undergoes a degenerate Hopf bifurcation and two limit cycles (the inner one is stable and the outer one is unstable) appear around $\overline{E}_{2}$ (see Figure 4(c)). There exists bistable phenomena, that is the stable manifold of the saddle or unstable limit cycle can be treated as a separation curve between the basins of attraction of the origin and stable limit cycle or $\overline{E}_{2}$. Hence, the main results of this paper complement and improve the previous paper ^[30]. In summary, the analysis about bifurcation enriches the dynamical behaviors of a Holling-Tanner predator-prey model with strong Allee effect. For the Holling-Tanner predator-prey model with strong Allee effect and Beddington-DeAngelis functional response, the dynamic properties of the system are worth further study by the method presented in this paper.

Acknowledgments

This work was supported by the Natural Science Foundation of Fujian Province (2021J01613, 2021J011032) and the Scientific Research Foundation of Minjiang University (MJY22027).

Conflict of interest

The authors declare there is no conflict of interest.

Appendix A.   Coefficients in the proof of Theorem 3.1

Appendix B.   Coefficients in the proof of Theorem 3.4

Appendix C.   Coefficients in the proof of Theorem 3.5

Appendix D.   Coefficients in the proof of Theorem 3.7