A Seasonally Perturbed Prey-Predator Model with Smith Growth Rate

Qinggui Guo

College of Science, University of Shanghai for Science and Technology, Shanghai, China.

DOI: 10.4236/ojapps.2026.164071                

Abstract

This paper investigates the dynamic behaviors of a predator-prey model with seasonal perturbations, incorporating the Smith growth pattern. First, for the unperturbed autonomous system, we conduct a detailed analysis of the existence and stability of equilibrium points, and derive the conditions under which a Hopf bifurcation occurs. Subsequently, taking into account seasonal variations in real ecosystems, we introduce a periodic perturbation by considering the prey growth rate as a time-dependent periodic function, thereby establishing a non-autonomous model. Moreover, we provide sufficient conditions for the existence of asymptotically stable periodic solutions. Numerical simulations are presented to validate the theoretical results.

Keywords

Smith Growth, Seasonal Perturbation, Hopf Bifurcation, Periodic Solution

1. Introduction

Since the formulation of the classical model by Lotka [1] and Volterra [2] in the 1920s, the study of predator-prey models has become one of the most prominent topics in ecology, leading to numerous significant findings. The general form of the predator-prey model is given by:

$\begin{array}{l}\dot{u}=ug\left(u\right)-vp\left(u,v\right),\\ \dot{v}=ep\left(u,v\right)-dv,\end{array}$ (1)

where $u\left(t\right)$ , $v\left(t\right)$ denote the population densities of prey and predator at time $t$ respectively. $g\left(u\right)$ represents the intrinsic growth rate of the prey in the absence of predators, $p\left(u,v\right)$ is the functional response function, and $e$ denotes the conversion efficiency representing the predator’s consumption of prey.

The classical predator-prey model employs the logistic growth pattern [3], defined as $g\left(u\right)=r\left(1-u\right)$ . This assumes a strictly negative correlation between the per capita growth rate and population density. Although this assumption is concise, it fails to accurately characterize population dynamics in resource-limited environments, particularly at low population densities. Consequently, in 1963, Smith [4] proposed the renowned Smith growth model:

$\dot{u}=ru\left(\frac{K-u}{K+\left(\frac{r}{c}\right)u}\right).$ (2)

Smith pointed out that even when population growth ceases, individuals still consume energy to sustain basal metabolism and replace deceased members; this portion of consumption is termed the “subsistence rate ($c$ )”. Furthermore, to describe the impact of seasonal variations on population dynamics, and incorporating the Smith growth model with reference to [5]-[8], we establish the following dynamical system:

$\begin{array}{l}\dot{u}=r\left(t\right)u\left(\frac{1-\frac{u}{k}}{1+\frac{\theta u}{k}}\right)-\frac{\beta uv}{\alpha +u},\\ \dot{v}=sv\left(1-\frac{hv}{u}\right).\end{array}$ (3)

All parameters are positive constants. Specifically, $r\left(t\right)=r\left(1+\gamma \sin \left(2\pi t\right)\right)$ denotes the intrinsic growth rate of the prey, $k$ is the environmental carrying capacity, $\theta$ represents the resource limitation parameter, $\beta$ is the maximum predation rate, and $\alpha$ is the half-saturation constant. For convenience, we introduce the following nondimensional transformation for system (3):

$x=\frac{1}{k}u,y=\frac{h}{k}v,\tau =st,$ (4)

then system (3) becomes (still denote $\tau$ by $t$ ):

$\begin{array}{l}\dot{x}=\frac{\overline{r}x\left(1-x\right)}{1+ax}-\frac{cxy}{x+b},\\ \dot{y}=y\left(1-\frac{y}{x}\right),\end{array}$ (5)

where

$\overline{r}=\frac{r\left(t\right)}{s},a=\theta ,b=\frac{\alpha }{k},c=\frac{\beta }{hs},$

of course, such a transformation is not unique. The transformed system takes a more concise form, which simplifies our analysis, and the rescaled parameters retain biological interpretability—for instance, $\overline{r}$ represents the efficiency of energy transfer from prey to predator.

2. Autonomous System

We first investigate the autonomous system of (5) by setting $\gamma =0$ and replacing $\overline{r}$ with $r$ :

$\begin{array}{l}\dot{x}=\frac{rx\left(1-x\right)}{1+ax}-\frac{cxy}{x+b},\\ \dot{y}=y\left(1-\frac{y}{x}\right).\end{array}$ (6)

For practical considerations, we restrict our study to the positive solutions of system (6), i.e.,${ℝ}_{+}^2=\left\{\left(x,y\right)|x\ge 0,y\ge 0\right\}$ . Furthermore, it is straightforward to obtain that the positively invariant region of system (6) is given by:

$\Omega =\left\{\left(x,y\right)|0\le x\le 1,0\le y\le 1\right\}.$ (7)

2.1. Equilibria

Since system (6) is undefined at $x=0$ , it is evident that the system possesses a unique boundary equilibrium point $E_1=\left(1,0\right)$ . Furthermore, any positive equilibrium point of system (6) must satisfy:

$\begin{array}{l}f\left(x,y\right)=\frac{r\left(1-x\right)}{1+ax}-\frac{cy}{x+b}=0,\\ g\left(x,y\right)=y\left(1-\frac{y}{x}\right)=0.\end{array}$ (8)

From the second equation, we get $y=x$ , substituting it into the first equation, we get:

$\left(-r-ac\right)x^2+\left(r-c-br\right)x+br=0,$ (9)

this is a quadratic equation, by applying the quadratic formula, we can readily compute that $\text{Δ}={\left(c-r+br\right)}^2+br\left(4r+4ac\right)>0$ ; hence, a real root necessarily exists. Moreover, by examining the relationships among the coefficients, we conclude that this root is unique; therefore, system (6) has exactly one positive equilibrium point $E_0=\left(x_0,y_0\right)$ where

$x_0=-\frac{c-r+br-\sqrt{c^2+\left(2b+4ab-2\right)cr+\left(4b+{\left(b-1\right)}^2\right)r^2}}{2r+2ac},y_0=x_0.$ (10)

Next, we will discuss the stability of the equilibrium point $E_0$ . After some simple calculations, we obtain the Jacobi matrix of system (6) at the equilibrium point $E_0$ .

$J\left(E_0\right)=\left(\begin{array}{cc}T\left(x_0\right)& -\frac{c{x}_0}{b+x_0}\\ 1& -1\end{array}\right),$ (11)

where $T\left(x_0\right)=-\frac{r\left(a{x}_0^2+2x_0-1\right)}{{\left(1+a{x}_0\right)}^2}-\frac{c{x}_{0}b}{{\left(x_0+b\right)}^2}$ , accordingly, we obtain its characteristic equation:

${\lambda }^2-tr\left(J\left(E_0\right)\right)\lambda +det\left(J\left(E_0\right)\right)=0,$ (12)

$\begin{array}{l}tr\left(J\left(E_0\right)\right)=T\left(x_0\right)-1,\\ det\left(J\left(E_0\right)\right)=-T\left(x_0\right)+\frac{c{x}_0}{b+x_0}.\end{array}$ (13)

Define

$a^{\text{*}}=\frac{1-2x_0}{x_0^2},$ (14)

$\begin{array}{l}{r}_1=-\frac{{\left(a{x}_0+1\right)}^2\left(2b{x}_0+b^2+x_0^2+bc{x}_0\right)}{{\left(b+x_0\right)}^2\left(a{x}_0^2+2x_0-1\right)},\\ r_2=-\frac{c{x}_0{\left(a{x}_0+1\right)}^2\left(2b+x_0\right)}{{\left(b+x_0\right)}^2\left(a{x}_0^2+2x_0-1\right)}.\end{array}$ (15)

If the real parts of both eigenvalues of (12) are negative, then the equilibrium point $E_0$ is asymptotically stable; therefore, we have the following conclusion:

Theorem 1. The system (8) has exactly one positive equilibrium point $E_0\left(x_0,y_0\right)$ ; furthermore,

1) If $a\ge a^{\text{*}}$ , then equilibrium point $E_0$ is locally asymptotically stable.

2) If $0<a<a^{\text{*}}$ , then when $0<r<\text{min}\left\{r_1,r_2\right\}$ , the equilibrium point $E_0$ is locally asymptotically stable; otherwise, it is unstable.

Proof. From (12), we know that the eigenvalues of the equilibrium point $E_0$ all have negative real parts if and only if

$\begin{array}{l}tr\left(J\left(E_0\right)\right)<0,\\ det\left(J\left(E_0\right)\right)>0,\end{array}$ (16)

if $a\ge a^{\text{*}}$ , then $T\left(x_0\right)<0$ , here must exist $tr\left(J\left(E_0\right)\right)<0$ and $det\left(J\left(E_0\right)\right)>0$ . If $0<a<a^{\text{*}}$ , then $T\left(x_0\right)>0$ , after a simple calculation, it follows that for the above equation to hold, the following must be true:

$0<r<\text{min}\left\{r_1,r_2\right\},$

the values of $r_1$ and $r_2$ depend on $x_0$ , this concludes the proof.

2.2. Hopf Bifurcation

The analysis in the previous section shows that system (6) possesses a unique equilibrium point. Consequently, neither a saddle-node bifurcation nor a Bogdanov-Takens (B-T) bifurcation can occur. In this section, we prove that under certain parameter conditions, the equilibrium $E_0$ can be a central equilibrium, and system (6) may undergo a Hopf bifurcation at $E_0$ .

From (9) and $tr\left(J\left(E_0\right)\right)=0$ , we can derive expressions for $r$ and $c$ in terms of $a$ , $b$ , and $x_0$ :

$\begin{array}{l}r=r_0=-\frac{\left(b+x_0\right){\left(a{x}_0+1\right)}^2}{x_0\left(a{x}_0^2+2x_0+b+ab-1\right)},\\ c=c_0=\frac{{\left(b+x_0\right)}^2\left(a{x}_0+1\right)\left(x_0-1\right)}{x_0^2\left(a{x}_0^2+2x_0+b+ab-1\right)}.\end{array}$ (17)

To ensure all parameters are positive, the following assumptions must be made:

${\tau }_H:=\left\{\left(a,b,x_0\right)\in {ℝ}_{+}^3:0<a<\frac{1-2x_0-b}{x_0^2+b},0<b<1-2x_0,0<x_0<\frac{1}{2}\right\}.$ (18)

Next, we will discuss the Hopf bifurcation of system (6) in the neighbourhood of $E_0$ , and the formation of the limit cycle when the parameters $\left(a,b,x_0\right)\in {\tau }_H$ lie in the Hopf set and $r=r_0$ , $c=c_0$ . First, apply the following transformation to system (6):

$x=x_1+x_0,y=y_1+y_0,$ (19)

then system (6) becomes;

$\begin{array}{l}{\dot{x}}_1=a_{10}{x}_1+a_{01}{y}_1+a_{20}{x}_1^2+a_{11}{x}_1{y}_1+a_{30}{x}_1^3+a_{21}{x}_1^2{y}_1+o\left({\left|x_1,y_1\right|}^3\right),\\ {\dot{y}}_1=b_{10}{x}_1+b_{01}{y}_1+b_{20}{x}_1^2+b_{11}{x}_1{y}_1+b_{02}{y}_1^2+b_{30}{x}_1^3+b_{12}{x}_1{y}_1^2+b_{21}{x}_1^2{y}_1+o\left({\left|x_1,y_1\right|}^3\right),\end{array}$ (20)

where $a_{ij}$ , $b_{ij}$ are given in Appendix. let $\delta =a_{10}{b}_{01}-a_{01}{b}_{10}$ . Clearly, when the Hopf condition is satisfied,

$\delta =-\frac{b+a{x}_0{}^2+x_0^2+2ab{x}_0-ab{x}_0^2}{x_0\left(a{x}_0^2+2x_0+b+ab-1\right)}>0,$

then let

$x_1=-\frac{1}{b_{10}}u-\frac{a_{10}}{\sqrt{\delta }}v,t=\frac{1}{\sqrt{\delta }}\tau ,$ (21)

then system (20) is transformed into

$\begin{array}{l}{\dot{x}}_1=y_1+c_{20}{x}_1^2+c_{11}{x}_1{y}_1+c_{02}{y}_1^2+c_{30}{x}_1^3+c_{12}{x}_1{y}_1^2+c_{21}{x}_1^2{y}_1+c_{03}{y}_1^3+o\left({\left|x_1,y_1\right|}^3\right),\\ {\dot{y}}_1=-x_1+d_{20}{x}_1^2+d_{30}{x}_1^3+d_{21}{x}_1^2{y}_1+o\left({\left|x_1,y_1\right|}^3\right),\end{array}$ (22)

where $c_{ij}$ , $d_{ij}$ are given in Appendix. Using the formal series method in Zhang [9], we can obtain the first focal values (or Lyapunov coefficients) as follows:

$W_1=\frac{\left(x_0-1\right)W_{11}}{4x_0^2\left(b+x_0\right)\left(a{x}_0+1\right)\sqrt{-\frac{{\sigma }_1}{x_0\left(a{x}_0^2+2x_0+b+ab-1\right)}}\left(a{x}_0^2+2x_0+b+ab-1\right){\sigma }_1^2},$ (23)

where ${\sigma }_1$ , $W_{11}$ are given in Appendix. According to Poincare-Andronov-Hopf bifurcation theorem, we can draw the following theorem

Theorem 2. If the parameters of system (6) satisfy $r=r_0$ , $c=c_0$ , and $\left(a,b,x_0\right)\in {\tau }_H$ , then,

1) If $W_1<0$ , system (6) undergoes a supercritical Hopf bifurcation, giving rise to a stable limit cycle near the equilibrium $E_0$ .

2) If $W_1>0$ , system (6) undergoes a subcritical Hopf bifurcation, giving rise to an unstable limit cycle near the equilibrium $E_0$ .

3) If $W_1=0$ , system (6) undergoes a degenerate Hopf bifurcation.

3. Perturbation System

In this section we consider model (5) and assume that the growth rate of prey is a periodic function with respect to t, $\overline{r}=\frac{r\left(t\right)}{s}=r^{\prime }\left(1+\gamma \sin \left(2\pi t\right)\right)$ , where $\gamma \in \left(0,1\right)$ , stands for Seasonal intensity. The growth rate reaches a maximum rate $r^{\prime }\left(1+\gamma \right)$ at time $t=\frac{1}{4}+i$ , where $i$ is an integer (representing the year), and a minimum value $r^{\prime }\left(1+\gamma \right)$ when $t=\frac{3}{4}+i$ . Nest, we will theoretically analyze the existence of stable periodic solutions of system (5), We rewrite system (5) as

$\dot{Y}=f\left(Y\right)+\gamma g\left(t,Y\right),$ (24)

where

$Y={\left(x,y\right)}^{\text{T}},$

$f\left(Y\right)={\left(r^{\prime }x\left(\frac{1-x}{1+ax}\right)-\frac{cxy}{b+x},y\left(1-\frac{y}{x}\right)\right)}^{\text{T}},$

$g\left(t,Y\right)={\left(r^{\prime }\sin \left(2\pi t\right)x\left(\frac{1-x}{1+ax}\right),0\right)}^{\text{T}}.$

Before entering the main content of this section, we need to use the following lemma, which is from Brauer [10].

Lemma 3. Let $f\left(Y\right)$ and $g\left(t,Y\right)$ be continuously differentiable with respect to the components of $Y$ and let $g\left(t,Y\right)$ be periodic in $t$ with period $w$ for each $Y$ . Let $x_{\infty }$ be a critical point of undisturbed system, i.e. a solution of the equation $f\left(x_{\infty }\right)=0$ , which is asymptotically stable in the strong sense that the eigenvalues of the matrix $f_Y\left[x_{\infty }\right]$ all have negative real part. Then the perturbed system has an asymptotically stable periodic solution $p\left(t,\gamma \right)$ of the same period $w$ for all sufficiently small $\gamma$ , with:

$\underset{\gamma \to 0}{\text{lim}}p\left(t,\gamma \right)=x_{\infty }.$ (25)

Applying Theorem 1 with the Lemma 3, it should be noted that Lemma 3 is applicable only within the region of asymptotically stable equilibria. Consequently, we have the following theorem:

Theorem 4. When $\gamma =0$ , then system (24) is equivalent to system (6). In this case, the system has a unique positive equilibrium, if this equilibrium is locally asymptotically stable, then, after introducing perturbations, for all sufficiently small perturbations $\gamma$ , the perturbed system (24) possesses an asymptotically stable periodic solution $p\left(t,Y\right)$ , and

$\underset{\gamma \to 0}{\text{lim}}p\left(t,\gamma \right)=E_0$ (26)

Proof. Since the local asymptotic stability of the equilibrium point of the unperturbed system is equivalent to all its eigenvalues having negative real parts, the condition in Lemma 3 is equivalent to that in Theorem 1. Based on the preceding analysis, the proof of Theorem 4 is complete.

Theorem 4 indicates that if the unperturbed system possesses a locally asymptotically stable equilibrium, the introduction of small perturbations leads to an asymptotically stable periodic solution.

4. Numerical Simulations

Next, we will validate our analysis through numerical simulations. Figure 1 illustrates the Hopf branching process for system (6). We selected a set of parameter values, the parameters for Figure 1(a) are: $r=2.5$ ,$a=0.16$ ,$b=0.2$ ,$c=3.9137$ ; the parameters for Figure 1(b) are: $r=3.4$ ,$a=0.16$ ,$b=0.1$ ,$c=3.9137$ . As can be seen, when the parameters change, a stable limit cycle appears near the equilibrium point $E_0$ .

The Poincaré map of system (24) in the (x, y)-plane and the corresponding attractor are shown in Figure 2(a) and Figure 2(b), respectively, with parameters $r=1$ ,$a=0.4$ ,$b=0.39$ ,$c=0.56$ , and $\gamma =0.1$ , and initial density $\left(x_0,y_0\right)=\left(0.56,0.56\right)$ , which lies within the basin of attraction of the stable focus $E_0$ . The attractor of the Poincaré map is a fixed point. Figure 3 displays the corresponding time series diagram of the prey population.

(a) A stable focal point $E_0$ (b) A stable limit cycle

Figure 1. Hopf bifurcation of system (6).

Figure 2. The Poincaré map of system (24).

Figure 3. The time series of the prey.

As shown in Figures 4-5 and Figures 6-7, we present the phase diagrams and time series for the two groups, respectively. The parameters are as above, but the initial value is $\left(x_0,y_0\right)=\left(0.689,0.689\right)$ . When $\gamma =0.01$ , a stable periodic solution is clearly observed, indicating oscillatory coexistence of the two species. When $\gamma =0$ , the periodic solution vanishes, and both species eventually approach a stable equilibrium state, the equilibrium point $E_0$ becomes a stable focus. Figure 5(a) illustrates the phase space trajectories, where the orbit originating from the initial value is observed to converge to the equilibrium point $E_0$ .

Figure 4. The phase portrait of system (24).

(a) (b)

Figure 5. $\gamma =0.01,\left(x_0,y_0\right)=\left(0.689,0.689\right)$ .

Figure 6. A curve-fitting equation.

(a) (b)

Figure 7. $\gamma =0,\left(x_0,y_0\right)=\left(0.689,0.689\right)$ .

5. Conclusions

This study investigates the dynamics of a predator-prey model incorporating Smith growth and seasonal perturbations, analyzing both the unperturbed and perturbed systems. Our analysis demonstrates that the system possesses a unique positive equilibrium point in all cases, with stable coexistence achievable under specific parameter conditions. Notably, parameters $r$ and $a$ play a critical role in determining system stability. There exists a threshold for parameter $a$ : when $a$ exceeds this threshold, the stability of the equilibrium is guaranteed, this implies that when resource limitations are sufficiently strong, the populations of both species will eventually converge to a stable state, thereby achieving long-term coexistence; conversely, when $a$ falls below this value, the dynamics change, and the equilibrium point may lose its stability. In this regime, the stability of the equilibrium is contingent upon the energy conversion efficiency $r$ .

Furthermore, by analyzing the distribution of characteristic equation roots and calculating the first Lyapunov coefficient, we verify that the system undergoes both supercritical and subcritical Hopf bifurcations near critical values, leading to the emergence of stable limit cycles. Under specific parametric conditions, the system undergoes a supercritical Hopf bifurcation. At this critical point, the originally stable equilibrium loses its stability and transitions smoothly into a stable periodic oscillation. Consequently, the populations of both species exhibit sustained, regular fluctuations, thereby maintaining long-term stable coexistence. Conversely, under alternative parametric conditions, the system undergoes a subcritical Hopf bifurcation. This scenario is far more precarious, potentially leading to a catastrophic collapse of the population.

Upon introducing seasonal perturbations, we also prove the existence of periodic solutions under appropriate conditions. It is imperative to note that in the limit where perturbation intensity approaches zero, a system initially stationary at an equilibrium point transitions into one characterized by infinitesimal periodic oscillations near said equilibrium, precluding any unbounded behavior. Numerical simulations are provided to validate our theoretical analysis, illustrating the transition from stable equilibria to periodic oscillations. Consequently, the prey growth rate, the Smith term, and seasonal perturbations exert significant influences on population persistence and stability.

Appendix

$a_{10}=1,a_{01}=-\frac{\left(b+x_0\right)\left(a{x}_0+1\right)\left(x_0-1\right)}{x_0\left(a{x}_0^2+2x_0+b+ab-1\right)},b_{10}=1,b_{01}=-1,$

$a_{20}=\frac{a^{2}b{x}_0^3-a^{2}b{x}_0^2+a{b}^2+2ab{x}_0^2+a{x}_0^2+b^2+3b{x}_0-b+x_0^2}{x_0\left(b+x_0\right)\left(a{x}_0+1\right)\left(a{x}_0^2+2x_0+b+ab-1\right)},$

$a_{11}=-\frac{b\left(a{x}_0+1\right)\left(x_0-1\right)}{x_0^2\left(a{x}_0^2+2x_0+b+ab-1\right)},$

$a_{30}=-\frac{a{b}^3+x_0{a}_{300}+a^2{b}^3}{x_0{\left(b+x_0\right)}^2{\left(a{x}_0+1\right)}^2\left(a{x}_0^2+2x_0+b+ab-1\right)},$

$a_{21}=\frac{b\left(a{x}_0+1\right)\left(x_0-1\right)}{x_0^2\left(b+x_0\right)\left(a{x}_0^2+2x_0+b+ab-1\right)},b_{20}=-\frac{1}{x_0},b_{11}=\frac{2}{x_0},b_{02}=-\frac{1}{x_0},$

$b_{30}=\frac{1}{x_0^2},b_{12}=\frac{1}{x_0^2},b_{21}=-\frac{2}{x_0^2},$

$a_{300}=b-3ab+3a{b}^2+x_0\left(6ab+x_0\left(a+3a^{2}b-a^{3}b+a^2+a^{3}b{x}_0\right)\right),$

$c_{20}=\frac{c_{200}}{x_0\left(b+x_0\right)\left(a{x}_0+1\right)\sqrt{-\frac{b+a{x}_0^2+x_0^2+2ab{x}_0-ab{x}_0^2}{x_0{c}_{201}}}{c}_{201}},$

$c_{11}=-\frac{c_{111}}{x_0\left(b+x_0\right)\left(a{x}_0+1\right)\left(b+a{x}_0^2+x_0^2+2ab{x}_0-ab{x}_0^2\right)},$

$c_{02}=-\frac{x_0\left(3a{b}^2+x_0\left(2b+2ab+x_0\left(-a^2{b}^2+a+1\right)-2a{b}^2+a^2{b}^2\right)\right)+b^2}{x_0\left(b+x_0\right)\left(a{x}_0+1\right)\sqrt{-\frac{c_{021}}{x_0\left(a{x}_0{}^2+2x_0+b+ab-1\right)}}{c}_{021}},$

$c_{30}=\frac{\frac{1}{x_0^2}+\frac{a\left(b+x_0\right)\left(a+1\right)}{x_0{\left(a{x}_0+1\right)}^2{c}_{300}}+\frac{b\left(a{x}_0+1\right)\left(x_0-1\right)}{x_0{\left(b+x_0\right)}^2{c}_{300}}}{\sqrt{-\frac{b+a{x}_0^2+x_0^2+2ab{x}_0-ab{x}_0^2}{x_0{c}_{300}}}},$

$c_{12}=-\frac{c_{121}}{x_0{\left(b+x_0\right)}^2{\left(a{x}_0+1\right)}^2\sqrt{-\frac{c_{122}}{x_0\left(a{x}_0^2+2x_0+b+ab-1\right)}}{c}_{122}},$

$c_{21}=-\frac{c_{210}}{x_0{\left(b+x_0\right)}^2{\left(a{x}_0+1\right)}^2\left(b+a{x}_0^2+x_0^2+2ab{x}_0-ab{x}_0^2\right)},$

$c_{03}=\frac{c_{031}}{x_0{\left(b+x_0\right)}^2{\left(a{x}_0+1\right)}^2\frac{c_{032}}{x_0\left(a{x}_0^2+2x_0+b+ab-1\right)}{c}_{032}},$

$d_{20}=-\frac{1}{x_0},d_{30}=-\frac{1}{x_0^2},d_{21}=-\frac{1}{x_0^2\sqrt{\frac{\left(b+x_0\right)\left(a{x}_0+1\right)\left(x_0-1\right)}{x_0\left(a{x}_0^2+2x_0+b+ab-1\right)}-1}},$

$\begin{array}{c}{c}_{200}=a^2{b}^2{x}_0+2a^{2}b{x}_0^3+a^2{x}_0^4+a{b}^2{x}_0+2a{b}^2+6ab{x}_0^2\\ +3a{x}_0^3+2b^2+6b{x}_0-2b+3x_0^2-x_0,\end{array}$

$c_{201}=a{x}_0^2+2x_0+b+ab-1,$

$\begin{array}{c}{c}_{111}=-a^2{b}^2{x}_0^3+a^2{b}^2{x}_0^2+a^{2}b{x}_0^4-a^{2}b{x}_0^3-2a{b}^2{x}_0^2+4a{b}^2{x}_0+2ab{x}_0^3\\ +2ab{x}_0^2+2a{x}_0^3+b^2{x}_0+b^2+5b{x}_0^2-b{x}_0+2x_0^3,\end{array}$

$c_{021}=b+a{x}_0^2+x_0^2+2ab{x}_0-ab{x}_0^2,$

$c_{300}=a{x}_0^2+2x_0+b+ab-1,$

$\begin{array}{c}{c}_{121}=-2a^3{b}^2{x}_0^4+2a^3{b}^2{x}_0^3+a^{3}b{x}_0^5-a^{3}b{x}_0^4+3a^2{b}^3{x}_0-6a^2{b}^2{x}_0^3\\ +15a^2{b}^2{x}_0^2+3a^{2}b{x}_0^4+6a^{2}b{x}_0^3+3a^2{x}_0^4+3a{b}^3{x}_0+3a{b}^2{x}_0^2\\ +6a{b}^2{x}_0+12ab{x}_0^3-3ab{x}_0^2+3a{x}_0^4-2b^2{x}_0+2b^2+b{x}_0^2-b{x}_0,\end{array}$

$c_{122}=b+a{x}_0^2+x_0^2+2ab{x}_0-ab{x}_0^2,$

$\begin{array}{c}{c}_{210}=a^3{b}^3{x}_0^2+3a^3{b}^2{x}_0^3+4a^{3}b{x}_0^5-a^{3}b{x}_0^4+a^3{x}_0^6+a^2{b}^3{x}_0^2\\ +5a^2{b}^3{x}_0+3a^2{b}^2{x}_0^3+15a^2{b}^2{x}_0^2+15a^{2}b{x}_0^4+3a^{2}b{x}_0^3+4a^2{x}_0^5\\ +2a^2{x}_0^4+5a{b}^3{x}_0+a{b}^3+15a{b}^2{x}_0^2+3a{b}^2{x}_0+27ab{x}_0^3-9ab{x}_0^2\\ +8a{x}_0^4-2a{x}_0^3+b^3+3b^2{x}_0+7b{x}_0^2-4b{x}_0+2x_0^3-x_0^2\mathrm{,}\end{array}$

$\begin{array}{c}{W}_{11}=-a^4{b}^4{x}_0^3-a^4{b}^3{x}_0^6+3a^4{b}^3{x}_0^5-3a^4{b}^3{x}_0^4+a^4{b}^2{x}_0^6+a^{4}b{x}_0^6\\ -3a^3{b}^4{x}_0^2-3a^3{b}^3{x}_0^5+12a^3{b}^3{x}_0^4-11a^3{b}^3{x}_0^3+a^3{b}^2{x}_0^6\\ +6a^3{b}^2{x}_0^5-3a^3{b}^2{x}_0^4+2a^{3}b{x}_0^6-a^3{x}_0^6+a^2{b}^4{x}_0^3-3a^2{b}^4{x}_0^2\\ +7a^2{b}^3{x}_0^3-3a^2{b}^3{x}_0^2+6a^2{b}^2{x}_0^5-3a^2{b}^2{x}_0^4+9a^2{b}^2{x}_0^3\\ +a^{2}b{x}_0^6-3a^{2}b{x}_0^5+6a^{2}b{x}_0^4-2a^2{x}_0^6+a{b}^4-2a{b}^3{x}_0^3+12a{b}^3{x}_0^2\\ +10a{b}^2{x}_0^3+6a{b}^2{x}_0^2-3ab{x}_0^5+6ab{x}_0{}^4+3ab{x}_0^3-a{x}_0^6\\ +b^4+6b^3{x}_0-b^3+b^2{x}_0^3+6b^2{x}_0^2+3b{x}_0^3\end{array}$

$\begin{array}{c}{c}_{031}=-a^3{b}^2{x}_0^4+a^3{b}^2{x}_0^3+a^2{b}^3{x}_0-3a^2{b}^2{x}_0^3+6a^2{b}^2{x}_0^2+3a^{2}b{x}_0^3\\ +a^2{x}_0^4+a{b}^3{x}_0+3a{b}^2{x}_0+3ab{x}_0^3+a{x}_0^4-b^2{x}_0+b^2,\end{array}$

$c_{032}=b+a{x}_0^2+x_0^2+2ab{x}_0-ab{x}_0^2.$

${\sigma }_1=b+a{x}_0^2+x_0^2+2ab{x}_0-ab{x}_0^2.$

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Lotka, A.J. (1925) Elements of Physical Biology. Williams & Wilkins.
[2]	Volterra, V. (1926) Variazioni e fluttuazioni del numero d’individui in specie animali conviventi. Società anonima tipografica “Leonardo da Vinci”.
[3]	Verhulst, P.F. (1838) Notice sur la loi que la population suit dansson accroissement. Correspondence mathematique et physique, 10, 113-129.
[4]	Smith, F.E. (1963) Population Dynamics in Daphnia Magna and a New Model for Population Growth. Ecology, 44, 651-663.[CrossRef]
[5]	Holling, C.S. (1965) The Functional Response of Predators to Prey Density and Its Role in Mimicry and Population Regulation. Memoirs of the Entomological Society of Canada, 97, 5-60.[CrossRef]
[6]	Du, Z.J., Chen, X. and Feng, Z.S. (2014) Multiple Positive Periodic Solutions to a Predator-Prey Model Withleslie-Gower Holling-Type II Functional Response and Harvesting Terms. Discrete and Continuous Dynamical Systems S, 7, 1203-1214.[CrossRef]
[7]	Mukherjee, D. (2003) Persistence and Bifurcation Analysis on a Predator-Prey System of Holling Type. ESAIM: Mathematical Modelling and Numerical Analysis, 37, 339-344.[CrossRef]
[8]	Leslie, P.H. and Gower, J.C. (1960) The Properties of a Stochastic Model for the Predator-Prey Type of Interaction between Two Species. Biometrika, 47, 219-234.[CrossRef]
[9]	Zhang, Z.F. (1992) Qualitative Theory of Differential Equations, Volume 101. American Mathematical Society.
[10]	Brauer, F. (1977) Periodic Solutions of Some Ecological Models. Journal of Theoretical Biology, 69, 143-152.[CrossRef] [PubMed]