The Effects of Harvesting on the Dynamics of a Leslie–Gower Model

with x1(0)≥ 0 and x2(0)≥ 0, where x1(τ) and x2(τ) are the prey and predator population densities, respectively, r, s, a1, a2, n, e1, e2 > 0, and τ is the time. Note that (a2x2/(n + x1)) is Leslie–Gower term in which the carrying capacity of the predator’s environment is a linear function of the prey size (x1/a2) + (n/a2). (a1x1x2/(n + x1)) is the number of prey consumed by the predator in unit time which shows that when the number of the prey x1 is severe scarcity, and the predators can switch over to other populations as food. Constants r and s are the intrinsic growth rate of the prey and predator, respectively, and e1 and e2 denote the harvesting efforts for the prey and predator, respectively. Since the first prey-predator dynamical models which is the Lotka–Volterra model was built in the 1920s by Mathematician Lotka and Volterra, more and more researchers are interested in such issues, and they start from different angles to think the problem and many important results have been obtained [1–10]. In particular, in 2003, Aziz-Alaoui and Daher Okiye [11] considered the following Leslie–Gower predator-prey model:


Introduction
In this paper, we consider Leslie-Gower predator-prey model with harvesting effect, with x 1 (0) ≥ 0 and x 2 (0) ≥ 0, where x 1 (τ) and x 2 (τ) are the prey and predator population densities, respectively, r, s, a 1 , a 2 , n, e 1 , e 2 > 0, and τ is the time. Note that (a 2 x 2 /(n + x 1 )) is Leslie-Gower term in which the carrying capacity of the predator's environment is a linear function of the prey size (x 1 /a 2 ) + (n/a 2 ). (a 1 x 1 x 2 /(n + x 1 )) is the number of prey consumed by the predator in unit time which shows that when the number of the prey x 1 is severe scarcity, and the predators can switch over to other populations as food. Constants r and s are the intrinsic growth rate of the prey and predator, respectively, and e 1 and e 2 denote the harvesting efforts for the prey and predator, respectively.
Since the first prey-predator dynamical models which is the Lotka-Volterra model was built in the 1920s by Mathematician Lotka and Volterra, more and more researchers are interested in such issues, and they start from different angles to think the problem and many important results have been obtained [1][2][3][4][5][6][7][8][9][10]. In particular, in 2003, Aziz-Alaoui and Daher Okiye [11] considered the following Leslie-Gower predator-prey model: where x 1 is the numbers of prey and x 2 is the numbers of predators. Existence and stability of the fixed points were studied by using the Lyapunov function. In 2006, Lin and Ho [12] discussed the local and global stability for system (2) by using Poincaré-Bendixson theorem and Dulac's criterion. Harvesting is an effective way for humans to control the size of predators and prey so that the population has continued to develop healthily and produced good economic benefits [13][14][15][16]. Academically, researchers often only consider the harvesting of prey in order to control the size of the population. In 2010, Zhu and Lan [17] investigated the Leslie-Gower predator-prey systems: In 2013, Gupta and Chandra [18] discussed the following Leslie-Gower predator-prey model with harvesting on the prey and the environment providing the same protection to both the predator and prey: For ecological balance and healthy economic development, for fisheries, wildlife resources, etc., we not only need to consider the harvesting of prey, but also the predator. erefore, in this paper, we study Leslie-Gower predatorprey model (1) with harvesting on the prey and predator.
where α is a positive constant, t ≥ 0, and x(0) > 0, we have Lemma 2 (see [20,21]). Consider system _ X � f(X, α) and suppose that f(X 0 , α 0 ) � 0, n × n Jacobian matrix (J ≡ Df(X 0 , α 0 )) has a simple eigenvalue s � 0 with eigenvector V, and the transpose of the Jacobian matrix J T has an eigenvector W to the eigenvalue s � 0. en, the system _ X � f(X, α) experiences a transcritical bifurcation at the equilibrium point X 0 as the control parameter α passes through the bifurcation value α � α 0 if the following conditions are satisfied: e rest of this paper is organized as follows. In Section 2, we study boundary of solutions. In Section 3, we discuss existence of equilibria points. In Section 4, we discuss stability of the equilibrium points.

Boundedness of Solutions
In this section, we prove that every solution of system (5) is positive and uniformly bounded with initial conditions Theorem 1. Consider system (5). For any given initial conditions (x 0 , y 0 ) ∈ R 2 + , the solution (x(t), y(t)) of system (5) exists and is unique and positive and ultimate bounded.
Obviously, function f(x(t), y(t)) is continuous differentiable on (x, y) ∈ R 2 + , so for any given the initial conditions (x 0 , y 0 ) ∈ R 2 + , the solution (x(t), y(t)) of the system (5) exists and is unique. Furthermore, x− axis and y− axis are the solutions of the system (5); by the uniqueness of the solution, the solutions (x(t), y(t)) of the system (5) with the initial value x 0 > 0, y 0 > 0 cannot cross with x− axis and y− axis.
Next, we show the solutions (x(t), y(t)) of system (5) with the initial value x 0 > 0 and y 0 > 0 which is ultimate bounded.
From system (5), we have Combining Lemma 1, we have So, the following inequality is established: (10) By Lemma 1, we have where

Existence of Equilibria
In order to find the equilibrium points of system (5), we let f(x(t), y(t)) � 0, i.e., It is clear that equation (12) has a trivial solution E 0 ≔ (0, 0). Furthermore, by calculation, we find other solutions of equation (12): (13) is shown in Figure 1.
erefore, we have the following results.where

Theorem 2. Consider system (5) admits x− axial only and y− axial equilibria under following conditions.
(i) e x− axial equilibrium, is a boundary equilibrium of system (5) if and only if (ii) e y− axial equilibrium, is a boundary equilibrium of system (5) if and only if Theorem 3. Consider system (5) admits a unique positive equilibrium, if only if Remark 2. Existence regions for equilibrium points of system (12) is shown in Figure 2. Equilibrium points E 1 and E 2 exist, but positive equilibrium point E 3 does not in region II and equilibrium points E 1 , E 2 , and E 3 coexist in region

Stability of of Equilibria
Proof. Firstly, we show the equilibria E 0 is a unstable equilibrium. e Jacobian matrix about E 0 is given by Secondly, we show the equilibria E 1 is a unstable equilibrium. e Jacobian matrix about E 1 is given by Discrete Dynamics in Nature and Society 3 Proof.
e Jacobian matrix about E 2 is given by erefore, E 2 is locally asymptotically stable.
) hold, then the positive equilibrium E 3 is locally asymptotically stable. Furthermore, assume that there is α < ρβ; then, the positive equilibrium E 3 is globally asymptotically stable.

Proof.
e Jacobian matrix about E 3 is given by with trace A direct calculation gives we obtain erefore, the positive equilibrium E 3 is unstable.

Local Bifurcation
hold, then system (5) undergoes a Hopf bifurcation with respect to bifurcation parameter m around the equilibrium point E 3 � (x ∞ , y ∞ ). Furthermore, the direction of the Hopf bifurcation is subcritical and the bifurcation periodic solutions are orbitally asymptotically stable if
e direction of the Hopf bifurcation is supercritical and the bifurcation periodic solutions are unstable if Proof. From (39), we have So, m H > 0. e Jacobian matrix of system (5) evaluated at the point E 3 is given by e trace z � Tr(J E 3 ) and the determinant D � Det(J E 3 ) of Jacobian matrix J E 3 are given by In addition, we have So, (zz/zm)| m H < 0. erefore, this guarantees the existence of Hopf bifurcation around E 3 .
We translate the equilibrium E 3 to the origin by the translation x � x − x ∞ and y � y − y ∞ . For the sake of convenience, we still denote x and y by x and y, respectively. So, the system (5) becomes Rewrite system (47) to where Denote the eigenvalues of J E 3 by φ + iω with φ � (z/2) and ω � ( Discrete Dynamics in Nature and Society where G � ((ω(m + x ∞ ))/αx ∞ ) and Obviously, and (52) By the transformation, where with In order to determine the stability of the periodic solution, we need to calculate the sign of the coefficient b(m H ), which is given by where all partial derivatives are evaluated at the bifurcation point (0, 0, m H ). ), ω 0 , G 0 , and N 0 , we have 8 Discrete Dynamics in Nature and Society where m � m H and Discrete Dynamics in Nature and Society 9 Combining (41), we have b(m H ) < 0.
erefore, according to Poincare-Andronow's Hopf bifurcation theory, we have the direction of the Hopf bifurcation is subcritical and the bifurcation periodic solutions are orbitally asymptotically stable.
In addition, combining (42), we have b(m H ) > 0. erefore, according to Poincare-Andronow's Hopf bifurcation theory, we have the direction of the Hopf bifurcation is supercritical and the bifurcation periodic solutions are unstable.

Numerical Illustrations
In this section, we perform numerical simulations about system (5). Figure 3 shows that E 0 is an unstable node point, E 1 is a saddle point, E 3 does not exist, and E 2 is asymptotically stable and every orbit tends to it. Figure 4 shows that E 0 is an unstable node point, E 1 is a saddle point, E 2 is unstable, E 3 is unstable, and there is a limit cycle around E 3 to which every orbit tends. Figure 5 shows that E 0 is an unstable node point, E 1 is unstable and E 2 is also unstable, but E 3 is asymptotically stable and every orbit approaches this equilibrium.

Data Availability
No data were used to support this study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.