ABSTRACT

In this paper we analytically and numerically consider the dynamical behavior of a certain predator-prey system with Holling type II functional response, including local and global stability analysis, existence of limit cycles, transcritical and Hopf bifurcations. Mathematical theory derivation mainly focuses on the existence and stability of equilibrium point as well as threshold conditions for transcritical and Hopf bifurcation, which can in turn provide a theoretical support for numerical simulation. Numerical analysis indicates that theoretical derivation results are correct and feasible. In addition, it is successful to show that the dynamical behavior of this predator-prey system mainly depends on some critical parameters and mathematical relationships. All these results are expected to be meaningful in the study of the dynamic complexity of predatory ecosystem.

Keywords:

Predator-Prey System, Holling Type II Functional Response, Global Stability, Bifurcation Analysis

1. Introduction

In 1965, C.S. Holling proposed three kinds of functional responses for different kinds of species to model the phenomena of predation, which made the traditional non-autonomous Lotka-Volterra predator-prey system more realistic [1] [2] [3] [4]. These functional responses describe how predators transform the harvested prey into the growth of itself and were discussed by numbers of researchers, including the stability of equilibrium points, existence of Hopf bifurcation, limit cycles, homoclinic loops, and even catastrophe [5].

For the Rosenzweig-MacArthur Model (R-M model) or the predator-prey model with Holling type II functional response, [6] studied stability of the R-M model by using graphical method, [7] studied global stability of the R-M model. In [8], conditions for an interior equilibrium are given, and the stability of this equilibrium is analyzed. Certain critical cases, some of which cannot occur in the usual model are also discussed. (If m = 1, q(y) = 0, the above system reduces to the “usual” system.) [9] investigated a predator-prey system for the global stability and existence of limit cycles; they proved that there exist at least two limit cycles by using qualitative analysis and the idea of Poincare-Bendixson theory.

In recent years, for a more complicated system with Holling type II functional response, Liu et al. [10] investigated a predator-prey model with Holling type II functional response incorporating a constant prey refuge. The authors studied the instability and global stability properties of the equilibrium points and the existence and uniqueness of limit cycle. Lv et al. [11] considered a model to describe the harvesting for the phytoplankton and zooplankton based on plausible toxic-phytoplankton-zooplankton systems. Shanbing Li et al. [12] studied a spatially heterogeneous predator-prey model where the interaction is governed by Holling II functional response. They showed that the degeneracy for the prey and predator has distinctly different effects on the coexistence states of the two species when intrinsic growth rate of the prey is above a certain critical value. P. D. N. Srinivasu et al. [13] considered a prey-predator model with Holling type II of predation and independent harvesting in either species. Their study showed that using the harvesting efforts as controls can break the cyclic behavior of the system and drive it to a required state. By introducing impulsive control strategy, Yongzhen Pei et al. [2] investigated the dynamics behaviors of one-prey multi-predator model with Holling type II functional response with the help of Floquet theorem and small amplitude perturbation method. They showed that multi-predator impulsive control strategy is more effective than the classical and single one.

Motivated by the above works, especially the references [8] and [9], we investigate a certain predator-prey system with density-dependent predator specific death rate and predator mutual interference incorporating a square term, which is described by the following nonlinear ordinary differential equations

$\stackrel{\dot{}}{x}=r_{1}x\left(1-\frac{x}{K_1}\right)-\frac{\alpha xy}{a+x}-m_{1}x\mathrm{,}$ (1a)

$\stackrel{\dot{}}{y}=\frac{\alpha exy}{a+x}-m_{2}y-d{y}^2$ (1b)

subject to initial conditions $x\left(0\right),y\left(0\right)>0$. Here, x and y are the prey and predator densities at time t, respectively. All the above positive constants have biological considerations. Parameter $r_1$ denotes the intrinsic growth rate of the prey; $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 species, respectively; e is

the biomass conversion; d is the intra-specific competition coefficient. The term $\frac{\alpha xy}{a+x}$ named Holling type II functional response describes the functional response of the predator. The specific growth term $r_{1}x\left(1-\frac{x}{K_1}\right)$ governs the increase

of the prey in the lack of predator. The square term $d{y}^2$ denotes intrinsic decrease of the predator. 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) [9].

The main purpose and conclusions of this paper are to analyze equilibrium stability and bifurcations of the above system by using qualitative analysis and bifurcation theories. This paper is organized as follows. Preliminaries, such as non-negativity, boundness, permanence and invariant sets are given in Section 2. In Section 3, we give existence and sufficient conditions for the stability analysis of equilibrium points. In critical cases of interior equilibrium, we describe sufficient criteria to illustrate focus-center problem and problem of higher-order equilibrium. Existence and non-existence conditions of limit cycle(s) are also presented in this section. In Section 4, analyses of transcritical and Hopf bifurcations are presented. Here we choose d as a Hopf bifurcation parameter to consider limit cycle(s) further and point out that the control parameter d how to affect the complexity of this system. In Section 5, we carry out numerical simulations and conclusions of our system.

2. Preliminaries

2.1. Invariant Sets

Denotes the first quadrant as $R^{+2}$, and its closure is $R_{+}^2=\overline{R^{+2}}$. For practical biological consideration, system (1) is defined on the domain $R^{+2}$ and all the solutions are positive. It is easy to prove lemma 1. Thus, any trajectory starting from $R_{+2}$ cannot cross the coordinate axes.

Lemma 1. All the solutions of system (1) are non-negative with the non-negative initial value $\left(x\left(t_0\right)\mathrm{,}y\left(t_0\right)\right)$, i.e. $R_{+}^2$ is an invariant set.

We will give some cases about invariant sets for system (1) in supplementary to lemma 1.

Remark. a) If $A\ge \max \left\{x_2,0\right\},x_2:=K_1\left(1-\frac{m_1}{r_1}\right)$ and $B\ge max\left\{\frac{\alpha e-m_2}{d}\mathrm{,0}\right\}$, domain $\left[\mathrm{0,}A\right]\times \left[\mathrm{0,}B\right]$ is an invariant set; b) Domain $\left[A_1\mathrm{,}{A}_2\right]\times \left[\mathrm{0,}+\infty \right]$ and $\left[\mathrm{0,}+\infty \right]\times \left[B_1\mathrm{,}{B}_2\right]$ is not an invariant set, if $A_1,B_1>0$ ; c) If $r_1\le m_1$ or $\alpha e\le m_2$, domain $\left[A_1\mathrm{,}{A}_2\right]\times \left[B_1\mathrm{,}{B}_2\right]$ ( $A_1,B_1>0$ ) is not an invariant set.

2.2. Boundness

Theorem 1 All the solutions of system (1) are bounded with non-negative initial conditions $\left(x\left(t_0\right)\mathrm{,}y\left(t_0\right)\right)$.

Proof. From the equation (1a) we have inequality $\stackrel{\dot{}}{x}\le x\left(r_1-r_{1}x/K_1\right)$. With the help of [14], this implies that $x\left(t\right)\le max\left\{x\left(t_0\right)\mathrm{,}{K}_1\right\}$. We denote the positive upper bound as $M_1$. A similar argument, from the Equation (1b) we have inequality

$\stackrel{\dot{}}{y}\le y\left(\frac{\alpha ex}{a+x}-dy\right)\le y\left(\frac{\alpha e{M}_1}{a+M_1}-dy\right)\mathrm{.}$ (2)

Similarly, this implies that $y\left(t\right)\le max\left\{y\left(t_0\right)\mathrm{,}\frac{\alpha e{M}_1}{a+M_1}\right\}\mathrm{<mo>\; :\; </mo>}=M_2$. This completes the proof and the system under consideration is dissipative.

Lemma 2. All the solutions of system (1) satisfy

$\underset{t\to \infty }{\lim \sup }x\left(t\right)\le \frac{r_1{K}_1}{4},\underset{t\to \infty }{\lim \sup }y\left(t\right)\le \frac{r_1{K}_1}{4e_1}.$ (3)

Proof. Taking an auxiliary function $z=ex+y$, it is obvious to see that

$\stackrel{\dot{}}{z}\le e{r}_{1}x\left(1-x/K_1\right)-mz\le e{r}_1{K}_1/4-mz,$ (4)

where $m=min\left\{m_1\mathrm{,}{m}_2\right\}$. By applying the theory of differential inequality [15], we obtain

$z\left(t\right)\le \frac{e{r}_1{K}_1}{4}\left(1-{\text{e}}^{-m\left(t-t_0\right)}\right)+z\left(t_0\right){\text{e}}^{-m\left(t-t_0\right)}\mathrm{,}$ (5)

which, upon letting $t\to \infty$, yields $\underset{t\to \infty }{\lim \sup }z\left(t\right)\le e{r}_1{K}_1/4$. This completes the proof.

Remark. It is clear to see that all the solutions of system (0) are uniformly bounded.

Theorem 2. If $r_1\le m_1$, then the system (1) is Lagrange stable.

Proof. Taking the positive definite Lyapunov function $V\left(x\mathrm{,}y\right)=ex+y$ with infinitesimally small upper and large lower bounds, computing its Dini derivative along the trajectories of Equations (1), we have

$D^{+}V\le -\frac{e{r}_1{x}^2}{K_1}-m_{2}y-d{y}^2.$ (6)

Thus $D^{+}V$ is negative definite. This completes the proof.

2.3. Permanence

Theorem 3. If parameters of system (1) satisfy

${\omega }_1:=\frac{K_1}{r_1}\left(r_1-m_1-\frac{\alpha M_2}{a}\right)>0,$ (7a)

${\omega }_2:=\left[\frac{\alpha e\left(1-\lambda \right){\omega }_1}{a+e\left(1-\lambda \right){\omega }_1}-m_2\right]/d>0,$ (7b)

where $\lambda \in \left(\mathrm{0,1}\right)$, then system (1) is permanent.

Proof. By theorem (2) we have a positive upper bound ${\xi }_1$ such that

$max\left\{\underset{t\to \infty }{\lim \sup }x\left(t\right)\mathrm{,}\underset{t\to \infty }{\lim \sup }y\left(t\right)\right\}\le {\xi }_1\mathrm{.}$ (8)

From Equation (1a) we have

$\stackrel{\dot{}}{x}\ge x\left(r_1-m_1-\alpha M_2/a-r_{1}x/K_1\right)\mathrm{,}$

here $M_2$ is defined in theorem 2. With the help of [14], this implies that $\underset{t\to \infty }{\lim \inf }x\left(t\right)\ge {\omega }_1$. For ${\epsilon }_0=\lambda {\omega }_1$, there must exit T, such that for any $t\ge T$ we have $x\left(t\right)\ge \left(1-\lambda \right){\omega }_1$, thus

$\stackrel{\dot{}}{y}\ge \left[\frac{\alpha e\left(1-\lambda \right){\omega }_1}{a+e\left(1-\lambda \right){\omega }_1}-m_2\right]y-d{y}^2\mathrm{.}$

Similarly, this implies that $\underset{t\to \infty }{\lim \inf }y\left(t\right)\ge {\omega }_2$. This completes the proof and we get the permanence of the system.

3. Equilibria

In this section we will discuss equilibria of the system (1) with their sufficient existence conditions and stability analysis.

It is obvious to see that our system has following trivial equilibria: $E_0=\left(\mathrm{0,0}\right)$,$E_1=\left(\mathrm{0,}{y}_1\right)$,$E_2=\left(x_2\mathrm{,0}\right)$, where $y_1=-m_2/d$. For practical or biological considerations, we omit the singular point $E_1$. The point $E_2$ is a desired equilibrium only if $r_1>m_1$. When $r_1=m_1$, the critical point $E_2$ becomes $E_0$.

Then we make a special effort to derive the existence conditions of an interior equilibrium denoted by $E_{\ast }=\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$, where it can be given from the following algebraic equations

$f\left(x\mathrm{,}y\right)\mathrm{<mo>\; :\; </mo>}=r_1\left(1-\frac{x}{K_1}\right)-\frac{\alpha y}{a+x}-m_1=\mathrm{0,}$ (9a)

$g\left(x,y\right):=\frac{\alpha ex}{a+x}-m_2-dy=0.$ (9b)

From Equation (9a) we know $y\to a\left(r_1-m_1\right)/\alpha$ as $x\to 0$ ; $x\to x_2$ as $y\to 0$ ; from Equation (9b) we know $y\to -m_2/d$ as $x\to 0$ ; $x\to x_0$ as $y\to 0$,$x_0=a{m}_2/\left(\alpha e-m_2\right)$. Thus, if the following conditions $m_1<r_1$,$0<x_0<x_2$ satisfy, an interior equilibrium $E_{\mathrm{<mo>\; *\; </mo>}}$ exists and meanwhile we have

$x_0<x_{\ast }<x_2$,$0<y_{\ast }<\bar{y}$, where $\bar{y}=\frac{r_1}{4a{K}_1}{\left(a+x_2\right)}^2$.

3.1. Equilibria Type

Here we will use the Routh-Hurwitz criteria and the Perron’s theorems to analyze local stability and type of these equilibria by the nature of eigenvalues of Jacobian matrices around them, respectively. We observed that the Jacobian matrix of the system (1) is

$J=\left[\begin{array}{cc}{r}_1-\frac{2r_{1}x}{K_1}-\frac{\alpha ay}{{\left(a+x\right)}^2}-m_1& -\frac{\alpha x}{a+x}\\ \frac{\alpha eay}{{\left(a+x\right)}^2}& \frac{\alpha ex}{a+x}-m_2-2dy\end{array}\right]\mathrm{.}$ (10)

For zero point $E_0$, from its Jacobian matrix $J\left(E_0\right)$ we have: a) If $r_1<m_1$,$E_0$ is an asymptotically stable node; b) If $r_1>m_1$,$E_0$ is a saddle.

For the critical value $r_1=m_1$,$E_0$ is a higher order singular point. According to Frommer’s method, by using its Jacobian matrix and polar-coordinate-transformation: $x=\rho \cos \theta ,y=\rho \sin \theta$, we derive $R\left(\theta \right)=-m_{2}si{n}^2\theta$ and characteristic equation $\Theta \left(\theta \right)=-m_2\sin \theta \cos \theta$. For the first simple real root ${\theta }_1\mathrm{<mo>\; =\; </mo>0}$ of the characteristic equation, from the calculation $R\left({\theta }_1\right)=0$,$R^{\prime }\left({\theta }_1\right)=0$,$\Theta \left({\theta }_1\right)\ne 0$ and $k=1$, we know that $\theta ={\theta }_1$ is a characteristic direction. While for the second simple root ${\theta }_2=\pi /2$, it is clear to see an expansion $\Theta \left(\theta \right)=m_2\left(\theta -{\theta }_2\right)+\cdots$, i.e. $p=1$ (odd number) and $C_2^{\left(1\right)}\left(=m_2\right)$ has the opposite sign as $R\left({\theta }_2\right)$. Thus $\theta ={\theta }_2$ is also a characteristic direction. By the sign of $R\left({\theta }_2\right)\Theta \left(\theta \right)=m_2^{2}sin\theta cos\theta$, we know that there is only one orbit tending to the critical point $E_0$ along the direction ${\theta }_2$ in the normal regions of second type.

Taking the transformation $t=\frac{-\tau }{m_2}$, we have following normal form (retain the symbol “ $\cdot$ ” for simplicity)

$\stackrel{\dot{}}{x}=x+\Phi \left(x,y\right),\stackrel{\dot{}}{y}=\Psi \left(x,y\right),$ (11)

where $\Phi \left(x\mathrm{,}y\right)$,$\Psi \left(x\mathrm{,}y\right)=o\left(\left|x\mathrm{,}y\right|\right)$ in the right hand side are analytical functions, with terms starting from second order. From center manifold $x=h\left(y\right)$,

the series of function $\psi \left(y\right)=\Psi \left(h\left(y\right),y\right)=\frac{r_1}{m_2{K}_1}{y}^2+\cdots$ with respect to y, and

$m=2$ (even number), the critical point $E_0$ is a saddle node and the parabolic sector is on the right half (x,y)-plane.

For point $E_2$, from its Jacobian matrix $J\left(E_2\right)$ and existence condition we have: a) If $\frac{\alpha e{x}_2}{a+x_2}<m_2$,$E_2$ is an asymptotically stable node; b) If $\frac{\alpha e{x}_2}{a+x_2}>m_2$,$E_2$ is a saddle.

For the critical value $\frac{\alpha e{x}_2}{a+x_2}=m_2$,$E_2$ is a higher order singular point. By

using its Jacobian matrix and polar-coordinate-transformation, we derive $R\left(\theta \right)=cos\theta \chi \left(\theta \right)$ and characteristic equation $\Theta \left(\theta \right)=-sin\theta \chi \left(\theta \right)$, where

$\chi \left(\theta \right)=\left(m_1-r_1\right)cos\theta -\frac{m_2}{e}sin\theta \mathrm{.}$

The real simple root ${\theta }_1=0$ is in $\left[\mathrm{0,}\pi /2\right]$, it is clear to see an expansion $\Theta \left(\theta \right)=-\left(m_1-r_1\right)\theta +\cdots$ and $R\left({\theta }_1\right)\ne 0$, i.e. $p=1$ (odd number) and $C_1^{\left(1\right)}\left(=r_1-m_1\right)$ has contrary sign with $R\left({\theta }_1\right)$. Thus $\theta ={\theta }_1$ is indeed a characteristic direction. While $R\left({\theta }_1\right)\Theta \left(\theta \right)=\left(r_1-m_1\right)\chi \left(\theta \right)\sin \theta$ changes its sign in the small neighbourhood of ${\theta }_1$, thus we know that there is only one orbit tending to the critical point $E_2$ along the direction ${\theta }_1$ in the normal regions of second type. Another real simple root ${\theta }_2$ is in $\left(\pi /2\mathrm{,}\pi \right)$, it is clear to see that $R\left({\theta }_2\right)=0$ and

${\Theta }^{\prime }\left({\theta }_2\right)=\left(m_1-r_1\right)\left[1+\frac{{\left(m_1-r_1\right)}^2{e}^2}{m_2^2}\right]co{s}^2{\theta }_2<\mathrm{0,}$

i.e. $k=1$, thus $\theta ={\theta }_1$ is also a characteristic direction.

Taking transformations $x\to x+x_2$ ; $x=u+\frac{\alpha x_2}{a+x_2}v$,$y=\left(m_1-r_1\right)v$ and $t=\frac{\tau }{m_1-r_1}$, we have likewise real Jordan normal form

$\stackrel{\dot{}}{u}=u+\Phi \left(u,v\right),\stackrel{\dot{}}{v}=\Psi \left(u,v\right),$ (12)

in here $\Phi \left(u,v\right)$,$\Psi \left(u\mathrm{,}v\right)=o\left(\left|u\mathrm{,}v\right|\right)$ in the right hand side are analytical functions, with terms starting from second order. From the center manifold $u=h\left(v\right)$, the series of function $\psi \left(v\right)=\Psi \left(h\left(v\right),v\right)=a_2{v}^2+\cdots$ with respect to v is derived, where

$a_2=-{\left(m_1-r_1\right)}^{2}d+\frac{\left(m_1-r_1\right){\alpha }^{2}ae{x}_2}{{\left(a+x_2\right)}^3}<0,$ (13)

Combining with $m=2$ (even number), the critical point $E_{\mathrm{<mo>\; *\; </mo>}}$ is a saddle node and the parabolic sector is on the left half (u, v)-plane.

For the interior equilibrium $E_{\mathrm{<mo>\; *\; </mo>}}$, from its Jacobian matrix $J\left(E_{\mathrm{<mo>\; *\; </mo>}}\right)$ and existence conditions, denoting discriminant $\Delta =A_1^2-4A_2$ with notations

$A_1:=trJ\left(E_{*}\right)=\frac{\alpha x_{*}{y}_{*}}{{\left(a+x_{*}\right)}^2}-\frac{r_1{x}_{*}}{K_1}-d{y}_{*},$ (14a)

$A_2:=\det J\left(E_{*}\right)=\frac{{\alpha }^{2}ea{x}_{*}{y}_{*}}{{\left(a+x_{*}\right)}^3}-d{y}_{*}\left(A_1+d{y}_{*}\right),$ (14b)

we have:

(a) If $A_1<0$ and

(a1) $A_2>0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is an asymptotically stable node; (a2) $A_2<0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is a saddle;

(b) If $A_1=0$ and

(b1) $A_2>0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is a center or a focus; (b2) $A_2<0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is a saddle;

(c) If $A_1>0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is unstable and

(c1) $\Delta =0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is a node; (c2) $\Delta <0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is a focus; (c3) $\Delta >0$ and $A_2>0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is a node; (c4) $\Delta >0$ and $A_2<0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is a saddle.

Remark. If $A_1+d{y}_{*}<0$,$E_{\mathrm{<mo>\; *\; </mo>}}$ is asymptotically stable.

3.2. The Special Case: $A_1=0,A_2>0$

If the interior equilibrium $E_{\mathrm{<mo>\; *\; </mo>}}$ exists and $A_1\mathrm{<mo>\; =\; </mo>0,}{A}_2\mathrm{<mo>\; >\; </mo>0}$, i.e. the Jacobian matrix $J\left(E_{\mathrm{<mo>\; *\; </mo>}}\right)$ has a pair of purely imaginary eigenvalues ${\lambda }_{\mathrm{1,2}}=\pm i\beta$, in here $\beta =\sqrt{A_2}$, or $E_{\mathrm{<mo>\; *\; </mo>}}$ is a center of local linear approximation of system (1). We consider this critical case or the center and focus problem by the H. Poincare’s formal series method. With the application of linear transforms

$x\to x+x_{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}y\to y+y_{\mathrm{<mo>\; *\; </mo>}}$ and $x=d{y}_{\mathrm{<mo>\; *\; </mo>}}u+\beta v\mathrm{,}y=\frac{\alpha ea{y}_{\mathrm{<mo>\; *\; </mo>}}}{{\left(a+x_{\mathrm{<mo>\; *\; </mo>}}\right)}^2}u$ (or $x=-\beta u+d{y}_{*}v,y=\frac{\alpha ea{y}_{*}}{{\left(a+x_{*}\right)}^2}v$ ), the above system (1) transforms into following normal form

$\stackrel{\dot{}}{u}=\beta v+P\left(u,v\right),\stackrel{\dot{}}{v}=-\beta u+Q\left(u,v\right),$ (15)

where higher order terms $P\left(u,v\right)$ and $Q\left(u,v\right)$ in the right hand side of the above system are all real analytical functions of u and v, with terms starting from second order.

Suppose that the above system has a first integral in the form of $F\left(u\mathrm{,}v\right)=u^2+v^2+F_3\left(u\mathrm{,}v\right)+F_4\left(u\mathrm{,}v\right)+\cdots$, where $F_k\left(u\mathrm{,}v\right)$ is a homogeneous polynomial with order k ( $k=\mathrm{3,4,}\cdots$ ). Computing its derivative along the above system, we have its total derivative

$\frac{\text{d}F}{\text{d}t}=\left(2u+F_{3u}+\cdots \right)\left(\beta v+P_2+\cdots \right)+\left(2v+F_{3v}+\cdots \right)\left(-\beta u+Q_2+\cdots \right)\mathrm{,}$ (16)

where $P_k\left(u\mathrm{,}v\right)$ and $Q_k\left(u\mathrm{,}v\right)$ are homogeneous polynomials in u, v of degree k ( $k=\mathrm{2,3,}\cdots$ ). In order to get Fourier series for further analysis, we firstly take the polar-coordinate-transformation $u=r\cos \theta ,v=r\sin \theta$ and assume $F_k=r^k{\Phi }_k\left(\theta \right)$ ( $k=\mathrm{3,4,}\cdots$ ), then following recurrence relations are derived by comparing coefficients:

$\begin{array}{l}\frac{\text{d}{\Phi }_3}{\text{d}\theta }={\varphi }_3\left(\theta \right),{\varphi }_3\left(\theta \right):=\frac{2}{\beta r^3}\left(u{P}_2+v{Q}_2\right),\\ \frac{\text{d}{\Phi }_4}{\text{d}\theta }={\varphi }_4\left(\theta \right),{\varphi }_4\left(\theta \right):=\frac{1}{\beta r^4}\left(2u{P}_3+P_2{F}_{3u}+2v{Q}_3+Q_2{F}_{3v}\right),\\ ⋮\end{array}$

From a complicated expression of function ${\varphi }_3\left(\theta \right)$ or the integral ${\displaystyle {\int }_0^{2\pi }}{\varphi }_3\left(\theta \right)\text{d}\theta =0$, we know that ${\varphi }_3\left(\theta \right)$ is precisely a periodic function of period $2\pi$. Thus the focus value $c_0^{\left(1\right)}$ of $E_{\mathrm{<mo>\; *\; </mo>}}$ reads

$\begin{array}{c}{c}_0^{\left(4\right)}=-\frac{e{\alpha }^2{s}_1^{2}a{s}_2}{2K_1^3{\bar{a}}^{11}{\beta }^3}[r_1^3{s}_1^8+6a{r}_1^3{s}_1^7+\left(a\alpha e{K}_1{r}_1^2+15a^2{r}_1^3-\alpha K_1{r}_1^2{s}_2\right)s_1^6\underset{}{\overset{}{}}\\ +20a\left(\frac{1}{4}ae\alpha K_1+a^2{r}_1-\frac{9}{40}{K}_1{s}_2\alpha \right)r_1^2{s}_1^5\\ +\left(-\frac{1}{2}ae{\alpha }^2{K}_1^2{r}_1{s}_2+10a^{3}e\alpha K_1{r}_1^2-8a^2\alpha K_1{r}_1^2{s}_2+15a^4{r}_1^3\right)s_1^4\\ +6a{r}_1\left(\frac{5}{3}{a}^{3}e\alpha K_1{r}_1-\frac{1}{12}{\alpha }^{2}ea{s}_2{K}_1^2+a^4{r}_1^2-\frac{7}{6}{a}^2\alpha K_1{r}_1{s}_2+\frac{1}{6}{\alpha }^2{K}_1^2{s}_2^2\right)s_1^3\end{array}$

$\begin{array}{c}+a(\frac{3}{2}{a}^{2}e{\alpha }^2{K}_1^2{r}_1{s}_2+5a^{4}e\alpha K_1{r}_1^2+\frac{1}{2}{\alpha }^{3}e{K}_1^3{s}_2^2-3a^3\alpha K_1{r}_1^2{s}_2\\ \underset{}{\overset{}{}}+2a{\alpha }^2{K}_1^2{r}_1{s}_2^2+a^5{r}_1^3)s_1^2+a\alpha K_1(\frac{5}{2}{a}^{3}e\alpha K_1{r}_1{s}_2+a^{5}e{r}_1^2+\frac{1}{2}ae{\alpha }^2{K}_1^2{s}_2^2\\ -\frac{1}{2}{a}^4{r}_1^2{s}_2+a^2\alpha K_1{r}_1{s}_2^2-\frac{1}{2}{\alpha }^2{K}_1^2{s}_2^3)s_1+a^5{\alpha }^{2}e{K}_1^2{r}_1{s}_2],\end{array}$ (17)

where $\bar{a}\mathrm{<mo>\; :\; </mo>}=a+x_{\mathrm{<mo>\; *\; </mo>}}$,$s_1\mathrm{<mo>\; :\; </mo>}=x_{\mathrm{<mo>\; *\; </mo>}}$ and $s_2\mathrm{<mo>\; :\; </mo>}=y_{\mathrm{<mo>\; *\; </mo>}}$ for later convenience. When $c_0^{\left(4\right)}\ne 0$, the point $E_{\mathrm{<mo>\; *\; </mo>}}$ is a fine focus of order 1. If $c_0^{\left(4\right)}>0$,$E_{\mathrm{<mo>\; *\; </mo>}}$ is an unstable fine focus; If $c_0^{\left(4\right)}<0$,$E_{\mathrm{<mo>\; *\; </mo>}}$ is a stable fine focus. When $c_0^{\left(4\right)}=0$, continuing this procedure to obtain $c_0^{\mathrm{<mo\; stretchy="false">(\; </mo>6\; <mo\; stretchy="false">)\; </mo>}}$.

3.3. The Special Case: $A_2=0,A_1\ne 0$

In this section, we consider a critical case: $A_2=\mathrm{0,}{A}_1\ne 0$, i.e. $E_{\mathrm{<mo>\; *\; </mo>}}$ is a critical point of higher order and its Jacobian matrix with zero eigenvalue can be rewritten as

$J\left(E_{\mathrm{<mo>\; *\; </mo>}}\right)=\left[\begin{array}{cc}\frac{\alpha \lambda x_{\mathrm{<mo>\; *\; </mo>}}}{a+x_{\mathrm{<mo>\; *\; </mo>}}}& -\frac{\alpha x_{\mathrm{<mo>\; *\; </mo>}}}{a+x_{\mathrm{<mo>\; *\; </mo>}}}\\ \lambda d{y}_{\mathrm{<mo>\; *\; </mo>}}& -d{y}_{\mathrm{<mo>\; *\; </mo>}}\end{array}\right]\mathrm{,}$ (18)

where positive parameter $\lambda \mathrm{<mo>\; :\; </mo>}=\frac{\alpha ea}{d\left(a+x_{\mathrm{<mo>\; *\; </mo>}}\right)}$. By using its Jacobian matrix and the polar-coordinate-transformation, we derive function

$R\left(\theta \right)=\left(\frac{\alpha x_{\mathrm{<mo>\; *\; </mo>}}}{a+x_{\mathrm{<mo>\; *\; </mo>}}}cos\theta +d{y}_{\mathrm{<mo>\; *\; </mo>}}sin\theta \right)\left(\lambda cos\theta -sin\theta \right)$ (19)

and its characteristic equation

$\Theta \left(\theta \right)=\left(d{y}_{*}\cos \theta -\frac{\alpha x_{*}}{a+x_{*}}\sin \theta \right)\left(\lambda \cos \theta -\sin \theta \right)=0.$ (20)

It is clear to see that the real roots ${\theta }_{\mathrm{1,2}}$ in $\left(\mathrm{0,}\pi /2\right)\cup \left(\pi \mathrm{,}3\pi /2\right)$ respectively satisfy

$\tan {\theta }_1=\frac{d{y}_{*}\left(a+x_{*}\right)}{\alpha x_{*}},\tan {\theta }_2=\lambda .$

Case I: $\theta ={\theta }_1$

In this case, for $\theta ={\theta }_1$,$R\left({\theta }_1\right)\ne 0$ if and only if $tan{\theta }_1\ne \lambda$. With some calculation, we derive $R\left({\theta }_1\right)=-c_1^{\left(1\right)}$ and series of function $\Theta \left(\theta \right)$ at ${\theta }_1$ : $\Theta \left(\theta \right)=c_1^{\left(1\right)}\left(\theta -{\theta }_1\right)+\cdots$, where

$c_1^{\left(1\right)}=\frac{\left[d^2{y}_{\mathrm{<mo>\; *\; </mo>}}^2{\left(a+x_{\mathrm{<mo>\; *\; </mo>}}\right)}^2+{\alpha }^2{x}_{\mathrm{<mo>\; *\; </mo>}}^2\right]\left(ad{y}_{\mathrm{<mo>\; *\; </mo>}}-\alpha \lambda x_{\mathrm{<mo>\; *\; </mo>}}+d{x}_{\mathrm{<mo>\; *\; </mo>}}{y}_{\mathrm{<mo>\; *\; </mo>}}\right)sin{\left({\theta }_1\right)}^2}{{\left(a+x_{\mathrm{<mo>\; *\; </mo>}}\right)}^3{d}^2{y}_{\mathrm{<mo>\; *\; </mo>}}^2}\mathrm{,}$ (21)

i.e. $p=1$ (odd number), $R\left({\theta }_1\right)$ and $c_1^{\left(1\right)}$ have opposite sign. Thus $\theta ={\theta }_1$ is a characteristic direction. While $R\left({\theta }_1\right)\Theta \left(\theta \right)$ changes its sign in the small neighbourhood of ${\theta }_1$, thus we know that there is only one orbit tending to the critical point $E_2$ along the direction ${\theta }_1$ in the normal regions of second type.

If $tan{\theta }_1=\lambda$, then $R\left({\theta }_1\right)=c_1^{\left(1\right)}=0$, while

$R^{\prime }({\theta }_1)=-\frac{\left[\left(d^2{y}_{*}^2+{\alpha }^2\right)x_{*}^2+2a{d}^2{x}_{*}{y}_{*}^2+a^2{d}^2{y}_{*}^2\right]\left[\left(d\lambda y_{*}+\alpha \right)x_{*}+\lambda d{y}_{*}a\right]\sin {\left({\theta }_1\right)}^2}{{\left(a+x_{*}\right)}^3{d}^2{y}_{*}^2}<0.$

Thus $\theta ={\theta }_1$ is also a characteristic direction.

Case II: $\theta ={\theta }_2$

In this case, for $\theta ={\theta }_2$, we have $R\left({\theta }_2\right)=0$,

$R^{\prime }\left({\theta }_1\right)=-\frac{\left({\lambda }^2+1\right)\left[a{x}_{*}+\lambda d{y}_{*}\left(a+x_{*}\right)\right]\sin {\left({\theta }_2\right)}^2}{{\lambda }^2\left(a+x_{*}\right)}<0$

and expression of function $\Theta \left(\theta \right)$ at ${\theta }_2$ : $\Theta \left(\theta \right)=c_2^{\left(1\right)}\left(\theta -{\theta }_2\right)+\cdots$, where

$c_2^{\left(1\right)}=-\frac{\left({\lambda }^2+1\right)\left(ad{y}_{\mathrm{<mo>\; *\; </mo>}}-\alpha \lambda x_{\mathrm{<mo>\; *\; </mo>}}+d{x}_{\mathrm{<mo>\; *\; </mo>}}{y}_{\mathrm{<mo>\; *\; </mo>}}\right)sin{\left({\theta }_2\right)}^2}{{\lambda }^2\left(a+x_{\mathrm{<mo>\; *\; </mo>}}\right)}\mathrm{.}$ (22)

If $tan{\theta }_1\ne \lambda$, then $c_2^{\left(1\right)}\ne 0$ and $\theta ={\theta }_1$ is not a special direction; If $tan{\theta }_1=\lambda$, then $\theta ={\theta }_1$ is a characteristic direction.

From now on, making a change of variables $x\to x+x_{\mathrm{<mo>\; *\; </mo>}}$,$y\to y+y_{\mathrm{<mo>\; *\; </mo>}}$, and taking transformations $x=\frac{\alpha x_{\mathrm{<mo>\; *\; </mo>}}u}{a+x_{\mathrm{<mo>\; *\; </mo>}}}+v$,$y=d{y}_{\mathrm{<mo>\; *\; </mo>}}u+\lambda v$ and $t=\frac{\tau }{A_1}$. Substituting them into the system (1), we have following normal form

$\stackrel{\dot{}}{u}=u+\Phi \left(u,v\right),\stackrel{\dot{}}{v}=\Psi \left(u,v\right),$ (23)

where $\Phi \left(u,v\right)$,$\Psi \left(u\mathrm{,}v\right)=o\left(\left|u\mathrm{,}v\right|\right)$ in the right hand side are analytical functions, with terms starting from second order. From center manifold $u=h\left(v\right)$ and series of function $\psi \left(v\right)=\Psi \left(h\left(v\right),v\right)=a_2{v}^2+a_3{v}^3+\cdots$ with respect to v, where

$a_2=-\frac{d{y}_{\mathrm{<mo>\; *\; </mo>}}\left(r_{1}d{\bar{a}}^4-a\alpha K_1\left(\alpha e+d{y}_{\mathrm{<mo>\; *\; </mo>}}\right)\bar{a}+2{\alpha }^{2}e{a}^2{K}_1\right)}{\bar{a}{K}_1\left(d^2{\bar{a}}^3{y}_{\mathrm{<mo>\; *\; </mo>}}-a{\alpha }^{2}e{x}_{\mathrm{<mo>\; *\; </mo>}}\right)A_1}\mathrm{,}$ (24)

one can judge the type of critical point $E_{\mathrm{<mo>\; *\; </mo>}}$ in this special case. When $a_2\ne 0$, i.e. $r_{1}d{\bar{a}}^4+2{\alpha }^{2}e{a}^2{K}_1\ne a\alpha K_1\left(\alpha e+d{y}_{\mathrm{<mo>\; *\; </mo>}}\right)\bar{a}$, combining with $m=2$ (even number), the critical point $E_{\mathrm{<mo>\; *\; </mo>}}$ is a saddle node. Furthermore, if the coefficient $a_2>0$ (or < 0), then the parabolic sector is on the right (or left) half (u,v)-plane. If $a_2=0$, continuing the above series to obtain

$\begin{array}{c}{a}_3=\frac{a\alpha }{d{K}_1^2{A}_1^2{\bar{a}}^6{\left(-s_{1}ae{\alpha }^2+d^2{s}_2{\bar{a}}^3\right)}^2}[a^4{e}^4{K}_1^2{s}_1^2{\alpha }^6+d{K}_1^2\bar{a}{s}_1{s}_2{a}^3\left(4a-3\bar{a}\right)e^3{\alpha }^5\\ -d{K}_1{\bar{a}}^2{s}_1{a}^2(-r_{1}e{\bar{a}}^3+\left(\left(d{s}_2+A_1\right)K_1+r_{1}a\right)e{\bar{a}}^2-K_1\left(ae{A}_1+2d{s}_2^2\right)\bar{a}\\ +5a{s}_2^{2}d{K}_1)e^2{\alpha }^4+d^2{s}_2{K}_1{\bar{a}}^3(-2r_{1}e{\bar{a}}^4+\left(\left(2d{s}_2+A_1\right)K_1+7r_{1}a\right)e{\bar{a}}^3\\ -\left(\left(4d{s}_2+3A_1\right)K_1+5r_{1}a\right)ae{\bar{a}}^2+\left(a\left(3d{s}_2+2A_1\right)e+2d{s}_2^2\right)K_{1}a\bar{a}\end{array}$

$\begin{array}{c}-2a^2{s}_2^{2}d{K}_1)ae{\alpha }^3-d^3{s}_2{K}_1{\bar{a}}^{5}a(r_{1}e{\bar{a}}^3+\left(\left(\left(-d{s}_2-A_1\right)K_1-r_{1}a\right)e+4r_1{s}_2\right){\bar{a}}^2\\ +\left(a{K}_1\left(d{s}_2+A_1\right)e+\left(-2d{s}_2^2-A_1{s}_2\right)K_1-4r_1{s}_{2}a\right)\bar{a}+K_1{s}_{2}a\left(3d{s}_2+A_1\right))e{\alpha }^2\\ -d^3{s}_2{\bar{a}}^7(-2r_1^2{\bar{a}}^3+\left(2d{K}_1{r}_1{s}_2+2a{r}_1^2\right){\bar{a}}^2-d{K}_1{s}_2\left(3r_{1}a-A_1{K}_1\right)\bar{a}\\ +a{s}_{2}d{K}_1^2\left(d{s}_2-2A_1\right))e\alpha -s_2^3{d}^5{\bar{a}}^8{A}_1{K}_1^2].\end{array}$ (25)

If $a_3>0$,$E_{\mathrm{<mo>\; *\; </mo>}}$ is an unstable node; If $a_3<0$,$E_{\mathrm{<mo>\; *\; </mo>}}$ is a saddle. When $a_3=0$, continuing the above procedure to obtain $a_4$ and using the same criteria of $a_2$.

3.4. The Special Case: $A_{\text{1}}=A_{\text{2}}=0$

In this section, we consider a critical case: $A_1=A_2=0$, i.e. $E_{\mathrm{<mo>\; *\; </mo>}}$ is a critical point of higher order and its Jacobian matrix with two zero eigenvalues can be rewritten as

$J\left(E_{\mathrm{<mo>\; *\; </mo>}}\right)=d{y}_{\mathrm{<mo>\; *\; </mo>}}\left[\begin{array}{cc}1& -\frac{1}{\lambda }\\ \lambda & -1\end{array}\right]\mathrm{.}$ (26)

From the above section we know that its characteristic equation is $\Theta \left(\theta \right)=\lambda \cos \theta -\sin \theta =0$, and function

$R\left(\theta \right)=d{y}_{*}\left(\frac{1}{\lambda }\cos \theta +\sin \theta \right)\left(\lambda \cos \theta -\sin \theta \right).$

Similarly, $\theta ={\theta }_1$ is a special direction, where $tan{\theta }_1=\lambda$.

Making a change of variables and taking transformations $x=u+v,y=\lambda u$ ; $t=\frac{\tau }{d{y}_{\mathrm{<mo>\; *\; </mo>}}}$, we have following normal form

$\stackrel{\dot{}}{u}=v+\Phi \left(u,v\right),\stackrel{\dot{}}{v}=\Psi \left(u,v\right).$ (27)

From center manifold $v=h\left(u\right)$, we have series of function $\psi \left(v\right)=\Psi \left(u\mathrm{,}h\left(u\right)\right)$ : $\psi \left(v\right)=a_2{u}^2+a_3{u}^3+\cdots$ and series of following function:

${\left[{\Phi }_u\left(u,v\right)+{\Psi }_v\left(u,v\right)\right]|}_{v=h\left(u\right)}=b_{1}u+b_2{u}^2+\cdots .$

where coefficients

$\begin{array}{l}{a}_2=\frac{2a\alpha \left(e+\lambda \right)\left[\left(a+s_1\right)\lambda -\frac{3}{2}{s}_2\right]}{{\left(a+s_1\right)}^4\lambda d{s}_2}\mathrm{,}\\ b_1=-\frac{\alpha a\left(-e\bar{a}-\lambda \bar{a}+2s_2\right)K_1-2{\bar{a}}^3{r}_1}{K_1{\bar{a}}^{3}d{s}_2}\mathrm{,}\\ ⋮\end{array}$ (28)

If $a_2\mathrm{,}{b}_1\ne 0$, i.e. $k=2$ (even number), $n=m=1$, or $a_2\ne 0,b_1=0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is a degenerate singular point. If $a_2=0,b_1\ne 0$, but

$a_3=\frac{-\left(e+\lambda \right)\alpha a{K}_1-6{\bar{a}}^2{r}_1}{2e\alpha a{\bar{a}}^2{K}_1}<\mathrm{0,}$ (29)

i.e. $k=3$ (odd number), $n=m=1$, and combining a discriminant $\mu =b_1^2+8a_3\ge 0$, where

$\begin{array}{c}\mu =\frac{1}{4{\bar{a}}^6{K}_1^2{d}^2{\lambda }^2}[(25\left({\text{e}}^2+\left(22/25\right)e\lambda +\left(1/25\right){\lambda }^2\right)a^2{\alpha }^2-68{\bar{a}}^2\lambda d\left(e+\left(11/17\right)\lambda \right)a\alpha \\ +36d^2{\bar{a}}^4{\lambda }^2)K_1^2+72{\bar{a}}^2\left(-\left(11/6\left(e+\left(1/11\right)\lambda \right)\right)a\alpha +{\bar{a}}^2\lambda d\right)r_1{K}_1+36{\bar{a}}^4{r}_1^2]\end{array}$ (30)

then $E_{\mathrm{<mo>\; *\; </mo>}}$ is a complicated singular point and its neighbourhood $U\left(E_{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\delta \right)$ consists of one hyperbolic sector and one elliptic sector, if $\mu <0$ and $b_1\ne 0$, or just $b_1=0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is a center or a focus.

3.5. Global Stability

In this section, we will give following theorems to illustrate the global stability by constructing proper Lyapunov functions.

Theorem 4. If $r_1\le m_1$, then the zero point $E_0$ is globally asymptotical stable.

Proof. Taking an unbounded positive definite Lyapunov function in lemma 2, we know that $D^{+}V$ is negative definite, similarly. This completes the proof.

Remark. Taking an unbounded positive definite Lyapunov function $V\left(x\mathrm{,}y\right)=\frac{1}{2}\left(u{x}^2+y^2\right)$, where positive constant $u>\frac{{\left(m_2-\alpha e\right)}^2}{4a\alpha d}$, and computing its Dini derivative along orbits of Equations (1), we have

$D^{+}V=-\frac{K_{1}g\left(x\mathrm{,}y\right)+u{r}_1{x}^3\left(a+x\right)+u{K}_1\left(m_1-r_1\right)x^3}{K_1\left(a+x\right)}\mathrm{,}$ (31)

where auxiliary function

$g\left(x\mathrm{,}y\right)=\left[ua\left(m_1-r_1\right)x^2+dx{y}^3+a{m}_2{y}^2\right]+\left[u\alpha x^2+\left(m_2-\alpha e\right)xy+ad{y}^2\right]y\mathrm{.}$

On account of the discriminant $\Delta ={\left(m_2-\alpha e\right)}^2-4u\alpha ad<0$ of a quadratic function, it is obvious that $g\left(x\mathrm{,}y\right)\ge 0$ or $D^{+}V$ is negative definite. This also completes the proof.

Theorem 5. Assume that the equilibrium $E_2$ exists and $\frac{\alpha e{x}_2}{a}\le m_2$, then $E_2$ is globally asymptotical stable.

Proof. Taking an unbounded positive definite Lyapunov function

$V\left(x\mathrm{,}y\right)=x-x_2-x_{2}ln\left(x/x_2\right)+y/e\mathrm{,}$ (32)

and computing its Dini derivative along orbits of Equations (1), we have

$D^{+}V\le -\frac{r_1}{K_1}{\left(x-x_2\right)}^2+\frac{y}{e}\left(\frac{\alpha e{x}_2}{a}-m_2\right)-\frac{d}{e}{y}^2$ (33)

This completes the proof.

Theorem 6. Assume that the equilibrium $E_2$ exists and $\alpha e\le m_2,x_2>2a$, then $E_2$ is globally asymptotical stable.

Proof. Taking an unbounded positive definite Lyapunov function $V\left(x\mathrm{,}y\right)=\frac{1}{2}\left[u{\left(x-x_2\right)}^2+y^2\right]$ with positive constant

$u\le \frac{4a{m}_2{r}_1\left(x_2-2a\right)}{{\alpha }^2{x}_2{K}_1}\mathrm{,}$ (34)

and computing its Dini derivative along orbits of Equations (1), we have $D^{+}V=\frac{g\left(x\mathrm{,}y\right)}{K_1\left(a+x\right)}$,$g\left(x\mathrm{,}y\right)=g_1\left(x\right)+g_2\left(x\mathrm{,}y\right)+g_3\left(x\mathrm{,}y\right)$, where auxiliary functions

$\begin{array}{l}{g}_1\left(x\right)=-u{r}_1{x}^3\left[x-\left(2x_2-a\right)\right],\\ g_2\left(x,y\right)=u{x}_2{r}_1\left(2a-x_2\right)x^2+u\alpha x_2{K}_{1}xy-a{K}_1{m}_2{y}^2,\\ g_3\left(x,y\right)=-d{K}_{1}x{y}^3+K_1\left(\alpha e-m_2\right)x{y}^2-ua{x}_2^2{r}_{1}x-ad{K}_1{y}^3-u\alpha K_1{x}^{2}y.\end{array}$