Stability and bifurcation analysis of the Bazykin's predator-prey ecosystem with Holling type Ⅱ functional response

1.   Introduction

As we all know, the classic Lotka-Volterra predator-prey system has been used to simulate predatory phenomena in nature, and its impact on the researches of mathematical biology and ecology will be roughly equivalent to the atomic bomb effect. In 1965, Holling proposed several different functional responses to characterize the dynamic predator-prey relationship between populations, which could describe the predator population how to transform the captured prey population into its own growth ability. Generally speaking, functional responses in predator-prey systems mainly depend on many internal and external key factors, such as the densities of prey and predator. On the other hand, prey-dependent functional responses are an important role in mathematical ecology, especially the dynamics of predator-prey systems. In recent decades, many scholars have done a lot of research on the predator-prey function response, and have made some excellent research results. Now, it is briefly summarized to enrich the population dynamics modeling system.

(ⅰ) The Holling type Ⅰ functional response ^[1] is

where $c$ is a positive constant.

(ⅱ) The Holling type Ⅱ or the Michaelis-Menten functional response ^[1,2,3,4,5] is

where $\alpha$, $a$ and $w$ are positive constants.

(ⅲ) The Holling type Ⅲ or the $p = 2$ S-type functional response ^[1,4,5] is

where $\alpha$ and $\beta$ are positive constants. The generalized Holling type Ⅲ or sigmoidal is $\Phi(x) = \frac{m x^2}{a x^2 +bx+1}$ with $b < -2 \sqrt{a}$, where $m$ and $a$ are positive constants, $b$ is a constant. When $b = 0$, the function $\Phi(x)$ can reduce to above Holling type Ⅲ functional response. The S-type functional response with index $p$ is $\Phi(x) = \frac{\alpha x^p}{\beta + x^p}$, where $\alpha$, $\beta$ and $p$ are positive constants.

(ⅳ) The generalized Holling Ⅳ or the Monod-Haldane functional response ^[1,4,5,6,7] is

where $m$ and $a$ are all positive constants, but $b$ is a constant. When $b = 0$, it is called the Holling type Ⅳ functional response or the simplified Holling type Ⅳ functional response ^[8,9].

(ⅴ) The Beddington-DeAngelis type functional response ^{[4,10,11,12,13]} is

where $\alpha$, $a$, $b$ and $c$ are positive constants, which is originally and independently introduced by Beddington and DeAngelis ^[10,11]. At the same time, it is similar to the Holling type Ⅱ functional response incorporating an extra term $cy$ in denominator, which can describe mutual interference among predators ^[12,13]. This functional response has some same qualitative features as the ratio-dependent form, but can avoid some singular behaviors of ratio-dependent models at low densities ^[12].

(ⅵ) The Hassell-Varley type functional response ^{[14,15,16,17,18]} is

where $A$ and $m$ are positive constants, $\gamma$ is the Hassell-Varley constant. When $\gamma = 0$ or $\gamma = 1$, it can be viewed as limiting cases mathematically. Specially, when $\gamma = 0$, it is the Holling type Ⅱ functional response regardless of constants.

(ⅶ) The Crowley-Martin type functional response ^[19] is

where $\alpha$, $a$, $b$, $c$ and $e$ are positive constants. It involves the interference among individual predator engaged in searching and handling the preys. While in ^[20,21,22], the authors particularly take the Crowley-Martin functional response in the new type of $\Phi(x, y) = \frac{\alpha x}{(1+ax)(1+by)}$, which is proposed by Bazykin and is an immense breakthrough of the Holling type Ⅱ and Beddington-DeAngelis functional responses.

In this paper, we will continually concentrate on a detailed discussion in the well-known Bazykin's predator-prey ecosystem with the Holling type Ⅱ functional response and interspecific density-restricted effect on the predator, which is also a variation of Volterra's classical predator-prey model and is expressed in the form of following ordinary differential equations (ODEs) ^[23,24]:

Here the functions $x = x(t)$ and $y = y(t)$ represent densities of the prey population and predator population at time $t$, respectively. All positive parameters have practically biological meanings: $r_1$ denotes the intrinsic growth rate of the prey population, $K_1$ represents the carrying capacity of the environment, $a$ is the half-saturation constant; $\alpha$ is the search efficiency of predator for prey, $m_1$ and $m_2$ are the mortality rate of the prey and predator population, $e$ is the biomass conversion, $d$ is the intra-specific competition coefficient. The specific growth term $r_1 x(1-\frac{x}{K_1})$ governs the increase of the prey in the lack of predator. The square non-linear term term $d y^2$, denotes intrinsic decrease of the predator, and represents interspecific density-restricted effect on the predator. Excluding the $d y^2$, the system (1.1) is based on the classical Gause type predator-prey system, which expresses the following form ^[25]:

where $g(x, k)$ is a continues function for $x > 0$, $p(x)$ is a functional response of predators, which is called Holling-type-Ⅱ predator-prey model as well. Such ordinary differential system of predator-prey populations is familiar to the Lotka-Volterra system, in which populations have the addition of damping terms(or self inhabit). The positivity of solutions with respect to initial condition $x(0) > 0$ and $y(0) > 0$ is easy to prove and we omit its proof. This system with positive and bounded solutions is also well behaved as we intuits from the biological significance.

Bazykin ^[24] wholly discussed the stability of equilibria, global existence of limit cycles, global attractivity of such equilibria, Hopf and codimension $2$ bifurcations. In ^[26], analytical description and alteration of local stability were given. Here the authors investigated a familiar Lotka-Volterra system, in which the populations have self-inhibit(the addition of damping term) for global stability and existence of limit cycles ^[27]:

where specific growth rate governs the growth of the prey in the absence of predator and it has an increase (or decrease) intrinsic rate on the predator. They already proved the existence of two limit cycles with the help of idea from the Poincare-Bendixson theory. Obviously, at special case $k_2 = 0$, we note that the system (1.3) can reduce to above system (1.1), which was analyzed by ^[28] and ^[29] for stability of equilibria, Hopf bifurcations, global attraction and codimension two bifurcations as well. While in the respect of global behavior, the system (1.1) was investigated by ^[30]. In ^[31], for a particular form of the system (1.1) with a modified Holling type Ⅱ functional response $\frac{\beta (x-m)}{1+ \alpha (x-m)}$ incorporating a constant prey refuge $m$, the authors therein gave sufficient conditions to guarantee the global stability of the positive equilibrium and uniqueness of a stable limit cycle. In ^[32], the authors revealed a rich bifurcation structure, including supercritical and subcritical Bogdanov-Takens bifurcation.

Based on classical biological manipulation theory, the Bazykin's predator-prey ecosystem can be used to explore the dynamic relationship between Microcystis aeruginosa and filter-feeding fish from the perspective of population dynamics, where $x(t)$ and $y(t)$ represent respectively the density of Microcystis aeruginosa and filter-feeding fish (bighead carp and silver carp), the growth kinetics function of Microcystis aeruginosa $x$ is $r_1 x(1-\frac{x}{K_1})$ with intrinsic growth rate $r_1$ and maximum environmental capacity $K_1$. The grazing function of filter-feeding fish $y$ is $\frac{\alpha xy}{1+ax}$ with capture coefficient $\alpha$ and density restriction coefficient $a$. Furthermore, the parameter $m_1$ and $m_{2}$ are natural mortality of Microcystis aeruginosa and filter-feeding fish, the parameter $d$ and $e$ are internal competition coefficient and energy conversion rate of filter feeding fish. In order to deeply explore the dynamic relationship between Microcystis aeruginosa and filter-feeding fish, it is necessary to investigate some bifurcation dynamic behaviors of the Bazykin's predator-prey ecosystem, thus we mainly focus on the stability and bifurcation of the Bazykin's predator-prey ecosystem in this paper. Firstly, we investigate the existence and stability of hyperbolic equilibrium point and non-hyperbolic equilibrium point, and conveniently study the cusp of condimension 3. Secondly, we explore Hopf bifurcation and Bogdanov-Takens bifurcation in detail, and give some sufficient threshold conditions. Finally, for Hopf bifurcation dynamics, we especially analyze the limit cycle via a perturbation procedure and canonical transformation. Moreover, we think that these mathematical analysis results can provide a theoretical basis for numerical simulation, which can give some biological interpretation for Hopf bifurcation and Bogdanov-Takens bifurcation, hence, we mainly study stability and bifurcation dynamics of the Bazykin's predator-prey ecosystem from the perspective of mathematical analysis, other biological significance issues will be completed in the follow-up work.

2.   Existence and stability analysis of equilibrium point

All solutions of the system (1.1) are non-negative and bounded with initial conditions $x(0) \geqslant 0$, $y(0) \geqslant 0$, thus it is namely dissipative in the first quadrant $\mathbb{R}^{+2}$ of $2$ dimensional space $\mathbb{R}^{2}$ and well-defined on the closed domain $\mathbb{R}_+^2 = \overline{\mathbb{R}^{+2}}$. Furthermore, the system is uniformly bounded with $\limsup \limits _{t \rightarrow + \infty} x(t) \leqslant M_1$, $\limsup \limits _{t \rightarrow + \infty} y(t) \leqslant M_2$, in which two positive constants $M_1$ and $M_2$ only depend on parameters $r_1$, $K_1$, $\alpha$, $a$, $m_1$, $e$, $m_2$ and $d$ ^[33]. In other words, the system (1.1) is confined in the domain

with $z = ex+y$ and a constant $M > 0$. Moreover, the system (1.1) is permanent if the value of all parameters can satisfy

with $\lambda \in (0, 1)$ ^[34].

2.1.   Boundary equilibrium point

It is clear that the system (1.1) admits two biological boundary equilibria $E_0 : = (0, 0)$ and $E_2 : = (x_2, 0)$(if $r_1 > m_1$), where $x_2 = K_1 (1-\frac{m_1}{r_1})$. For the Jacobian matrix $J(E_0) = {\rm{diag}} \{r_1-m_1, -m_2\}$ at point $E_0$, it shows that $E_0$ is a hyperbolic asymptotically stable node(unstable node) when $r_1 < m_1$($r_1 > m_1$), while it is a stable node(non-hyperbolic attractor) with only one zero eigenvalue if $r_1 = m_1$ (see the Theorem 7.1 in ^[35]). The Jacobian matrix at point $E_2$ is

then, the hyperbolic point $E_2$ is an asymptotically stable node(unstable node) when $\frac{\alpha e x_2}{a+x_2} < m_2$($\frac{\alpha e x_2}{a+x_2} > m_2$). When $\frac{\alpha e x_2}{a+x_2} = m_2$, $E_2$ is also a stable node(non-hyperbolic attractor) with characteristic direction $\tan(\theta) = \frac{(m_1-r_1) e}{m_2}$ under the polar coordinate transformation.

Now, we use an example to verify the stability of the equilibrium point $E_0$ and $E_2$ with $r_1 = 0.6$, $a = 1.5$, $\alpha = 0.5$, $e = 0.6$, $K_1 = 20$ and $d = 0.1$. It is easy to find from Figure 1(a) that $E_0$ is a stable node with characteristic direction $\theta = 0$ when $r_1 = m_1$ and $m_2 = 0.06$. Furthermore, it is obvious to see from Figure 1(b) that $E_2$ is a stable node with characteristic direction $\tan(\theta) = \frac{(m_1-r_1) e}{m_2}$ when $\frac{\alpha e x_2}{a+x_2} = m_2$ and $m_1 = 0.3$.

2.2.   Interior equilibrium point

2.2.1.   Hyperbolic equilibrium point

At first, an interior equilibrium $E_* = (x_*, y_*)$ of the system (1.1) always satisfies following algebraic equations

Furthermore, from Sylvester's resultants in polynomial Equation (2.4), components $x_*$ and $y_*$ must be positive roots of third-order polynomial equations(cubic equations) $p(x) = \sum \limits_{i = 0} ^3 a_i x^i = 0$ and $q(y) = \sum \limits_{i = 0} ^3 b_i y^i = 0$, respectively. Here the coefficients are listed as follows:

Then, we define these complicated expressions

and

which is a discriminant of above cubic equations $p(x) = 0$ and $q(y) = 0$ for later use respectively ^[26]. The Eq (2.4) also implies that such interior equilibrium $E_*$ does not exist when one of conditions holds: (ⅰ) $r_1 \leqslant m_1$; (ⅱ) $\alpha e \leqslant m_2$. The rest of our paper always assume the necessary existence condition $r_1 > m_1$ and $\alpha e > m_2$. If condition $0 < \frac{a m_2}{\alpha e-m_2} < x_2$ holds, it is cleat that $E_*$ always exists.

Thus we define the Jacobian matrix of the system (1.1) at an interior equilibrium $E_* = (x_*, y_*)$ as

The trace, determinant and discriminant of the matrix $J(E_*)$ are denoted as $A_1 = A_1(E_*) : = tr[J(E_*)]$, $A_2 = A_2(E_*) : = \det[J(E_*)]$ and $\Delta_* = \Delta_*(E_*) : = A_1^2 - 4 A_2$, respectively. Then, by using the Perron's theorems and the Routh-Hurwitz's criteria, we have following local stability of a hyperbolic equilibrium $E_*$ in general cases:

(a) If $A_1 < 0$ and (a1) $A_2 > 0, \Delta_* \geqslant 0$, then $E_*$ is an asymptotically stable node; (a2) $A_2 > 0, \Delta_* < 0$, then $E_*$ is an asymptotically stable focus; (a3) $A_2 < 0$, then $E_*$ is a saddle point;

(b) If $A_1 = 0$ and (b1) $A_2 > 0$, then $E_*$ is a center or a focus; (b2) $A_2 < 0$, then $E_*$ is a saddle point;

(c) If $A_1 > 0$, then $E_*$ is unstable and (c1) $\Delta_* = 0$, then $E_*$ is a node; (c2) $\Delta_* < 0$, then $E_*$ is a focus; (c3) $\Delta_* > 0$ and $A_2 > 0$, then $E_*$ is a node; (c4) $\Delta_* > 0$ and $A_2 < 0$, then $E_*$ is a saddle point.

A non-hyperbolic point $E_*$ is a stable(unstable) node if $A_1 < 0$($A_1 > 0$) and $A_2 = 0$. The nilpotent $E_*$ is probable a cusp of codimension at least $2$, which can ensure potential Bogdanov-Takens bifurcation when $A_1 = A_2 = 0$.

Obviously, when $r_1 \leqslant m_1$, the equilibrium point $E_0$ is globally asymptotically stable, which can be proved by using a Lyapunov function $V = ex+y$. Similarly, when $r_1 > m_1$ and $\alpha e < m_2$, the point $E_0$ is unstable and $E_*$ does not exist, thus $E_2$ is globally asymptotically stable. For equilibrium point $E_*$, we define a positive definite Lyapunov function

Now, along solutions of the system (1.1), differentiate $V$ with regard to time $t$ to obtain

Substituting the value of $\frac{dx}{dt}$ and $\frac{dy}{dt}$ from the system (1.1), we can get

Thus, it is obvious that if $\frac{\alpha y_*}{a(a+x_*)} < \frac{r_1}{K_1}$, then $\frac{dV}{dt} \leqslant 0$. This equality holds if and only if $(x, y) = E_{*}$, i.e., the equilibrium point $E_{*}$ is globally asymptotically stable.

Under a generalized condition $\frac{\alpha y_*}{a(a+x_*)} \leqslant \frac{r_1}{K_1}$, the hyperbolic point $E_*$ is a locally asymptotically stable focus or node since $A_1 < 0$ and $A_2 > 0$. Hence we assume $y_* = \frac{a(a+x_*) r_1 \rho}{\alpha K_1}$ with the introducing of a control variable $\rho \in (0, 1]$, components of the corresponding equilibrium point $E_* = (x_*, y_*)$ and restricted parameter $d$ are

Thus, we can obtain Theorem 2.1, which can guarantee that the equilibrium point $E_*$ is globally asymptotically stable.

Theorem 2.1. If the system (1.1) has a unique interior equilibrium point $E_*$ with $\frac{\alpha y_*}{a(a+x_*)} < \frac{r_1}{K_1}$, then the equilibrium point $E_*$ is globally asymptotically stable (hyperbolic focus or node).

At the same time, by using the Bendixson-Dulac criteria and a proper Dulac function $B(x, y) = \frac{1}{xy}$ or $B(x, y) = \frac{1}{(x+c)y}$ with a positive constant $c > 0$, then theorem 2.2 can also guarantee that a unique equilibrium point $E_*$ is globally asymptotically stable on account of the non-existence of closed orbits and limit cycles.

Theorem 2.2. If the system (1.1) has a unique interior equilibrium $E_*$, and $0 < \frac{a m_2}{\alpha e-m_2} < x_2 < a$, then the equilibrium $E_*$ is globally asymptotically stable.

In order to verify feasibility of Theorem 2.1 and 2.2, we will give some numerical simulations. If we take $r_1 = 0.6$, $m_1 = 0.2$, $m_2 = 0.1$, $\alpha = 0.5$, $a = 1.5$, $K_1 = 20$, $e = 0.6$ and $\rho = 0.9$, the calculation shows that the values of these parameters can satisfy with the condition of Theorem $2.1$, the equilibrium point $E_*$ is a globally asymptotically stable node, which can be seen in Figure 2(a). If we take $r_1 = 0.6$, $m_1 = 0.2$, $m_2 = 0.1$, $\alpha = 0.5$, $a = 1.5$, $e = 0.6$, $d = 0.05$ and $x_2 = a \rho(\rho = 0.9)$, then Theorem $2.2$ is true, thus the equilibrium point $E_*$ is a globally asymptotically stable point, which can be found from Figure 2(b). In a word, the equilibrium point $E_*$ is globally asymptotically stable under certain conditions.

2.2.2.   Non-hyperbolic equilibrium point: case (C1)

In this section, we mainly consider existence and stability of interior equilibrium in special cases, which can ensure the occurrence of Hopf bifurcation and Bogdanov-Takens bifurcation. In order to simplify $A_1$ as zero or get potential Hopf bifurcation, we take some parameters in a special case as

where $\lambda \in (0, 1)$ and $\mu = \sqrt{1+ 8 \lambda \bar{\mu}} > 1$ are two control variables for later use to scale parameters. Then the parameter $d$ is constrained by $d = \frac{4 \alpha (\mu+1) \varphi_d(\lambda, \mu)}{a (\mu+3)^2 \psi_d(\lambda, \mu)}$, where auxiliary functions are

It is quite clear that we can derive a little complicated expressions of a required interior equilibrium point $E_*^{(2)} = (x_*^{(2)}, y_*^{(2)})$ with $x_*^{(2)} = \frac{1}{4}a(\mu -1)$, $y_*^{(2)} = - \frac{ae(\mu+3)}{4(\mu+1)} \psi_d(\lambda, \mu)$, $A_1(E_*^{(2)}) = 0$ and $A_2(E_*^{(2)}) = \frac{- \alpha^2 e^2}{(\mu+3)^2 (\mu+1)} \varphi_{_{A_2}}(\lambda, \mu)$, in which $\varphi_{_{A_2}}(\lambda, \mu) = (\lambda-1)^2 \mu^3+(\lambda-1)(7 \lambda+5) \mu^2 +(15 \lambda^2+10 \lambda-1) \mu +9 \lambda^2-6 \lambda+5$ is also an auxiliary function.

Obviously, the expressions of $y_*^{(2)}$ and positive parameter $d$ indicate that $\varphi_d(\lambda, \mu) < 0$ and $\psi_d(\lambda, \mu) < 0$, i.e., $\mu > \mu_m : = \frac{3 \lambda+1}{1- \lambda}$. The unique positive root(stationary point) of the following equation

which can satisfy above condition, is $\mu_*^{(2)} = \frac{7 \lambda +5 + 2 \sqrt{\lambda^2+10 \lambda+7}}{3(1- \lambda)}$. Substituting it into $\varphi_{_{A_2}}(\lambda, \mu)$, we have $\varphi_{_{A_2}}(\lambda, \mu_*^{(2)}) = \frac{16}{27(\lambda-1)} \varphi_{_{\varphi A_2 \mu_*^{(2)} }}(\lambda)$, where $\varphi_{_{\varphi A_2 \mu_*^{(2)}}}(\lambda) = (\lambda^2+10 \lambda+7)^{\frac{3}{2}} - \lambda^3 -15 \lambda^2 + 60 \lambda +10$ is an auxiliary function of $\lambda$. Letting $\frac{d}{d \lambda} \varphi_{_{\varphi A_2 \mu_*^{(2)}}}(\lambda) = 0$, we have a unique negative root $\lambda = -5 - \frac{15}{4} \sqrt{2}$. Combining $\varphi_{_{\varphi A_2 \mu_*^{(2)}}}(0) = 10+7 \sqrt{7}$ and $\varphi_{_{\varphi A_2 \mu_*^{(2)}}}(\frac{1}{2}) = 79$, we know that the function $\varphi_{_{\varphi A_2 \mu_*^{(2)}}}(\lambda)$ must be monotonically increasing and has a positive minimum $\varphi_{_{\varphi A_2 \mu_*^{(2)}}}(0)$ on the interval $[0, 1]$, or this function is always positive on $(0, 1)$, while function $\varphi_{_{A_2}}(\lambda, \mu_*^{(2)})$ is always negative on $(0, 1)$. The second order partial derivative with respect to $\mu$ at the point $\mu_*^{(2)}$ is

which can ensure that the function $\varphi_{_{A_2}}(\lambda, \mu)$ owns a negative local minimum at the point $\mu_*^{(2)}$.

At this time, combining $\varphi_{_{A_2}}(\lambda, \mu_m) = \frac{32 \lambda}{\lambda -1} < 0$, we have a rough estimation of this function $\varphi_{_{A_2}}(\lambda, \mu)$ when $\mu > \mu_m$ and $\lambda \in (0, 1)$:

(ⅰ) When $\mu_m < \mu < \mu_*^{(2)}$, the function $\varphi_{_{A_2}}(\lambda, \mu)$ is monotonically decreasing and negative with respect to variable $\mu$.

(ⅱ) When $\mu > \mu_*^{(2)}$, the function $\varphi_{_{A_2}}(\lambda, \mu)$ is monotonically increasing with respect to variable $\mu$.

(ⅲ) When $\mu = \mu_*^{(2)}$, the function $\varphi_{_{A_2}}(\lambda, \mu)$ has a negative (local) minimum.

With the positive coefficient of leading order term $\mu^3$ in polynomial $\varphi_{_{A_2}}(\lambda, \mu)$ at hand, we have a positive value of variable $\mu$, which is also sufficiently large and is larger than $\mu_*^{(2)}$, such that $\varphi_{_{A_2}}(\lambda, \mu) > 0$. The zero theorem implies that the equation $\varphi_{_{A_2}}(\lambda, \mu) = 0$ must have a unique positive root, which is denoted as $\mu_1$ and on the right hand side of the point $\mu_*^{(2)}$.

All in all, we can confirm the classification of the interior equilibrium point $E_*^{(2)}$:

(ⅰ) When $\mu_m < \mu < \mu_1$, we have $\varphi_{_{A_2}}(\lambda, \mu) < 0$ and $A_2(E_*^{(2)}) > 0$, thus the equilibrium point $E_*^{(2)}$ is a focus or center;

(ⅱ) When $\mu = \mu_1$, namely $\varphi_{_{A_2}}(\lambda, \mu) = 0$ or $\mu$ satisfies the following cubic equation

we have a degenerate interior equilibrium $E_*^{(2)}$ with $A_1(E_*^{(2)}) = A_2(E_*^{(2)}) = 0$;

(ⅲ) When $\mu > \mu_1$, we have $\varphi_{_{A_2}}(\lambda, \mu) > 0$ and $A_2(E_*^{(2)}) < 0$, thus the equilibrium point $E_*^{(2)}$ is a saddle point.

Based on these analysis, we summarize two cases for consideration of stability and type of the interior equilibrium point $E_*^{(2)}$.

Case Ⅰ: $\mu_m < \mu < \mu_1$

When $\mu \in (\mu_m, \mu_1)$, making a linear transformation (Ⅰ): $x = u+x_*^{(2)}$, $y = v+y_*^{(2)}$, we can transfer the equilibrium point $E_*^{(2)}$ to the origin $O$. Then we can construct a transformation (Ⅱ): $u = -d y_*^{(2)} X+ \beta Y$, $v = - \frac{\alpha ea y_*^{(2)}}{(a+x_*^{(2)})^2} X$ to obtain the Jordan form of the system (1.1), thus the new system is

where $\beta = \sqrt{A_2(E_*^{(2)})} > 0$. Following Perko's book ^[36] or ^[37], the first Lyapunov number of the system (2.10) under the condition (2.8) is

where the auxiliary function is $\varphi_{\sigma}(\lambda, \mu) = (\lambda-1)^2 \mu^3+(\lambda-1)(9 \lambda-1) \mu^2+(27 \lambda^2-18 \lambda+15)(\mu+1)$. The partial derivative of function $\varphi_{\sigma}(\lambda, \mu)$ with respect to variable $\mu$ is

which is always positive since its discriminant $\Delta(\frac{\partial}{\partial \mu} \varphi_{\sigma}) = 144(\lambda-1)^2 (\lambda- \frac{11}{9}) < 0$. The special value $\varphi_{\sigma}(\lambda, \mu_m) = \frac{32}{1- \lambda} > 0$ can ensure that $\varphi_{\sigma}(\lambda, \mu)$ is a positive function in this case, i.e., $\sigma < 0$ or the equilibrium point $E_*^{(2)}$ is a stable multiple focus with multiplicity one.

Case Ⅱ: $\mu = \mu_1$

When $\mu = \mu_1$, we will show that the nilpotent(double-zero eigenvalue) $E_*^{(2)}$ is a BT cusp of codimension 2. Firstly, by using transformations (Ⅰ): $x = X+x_*^{(2)}$, $y = Y+y_*^{(2)}$ and

in original system

we can derive a new system

With the help of the Lemma $1$ in ^[38], such system (2.13) is equivalent to system

where the discriminants are

and $\varphi_{d_2}(\lambda, \mu) = (\lambda -1)(\lambda +7) \mu^2 + (6 \lambda^2 +16 \lambda +2) \mu +(9 \lambda^2 -6 \lambda +5)$ is an auxiliary function.

From the equations $\varphi_{_{A_2}}(\lambda, \mu) = 0$ and $\varphi_{d_2}(\lambda, \mu) = 0$, the Sylvester's resultant with respect to variable $\mu$ is

which implies $d_2 \neq 0$. On the other hand, for the quadratic function $\varphi_{d_2}(\lambda, \mu)$, which is also an downward opening parabola since the negative coefficient $(\lambda -1)(\lambda +7)$ is in the highest order term, its discriminant is $\Delta(\varphi_{d_2}) = 656 \lambda^2 -224 \lambda +144 > 0$, and the symmetry axis $\mu = \frac{6 \lambda^2 +16 \lambda +2}{2(1-\lambda)(\lambda +7)}$ is between the longitudinal axis $\mu = 0$ and the vertical line $\mu = \mu_m$. By using the special values $\varphi_{d_2}(\lambda, 0) > 0$ and $\varphi_{d_2}(\lambda, \mu_m) = \frac{32 \lambda}{\lambda -1} < 0$, we have $\varphi_{d_2}(\lambda, \mu_1) < 0$ and $d_2 > 0$. Therefore it completes the non-degeneracy condition $d_1 d_2 \neq 0$ (actually $d_1 d_2 < 0$) and the classification work of codimension $2$ cusps in this paper, which is meaningful to Bogdanov-Takens bifurcation analysis of codimension $2$.

Now, we will continually use transformations (Ⅲ): $u = p+ a_{02} pq$, $v = q- a_{20} p^2$; (Ⅳ): $p = w$, $q = z-c_{11} wz$ and (Ⅴ): $w = x_1+ \frac{1}{2} f_{02} x_1^2$, $z = y_1 + f_{02} x_1 y_1$ to derive a standard form with discriminants $d_1$ and $d_2$:

hence it also support above conclusions and we can obtain the Theorem $2.3$.

Theorem 2.3. As we take the value of parameters under the condition (2.8), the system (1.1) admits an interior equilibrium point $E_*^{(2)}$ with zero trace.

(ⅰ) When $\mu_m < \mu < \mu_1$, the equilibrium point $E_*^{(2)}$ is a stable multiple focus with multiplicity one.

(ⅱ) When $\mu = \mu_1$, the equilibrium point $E_*^{(2)}$ is a cusp of codimension $2$ (BT bifurcation point).

(ⅲ) When $\mu > \mu_1$, the equilibrium point $E_*^{(2)}$ is a saddle point.

Here we can take the value $\lambda = \frac{2}{3}$ to get some interesting result, the first positive(meaningful) root of Eq (2.9) is $\mu_1 = \mu_{1, \frac{2}{3}} : = 12+3 \sqrt{17} \approx 24.369317$ and $\mu_m = 9$ by using identities

Solving out the Eq (2.4), we have three possible interior equilibrium point $E_4^{(i)} = (x_4^{(i)}, y_4^{(i)})$ ($i = 1, 2, 3$), where some components are

and auxiliary functions are $\Phi_s(\mu) = \sqrt{(\mu-1)(\mu-3)(\mu-9) \varphi_s(\mu)}$ and $\varphi_s(\mu) = \mu^3-13 \mu^2-153 \mu-603$. With the techniques in Calculus, we know:

(ⅰ) when $\mu \in (\mu_m, \mu_s^+)$, $\mu_s^+ = \frac{13+2 \sqrt{157}}{3} \approx 12.687$, the function $\varphi_s(\mu)$ is negative and monotonically decreasing.

(ⅱ) when $\mu = \mu_s^+$, the function $\varphi_s(\mu)$ owns a negative (local) minimum.

(ⅲ) when $\mu > \mu_s^+$, the function $\varphi_s(\mu)$ is monotonically increasing and has a unique positive root, which is denoted as $\mu_s$, where

Case 1. When $\mu_s < \mu < \mu_1$ or $\mu > \mu_1$, that is $\Phi_s(\mu) > 0$, the system (1.1) has three interior equilibrium point $E_4^{(i)} = (x_4^{(i)}, y_4^{(i)})$ ($i = 1, 2, 3$) due to the following inequalities

Here the equilibrium point $E_4^{(2)}$ is a stable multiple focus with multiplicity one when $\mu < \mu_1$, while it becomes a saddle point when $\mu > \mu_1$.

Case 2. When $\mu = \mu_1$, there exist two interior equilibrium points, including an asymptotically stable node $E_4^{(1)} = (x_4^{(1)}, y_4^{(1)}) = (a(11+3 \sqrt{17}), 3ae)$ since

The second degenerate interior equilibrium point $E_4^{(2)} = (x_4^{(2)}, y_4^{(2)}) = (\frac{1}{4} x_4^{(1)}, \frac{3}{8} (1+ \sqrt{17})ae)$ is a cusp of codimension $2$. Furthermore, it is worth noting that the two equilibrium point $E_4^{(2)}$ and $E_4^{(3)}$ can coincide with each other.

Case 3. When $\mu = \mu_s$ or $\Phi_s(\mu) = 0$, there exist two interior equilibria, that is to say, above two interior equilibrium point $E_4^{(1)}$ and $E_4^{(3)}$ can coincide with each other (but we still denote it as $E_4^{(1)}$), where

It is quite clear that this point $E_4^{(1)}$ is a stable node when $A_1(E_4^{(1)}) = - \frac{128 \alpha e(55 \mu_s^2+324 \mu_s+981)}{9(\mu_s-9)(\mu_s^2-1)(\mu_s^2-9)} < 0$ and $A_2(E_4^{(1)}) = 0$. Seeing ^[35] in detail, the interior equilibrium $E_4^{(2)}$ is still a stable multiple focus with multiplicity one.

Case 4. When $\mu_m < \mu < \mu_s$ or $\Phi_s(\mu) < 0$, there exists a unique stable multiple focus $E_4^{(2)}$ with multiplicity one.

In order to verify the feasibility of theoretical derivation, we take $r_1 = 0.6$ $\alpha = 0.5$, $a = 1.5$ and $e = 0.6$, then some numerical simulations are implemented. Figure 3 depicts the curves of functions $\varphi_{_{A_2}}(\lambda, \mu)$, $\varphi_d(\lambda, \mu)$ and $\psi_d(\lambda, \mu)$, which can show the existence of key values. When $\mu = 20 < \mu_1$, the unique equilibrium point $E_4^{(2)}$ is a stable multiple focus with multiplicity one(see Figure 4(a)). When $\mu = 25 > \mu_1$, there exist three interior equilibrium points including a stable node $E_4^{(1)}$, a saddle point $E_4^{(2)}$ and an unstable node $E_4^{(3)}$ (see Figure 4(b)). When $\mu = \mu_1$, a cusp of codimension $2$ $E_4^{(2)}$ and a stable node $E_4^{(1)}$ will occur (see Figure 4(c)). When $\mu = \mu_s$, there exist two equilibrium points including a stable multiple focus $E_4^{(2)}$ with multiplicity one and a stable node $E_4^{(1)}$ (see Figure 5(a)). When $\mu = 21.5$, there exist three equilibrium points including a stable multiple focus $E_4^{(2)}$ with multiplicity one, a saddle point $E_4^{(3)}$ and a stable node $E_4^{(1)}$ (see Figure 5(b)). In a word, there are several kinds of internal equilibrium points with different characteristics in the system (1.1).

2.2.3.   Non-hyperbolic equilibrium point: Case (C2)

This subsection will show the existence and stability of interior equilibrium point in another special case, which can also ensure potential Hopf bifurcation and Bogdanov-Takens bifurcation. As we take $m_1 = r_1 - \alpha e - 3 m_2$, $m_2 = \lambda \alpha e$ and $K_1 = \frac{a r_1 s}{4 \lambda \alpha e}$ with dimensionless "control variables" $\lambda \in (0, 1)$, $\mu > 1$ and $s = \sqrt{\mu^2-1} > 0$ for later use, then an interior equilibrium point $E_*$, the determinant $A_2(E_*)$ and the confined positive parameter $d$ are listed as follows:

where the mentioned auxiliary functions $\varphi_{y_*} = \varphi_{y_*}(\lambda, \mu)$, $\varphi_d = \varphi_d(\lambda, \mu)$, $\psi_d = \psi_d(\lambda, \mu)$, $\varphi_{A_2} = \varphi_{A_2}(\lambda, \mu)$ are:

The inequality $d > 0$ or $\varphi_d < 0$ deduces a lower bound $\mu_m = \frac{5 \lambda^2+2 \lambda+1}{(1- \lambda)(3 \lambda+1)}$. The generalized expression of the first positive root $\mu_1 = \mu_1(\lambda)$ in equation $\varphi_{A_2}(\lambda, \mu) = 0$ is

where

To reveal complexity in this special case, taking a fascinating value $\lambda = \frac{5}{9}$ with identities

and denoting this case as (C2), we have a threshold $\mu_1 = \mu_{1, \frac{5}{9}} = \frac{2419}{240} + \frac{61}{240} \sqrt{1921}$ from equation $\varphi_{A_2}(\frac{5}{9}, \mu) = 0$. Therefore, we conversely discuss following two cases and the above mentioned interior equilibrium point $E_*$(denoted as $E_5^{(2)}$) could be a multiple focus or center with multiplicity one when $\mu \in (\mu_m, \mu_1)$. Following the Eq (2.11) and steps in above subsection, the first Lyapunov coefficient is

where $\psi_{\sigma}(\mu) = 14 \mu^3+(14s+47) \mu^2+(47s-24) \mu+91s-37 > 0$, and all coefficients in auxiliary function $\varphi_{\sigma}(\mu)$ is listed in the Appendix for completeness. If $\mu < \mu_{\sigma}(>\mu_{\sigma})$ or $\sigma < 0(>0)$, the equilibrium point $E_5^{(2)}$ is a stable(unstable) multiple focus with multiplicity one, the system (1.1) undergoes a non-degenerate supercritical(subcritical) Hopf bifurcation around $E_5^{(2)}$, and limit cycles generated by this critical point are stable(unstable). On occasion, there may exist some parameter values such that $\sigma = 0$ or the system (1.1) may undergoes a degenerate Hopf bifurcation for some values of parameters ^[36]. Accompanying with the Calculus, Figure 6(a) is the curve of function $\varphi_{\sigma}(\mu)$, which is used to guess the unique positive root of equation $\varphi_{\sigma}(\mu) = 0$:

At $\mu = \mu_{\sigma}$, by using successor function method, the second focal(Lyapunov) quantity $g_5 \approx - \frac{0.00003569 \alpha^4 e^4}{a^4} < 0$ ensures that the equilibrium point $E_5^{(2)}$ is a stable weak focus of order 2 ^[36]. The Bautin bifurcation(generalized Hopf bifurcation) may occur.

When $\mu = \mu_{1, \frac{5}{9}}$, the system (1.1) owns two separate interior equilibrium points $E_5^{(1)} = (\frac{(151+9 \sqrt{1921})a}{200}, \frac{4(56+ \sqrt{1921})ae}{45})$ and $E_5^{(2)} = (\frac{(39+\sqrt{1921})a}{8}, \frac{(709+17 \sqrt{1921})ae}{81})$. The equilibrium point $E_5^{(1)}$ is an unstable node since

Now we construct a linear transformation (Ⅱ): $u = \frac{-36 X}{(\sqrt{1921}-11) \alpha}$, $v =$ $\frac{-\sqrt{1921}+31}{\sqrt{1921}-11}eX +Y$, then the system (1.1) is apparently equivalent to a new system in the standard form of Eq (2.14) with discriminants

The equilibrium point $E_5^{(2)}$ is just a cusp of codimension 2 due to the non-degeneracy condition $d_1 d_2 \neq 0$. Hence we can obtain Theorem 2.4.

Theorem 2.4. Under the conditions of case (C2), (ⅰ) when $\mu \in (\mu_m, \mu_{1, \frac{5}{9}}) \backslash \{ \mu_{\sigma} \}$, the equilibrium point $E_5^{(2)}$ is a multiply stable(unstable) focus with multiplicity one if $\mu < \mu_{\sigma}(>\mu_{\sigma})$; (ⅱ) when $\mu = \mu_{\sigma}$, the equilibrium point $E_5^{(2)}$ is a stable weak focus of order $2$; (ⅲ) when $\mu = \mu_{1, \frac{5}{9}}$, the equilibrium point $E_5^{(2)}$ is a cusp of codimension 2, and the point $E_5^{(1)}$ is an unstable node; (iv) when $\mu > \mu_{1, \frac{5}{9}}$, the equilibrium point $E_5^{(2)}$ is a saddle point; (v) the equilibrium point $E_5^{(1)}$ is just an unstable node.

In order to verify the feasibility of theoretical derivation, we will give some numerical simulations. Figure 6(b) is the curves of functions $\varphi_{A_2}$(red) and $\varphi_d$(blue) with $\lambda = \frac{5}{9}$, which mainly displays the threshold value of control parameters. And then, we take $r_1 = 1$, $a = 1.5$, $\alpha = 0.5$ and $e = 0.6$, the system (1.1) exists a multiple stable focus with multiplicity one (see Figure 7(a)), a unique multiple unstable focus with multiplicity one and a limit cycle (see Figure 7(b)), a cusp of codimension $2$ and an unstable node (see Figure 7(c)), a unique stable weak focus of order $2$ (see Figure 7(d)). In a word, the system (1.1) has different internal equilibrium points with the value change of key parameters.

2.3.   Cusp of codimension 3

In the following, we will investigate a cusp of codimension $3$ in the system (1.1). First of all, translating the equilibrium point $E_* = (x_*, y_*)$ to the origin $O$ via a transformation (Ⅰ): $x = X+x_*$, $y = Y+y_*$, we obtain

Next, following the technique in above subsection and making a linear transformation (Ⅱ):

the above system can be rewritten as the form

From the Lemma $1$ in ^[38], the system (2.24) is equivalent to the system in standard form

after some nonsingular transformations in the neighborhood of $O$, where $d_1 = b_{20}$, $d_2 = b_{11}+2 a_{20}$ are discriminates. Solving out an degenerate equilibrium $E_3 = (\frac{a(2 \alpha e+m_2)}{\alpha e-m_2}, \frac{12 e^2 \alpha a}{2 \alpha e+m_2})$ from the condition $d_1 = 0$, which also satisfies $A_1 = 0$, $A_2 = 0$ and $p_x = \Delta_x = p_y = \Delta_y = 0$, thus we have parameters $r_1$, $K_1$ and $d$ (suppose they are positive) with restrictions

These restrictions (2.26) can deduce another discriminate $d_2 = \frac{(8 \alpha e+7 m_2)(m_2 - \alpha e)}{18a \alpha e^2} \neq 0$, i.e., the degeneracy condition $d_1 d_2 = 0$, which also yields that the equilibrium point $E_3$ is a cusp of codimension at least 3. Indeed, the degenerate equilibrium point $E_3$ is a codimension 3 Bogdanov-Takens singularity(focus or center) after some nonsingular transformations. Finally, the existence of the equilibrium point $E_3$ can be seen from Figure 8 in details with parameters $a = 1.5$, $\alpha = 0.5$, $m_1 = 0.6$, $m_2 = 0.06$ and $e = 0.6$. At the same time, we can obtain the Theorem 2.5.

Theorem 2.5. The the degenerate equilibrium point $E_3$ with conditions 2.26 is a codimension 3 Bogdanov-Takens singularity(focus or center).

3.   Bifurcations analysis

In this section, for the interior equilibrium point $E_4^{(2)}$ and $E_5^{(2)}$, we mainly concentrate on Hopf bifurcation curve when $\mu \in (\mu_m, \mu_1)$ and codimension 2 Bogdanov-Takens bifurcation when $\mu = \mu_1$, respectively.

3.1.   Hopf bifurcation

In the case of (C1), we firstly discuss the existence of Hopf bifurcation curve when $\mu_m < \mu < \mu_{1, \frac{2}{3}}$. We will choose $m_2$ and $d$ as Hopf bifurcation controlling parameters and consider the following perturbed system

where $\delta = (\delta_1, \delta_2)$ is a sufficiently small parameter vector in a neighbourhood of the origin $O = (0, 0)$ in the parameter plane. Letting $\delta \neq (0, 0)$, we suppose the equilibrium point $E_*$ as $(x_*, y_*)$, where $x_* = x_4^{(2)}+w$, $w$ is a sufficiently small variable and component

Substituting it into $A_1$ and $A_2$, we have the solution $\delta_1 = \delta_1(w)$, $\delta_2 = \delta_2(w)$, where

At this time, the approximation of the required Hopf bifurcation curve $Hp$ in a small neighbourhood of the origin is a straight line with slope

where $\varphi_k(\mu) = \mu^4-12 \mu^3+82 \mu^2+12 \mu+45$. Noticing that

so the function $\varphi_k(\mu)$ is monotonically increasing and positive when $\mu > \mu_m$, or $k_{Hp} < 0$. Hence we can rewrite the determinant as $A_2 = - \frac{\alpha^2 e^2 [a(\mu-1)+4w]}{9 a^2 (\mu^2-1)^2 (a \mu+3a+4 w)^2} \phi_{_{A_2}}(\mu, w)$, where an auxiliary function is

Thus, the bifurcation curve $Hp$ of the system (3.1) at the equilibrium point $E_4^{(2)}$ is analytically defined by the solution (3.2), the variables $\mu$ and $w$ can ensure the existence of the equilibrium point $E_4^{(2)}$ and $A_2 > 0$, and we can obtain the Theorem 3.1.

Theorem 3.1. (Hopf bifurcation curve) For the equilibrium point $E_4^{(2)}$ with $\mu_m < \mu < \mu_{1, \frac{2}{3}}$($A_1 = 0$ and $A_2 > 0$), when parameter $\delta$ varies in a small neighbourhood of the origin in parameter plane, the Hopf bifurcation curve of the system (3.1) is defined by (3.2) (notice the range of parameter $w$) and the approximation is a straight line with slope $k$ in a small neighbourhood of the origin. Furthermore, the curve divides a small neighbourhood of the origin in the parameter plane into two regions Ⅰ and Ⅱ, in which dynamical behaviors of the system (3.1) can be exhibited.

For (C2), starting from a perturbed system (3.1) with $\mu \in (\mu_m, \mu_{1, \frac{5}{9}}) \backslash \{ \mu_{\sigma} \}$ and $y_* = \frac{4e(6as-5 x_*)(a+x_*)}{9as}$. Solution $\delta_1 = \delta_1(w)$, $\delta_2 = \delta_2(w)$ are

Therefore, the slope $k_{Hp}(\mu)$ of the approximate straight line of the Hopf bifurcation curve Hp at $O$ is

At the same time, we can obtain the Theorem 3.2.

Theorem 3.2. (Hopf bifurcation curve) For the equilibrium point $E_5^{(2)}$ with $\mu \in (\mu_m, \mu_{1, \frac{5}{9}}) \backslash \{ \mu_{\sigma} \}$, when parameter $\delta$ varies in a small neighbourhood of the origin in parameter plane, the Hopf bifurcation curve of the system (3.1) is defined by (3.4) (notice the range of parameter $w$) and the approximation is a straight line with slope $k$ in a small neighbourhood of the origin. Furthermore, the curve divides a small neighbourhood of the origin in the parameter plane into two regions Ⅰ and Ⅱ, in which dynamical behaviors of the system (3.1) can be exhibited.

In order to verify the feasibility of the Theorems 3.1 and 3.2, we will give some numerical simulations. For the equilibrium point $E_4^{(2)}$ with $r_1 = 0.6$ $\alpha = 0.5$, $a = 1.5$, $e = 0.6$ and $\mu = 10$, the Hopf bifurcation curve with $w$ in parameter plane is

and the straight line of the approximate representation of bifurcation curve $Hp$ is $\delta_2 \approx -0.704575 \delta_1$. Figure 9(a) depicts the Hopf bifurcation curve, Figure 10(a), (b) depict phase diagrams of asymptotically stable focus and unstable focus(with a limit cycle) corresponding to $\delta_1 = \delta_2 = 0.001$ (in region Ⅰ) and $\delta_1 = \delta_2 = -0.001$(in region Ⅱ) respectively. For the equilibrium point $E_5^{(2)}$ with $r_1 = 0.6$ $\alpha = 0.5$, $a = 1.5$, $e = 0.6$ and $\mu = 8$ in the case (C2), the bifurcation curve Hp is analytically formulated by such solutions of the system (3.1), i.e., $H_{p} = \{\delta \mid \delta \, {\rm{ satisfies }} \, (3.1) \}$, which can be seen from Figure 9(b), and a Hopf bifurcation curve Hp(red) is defined by

and accompanied by its corresponding slope $k_{Hp}$(dashed blue line) at $O$. Furthermore, we will note

then the curve Hp divides a parameter plane into separate two regions. Figure 11(a), (b) present phase diagrams of a stable node and an unstable focus with respect to $\delta_1 = 0.0001$, $\delta_2 = 0.00001$ and $\delta_1 = -0.0001$, $\delta_2 = -0.00001$ when $\mu = 8$, respectively. Figure 11(c), (d) present phase diagrams with same values of parameters $\delta_1$ and $\delta_2$ when $\mu = 10$. In a word, it is obvious to see from the numerical simulation works that the Hopf bifurcation can occur for the equilibrium point $E_4^{(2)}$ and $E_5^{(2)}$, which also indirectly proves the validity of the theoretical derivation.

3.2.   Bogdanov-Takens bifurcation of codimension 2

We firstly recall the system (1.1) with a cusp $E_4^{(2)}$ of codimension 2 when parameter conditions $\lambda = \frac{2}{3}$ and $\mu = \mu_{1, \frac{2}{3}}$ hold, in other words, we can start with an unfolding system

by introducing bifurcation parameters $e$ and $m_2$. Naturally, a parameter vector $\delta = (\delta_1, \delta_2)$ is in a sufficiently small neighbourhood of the origin $O$ in the parameter plane. By using the transformation (Ⅰ): $x = X+ x_4^{(2)}$, $y = Y+ y_4^{(2)}$ and expanding such system in a power series around the origin, it can be rewritten as

where $b_{00}(0, 0) = 0$. Secondly, we will use a transformation (Ⅱ): $u = X$, $v = F_1(X, Y)$, the system (3.9) can be reduced to a new system

where $d_{ij}(\delta)$ can be expressed by coefficients $a_{ij}$ and $b_{ij}(\delta)$, $d_{ij}(0, 0) = 0$ ($i+j \leqslant 1$). Thirdly, making a transformation (Ⅲ): $p = u+ \frac{d_{01}(\delta)}{d_{11}(\delta)}$, $q = v$ since $d_{11}(0, 0) = \frac{(5 \sqrt{17}-21) \alpha e}{12a} \neq 0$, we have

where $f_{ij}(\delta)$ can be expressed by $d_{ij}(\delta)$, but we will omit them here. At the same tome, there exist $f_{01}(\delta) = 0$ and $f_{20}(0, 0) = \frac{(161-39 \sqrt{17}) \alpha^2 e^2}{54a} > 0$. Then we will construct a transformation (Ⅳ): $w = p$, $z = (1-f_{02}(\delta) p)q$, $dt = (1-f_{02}(\delta) p) d \tau$, one can rewrite above system as(the symbol $\tau$ is denoted as $t$)

Similarly, we omit the expressions of coefficients $h_{ij}(\delta)$ although they are expressed iteratively by $f_{ij}(\delta)$, and $h_{20}(0, 0) = f_{20}(0, 0) > 0$, $h_{11}(0, 0) = d_{11}(0, 0) < 0$. That is to say, there is a small neighbourhood of the origin such that $h_{20}(\delta)$ is positive and $h_{11}(\delta)$ is negative when $\delta$ falls in this neighbourhood. Finally, the transformation (Ⅴ): $m = \frac{h_{11}(\delta)^2}{h_{20}(\delta)} w$, $n = \frac{h_{11}(\delta)^3}{h_{20}(\delta)^2} z$, $dt = \frac{h_{11}(\delta)}{h_{20}(\delta)} d \tau$ converts from above system to a generic normal form of Bogdanov-Takens bifurcation

where the symbol $\tau$ is still denoted as $t$, and two discriminants are

This system (3.8) is indeed a generic family unfolding at the codimension 2 cusp $E_4^{(2)}$ according to the rank of a Jacobian matrix or the nonzero Jacobian determinant

Therefore we obtain local approximated representations of saddle-node (SN), Hopf (H) and homoclinic (HL) bifurcation curves up to second-order with slope $k_{BT} = \frac{(1+\sqrt{17}) \alpha}{6} \approx 0.853851 \alpha > 0$ around the origin for the system (3.8) ^[39]. These bifurcation curves can divide the parameter plane into several regions, which can exhibit separately dynamical behaviors.

(ⅰ) The saddle-node bifurcation curve is formulated by

(ⅱ) The Hopf bifurcation curve is formulated by

(ⅲ) The homoclinic bifurcation curve is formulated by

Thus, we can obtain the Theorem 3.3.

Theorem 3.3. (Bogdanov-Takens bifurcation of codimension 2). For the unfolding system (3.8) with bifurcation parameters $e$ and $m_2$, in a small neighbourhood of the equilibrium point $E_4^{(2)}$, the system undergoes an attracting Bogdanov-Takens bifurcation of codimension 2 when the value of parameter $\delta$ varies in such sufficiently small neighbourhood of the origin. Furthermore, this system is a generic family unfolding at the cusp $E_4^{(2)}$ of codimension 2.

Here we will take $r_1 = 0.6$, $a = 1.5$, $\alpha = 0.5$ and $e = 0.6$, then Figure 12 depicts the saddle-node, Hopf and homoclinic bifurcation curves in different colors, which can show the existence of critical thresholds.

(ⅰ) When $\delta_1 = \delta_2 = 0$, it is evident from the Theorem 2.3 that there exist two interior equilibrium points, including an asymptotically stable node $E_4^{(1)}$ and a Bogdanov-Takens cusp of codimension 2 $E_4^{(2)}$.

(ⅱ) When $\delta_1 = 0.01$, $\delta_2 = 0$ ($\delta$ lies in positive $\delta_1$ axis) or $\delta$ falls in region Ⅰ (the region between saddle-node bifurcation curve $SN_2$ and homoclinic bifurcation curve), there exist three interior equilibrium points, where two interior equilibrium points are bifurcated from a stable node in (viii), which can be seen from Figure 13(a), (b).

(ⅲ) When $(\delta_1, \delta_2) \approx (0.01, 4.265486 \times 10^{-3})$ or $\delta$ lies in homoclinic bifurcation curve, there exist three interior equilibrium points and a homoclinic loop.

(ⅳ) When $(\delta_1, \delta_2) \approx (0.01, 4.267370 \times 10^{-3})$ or $\delta$ falls in region Ⅱ(the region between homoclinic bifurcation curve and Hopf bifurcation curve), there exist three interior equilibrium points, including an unstable focus, a saddle point and a stable node, which can be seen from Figure 14(a), (b).

(ⅴ) When $(\delta_1, \delta_2) \approx (0.01, 4.269255 \times 10^{-3})$ or $\delta$ lies in Hopf bifurcation curve, there exist three interior equilibrium points, including a saddle point and a stable node.

(ⅵ) When $(\delta_1, \delta_2) \approx (0.01, 4.271301 \times 10^{-3})$ or $\delta$ falls in region Ⅲ(the region between Hopf bifurcation curve and saddle-node bifurcation curve $SN_1$), there exist three interior equilibrium points, including a stable focus which is unstable in case (ⅳ), which can be seen from Figure 15(a). However, by combining the case (ⅳ), it can ensure the potential Hopf bifurcation, but the homoclinic loop is broken.

(ⅶ) When $(\delta_1, \delta_2)$ lies in saddle-node bifurcation curve, there exist two interior equilibrium points.

(ⅷ) When $(\delta_1, \delta_2) \approx (0.01, 8.546695 \times 10^{-3})$ or $\delta$ falls in region Ⅳ(the region on the left hand side of saddle-node bifurcation curve), there exists a unique stable node, which can be seen from Figure 15(b).

However, if we choose $m_2$ and $d$ as Bogdanov-Takens bifurcation parameters and rewrite the original system in the unfolding form (3.1) with $\mu = \mu_{1, \frac{2}{3}}$:

where $\delta_1$ and $\delta_2$ are sufficiently small parameters and the vector $\delta = (\delta_1, \delta_2)$ is in a small neighbourhood of the origin $O$ as well. Following the procedures above and the values of parameters, the unfolding system (3.1) is also a generic family unfolding at the codimension 2 Bogdanov-Takens cusp $E_4^{(2)}$ according to the nonzero Jacobian:

and we have representations of saddle-node, Hopf and homoclinic bifurcation curves with slope $k_{BT} = \frac{1- \sqrt{17}}{6ae} \approx - \frac{0.520518}{ae} < 0$ around the origin. Furthermore, it should be noticed that the slope $k$ can be viewed as the limiting case of the slope (3.3) when $\mu \rightarrow \mu_{1, \frac{2}{3}}$. At the same time, when $\delta$ lies on the Hopf bifurcation curve $H$, for instance, $(\delta_1, \delta_2) \approx (1 \times 10^{-5}, -5.783264 \times 10^{-6})$, the Hopf bifurcation can not undergo; when $\delta$ lies on the homoclinic bifurcation curve $HL$, for instance, $(\delta_1, \delta_2) \approx (1 \times 10^{-5}, -5.783278 \times 10^{-6})$, the homoclinic loop does not exist. Moreover, it should be noticed more that these two cases both occur owing to $d_2 > 0$ or the minus of (3.2). While the saddle-node bifurcation curve up to second-order can be formulated by

For the saddle-node bifurcation curve, when $\delta$ lies on the region Ⅰ (the left hande side of the SN curve), there exist three interior equilibrium points. When $\delta$ lies on the region Ⅱ (the right hand side of the SN curve), there exists a unique interior equilibrium point. When $\delta$ lies on the saddle-node bifurcation curve $SN_1$, there exist three interior equilibrium points. When $\delta$ lies on the saddle-node bifurcation curve $SN_2$, there exists a unique interior equilibrium point. All the detailed results can be seen in the Figure 16 for saddle-node bifurcation curve in this novel phenomenon with $r_1 = 0.6$, $a = 1.5$, $\alpha = 0.5$ and $e = 0.6$.

Similarly, for the case (C2) with $\lambda = \frac{5}{9}$ and Bogdanov-Takens bifurcation parameters $m_2$ and $d$, according to a Jacobian matrix with rank 2 and

local approximated representations of saddle-node (SN), Hopf (H) and homoclinic (HL) bifurcation curves up to second order with slope $k_{BT} = \frac{709 - 17 \sqrt{1921}}{648 ae} < 0$ at $O$ are obtained rapidly. It is also the limitation $\lim \limits _{\mu \rightarrow \mu_{1, \frac{5}{9}}} k_{Hp}(\mu)$.

(ⅰ) The saddle-node bifurcation curve is formulated by

where the half curves $SN_1$($SN_2$) is the "right" ("left") part of curve $SN$ in the forth(second) quadrant, respectively.

(ⅱ) The Hopf bifurcation curve is formulated by

(ⅲ) The homoclinic bifurcation curve is formulated by

Thus, we can obtain the Theorem 3.4.

Theorem 3.4. (Bogdanov-Takens bifurcation of codimension 2) From the unfolding system (3.17) with bifurcation parameters $m_2$ and $d$ in the case (C2), this system undergoes an Bogdanov-Takens bifurcation of codimension 2 when $\delta$ varies in a sufficiently small neighbourhood of the origin. Furthermore, the system is a generic family unfolding at the cusp $E_5^{(2)}$ of codimension 2 as well.

For numerical simulation, we take $r_1 = 1$, $\alpha = 0.5$, $a = 1.5$ and $e = 0.6$, then Figure 17 depicts saddle-node (red), Hopf (green) and homoclinic (blue) bifurcation curves in different colors, which can show the existence of critical thresholds..

(ⅰ) When $\delta_1 = \delta_2 = 0$, it is evident that there exist two interior equilibrium points, including an unstable node $E_5^{(1)}$ and a codimension 2 cusp $E_5^{(2)}$.

(ⅱ) When $(\delta_1, \delta_2) \approx (0.001, -3.095598 \times 10^{-5})$ or $\delta$ falls in region Ⅰ (the region below saddle-node bifurcation curve SN), there exists a unique unstable node.

(ⅲ) When $(\delta_1, \delta_2) \approx (0.001, -6.188787 \times 10^{-5})$ or $\delta$ falls in region Ⅱ (the region between Hopf bifurcation curve and saddle-node bifurcation curve $SN_1$), there exist three interior equilibrium points, including an unstable node, a saddle point and an unstable focus.

(ⅳ) When $\delta$ lies in Hopf bifurcation curve, there exist three interior equilibrium points, including an unstable node, a saddle point, and a non-hyperbolic equilibrium (focus or center) with zero trace and positive determinant, which can ensure potential Hopf bifurcation.

(ⅴ) When $(\delta_1, \delta_2) \approx (0.001, -6.184072 \times 10^{-5})$ or $\delta$ falls in region Ⅲ (the region between homoclinic bifurcation curve and Hopf bifurcation curve), there exist three interior equilibrium points, including a stable focus which is unstable in case (ⅲ).

(ⅵ) When $\delta$ lies in homoclinic bifurcation curve, there exist three interior equilibrium points and a homoclinic loop, including an unstable node, a saddle point and a stable focus.

(ⅶ) When $(\delta_1, \delta_2) = (0.001, -3.090882 \times 10^{-5})$ or $\delta$ falls in region Ⅳ (the region between saddle-node bifurcation curve $SN_2$ and homoclinic bifurcation curve), there still exist three interior equilibrium points, including an unstable node, a saddle point and a stable focus.

(ⅷ) When $(\delta_1, \delta_2)$ lies in saddle-node bifurcation curve $SN_1$, there exists a unique unstable node.

(ⅸ) When $(\delta_1, \delta_2)$ lies in saddle-node bifurcation curve $SN_2$, there exist three interior equilibrium points.

4.   Limit cycles via perturbation procedure

Now we especially study the limit cycle generated by Hopf bifurcation, they showed the efficiency of the perturbation method by using a topological polynomial version of the classical Rosenzweig-MacArthur (R-M) predator-prey model in the paper ^[40]. In this section, we focus on the approximate calculation of limit cycles in the original predator-prey system (3.1) via a perturbation procedure and canonical transformation, which can be used to determine the limit cycles and their associated frequencies in general two-dimensional systems. Comparing it with the Lindstedt-Poincare (LP) method, the method can give accurate results, while the LP method is limited to weakly nonlinear systems, although it is simple and is frequently used as an algorithm to approximate steady-state periodic solutions in nonlinear oscillators ^[41]. Recalling an unfolding system mentioned in Subsection 3.1 with $\mu \in (\mu_m, \mu_{1, \frac{5}{9}})$, $\sigma < 0$ and a sufficiently small parameter vector $(\delta_1, \delta_2)$ in a neighbourhood of the origin $O$ in the parameter plane:

when $\mu = 8$, the equilibrium points $E_*$ is a multiple stable focus with multiplicity one, and the corresponding non-degenerate Hopf bifurcation is supercritical. According to the Figure 10(b), a limit cycle exists when sufficiently small $(\delta_1, \delta_2)$ falls in region Ⅱ. In this perturbation procedure, we firstly transfer the equilibrium point $E_*$ to the origin $O$ by using a linear transformation (Ⅰ): $x = X+x_*$, $y = Y+y_*$ and obtain a new system. Secondly, we construct a nonsingular transformation (Ⅱ): $X = \eta$, $Y = \frac{\sqrt{-J_{11}^2+2 J_{11} J_{22}-4 J_{12} J_{21}-J_{22}^2} \xi -(J_{11}-J_{22}) \eta}{2 J_{12}}$ such that the system (3.1) has its real Jordan's form:

where $\mathcal{A} = \frac{1}{2} A_1$, $\mathcal{B} = \frac{1}{2} \sqrt{- \Delta_*}$, $P_1(\xi, \eta) = \sum \limits _{i+j = 2} ^{\infty} \widetilde{a}_{ij} \xi^i \eta^j$ and $Q_1(\xi, \eta) = \sum \limits _{i+j = 2} ^{\infty} \widetilde{b}_{ij} \xi^i \eta^j$ are analytical functions. Obviously, the limit cycle is enclosing an unstable hyperbolic focus or non-hyperbolic node. Introducing a dimensionless time scale transformation $\tau = \omega t$ with frequency $\omega$ of a limit cycle, above system can be written as a canonical system:

where $\Omega = \frac{\omega}{\mathcal{B}}$, $\delta = \frac{\mathcal{A}}{\mathcal{B}}$, $P_2 = \frac{P_1}{\mathcal{B}}$ and $Q_2 = \frac{Q_1}{\mathcal{B}}$. Now we suppose that there exist series

in which $\epsilon = \sqrt{\frac{\delta \Omega_1}{1+ \delta \Omega_1}}$. Substituting them into Eq (4.3) and noticing the series expansions of functions in the right hand side, we recursively derive following coupled first-order differential equations of $\xi_n$ and $\eta_n$ in all orders of $\epsilon$:

and so on, where all involved coefficients are defined by $a_{ij} = \frac{1}{\mathcal{B}} \widetilde{a}_{ij}$, $b_{ij} = \frac{1}{\mathcal{B}} \widetilde{b}_{ij}$.

To illustrate the procedure process, we mainly concentrate on the values of parameters from Figure 10 with small parameters $\delta_1 = -0.0001$ and $\delta_2 = 0$. The steady-state solutions in order $\epsilon^1$ are $\xi_1 = -A \sin(\tau)$, $\eta_1 = A \cos(\tau)$, where $A$ is a constant to be determined. Then the straightforward solutions of $\xi_2$ and $\eta_2$ are

From the Eq (4.4), we have a second-order differential equation of $\eta_3$:

where $C_{31}$, $S_{31}$, $C_{33}$, $S_{33}$ are some constants. Letting $C_{31} = S_{31} = 0$ or eliminating the secular terms in this equation, we know $\gamma_2 \approx -0.011293 A^2-1$, $\Omega_1 \approx \frac{274.507342}{A^2}$. Hence we accordingly derive solutions

Similarly, based on the previous process, we obtain a second-order differential equation of $\eta_4$ in the form of

Here the coefficient $C_{41}$ yields $\gamma_3 = 0$, thus the solutions of $\eta_4$ and $\xi_4$ read

Repeating above mentioned steps, we iteratively and formally derive required constants $A \approx 4.202827$, $\gamma_4 \approx 0.206668$, $\gamma_5 = 0$, $\gamma_6 \approx -0.014240$, $\gamma_7 = 0$, $\cdots$ from coefficients $C_{51}$, $S_{51}$, $C_{61}$, $C_{71}$, $C_{81}$, $\cdots$ and solutions $\xi_n(\tau)$, $\eta_n(\tau)$($n = 5, 6, 7, 8$), which are listed in the Appendix. A generalized second-order ODE of $\eta_n$ reads

and the corresponding formal solutions via Leonhard Euler's method is

Here $C_{n, k}(\tau)$, $S_{n, k}(\tau)$, $\widetilde{C}_{n, k}(\tau)$, $\widetilde{S}_{n, k}(\tau)$, $\overline{C}_{n, k}(\tau)$ and $\overline{S}_{n, k}(\tau)$ are undetermined polynomials of variable $\tau$ with $C_{n, k}(0) : = C_{n, k}$, $S_{n, k}(0) : = S_{n, k}$, $S_{n, 0}(\tau) = \widetilde{S}_{n, 0}(\tau) = \overline{S}_{n, 0}(\tau) : = 0$; $n = 0, 1, 2, \cdots$, $k = 0, 1, 2, \cdots, n$. Finally, a $N$-order approximate solution to the limit cycle of the Eq (3.1) is

where

Furthermore, up to eight-order approximation of solutions $\xi(t)$ and $\eta(t)$, for comparison, Figure 18(a), (b) depict limit cycles via the Runge-Kutta 45 method(red) and perturbation procedure(blue) with values $\mu = 8$, $\delta_1 = -0.0001$, $\delta_2 = 0$ and $\mu = 5$, $\delta_1 = -0.001$, $\delta_2 = 0$, respectively. In the first figure, the red and blue curves almost coincide with each other. For the latter option, some main invariants are calculated as $A \approx 10.010405$, $\Omega_1 \approx 0.5067015052$, $\gamma_2 \approx -3.031960310$, $\gamma_3 = 0$, $\gamma_4 \approx -6.244350476$, $\gamma_5 = 0$, $\gamma_6 \approx -29.58780360$, $\gamma_7 \approx 183.5476191$, $\gamma_8 \approx 159.9270824$.

5.   Conclusions

Within the framework of predator-prey ecological dynamics, this paper mainly discusses the dynamic properties of the Bazykin's predator-prey ecosystem with Holling type Ⅱ functional response and interspecific density-restricted effects on the predators, including the existence and stability of all possible equilibrium points, Hopf bifurcation and Bogdanov-Takens bifurcation.

Aiming at all possible equilibrium points of the Bazykin's predator-prey ecosystem, mathematical theory works mainly investigate the existence and stability of boundary equilibrium point, hyperbolic equilibrium point, non-hyperbolic equilibrium point and cusp of condimension 3, and then give some corresponding threshold conditions of some key parameters, for example, by introducing control variables $\lambda$ and $\mu$, the stability analysis and type of the interior equilibrium point $E_*^{(2)}$($\ast = 4, 5$) are ascertainable and clear in detail, and this equilibrium point $E_*^{(2)}$ is a multiple focus with multiplicity one(or cusp of codimension 2 or saddle point) when $\mu_m < \mu < \mu_1$(or $\mu = \mu_1$ or $\mu > \mu_1$). At the same time, it is easy to find from numerical simulation works that all the equilibrium points derived from theoretical derivation always exist and have corresponding stability states, which indirectly prove that the Bazykin's predator-prey ecosystem has different intrinsic dynamic properties with the value change of critical parameters, which also represents that predator population and prey population have different coexistence modes.

For Hopf bifurcation and Bogdanov-Takens bifurcation, mathematical theory works mainly investigate Hopf bifurcation and Bogdanov-Takens bifurcation about the equilibrium point $E_{4}^{(2)}$ and $E_{5}^{(2)}$, and give some threshold conditions of some key parameters to ensure the occurrence of Hopf bifurcation and Bogdanov-Takens bifurcation. With the help of numerical simulation, the formulated Hopf bifurcation curve when $\mu_m < \mu < \mu_1$ and Bogdanov-Takens bifurcation of codimension 2 when $\mu = \mu_1$ are both presented, which can directly verify the validity and feasibility of theoretical derivation, and indirectly explain that the Bazykin's predator-prey ecosystem has complex bifurcation dynamic evolution process with the value change of critical parameters, such as saddle-node, Hopf and homoclinic bifurcation. Furthermore, it is obvious to find that the Bazykin's predator-prey ecosystem can exhibit different dynamical behaviours when parameter vector $(\delta_1, \delta_2)$ varies in different regions in a small neighbourhood of the origin $O$ in the parameter plane. However, it is also worth noting that the Hopf and homoclinic bifurcations do not exist if we choose $m_2$ and $d$ as bifurcation parameters in (C1). Moreover, we specifically investigate the limit problem by comparing the Runge-Kutta 45 method with perturbation procedure.

In the follow-up research works, based on the research results of this paper and biological manipulation theory, we will further explore the dynamic relationship between Microcystis aeruginosa and filter-feeding fish by using bifurcation dynamic analysis and explain the biological significance of the Bazykin's predator-prey ecosystem. At the same time, due to the difference of nutrient load, biological composition and hydrodynamic conditions in different water bodies, the dynamic relationship between filter feeding fish and Microcystis aeruginosa on the basic of the Bazykin's predator-prey ecosystem still needs to be further studied in the experiment.

In summary, all results in this paper further show that the system (1.1) proposed by the paper ^[24] has more abundant dynamic behavior, and vigorously develop the dynamic properties of the Bazykin's predator-prey ecosystem. Furthermore, the results of bifurcation and stability can be helpful to better understand the interaction mechanism between prey population and predator population in natural real ecosystem.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant no.31570364 and no.61871293) and the National Key Research and Development Program of China (Grant no. 2018YFE0103700). We thank teachers Min Zhao and Chuanjun Dai in Wenzhou University. Author S. T. Wang thank teachers Maotong Zhang and Shujing Zhao in Liushi No.3 Middle School. We also thank the Editor and Referee(s).

Conflict of interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Appendix

All coefficients in function $\varphi_{\sigma}(\mu) = \sum \limits _{i = 0} ^{26} a_i \mu^i$ are:

The required approximate expression of solutions $\xi_n(\tau)$, $\eta_n(\tau)$($n = 5, 6, 7, 8$) in Section 4 are: