Analysis of modified Holling-Tanner model with strong Allee effect

1.   Introduction

Predator-prey is one of three major types of interactions between species besides symbiosis and competition. Understanding the interactions between predators and their prey has been one of the leading research interests in population dynamics ^[1]. A predator-prey system is dominated by two important factors: the population growth function and the functional response. Tanner ^[2] presented a predator-prey system in which the environmental carrying capacity of the predator is proportional to the prey population size, and the reduction rate of the prey is proportional to the predator size. It takes the following form:

where $x(t)$ is the prey population and $y(t)$ is the predator population at time $t$; $r$ and $s$ are the intrinsic growth rate of the prey and predator respectively; $h, k, K$ and $D$ are all positive parameters. The term $\frac{kx}{x+D}$ is called Holling-Ⅱ functional response, generally reflecting the reduction rate of the prey caused by per capita of the predator.

By variable changes:

system (1) can be rewritten in a nondimensional form:

System (2) has been extensively studied ^[3,4,5,6,7]. Aziz-Alaoui and Okiye ^[8] considered the following system with alternative food sources for predators, which is called the modified Holling-Tanner model:

where $K_2 > 0$, $hx+K_2$ in system (3) is the new carrying capacity for predators. $K_2$ can be seen as an extra constant carrying capacity from all other food sources for predators. Several researchers ^{[9,10,11,12,13,14,15]} studied the existence of periodic solutions and bifurcation phenomena of system (3).

Allee effect refers to the phenomenon that low population density inhibits growth. Bioresearch indicates that clustering benefits the growth and survival of species. However, extreme sparsity and overcrowding will prevent population growth and negatively affect reproduction ^{[16,17,18,19]}. Every species has its optimal density. Species with small population densities are generally vulnerable. Once the population density falls below a critical level, interactions within the species will diminish. Ye et al. ^[20] considered a predator-prey model with a strong Allee effect and a nonconstant mortality rate. They found that a strong Allee effect may guarantee the coexistence of the species. Hu and Cao ^[21] considered a predator-prey model with Michaelis-Menten type predator harvesting. Their model exhibits a Bogdanov-Takens bifurcation of codimension 2. Xiang et al. ^[22] studied the Holling-Tanner model with constant prey harvesting. They found a degenerate Bogdanov-Takens bifurcation of codimension 4 and at least three limit cycles. In ^[23], Xiang et al. considered the Holling-Tanner model with predator and prey refuge and proved that this model undergoes a Bogdanov-Takens bifurcation of codimension 3. Arancibia et al. ^[24] adjusted the Holling-Tanner model by adding a strong Allee effect to prey. They found a Bogdanov-Takens bifurcation of codimension 2 and a heteroclinic bifurcation. It seems that a strong Allee effect gives rise to the heteroclinic bifurcation. Jia et al. ^[25] studied a modified Leslie-Gower model with a weak Allee effect on prey. This model undergoes a degenerate Bogdanov-Takens bifurcation of codimension 3 and has at least two limit cycles. Zhang and Qiao ^[26] analyzed the SIR model with vaccination and proved that this model undergoes a Bogdanov-Takens bifurcation of codimension 3 in some specific cases.

In this paper, we will analyze the following predator-prey model with the parameter $M$ indicating strong Allee effect:

where $\alpha > 0$ and $\beta > 0$ are real numbers. We shall assume $M\ll K$ from now on. If the original population of the prey is less than $M$, the prey will have a negative growth rate and become extinct ultimately. In system (4), $m > 0$ measures the extent of protection to which the environment provides to prey (or, to predators). It means that the growth rate of prey and predators are not negative infinity when $x$ reduces to zero.

As a first step in analyzing system (4), we nondimensionalize system (4) by writing

and it becomes

where $a = \frac{M}{K} < 1, \ b = \frac{m}{K}, \ c = \frac{s}{rK}$ and $\ d = \frac{\beta rK}{\alpha}.$

This paper is organized as follows. Section 2 discusses the stability of the equilibria. Section 3 deals with possible bifurcations that system (5) undergoes, such as Hopf bifurcation and Bogdanov-Takens bifurcation. Section 4 summarizes our conclusions.

2.   Equilibria and their stability

2.1.   Boundary equilibria and their stability

The equilibria of system (5) satisfy the following equations

We get four boundary equilibria on the axes: $E_0 = (0, 0), E_1 = (1, 0), E_2 = (a, 0)$ and $E_3 = (0, \frac{b}{d}).$ These equilibria are all hyperbolic because corresponding linearized matrices at the equilibria $E_i\ (i = 0, 1, 2, 3)$ are

Because we have assumed $0 < a < 1$, we get the stability of the boundary equilibria.

Theorem 2.1. System (5) has four equilibria $E_i\ (i = 0, 1, 2, 3)$ on the boundary. $E_0(0, 0)$ and $E_1(1, 0)$ are both hyperbolic saddles, $E_2(a, 0)$ is a hyperbolic unstable node and $E_3(0, \frac{b}{d})$ is a hyperbolic stable node (see Figure 1(c)).

2.2.   Positive equilibria and their stability

Theorem 2.2. If $\Delta = (a-1)^2-\frac{4}{d} > 0$, system (5) has two positive equilibria $E_4(x_4, \frac{x_4+b}{d})$ and $E_5(x_5, \frac{x_5+b}{d})$, where $x_4 = \frac{1+a-\sqrt{\Delta}}{2}$ and $x_5 = \frac{1+a+\sqrt{\Delta}}{2}$. $E_4$ is always a hyperbolic saddle. $E_5$ may be a source, a sink or a center depending on the parametric values (see Figure 1(a)).

Proof. The positive equilibria of system (5) satisfy the equations

If $\Delta = (a-1)^2-\frac{4}{d} > 0$, the equation $x^2-(1+a)x+a+\frac{1}{d} = 0$ has two roots $x_i\ (i = 4, 5)$.

The linearized matrix of system (5) at $E_i = (x_i, \frac{x_i+b}{d})\ (i = 4, 5)$ is

The determinant of $E_4$ is $detJ_{E_4} = -cx_4\sqrt{\Delta} < 0,$ from which we get that $E_4$ is a hyperbolic saddle.

Because $detJ_{E_5} = cx_5\sqrt{\Delta} > 0$ and $trJ_{E_5} = -x_5\sqrt{\Delta}-c+\frac{x_5}{d\left(x_5+b\right)},$ the stability of $E_5$ depends on the parameters $a, b, c$ and $d$. That is, $E_5$ is a sink if $trJ_{E_5} < 0$, a source if $trJ_{E_5} > 0$.

Theorem 2.3. If $\Delta = (a-1)^2-\frac{4}{d} = 0$, i.e., $d = \frac{4}{(a-1)^2}$, system (5) has a unique positive equilibrium $E_6(x_6, y_6)$, where $x_6 = \frac{1+a}{2}$ and $y_6 = \frac{1+a+2b}{2d}$ (see Figure 1(b)).

1) If $\ c > \frac{1+a}{d(1+a+2b)}\ (c < \frac{1+a}{d(1+a+2b)})$, $E_6$ is an attracting (a repelling) saddle-node;

2) If $\ c = \frac{1+a}{d(1+a+2b)}$, $E_6$ is a nilpotent cusp of codimension 2 (i.e., the Bogdanov-Takens singularity).

The phase portraits are given in Figure 3.

Proof. The linearized matrix at the equilibrium $E_6$ is

$E_6$ is not hyperbolic because $detJ_{E_6} = 0$. By a shift transformation $x-x_6\rightarrow x$ and $y-y_6\rightarrow y$, system (5) becomes

The eigenvalues of $J_{E_6}$ are $\lambda_1 = 0$ and $\lambda_2 = \frac{1+a}{d(1+a+2b)}-c$.

Case 1. $\lambda_2\neq0$, i.e., $c\neq \frac{1+a}{d(1+a+2b)}$.

By variable changes

system (6) can be rewritten as

By Theorem 7.1 in chapter 2 of ^[27], we obtain that $E_6$ is a saddle-node. There is a parabolic sector neighborhood in which all trajectories approach to $E_6$ when $\lambda_2 < 0$, and leave it when $\lambda_2 > 0$.

Case 2. $\lambda_2 = 0$, i.e., $c = \frac{1+a}{d(1+a+2b)}$.

According to ^[28], the cusp is a kind of nonhyperbolic critical point. The cusp can be illustrated by the following example:

The phase portrait for this system is shown in Figure 2. The neighborhood of the origin consists of two hyperbolic sectors and two separatrices.

Now consider the case when the linearized matrix (denoted as $A$) has two zero eigenvalues, i.e., $det\ A = 0, tr\ A = 0$, but $A\neq 0$. In this case it is shown in ^[28], that the system can be put in the normal form:

where $h(\xi), g(\xi)$ and $R(\xi, \eta)$ are analytic in a neighborhood of the origin, $h(0) = g(0) = 0, k\ge 2, a_k\neq 0$ and $n\ge 1$. The following lemma is given in ^[28].

Lemma 2.1. Let $k = 2m$ with $m\ge 1$ in system (8). Then the type of the origin is given by Table 1.

Next, we transform system (6) into a normal form by coordinate transformations. The corresponding linearized matrix of system (6) at $E_6$ is

By variable changes

system (6) can be rewritten as

Lemma 2.2. (^[28]) System (10)

is equivalent to system (11) near the origin, where system (11) is

Therefore, by Lemma 2.2, we can transform system (9) into the following form:

Using the notations of Lemma 2.1, we obtain $k = 2m, \ m = n = 1, \ a_{2m} = -\frac{(1+a)^2}{2(1+a+2b)}$ and $b_n = -\frac{d\left[(1 + a)^3 + (3 + a) (1 + 3 a) b + 4 (1 + a) b^2\right]}{(1+a+2b)^2}$. Consequently, by Lemma 2.1, we find that $E_6(x_6, y_6)$ is a degenerate critical point (cusp). More exactly, by the results in ^[28], $E_6$ is a cusp of codimension 2.

For example, we take $a = 0.1, b = 0.1$ and $d = 4.9383$, which satisfy $(a-1)^2-\frac{4}{d} = 0$. $E_6$ is a cusp for $c = \frac{1+a}{d(1+a+2b)} = 0.1713$. $E_6$ is a saddle-node with parabolic sector approaching it for $c = 0.3$. $E_6$ is a saddle-node with parabolic sector repelling it for $c = 0.05$ (see Figure 3).

3.   Bifurcation analysis

3.1.   Hopf bifurcation

As stated in Theorem 2.2, the positive equilibrium $E_5$ may be a center, source, or sink because $detJ_{E_5} > 0$. Considering $c$ as the bifurcation parameter, Hopf bifurcation occurs when $c = c_H$, where $c_H$ satisfies $trJ_{E_5}|_{c = c_H} = 0$. The local stability of $E_5$ changes when $c$ passes through $c = c_H$. We summarize our results in the following theorem.

Theorem 3.1. Hopf bifurcation occurs at $E_5$ in system (5) when $c = c_H > 0$.

Proof. Take $c$ as the bifurcation parameter. By $trJ_{E_5} = -x_5\sqrt{\Delta}+\frac{x_5}{d(x_5+b)}-c$, we have

Hopf bifurcation may occur when $c$ crosses $c_H$.

Next, we discuss the stability of $E_5$ as $c = c_H$. Moving $E_5$ to $(0, 0)$ by $X = x+x_5$ and $Y = y+y_5$, the Taylor expansion of system (5) at $E_5$ takes the form

where

With $\omega = \sqrt{-A_{10}^2-A_{01}B_{10}} > 0$ and

we transform system (13) into the following system

where $F(U, V)$ and $G(U, V)$ are the sum of those terms with orders not less than 2.

The stability of $O(0, 0)$ relies on the number

The following simplified expression of $\mathcal{K}$ can be obtained by Maple:

It is obvious that $I_3 > 0$ because all the elements are positive. Therefore, $E_5$ undergoes a subcritical bifurcation when $c = c_H$. When $c > c_H$ and $\left|c-c_H\right|\ll \varepsilon$, there is an unstable limit cycle (see Figure 4(b)).

3.2.   Bogdanov-Takens bifurcation

Theorem 3.2. When $c$ and $d$ are selected as two bifurcation parameters, $detJ_{E_5}|_{(c, d) = (c_{BT}, d_{BT})} = 0$ and $trJ_{E_5}|_{(c, d) = (c_{BT}, d_{BT})} = 0$, system (5) undergoes a Bogdanov-Takens bifurcation of codimension 2 in a small neighborhood of $E_5$ as $(c, d)$ varies near $(c_{BT}, d_{BT}) = (\frac{(a-1)^2}{4}\frac{1+a}{1+a+2b}, \frac{4}{(a-1)^2})$.

Proof. Perturb parameters $c$ and $d$ by $c = c_{BT}+\varepsilon_1$ and $d = d_{BT}+\varepsilon_2$, where $(\varepsilon_1, \varepsilon_2)$ is sufficiently small, and system (5) takes the following form:

The Taylor expansion of system (15) at $E_5(x_5, y_5)$ is

where

Take a $C^\infty$ change of coordinates around $(0, 0)$

then system (16) is changed into the following form

where

After rescaling the time by $(1-n_{02}u_1)t\rightarrow t$, system (17) is rewritten as

Letting $u_2 = u_1$ and $v_2 = v_1(1-n_{02}u_1)$, we get system (19) as follows

where

Case 1: For small $\varepsilon_i\ (i = 1, 2)$, if $\theta_{20} > 0$, by the following change of variables

system (19) becomes

where

To eliminate the $u_3$ term, letting $u_4 = u_3+\frac{s_{10}}{2}$ and $v_4 = v_3$, we get system (21) as follows

where

Clearly, $r_{11} = s_{11} = \frac{\theta_{11}}{\sqrt{\theta_{20}}}\neq 0$ if $\theta_{11}\neq 0$.

Setting $u_5 = r_{11}^2u_4, \ v_5 = r_{11}^3v_4$ and $\tau = \frac{1}{r_{11}}t$, we obtain the universal unfolding of system (16)

where

Case 2: For small $\varepsilon_i(i = 1, 2)$, if $\theta_{20} < 0$, by the following change of variables

system (19) becomes

where

To eliminate the $u_3$ term, letting $u^{\prime}_4 = u^{\prime}_3-\frac{s^{\prime}_{10}}{2}$ and $v^{\prime}_4 = v^{\prime}_3$, we get system (25) as follows

where

Clearly, $r^{\prime}_{11} = s^{\prime}_{11} = \frac{\theta_{11}}{\sqrt{-\theta_{20}}}\neq 0$ if $\theta_{11}\neq 0$.

Setting $u^{\prime}_5 = -r^{\prime 2}_{11}u^{\prime}_4, \ \ v^{\prime}_5 = r^{\prime 3}_{11}v^{\prime}_4$ and $\tau = -\frac{1}{r^{\prime}_{11}}t$, we obtain the universal unfolding of system (16)

where

Retain $\mu_1$ and $\mu_2$ to denote $\mu^{\prime}_1$ and $\mu^{\prime}_2$. If the matrix $\left|\frac{\partial (\mu_1, \mu_2)}{\partial (\varepsilon_1, \varepsilon_2)}\right|_{\varepsilon_1 = \varepsilon_2 = 0}$ is nonsingular, the parameter transformations (23) and (27) are homeomorphisms in a small neighborhood of $(0, 0)$, and $\mu_1, \mu_2$ are independent parameters. Direct computation shows that $\theta_{20} = -\frac{(1+a)^2(a-1)^2\left(a^2+2 a b+2 b^2+2 b+1\right)}{4(1+a+2 b)^3} < 0$ when $\varepsilon_i = 0\ (i = 1, 2)$ and

By Perko ^[28], we know that systems (22) and (26) undergo the Bogdanov-Takens bifurcation when $\varepsilon = (\varepsilon_1, \varepsilon_2)$ is in a small neighborhood of the origin. The local representations of the unfolding bifurcation curves are as follows ("+" for $\theta_{20} > 0$, "-" for $\theta_{20} < 0$):

1) The saddle-node bifurcation curve SN $= \left\{\left(\varepsilon_1, \varepsilon_2\right): \mu_1\left(\varepsilon_1, \varepsilon_2\right) = 0, \mu_2\left(\varepsilon_1, \varepsilon_2\right) \neq 0\right\}$;

2) The Hopf bifurcation curve H $= \left\{\left(\varepsilon_1, \varepsilon_2\right): \mu_2\left(\varepsilon_1, \varepsilon_2\right) = \pm \sqrt{-\mu_1\left(\varepsilon_1, \varepsilon_2\right)}, \mu_1\left(\varepsilon_1, \varepsilon_2\right) < 0\right\}$;

3) The homoclinic bifurcation curve HOM $= \left\{\left(\varepsilon_1, \varepsilon_2\right): \mu_2\left(\varepsilon_1, \varepsilon_2\right) = \pm \frac{5}{7} \sqrt{-\mu_1\left(\varepsilon_1, \varepsilon_2\right)}, \mu_1\left(\varepsilon_1, \varepsilon_2\right) < 0\right\}$.

For example, in system (15), we can set $a = 0.1, \ b = 0.1$ and $d_{BT} = \frac{4}{(a-1)^2}\approx 4.93827$. Further computation yields $x_5 = \frac{1+a}{2} = 0.55, \ y_5 = \frac{x_5+b}{d} = 0.131625$ and $c_{BT} = -2x_5^2+(1+a)x_5+\frac{x_5}{d(x_5+b)}\approx 0.171346.$ Since

the parametric transformation (27) is nonsingular. Moreover, $\theta_{20} = -0.139409- 0.81361\varepsilon_1- 0.058426\varepsilon_2 - 0.340982\varepsilon_1\varepsilon_2 < 0$ and $\theta_{11} = -1.05207 + 1.81818\varepsilon_1+ 0.126173\varepsilon_2 + 8.97868\varepsilon_1^2 + 1.67098\varepsilon_1\varepsilon_2 + 0.021619\varepsilon_2^2\neq 0$ for small $\varepsilon_i \ (i = 1, 2).$ The local representations of the bifurcation curves of system (15) up to second-order approximations are as follows. The details of the computation are given in Appendix.

1) The saddle-node bifurcation curve SN, is expressed as $\left\{\left(\varepsilon_1, \varepsilon_2\right): \varepsilon_2 = 0, \varepsilon_1 < 0\right\};$

2) The Hopf bifurcation curve H, is expressed as

3) The homoclinic bifurcation curve HOM, is expressed as

4) The heteroclinic bifurcation curve HET, can be detected with MATCONT, the numerical bifurcation package.

(a) When $\left(\varepsilon_1, \varepsilon_2\right) = (0, 0)$, $E_5$ is a cusp of codimension 2 (see Figure 6(a)).

(b) When $\left(\varepsilon_1, \varepsilon_2\right)$ is below SN, there are no positive equilibria (see Figure 6(b)).

(c) When $\left(\varepsilon_1, \varepsilon_2\right)$ crosses the SN curve and locates in the area between H and SN, there are two positive equilibria $E_4$ and $E_5$. $E_4$ is a saddle and $E_5$ is unstable (see Figure 6(c)).

(d) When $\left(\varepsilon_1, \varepsilon_2\right)$ crosses H, an unstable limit cycle will appear (see Figure 6(d)).

(e) When $\left(\varepsilon_1, \varepsilon_2\right)$ is on the curve HOM, there is an unstable homoclinic orbit (see Figure 6(e)).

(f) When $\left(\varepsilon_1, \varepsilon_2\right)$ is between the curve HOM and HET, $E_4$ connects with $E_2$ (see Figure 6(f)).

(g) When $\left(\varepsilon_1, \varepsilon_2\right)$ falls on the HET curve, $E_4$ connects with $E_1$ (see Figure 6(g)).

(h) When $\left(\varepsilon_1, \varepsilon_2\right)$ crosses the HET curve, $E_1$ connects with $E_5$ (see Figure 6(h)).

4.   Conclusions

This paper considers the modified Holling-Tanner model with a strong Allee effect. The aim is to explore the dynamical behaviors occurring in the predator-prey model with a strong Allee effect and alternative food sources for predators. Section 2 considers the equilibria and their stability. There exist four equilibria on the boundary. $E_0(0, 0)$ and $E_1(1, 0)$ are saddles, $E_2(a, 0)$ is an unstable node and $E_3(0, \frac{b}{d})$ is a stable node. As for the positive equilibria, when $d < \frac{4}{(a-1)^2}$, there is no positive equilibrium. When $d = \frac{4}{(a-1)^2}$, saddle-node bifurcation occurs and there is a unique positive equilibrium $E_6$, which is a cusp when $c = \frac{(a-1)^2}{4}\frac{1+a}{1+a+2b}$ and a saddle-node when $c \neq \frac{(a-1)^2}{4}\frac{1+a}{1+a+2b}$. When $d > \frac{4}{(a-1)^2}$, the saddle-node $E_6$ separates into two positive equilibria, a saddle $E_4$ and a hyperbolic equilibrium $E_5$. We examine the local stability of $E_5$ and find that $E_5$ is stable if $c > c_H$ and unstable if $c < c_H$. In Section 3.1, we prove that system (5) undergoes a Hopf bifurcation near $E_5$ when $c = c_H$ by showing that the constant $\mathcal{K} > 0$. In Section 3.2, we prove that system (5) exhibits a Bogdanov-Takens bifurcation of codimension 2 by calculating the universal unfolding near the cusp $E_6$. Besides, we give the bifurcation diagram with a little perturbation $(\varepsilon_1, \varepsilon_2)$ added to $(c_{BT}, d_{BT})$.

Our main result is that, after adding a strong Allee effect and alternative food sources, system (5) allows the independent survival of predators. In addition, a strong Allee effect makes the system more stable. Since, in the original Holling-Tanner system (2), there are at least two limit cycles ^[29], while after adding a strong Allee effect, there seems to be a unique stable limit cycle ^[24].

From the ecological viewpoint, alternative food sources for predators help them survive without prey. Besides, a strong Allee effect makes prey extinct at low density. Codimension of the Bogdanov-Takens bifurcation is at most two and the Hopf bifurcation is nondegenerate. We compare Holling-Tanner models with and without a strong Allee effect. It is found that a strong Allee effect increases and changes the dynamics. The first is that a strong Allee effect introduces saddle-node bifurcation. Second, heteroclinic bifurcation is brought about by the Allee effect, which means that the prey and predators may take different paths to reach distinct ultimate states. These indicate that the ecosystem may be sensitive to disturbances. It is essential to be aware of such bifurcations and protect the environment to weaken the Allee effect. In the future, we plan to study the system with different environmental protection for prey and predators. We will introduce $m_1$ and $m_2$ in place of $b$ to system (5) and consider the following equation:

If the functional response depends on the time as well as the prey and predator population, the model could exhibit more interesting dynamical behavior even chaos. We will do some research in our future works.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by National Natural Science Foundation of China (Grant No.62227810).

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

For $a = 0.1, \ b = 0.1$, system (16) takes the form

where

Take a $C^\infty$ change of coordinates around $(0, 0)$

then system (A.1) is changed into the following form

where

Letting $u_2 = u_1, v_2 = v_1(1-n_{02}u_1)$, we get system (A.3) as follows

where

For small $\varepsilon_i\ (i = 1, 2)$, $\theta_{20} < 0$, by the following change of variables

system (A.3) becomes

where

To eliminate the $u_3$ term, letting $u_4 = u_3-\frac{s_{10}}{2}, \ \ v_4 = v_3$, we get system (A.5) as follows

where

Setting $u_5 = -r^{ 2}_{11}u_4, \ \ v_5 = r^{ 3}_{11}v_4, \ \ \tau = -\frac{1}{r_{11}}t$, we obtain the universal unfolding of system (A.6)

where

1) The saddle-node bifurcation curve SN = $\left\{\left(\varepsilon_1, \varepsilon_2\right): \mu_1\left(\varepsilon_1, \varepsilon_2\right) = 0, \mu_2\left(\varepsilon_1, \varepsilon_2\right) \neq 0\right\}$;

2) The Hopf bifurcation curve H = $\left\{\left(\varepsilon_1, \varepsilon_2\right): \mu_2\left(\varepsilon_1, \varepsilon_2\right) = \pm \sqrt{-\mu_1\left(\varepsilon_1, \varepsilon_2\right)}, \mu_1\left(\varepsilon_1, \varepsilon_2\right) < 0\right\}$;

3) The homoclinic bifurcation curve HOM = $\left\{\left(\varepsilon_1, \varepsilon_2\right): \mu_2\left(\varepsilon_1, \varepsilon_2\right) = \pm \frac{5}{7} \sqrt{-\mu_1\left(\varepsilon_1, \varepsilon_2\right)}, \mu_1\left(\varepsilon_1, \varepsilon_2\right) < 0\right\}$.

With the specific parameters, we have

1) The saddle-node bifurcation curve, denoted SN, is expressed as

2) The Hopf bifurcation curve, denoted H, is expressed as

3) The homoclinic bifurcation curve, denoted HOM, is expressed as