The Dynamical Behavior of a Certain Predator-Prey System with Holling Type II Functional Response

Shuangte Wang¹, Hengguo Yu^1,2, Chuanjun Dai^2,3, Min Zhao^2,3

¹College of Mathematics and Physics, Wenzhou University, Wenzhou, China.
²Key Laboratory for Subtropical Oceans & Lakes Environment & Biological Resources Utilization Technology of Zhejiang, Wenzhou University, Wenzhou, China.
³School of Life and Environmental Science, Wenzhou University, Wenzhou, China.

DOI: 10.4236/jamp.2020.83042                

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}$

Conditions in this theorem imply that discriminant of a quadratic function in $g_2\left(x,y\right)$ is

${\Delta }_2={\left(u\alpha x_2{K}_1\right)}^2+4u{x}_2\left(2a-x_2\right)a{K}_1{m}_2\le \mathrm{0,}$

thus $g_2\left(x\mathrm{,}y\right)\le 0$, i.e. $g\left(x\mathrm{,}y\right)\le 0$. By using the Krasovskii’s theorem, the proof is completed.

Theorem 7. Assume that the equilibrium $E_{\mathrm{<mo>\; *\; </mo>}}$ exists and $\eta :=\frac{\alpha y_{*}}{a\left(a+x_{*}\right)}-\frac{r_1}{K_1}\le 0$, then $E_{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotical stable.

Proof. Taking an unbounded positive definite Lyapunov function

$V\left(x,y\right)=x-x_{*}-x_{*}\ln \left(x/x_{*}\right)+u\left(y-y_2-y_{*}\ln \left(y/y_{*}\right)\right),u=\frac{a+x_{*}}{ea}\le 0,$ (35)

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

$\begin{array}{c}{D}^{+}V=\left[\frac{\alpha y_{\mathrm{<mo>\; *\; </mo>}}}{\left(a+x\right)\left(a+x_{\mathrm{<mo>\; *\; </mo>}}\right)}-\frac{r_1}{K_1}\right]{\left(x-x_{\mathrm{<mo>\; *\; </mo>}}\right)}^2-ud{\left(y-y_{\mathrm{<mo>\; *\; </mo>}}\right)}^2\\ \le \eta {\left(x-x_{\mathrm{<mo>\; *\; </mo>}}\right)}^2-ud{\left(y-y_{\mathrm{<mo>\; *\; </mo>}}\right)}^2\mathrm{.}\end{array}$ (36)

From the condition we know that $D^{+}V$ is negative definite. By using the Krasovskii’s theorem, the proof is completed.

3.6. Closed Orbits and Limit Cycles

As we see, if $r_1\le m_1$, i.e. the point $E_{\mathrm{<mo>\; *\; </mo>}}$ does not exist, there are no closed orbits in $R^{+2}$. The trace of Jacobian matrix $J\left(E_{\mathrm{<mo>\; *\; </mo>}}\right)$ is non-positive if $r_1-m_1\le m_2-\alpha e$ whether the existence condition $r_1-m_1\le 0$ holds or not. By using the Bendixson criteria, there are no closed orbits in $R^{+2}$. Such closed orbits and limit cycles surrounding the existed point $E_{\mathrm{<mo>\; *\; </mo>}}$ cannot cross the x and y coordinate axes and are confined in an invariant domain $\left[\mathrm{0,}A\right]\times \left[\mathrm{0,}B\right]$, where the upper bounds A, B are all sufficiently large. The point $E_{\mathrm{<mo>\; *\; </mo>}}$ is not a saddle point for index $j\left(E_{\mathrm{<mo>\; *\; </mo>}}\right)=1$ on this occasion. In this section, we will construct proper Dulac functions to further illustrate the existence problem of closed orbits and limit cycles.

Theorem 8. If $x_2<a$ or $r_1\le m_1$, then for system (1), there are no closed orbits in $R^{+2}$.

Proof. Taking the diffeomorphism $\phi :u=x,v=y,\text{d}\tau =\text{d}t/\left(a+x\right)$ which preserving the orientation of time, the system (1) is topologically equivalent to following system [16] [17]

$\{\begin{array}{l}\stackrel{\dot{}}{x}=P\left(x,y\right):=x\left[\left(a+x\right)\left(r_1-m_1-\frac{r_{1}x}{K_1}\right)-\alpha y\right],\\ \stackrel{\dot{}}{y}=Q\left(x,y\right):=y\left[\alpha eax-\left(a+x\right)\left(m_2+dy\right)\right],\end{array}$ (37)

For notation simplicity, we still retain the symbols x, y, t. Define a Dulac function $B\left(x\mathrm{,}y\right)=x^{-1}{y}^{-1}$, along the above system, we have

$\frac{\partial \left(BP\right)}{\partial x}+\frac{\partial \left(BQ\right)}{\partial y}\le \frac{g\left(x\right)}{y}\mathrm{,}$

where the numerator

$g\left(x\right)=\left(r_1-m_1-\frac{a{r}_1}{K_1}\right)-\frac{2r_{1}x}{K_1}\le \frac{r_1}{K_1}\left(x_2-a\right)\mathrm{.}$

By using the Bendixson-Dulac criteria, the proof is completed.

Remark. (a) If $x_2<a$, taking Dulac function $B\left(x,y\right)=\frac{1}{\left(x+c\right)y},c\ge 0$, then there are no closed orbits; (b) If $2x_2\le a$, taking Dulac function $B\left(x\mathrm{,}y\right)=\frac{a+x}{xy}$, then there are also no closed orbits.

Theorem 9. If $K_1\le 4a{m}_1\mathrm{,}{r}_1-m_1\le m_2$, and $2a{r}_1>K_1\left[2\left(r_1-m_1\right)+\alpha e-m_2\right]$, then for system (1), there are no closed orbits in $R^{+2}$.

Proof. We can take another type of Dulac function ${\text{e}}^{ux},u>0$ for our discussion. Obviously, along system (37), we have

$\frac{\partial \left(BP\right)}{\partial x}+\frac{\partial \left(BQ\right)}{\partial y}=\frac{1}{K_1}{\displaystyle \underset{i+j=0}{\overset{3}{\sum }}}{a}_{ij}{x}^i{y}^j,$ (38)

where these coefficients $a_{30},a_{20},a_{10},\cdots$ in the summation are

$\begin{array}{l}{a}_{30}=-u{r}_1,\\ a_{20}:=\left(-a{r}_1-K_1{m}_1+K_1{r}_1\right)u-3r_1,\\ a_{10}:=a{K}_1\left(r_1-m_1\right)u+\left(\alpha e-m_2\right)K_1-2a{r}_1+2\left(r_1-m_1\right)K_1,\\ ⋮\\ a_{00}=-K_{1}a\left(m_1+m_2-r_1\right),\end{array}$

For a quadratic function $a_{30}{x}^3+K_1{r}_{1}u{x}^2-a{K}_{1}u{m}_{1}x$ of variable x, its discriminant implies that it is non-positive. Since other coefficients preserved are all non-positive, this completes the proof.

Remark. For the system (0) and the above form Dulac function ${\text{e}}^{ux},u>0$, if following conditions hold: $K_1\le 2a\mathrm{,}{r}_1-m_1\le m_2$, and $4a{r}_1>K_1\left[\left(r_1-m_1\right)+\alpha e-m_2\right]$, then there are no closed orbits. This procedure can also derive our theorem while it is uncomplicated but tedious.

Corollary 1. If $E_{\mathrm{<mo>\; *\; </mo>}}$ is asymptotically stable, then combining with conditions in theorem 8 or 9, $E_{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotical stable.

Remark. See reference [18] for the above corollary.

Theorem 10. If $x_2>a$ and the equilibrium point $E_{\mathrm{<mo>\; *\; </mo>}}$ is unstable, then system (1) has limit cycle(s).

Proof. We see that ${\frac{\text{d}x}{\text{d}t}|}_{x=x_2}=-\frac{\alpha x_{2}y}{a+x_2}<0$ with $y>0$. Thus the straight line

$x=x_2$ is an untangent line of the system (37). Taking the Dulac function $w\left(x\mathrm{,}y\right)=ex+y-l$, where $l>0$ is sufficiently large, and computing $w=0$ along the orbits of Equations (37), we derive

$\frac{\text{d}w}{\text{d}t}=-m_2\left(l-ex\right)-d{\left(l-ex\right)}^2+e{r}_{1}x\left(1-\frac{x}{K_1}\right)-e{m}_{1}x<0,$

where $0<x<x_2$. For this system, we can construct a Bendixson ring including unstable singular point $E_{\mathrm{<mo>\; *\; </mo>}}$. By the Poincare-Bendixson theorem, system (37) has at least one limit cycle in the first quadrant. This completes the proof.

4. Bifurcation Analysis

4.1. Transcritical Bifurcation

Without loss of generality, denoting the system (1) as

$\left(\begin{array}{c}\stackrel{\dot{}}{x}\\ \stackrel{\dot{}}{y}\end{array}\right)=f\left(x\mathrm{,}y\right)\mathrm{.}$ (39)

Firstly, by fixing values of all the parameters and varying the bifurcation parameter $m_1$, we observed that the two equilibrium points $E_0$ and $E_2$ collide with each other as $m_1$ cross the critical value $m_1=r_1$. Thus, there is a chance of bifurcation around this point $E_0$.

Theorem 11. The system (39) undergoes a transcritical bifurcation around $E_0$ when the parameters satisfy $m_1=m_1^{\left[TC\right]}:=r_1$.

Proof. For the Jacobian matrix at $E_0$, we have eigenvectors v and w corresponding to matrix $J\left(E_0\right)$ and its transpose, respectively: $v=w=\left(\begin{array}{c}1\\ 0\end{array}\right)$. Then the following transversality conditions are hold:

$\begin{array}{l}{w}^{\text{T}}{f}_{m_1}\left(E_0,m_1^{\left[TC\right]}\right)=\left(\begin{array}{cc}1& 0\end{array}\right)\left(\begin{array}{c}0\\ 0\end{array}\right)=0,\\ w^{\text{T}}\left[D{f}_{m_1}\left(E_0,m_1^{\left[TC\right]}\right)v\right]=\left(\begin{array}{cc}1& 0\end{array}\right)\left(\begin{array}{c}-1\\ 0\end{array}\right)=-1\ne 0,\\ w^{\text{T}}\left[D^{2}f\left(E_0,m_1^{\left[TC\right]}\right)\left(v,v\right)\right]=\left(\begin{array}{cc}1& 0\end{array}\right)\left(\begin{array}{c}-\frac{\alpha }{a}\\ 0\end{array}\right)=-\frac{\alpha }{a}\ne 0.\end{array}$ (40)

Clearly, one of the eigenvalues of the Jacobian matrix $J\left(E_0\right)$ is zero and the other is negative. From the Sotomayor’s theorem( see [19] [20]), the system (39) undergoes a transcritical bifurcation around $E_0$ at $m_1=m_1^{\left[TC\right]}$. This completes the proof.

Here we take a as bifurcation parameter for the consideration of $J\left(E_2\right)$.

Theorem 12. Assume that the parameters satisfy conditions for the existence of $E_2$ and $m_1<r_1$, then the system (39) undergoes a transcritical bifurcation around $E_2$ when the parameters satisfy:

$a=a^{\left[TC\right]}:=\frac{\left(\alpha e-m_2\right)x_2}{m_2}.$ (41)

Proof. For the Jacobian matrix at $E_2$, we have eigenvectors v and w corresponding to matrix $J\left(E_2\right)$ and its transpose, respectively: $v=\left(\begin{array}{c}1\\ \frac{m_1-r_1}{m_2}\end{array}\right)$, $w=\left(\begin{array}{c}0\\ 1\end{array}\right)$. Then the following transversality conditions are hold:

$\begin{array}{l}{w}^{\text{T}}{f}_a\left(E_0,a^{\left[TC\right]}\right)=\left(\begin{array}{cc}0& 1\end{array}\right)\left(\begin{array}{c}0\\ 0\end{array}\right)=0,\\ w^{\text{T}}\left[D{f}_a\left(E_0,a^{\left[TC\right]}\right)v\right]=\left(\begin{array}{cc}0& 1\end{array}\right)\frac{m_1-r_1}{e\left(a^{\left[TC\right]}+x_2\right)}\left(\begin{array}{c}1\\ -e\end{array}\right)=\frac{m_1-r_1}{e\left(a^{\left[TC\right]}+x_2\right)}\ne 0,\\ w^{\text{T}}\left[D^{2}f\left(E_0,a^{\left[TC\right]}\right)\left(v,v\right)\right]=\frac{2\left(m_1-r_1\right)}{m_2}\left[\frac{\alpha ea}{{\left(a^{\left[TC\right]}+x_2\right)}^2}+\frac{d\left(r_1-m_1\right)}{m_2}\right]<0.\end{array}$ (42)

It is clear to see that one of the eigenvalues of the Jacobian matrix $J\left(E_2\right)$ is zero and the other is negative. From the Sotomayor’s theorem, the system (39) undergoes a transcritical bifurcation around $E_2$ at $a=a^{\left[TC\right]}$. This completes the proof.

4.2. Hopf Bifurcation

We choose parameter d as a bifurcation parameter. Assume that the parameters satisfy conditions for the existence of $E_{\mathrm{<mo>\; *\; </mo>}}$ and

$d=d^{\left[H\right]}:=\frac{\alpha x_{*}}{{\left(a+x_{*}\right)}^2}-\frac{r_1{x}_{*}}{K_1{y}_{*}}<\sqrt{\frac{{\alpha }^{2}ea{x}_{*}}{y_{*}{\left(a+x_{*}\right)}^3}}.$ (43)

Define $\lambda \left(d\right)=u\left(d\right)\pm iv\left(d\right)$ as a pair of purely imaginary eigenvalues of matrix $J\left(E_{\mathrm{<mo>\; *\; </mo>}}\right)$, where $u\left(d\right)=A_1/2$. It is clear to see that, if following conditions hold at $d=d^{\left[H\right]}$ : $A_1=0$, $A_2>0$, and transversality condition

$u^{\prime }\left(d\right)=\frac{y_{*}\left[a^3{r}_1+\left(K_1\alpha e+4r_1{x}_{*}\right)a^2+\left(5r_1{x}_{*}^2+K_1\alpha \left(e{x}_{*}-y_{*}\right)\right)a+2r_1{x}_{*}^3\right]\alpha }{2\left(a+x_{*}\right)\left[-\left(a+x_{*}\right)\left(-K_1\alpha y_{*}+a^2{r}_1+2a{r}_1{x}_{*}+r_1{x}_{*}^2\right)d-K_{1}a{\alpha }^{2}e\right]}\ne 0,$ (44)

for instance, $r_1{x}_{\mathrm{<mo>\; *\; </mo>}}^2\ge K_1\alpha y_{\mathrm{<mo>\; *\; </mo>}}$. Then $E_{\mathrm{<mo>\; *\; </mo>}}$ changes its stability through Hopf bifurcation threshold the critical value $d=d^{\left[H\right]}$.

Theorem 13. Assume that the parameters satisfy conditions for the existence of equilibrium $E_{\mathrm{<mo>\; *\; </mo>}}$, the condition (43) and the transversality condition (44), then the system (1) undergoes a Hopf bifurcation around equilibrium $E_{\mathrm{<mo>\; *\; </mo>}}$ as parameter d passes through the bifurcation value $d^{\left[H\right]}$.

We will calculate the first Lyapunov number $\sigma$ at the point $E_{\mathrm{<mo>\; *\; </mo>}}$ of the system which is used to determine the stability of limit cycles around Hopf bifurcation. Therefore we first translate the point $E_{\mathrm{<mo>\; *\; </mo>}}$ to the origin by a non-singular linear transformation $\hat{x}=x-x_{\mathrm{<mo>\; *\; </mo>}}$, $\hat{y}=y-y_{\mathrm{<mo>\; *\; </mo>}}$, then the system (39) in power series around the origin(drop the hats for the sake of convenience as usual) is:

$\{\begin{array}{l}\stackrel{\dot{}}{x}=a_{10}x+a_{01}y+a_{20}{x}^2+a_{11}xy+a_{02}{y}^2+a_{30}{x}^3+a_{21}{x}^{2}y+a_{12}x{y}^2+a_{03}{y}^3+F_1\left(x,y\right),\\ \stackrel{\dot{}}{y}=b_{10}x+b_{01}y+b_{20}{x}^2+b_{11}xy+b_{02}{y}^2+b_{30}{x}^3+b_{21}{x}^{2}y+b_{12}x{y}^2+b_{03}{y}^3+F_2\left(x,y\right),\end{array}$ (45)

with coefficients

$\begin{array}{l}{a}_{10}=-\frac{r_1{x}_{*}}{K_1}+\frac{\alpha x_{*}{y}_{*}}{{\left(a+x_{*}\right)}^2},a_{01}=-\frac{\alpha x_{*}}{a+x_{*}},\\ a_{20}=-\frac{r_1}{K_1}+\frac{\alpha y_{*}}{{\left(a+x_{*}\right)}^2}-\frac{\alpha x_{*}{y}_{*}}{{\left(a+x_{*}\right)}^3},a_{02}=0,a_{11}=-\frac{\alpha a}{{(a+x_{*})}^2},\\ a_{30}=-\frac{\alpha y_{*}a}{{\left(a+x_{*}\right)}^4},a_{03}=0,a_{12}=0,a_{21}=\frac{\alpha a}{{\left(a+x_{*}\right)}^3};\\ b_{10}=\frac{\alpha e{y}_{*}a}{{\left(a+x_{*}\right)}^2},b_{01}=-d{y}_{*},b_{20}=-\frac{\alpha e{y}_{*}a}{{\left(a+x_{*}\right)}^3},b_{02}=-d,b_{11}=\frac{\alpha ea}{{\left(a+x_{*}\right)}^2},\\ b_{30}=\frac{\alpha e{y}_{*}a}{{\left(a+x_{*}\right)}^4},b_{03}=0,b_{12}=0,b_{21}=-\frac{\alpha ea}{{\left(a+x_{*}\right)}^3}.\end{array}$ (46)

and smooth functions $F_1\left(x,y\right)={\displaystyle {\sum }_{i+j=4}^{\infty }}{a}_{ij}{x}^i{y}^j$ and $F_2\left(x,y\right)={\displaystyle {\sum }_{i+j=4}^{\infty }}{b}_{ij}{x}^i{y}^j$. From [19] [21], we know that the first Lyapunov number for a planar system is given by

$\begin{array}{c}\sigma =-\frac{3\pi }{2a_{01}{A}_2^{3/2}}\{[a_{10}{b}_{10}\left(a_{11}^2+a_{11}{b}_{02}+a_{02}{b}_{11}\right)+a_{10}{a}_{01}\left(b_{11}^2+a_{20}{b}_{11}+a_{11}{b}_{02}\right)\\ +b_{10}^2\left(a_{11}{a}_{02}+2a_{02}{b}_{02}\right)-2a_{10}{b}_{10}\left(b_{02}^2-a_{20}{a}_{02}\right)-2a_{10}{a}_{01}\left(a_{20}^2-b_{20}{b}_{02}\right)\\ -a_{01}^2\left(2a_{20}{b}_{20}+b_{11}{b}_{20}\right)+\left(a_{01}{b}_{01}-2a_{10}^2\right)\left(b_{11}{b}_{02}-a_{11}{a}_{20}\right)]\\ -\left(a_{10}^2+a_{01}{b}_{10}\right)\left[3\left(b_{10}{b}_{03}-a_{10}{a}_{30}\right)+2a_{10}\left(a_{21}+b_{12}\right)+\left(b_{10}{a}_{12}-a_{01}{b}_{21}\right)\right]\},\end{array}$ (47)

By computing the first Lyapunov number at the point $E_{\mathrm{<mo>\; *\; </mo>}}$, we will know that the Hopf bifurcation is supercritical and limit cycle(s) is(are) stable if $\sigma <0$ ; the Hopf bifurcation is subcritical and limit cycle(s) is(are) unstable if $\sigma >0$. Hence we have to give the following numerical simulations in the next section, since the above expression is much too complicated.

5. Numerical Simulations and Conclusions

In this section we numerically give simulations of dynamical behavior between the prey and predator. In order to verify the feasibility and correctness of the theoretical derivation results, we will give some numerical simulations with a control parameter d and take $r_1=0.6$, $K_1=20$, $\alpha =0.5$, $a=1.5$, $m_1=0.2$, $m_2=0.06$ and $e=0.6$. Thus, when $d=0.05$, the equilibrium point $E_{\mathrm{<mo>\; *\; </mo>}}$ is (2.898588, 2.753889), then we can verify that the values of these parameters can meet theorem 7, which means that the equilibrium point $E_{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotical stable. The time series diagram and phase diagram of the system (1) can be seen in Figure 1. It is easy to see from Figure 1 that the population x and y are permanent and can approach a fixed value respectively. These results demonstrate that the system (1) is permanent and the equilibrium point $E_{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotical stable.

To investigate in detail how the control parameter affects the dynamical behavior of the system (1), Figure 2 and Figure 3 depict the time series diagram and

Figure 1. (a) Time series diagram of population x; (b) Time series diagram of population y; (c) Phase diagram of population x and y.

Figure 2. (a) Time series diagram of population x; (b) Time series diagram of population y; (c) Phase diagram of population x and y; (d) Enlarged phase diagram of population x and y.

phase diagram of the system (1). On the basis of theoretical derivation and numerical calculus, the critical threshold of the control parameter d is approximately 0.04910. When the value of the control parameter d is greater than 0.04910, the other solutions of the system (1) can gradually converge to the

Figure 3. (a) Time series diagram of population x; (b) Time series diagram of population y; (c) Phase diagram of population x and y; (d) Enlarged phase diagram of population x and y.

equilibrium point $E_{\mathrm{<mo>\; *\; </mo>}}$, which means that the equilibrium point $E_{\mathrm{<mo>\; *\; </mo>}}$ is asymptotical stable; this result can be seen in Figure 2 with $d=0.04915$. When the value of the control parameter d is less than 0.04910, the other solutions of the system can gradually converge to a limit cycle, which implies that the Hopf bifurcation occurs in the system (1), and this result can be seen in Figure 3 with $d=0.04905$. At the bifurcation parameter $d^{\left[H\right]}\approx 0.0490796$, the first Lyapunov number $\sigma <0$, thus the Hopf bifurcation is supercritical and a limit cycle generated by the critical point is stable.

Based on mathematical theory derivation and numerical simulation analysis, it is successful to show that dynamical behavior of this certain predator-prey system mainly depends on some critical parameters and relationships. Then, it is particularly significant to point out how the control parameter d affects the complexity of the system (1), which can lead this predator-prey system to emerge Hopf bifurcation. Moreover, these results also show that some control parameters can directly or indirectly affect the dynamic complexity of our predator-prey system.

Acknowledgements

We thank the Editor and the referee for their comments. This research is funded by the National Natural Science Foundation of China (Grant nos.31570364) and the National Key Research and Development Program of China (Grant No.2018YFE0103700). This support is greatly appreciated.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Holling, C.S. (1965) The Functional Response of Predator Density and Its Role in Mimicry and Population Regulations. Memoirs of the Entomological Society of Canada, 45, 3-60. https://doi.org/10.4039/entm9745fv
[2]	Pei, Y.Z., Chen, L.S., Zhang, Q.R. and Li, C.G. (2005) Extinction and Performance of One-Prey Multi-Predators of Holling Type II Function Response System with Impulsive Biological Control. Journal of Theoretical Biology, 235, 495-503. https://doi.org/10.1016/j.jtbi.2005.02.003
[3]	Misha, P. and Raw, S.N. (2019) Dynamical Complexities in a Predator-Prey System Involving Teams of Two Prey and One Predator. Journal of Applied Mathematics and Computing, 61, 1-24. https://doi.org/10.1007/s12190-018-01236-9
[4]	Huang, J.C., Ruan, S.G. and Song, J. (2014) Bifurcation in a Predator-Prey System of Leslie Type with Generalized Holling Type III Functional Response. Journal of Differential Equations, 257, 1721-1752. https://doi.org/10.1016/j.jde.2014.04.024
[5]	Zhang, C.M., Chen, W.C. and Yang, Y. (2006) Periodic Solutions and Global Asymptotic Stability of a Delayed Discrete Predator-Prey System with Holling II Type Functional Response. Journal of Systems Science and Complexity, 19, 449-460. https://doi.org/10.1007/s11424-006-0449-x
[6]	Rosenzweig, M.L. and MacArthur, R. (1963) Graphical Representation and Stability Conditions of Predator-Prey Interactions. American National, 97, 209-223. https://doi.org/10.1086/282272
[7]	Hsu, S.B. and Huang, T.W. (1995) Global Stability for a Class of Predator-Prey Systems. SIAM Journal on Applied Mathematics, 55, 763-783. https://doi.org/10.1137/S0036139993253201
[8]	Freedman, H.I. (1979) Stability Analysis of a Predator Prey System with Mutual Interference and Density Dependent Death Rate. Bulletin of Mathematical Biology, 41, 67-78. https://doi.org/10.1016/S0092-8240(79)80054-3
[9]	Wang, Y.Q. and Jing, Z.J. (1999) Multiple Limit Cycles and Global Stability in Predator-Prey Model. ACTA Mathematicae Applicatae Sinica, 15, 206-219. https://doi.org/10.1007/BF02720497
[10]	Chen, L.J., Chen, F.D. and Chen, L.J. (2010) Qualitative Analysis of a Predator-Prey Model with Holling Type II Functional Response Incorporating a Constant Prey Refuge. Nonlinear Analysis, 11, 246-252. https://doi.org/10.1016/j.nonrwa.2008.10.056
[11]	Lv, Y.F., Pei, Y.Z., Gao, S.J. and Li, C.G. (2010) Harvesting of a Phytoplankton-Zooplankton Model. Nonlinear Analysis, 11, 3608-3619. https://doi.org/10.1016/j.nonrwa.2010.01.007
[12]	Li, S.B., Wu, J.H. and Dong, Y.Y. (2018) Effects of a Degeneracy in a Diffusive Predator-Prey Model with Holling Functional Response. Nonlinear Analysis, 43, 78-95. https://doi.org/10.1016/j.nonrwa.2018.02.003
[13]	Srinivasu, P.D.N., Ismail, S. and Naidu, C.R. (2001) Global Dynamics and Controllability of a Harvested Prey-Predator System. Journal of Biological Systems, 9, 67-79. https://doi.org/10.1142/S0218339001000311
[14]	Chen, F. (2005) On a Nonlinear Nonautonomous Predator-Prey Model with Diffusion and Distributed Delay. Journal of Computational and Applied Mathematics, 180, 33-49. https://doi.org/10.1016/j.cam.2004.10.001
[15]	Birkhoff, G. and Rota, G.C. (1962) Ordinary Differential Equations Introductions to Higher Mathematics. Ginn and Company, Boston, MA.
[16]	Chicone, C. (2006) Ordinary Differential Equations with Applications. In: Texts in Applied Mathematics, Volume 34, World Scientific, Springer-Verlag, New York.
[17]	Andronov, A. (1973) Qualitative Theory of Second-Order Dynamic Systems. Volume 22054. Halsted Press, Canberra.
[18]	Jiang, Q.H. and Wang, J.H. (2013) Qualitative Analysis of a Harvested Predator-Prey System with Holling Type III Functional Response. Advanced in Difference Equations, 249. https://doi.org/10.1186/1687-1847-2013-249
[19]	Perko, L. (2001) Differential Equations and Dynamical Systems. In: Texts in Applied Mathematics, Volume 7, 3rd Edition, Springer-Verlag, New York. https://doi.org/10.1007/978-1-4613-0003-8
[20]	Pirayesh, B., Pazirandeh, A. and Akbari, M. (2016) Local Bifurcation Analysis in Nuclear Reactor Dynamics by Sotomayor’s Theorem. Annals of Nuclear Energy, 94, 716-731. https://doi.org/10.1016/j.anucene.2016.04.021
[21]	Kuznetsov, Y. (1998) Elements of Applied Bifurcation Theory. 2nd Edition, Springer-Verlag, New York, 96-100.