The Dynamics of Modified Leslie-Gower Predator-Prey Model Under the Influence of Nonlinear Harvesting and Fear Effect

A modified Leslie-Gower predator-prey model with fear effect and nonlinear harvesting is developed and investigated in this study. The predator is supposed to feed on the prey using Holling type-II functional response. The goal is to see how fear of predation and presence of harvesting affect the model's dynamics. The system's positivity and boundlessness are demonstrated. All conceivable equilibria's existence and stability requirements are established. All sorts of local bifurcation occurrence conditions are presented. Extensive numerical simulations of the proposed model are shown in form of Phase portraits and direction fields. That is to guarantee the correctness of the theoretical results of the dynamic behavior of the system and to confirm the existence of various forms of bifurcations. The fear rate is observed to have a stabilizing effect up to a threshold value, after which it leads to prey extinction. The harvesting coefficients, on the other hand, serve as control parameters that, when exceeded, trigger the system to extinction.


1.
Introduction One of the most popular subjects in biomathematics is population dynamics. The study of the evolution of diverse populations has always been of special interest, beginning with populations of a single species and progressing to more realistic models in which several species exist and interact in the same ecosystem. Models that explore competitive interactions, symbiosis, commensalism, or predator-prey dynamics are some of them. The predator-prey model has been extensively studied by mathematical and biological researchers since its introduction made by Alfred J. Lotka in 1925 and Vito Volterra [1] in 1926. They described the interaction between two species combined with the predator-prey relationship, where they define the problem with a system of polynomial differential equations of degree two. The importance of this problem lies in understanding the dynamics between two species (a predator and prey) that live together in the same environment and looking for suitable conditions that allow both species to survive in equilibria. Later, applications of these systems began to increase. New applications on population dynamics had been developed, and these systems have also been utilized to represent a variety of other natural phenomena. Besides the basic relationship given by the Lotka-Volterra model; many factors may affect species growth. So, this model was developed by many researchers taking into consideration various environmental factors that affect the existence and stability of this system, such as prey refuge [2][3][4], disease [5,6], delay [7], harvesting [7][8][9], Allee effect [4,10,11], age structure [12], sex structure and sexual favoritism [13], seasonal variation [14], and many other factors. The functional response is an essential part of the predator-prey model, which describes the change in prey number killed per individual predator per unit of time as a consequence of changes in prey density. The most commonly used functional response in the existing literature is a function of prey's density only (Holling I-III) [2,4,15,16], in which interfering among predators is not utilized whereas this will be common when predators contest for food. To address this important factor, functional responses (ratio-dependent [17], Beddington-DeAngelis [18], and Crowley-Martin [7] have been developed which do not rely just on the density of the prey but rather on the density of both prey and predator. While many predator-prey models considered a logistic growth of predators, Leslie and Gower [19] assumed that the predator grows logistically, where its carrying capacity is proportional to the density of prey ℎ 1 − where and are the populations of prey and predator respectively. The term is called Leslie-Gower term. On the other hand, predators can devour other populations when food is scarce, but their growth will be limited since their primary prey is scarce. To consider this issue, Aziz-Alaoui and Okiye [20] suggested a modified Leslie-Gower model by introducing a constant in the denominator of Leslie-Gower term that measures environmental protection for the predator to avoid singularities when = 0. Since then, many researchers have examined the modified Leslie-Gower models with a variety of functional responses [2,[21][22][23], harvesting [7,22] Allee effect [24], etc. Moreover, from a financial income point of view, it is significant to consider the harvesting of species in predator-prey models. In the literature, several types of harvesting strategies have been utilized. Some of them used constant harvesting, ℎ( ) = ℎ [25], proportional harvesting ℎ( ) = [26] where is the population that presents the harvesting (prey or predator), ageselective harvesting [27], while others considered nonlinear harvesting [7]. Nonlinear harvesting is more relevant than other approaches from both a financial and biological point of view [23,28]. Many researchers consider Holling type II harvesting ℎ( ) = . For 261 example, Gupta et al worked with a model with Holling type II harvesting in prey [23] and Holling type II harvesting in predator in [29]. Another factor to consider is that in some environments, prey may be afraid of predators and respond appropriately, making predator hunting more difficult. Due to fear of predation risk, the prey population can change its feeding area to a safer place and sacrifice their highest intake rate areas, increase their vigilance, regulate their strategies for reproduction, etc. In recent years, many experts began to study the predator-prey model with fear effect; see [3,6] The dynamics and bifurcations of a modified Leslie-Gower predator-prey model with Holling-II functional response and nonlinear harvesting in both the prey and predator communities are investigated in this paper, as well as the influence of the fear factor.

Mathematical Model Formulation
The study considers a predator-prey problem, with ( ) and ( ), respectively, representing prey and predator population densities at time T. Resource-consumer, plant-herbivore, parasitehost, tumor cells (virus)-immune system, susceptible-infectious interactions, and so on are examples. In the proposed model, prey population ( ) is considered logistically growing in absence of predator ( ) with a birth rate and level of fear induced by predator population such that: (1) The morality density is represented by the term , where is the natural death rate of prey. Also, the term is added to consider competition between prey community members, where is the intraspecific competition. Moreover, the interaction between prey and predator is assumed to follow Holling-II functional response. According to these considerations, the change in the density of prey takes the following form in the presence of the predator: where represents the maximum attack rate, represents the half-saturation constant of predation and the parameter ∈ (0,1) represents the availability constant rate of prey for predation due to the assumption of the existence of a (1 − ) constant rate of prey's refuge in the environment. The density of predator population is assumed to follow the modified Leslie-Gower predation as follows: where the parameter represents the intrinsic growth rates of the predator, is the maximum value which per capita reduction rate of can attend, and is the carrying capacity of the predator in the absence of the prey. In the proposed model, prey and predator are assumed to follow nonlinear harvesting with the harvesting function of Holling-II. The harvesting of prey is presented by the term , while that of predator is represented by the term , where , ( = 1,2) are the catchability coefficient of prey and predator respectively, is the effort made to harvest individuals and , ( = 1,2,3,4): suitable constants. It's worth noting that the effort represented by in both equations is considered to be the same for both species, making this model more appropriate for aquatic environments. Combining all the above assumptions give the following set of dynamical differential equations: where all of the parameters are assumed to be positive and described as above. Note that, using the scaling variables = , = , and = in the system (4) reduces the number of parameters from 17 to 12 parameters and the system (4) takes the following dimensionless form: with the initial conditions: and the dimensionless parameters are given by: Note that, since the right-hand side of the interaction functions of the system (5) are continuous and have continuous partial derivatives, then system (5) has a unique solution that belongs to the positive quadrant ℝ .

3.
Positivity and Boundedness Theorem (1) and theorem (2) below prove that the model formulation is ecologically relevant by showing that solutions of system (5) together with the initial condition (6) are positive and uniformly bounded. Theorem 1: All solutions of system (5) with initial conditions (6) remain positive forever. Proof: The proof is direct and hence it is omitted. Theorem 2: All the solutions of system (5) with initial conditions (6) are uniformly bounded. Proof: From the system (5), and this shows that the solution of the system ( ) ≤ 1 − ≔ as → ∞, by lemma 2 in (23). Clearly, due to the survival condition of the prey in the absence of a predator, we have always that 1 − > 0. Now substituting the maximum value in the second equation of system (5) gives that: Then, by solving the above differential inequality, it is observed that: which proves the boundedness of all the solutions. □ 4.
Existence and the Stability of Equilibria The presence of equilibrium points of the dimensionless system, as well as a qualitative analysis of their stability, are investigated in this section. The number of equilibrium points of system (5) depends on the parameter values. For example, Figure 1 shows that for the set of parameter values given in Table 1, the system has one trivial, one predator-free, one prey-free, and one interior equilibrium point. While, for set of parameter 263 values given in Table 2, the system has one trivial, one predator-free, two prey free and two interior equilibrium points.
(a) (b) Figure 1-The number of equilibrium points of the system (5) (a) there are 4 equilibrium points for set #1 of parameter values given in Table 1 (b) there are 6 equilibrium points for set #2 of parameter values given in Table 2.

4.1
Trivial Equilibrium Point It is clear that the trivial equilibrium point (0,0) at the origin is always exists.

4.2
Predator Free Equilibrium Points The boundary equilibria on −axis are calculated by solving the following quadratic equation: The roots of Eq. (7) depend on the parameters , , , and , so according to Descartes's rule there are the following cases: Case 1. There is no equilibrium point with ≥ 1, since the prey survives if the natural mortality rate is lower than the birth rate.  7) with < and are given by Prey Free Equilibrium Points The boundary equilibria on −axis are calculated by solving the following quadratic equation: Interior Equilibrium Points The positive interior points are found by solving system (5) for > 0 and > 0. It is obtained that: It is clear that each one of them has a positive value. However, * is the positive root of the polynomial (10) below:

5.
Stability Analysis of Equilibria In this section, the nonlinear system (5) is linearized around each equilibrium point using the Jacobian matrix to investigate the local stability of various equilibrium points. The Jacobian matrix of system (5) about an arbitrary point ( , ) is determined by: Recall that, if all eigenvalues of the Jacobian matrix at an equilibrium point have negative real parts then this point is locally asymptotically stable. Accordingly, the following theorems present the local stability conditions for each of the above equilibria. If the condition (i) holds, then < 0 and it is clear that < 0 ⇒ | ( * , * )| < 0 ⇒ * is a saddle point. If the condition (ii) holds, then > 0 ⇒ | ( * , * )| > 0 ⇒ * is unstable node when ( * ) > 0. Similarly, * is stable node when ( * ) < 0. □

Bifurcation Analysis
This section is dedicated to study some potential bifurcation scenarios at the stable equilibrium points of the system (5)  Clearly, has a zero eigenvalue with another negative eigenvalue, and the corresponding eigenvector for the zero eigenvalue can be written as: While the eigenvector corresponding to the zero eigenvalue of is determined as: Differentiating with respect to gives: Therefore, straightforward computation shows that: 267 ( , * ) = 0. Consequently, by Sotomayor's theorem, system (5) has no saddle-node bifurcation near and = * . Moreover, direct computation gives that: ( , * ) = 1 ≠ 0. Also, using the form of given by equation (12a) Hence, a pitchfork bifurcation takes place, and the proof is complete. □ Theorem 8: Assume that < ( + ) near any of the predator-free equilibrium points ( , 0), ( = 1,2), then if the parameter passes through the value * = , then the system (5) at this equilibrium point has i.No saddle-node bifurcation.
iii.A pitchfork bifurcation otherwise. Proof: At ( , 0), ( = 1,2), the Jacobian matrix of system (5) with = * becomes: Clearly, has a zero eigenvalue with another negative eigenvalue, and the corresponding eigenvector for the zero eigenvalue can be written as: .
Clearly < 0, due to given condition. While the eigenvector corresponding to the zero eigenvalue of is determined as: And, ( , * ) = 0. Consequently, by Sotomayor's theorem, system (5) has no saddle-node bifurcation near and = * . Moreover, direct computation gives that: ( , * ) = 1 ≠ 0. Also, using the form of given by equation (12a), and the eigenvectors with gives that: Clearly has a zero eigenvalue with another negative eigenvalue, and the corresponding eigenvector for the zero eigenvalue can be written as: ) . Clearly > 0, due to given condition. While the eigenvector corresponding to the zero eigenvalue of is determined as: And, ( , * ) = 0. Consequently, by Sotomayor's theorem, system (5) has no saddle-node bifurcation near and = * . Moreover, direct computation gives that:  Straightforward computation shows that | * | = * − * = 0, and hence * has a zero eigenvalue, and the corresponding eigenvector for this eigenvalue can be written as: While the eigenvector corresponding to the zero eigenvalue of * is determined as: And, ( * , * * ) = * 1 − * * ≠ 0. Consequently, by Sotomayor's theorem, system (5) satisfies the first condition of saddle-node bifurcation near * and * * . Moreover, direct computation with the use of conditions (17) and (18) gives that: , represent the elements of ( * , * * )( , ). Hence, a saddle node bifurcation takes place, and the proof is complete. □ 7.
Numerical Analysis The overall dynamics of system (5) are numerically investigated in this section for various sets of initial values and parameter values. This study aims to demonstrate the impacts of changing parameter values on the system's dynamical behavior and to validate the theoretical results achieved. The numerical simulations are carried out using a combination of parameter values given in Table (1) and Table (  As previously mentioned in section 4, there is a various number of equilibrium points for different values of parameters. Figure 1 shows two examples of that. Now, using the parameters given in Table 1 with various sets of initial points, system (5) is solved numerically, and then the trajectories that have been obtained are drawn in form of direction field and phase portrait as shown in Figure 2. It is clear from this figure that there are 4 equilibrium points: one trivial, one predator-free, one prey-free, and one interior equilibrium point. The interior equilibrium point is stable (spiral sink), while the others are unstable.   Table   (1), then instead of four equilibrium points, system (5) will have just two: the trivial and the predator-free equilibrium points. Both of them are nonhyperbolic equilibrium points. This agrees with the results of both Theorem 7 and Theorem 8.  Table 1 then instead of four equilibrium points, system (5) will have three equilibrium points: the trivial, one predator-free, and one prey-free equilibrium point. Moreover, the prey-free equilibrium point will change its behavior from saddle-node to nodal sink. This agrees with the result of Theorem 9.
(a) (b) Figure 5-The dynamical behavior of system (5) with parameter values in Table- As mentioned previously, system (5) has a rich dynamic behavior, for example, Figure (6) shows the behavior of the system under influence of set #2 of parameter values which are given in Table 2. In this case, there are six equilibrium points: (i) one trivial saddle equilibrium point (ii) one predator-free nodal sink (iii) two prey-free equilibrium points, one of them is saddle point and the other is a nodal sink (iv) two interior equilibrium points, one spiral source and the other is a saddle point. The dynamic behavior near the 2 interior equilibrium points. Figure 7 shows that for set #1, as the value of the level of fear induced by predator population ( ) increases, the unique interior equilibrium point converges to the prey-free equilibrium point and preserves its behavior as a spiral sink. While Figure 8 shows that for set #2 with an increase in the level of fear induced by predator population ( ), the two interior equilibrium 274 points converge to each other and then unify and disappear. When these two points exist, they preserve their behavior, the right one is always saddle-node while the left one is always a spiral source.
Moreover, to explore the effect of harvesting, Figure 9 shows that for set #1, as the value of the catchability coefficient of prey ( ) increases, the unique interior equilibrium point converges to the prey-free equilibrium point and preserves its behavior as a spiral sink. While Figure 10 shows that for set #2, the two interior equilibrium points converge to each other and then unify and disappear. When these two points exist, they preserve their behavior, the right one is always saddle-node while the left one is always a spiral source. Also, Figure 11 shows that for set #1, as the value of the catchability coefficient of predator ( ) increases, the unique interior equilibrium point converges to predator-free equilibrium point and preserves its behavior as a spiral sink. At the same time, for set #2, the two interior equilibrium points follow a similar scenario to that appear in Figures 8 and 10, where the two points combine with each other and then disappear. On the other hand, decreasing the value of this parameter toward zero makes the right interior equilibrium point goes toward the predatorfree equilibrium point, while the left interior equilibrium point goes toward the prey-free equilibrium point and then toward the trivial equilibrium point as shown in Figure 12.

Discussion and Conclusions
In this study, a modified Leslie-Gower predator-prey model is proposed and discussed. Fear of predation and harvesting are both factored into the model's design and studied. According to Holling type-II functional response, which is also used to characterize the harvesting process, the predator consumes food. The solution's positivity and boundlessness have been demonstrated. All of the system's probable equilibrium positions are identified. These equilibrium points' local stability is explored, and their requirements are provided. The system shows it has several boundaries and interior equilibrium points, which complicates the study. The effect of changing the parameters on the model's dynamics is investigated by looking at the possibilities of local bifurcation types such as saddle-node, transcritical, and pitchfork bifurcation. Two sets of hypothetical data parameters are used to analyze the system's global 280 dynamics numerically. The existing equilibrium locations and the solution behavior surrounding them are depicted using phase portraits and direction fields. The system has a single globally asymptotically stable interior equilibrium point with two axial saddle-node equilibrium points and an unstable trivial equilibrium point for the data in set #1. The system has two unstable interior equilibrium points (one unstable spiral and the second is a saddle point) with three axial equilibrium points (two of which are locally asymptotically stable with their own basin of attraction for each of them and the third point is a saddle point) and a saddlenode trivial equilibrium point for data set #2. This ensures that the system's dynamics are rich, and the number of equilibrium points is determined by the data parameters. The impact of fear rate on the system's dynamics is also explored. For set # 1 data, increasing fear rate produces a gradual decrease in prey population until the system reaches prey extinction, and the interior equilibrium point coincides with the prey-free equilibrium point, which becomes globally asymptotically stable. In the set #2 data, however, raising the fear rate causes the two inner equilibrium points in the positive quadrant to be combined, resulting in a saddle-node point. According to the preceding discovery, fear rate acts as a bifurcation parameter, causing extinction in prey species when its value rises above a certain threshold. When the value of the catchability coefficient of the prey increases, similar observations are made as with the fear rate. Furthermore, when using set # 1 data, as the value of the predator's catchability coefficient increases, the system behaves similarly to that obtained when increasing the value of the prey's catchability coefficient, except that the interior equilibrium point combines with the predator-free equilibrium point, which then becomes globally asymptotically stable. The system behaves similarly when using set # 2 data. Finally, given set # 2 data, decreasing the predator species' catchability coefficient causes each of the interior equilibrium points to be combined with the nearest axial equilibrium point, and the system's solution approaches asymptotically to the prey-free equilibrium point.