The Dynamics of a Tritrophic Leslie-Gower Food-Web System with the Effect of Fear

The avoidance strategy of prey to predation and the predation strategy for predators are important topics in evolutionary biology. Both prey and predators adjust their behaviors in order to obtain the maximal benefits and to raise their biomass for each. Therefore, this paper is aimed at studying the impact of prey’s fear and group defense against predation on the dynamics of the food-web model. Consequently, in this paper, a mathematical model that describes a tritrophic LeslieGower food-web system is formulated. Sokol-Howell type of function response is adapted to describe the predation process due to the prey’s group defensive capability. The effects of fear due to the predation process are considered in the first two levels. It is assumed that the generalist predator grows logistically using the Leslie-Gower type of growth function. All the solution properties of the model are studied. Local dynamics behaviors are investigated. The basin of attraction for each equilibrium is determined using the Lyapunov function. The conditions of persistence of the model are specified. The study of local bifurcation in the model is done. Numerical simulations are implemented to show the obtained results. It is watched that the system is wealthy in its dynamics including chaos. The fear factor works as a stabilizing factor in the system up to a specific level; otherwise, it leads to the extinction of the predator. However, increasing the prey’s group defense leads to extinction in predator species.


Introduction
Mathematical modeling is used to understand the interaction of organisms with their surrounding environment, as well as species evolution that is well-known biological evolution. Nonlinear differential equation systems have always played an important role in simulating such mathematical models. Therefore, in the last few decades, there are increasing interests in the study of the dynamics of such systems. After the pioneering prey-predator model by Lotka-Volterra, the prey-predator models have been essential in mathematical ecology and have been extensively investigated. The Lotka-Volterra model was then adjusted by compiling a prey's logistic growth and a Holling-type II functional response to describe the predation [1].
It is completely recognized that a prey-predator is fundamental interaction for drawing the dynamics of food webs in the real-world environment, predation risks can negatively affect prey biomass and growth efficiency, and thus, predators affect the structure of natural communities. Prey can die of fear, and it can also reduce fertility, as a result of which fear can be just as important as outright killing by predators in affecting prey numbers [2][3][4]. Many prey-predator models have been suggested and investigated thoroughly in which either the predator kills the prey or the presence of the predator affects the behavior of the prey population due to the fear of the predation process [5]. Upadhyay and Mishra [6] modeled a prey-predator system taken into account the fear effect in prey growth. They watched that a relative rise in the scale of fear may minimize the size of the total population. Panday et al. [7] investigated the action of fear in a three-species food chain model, where the evolution rate of lower predator is minimized due to the fear's action from a higher predator, and the evolution rate of prey is reduced due to the fear's action from the lower predator. Das and Samanta [8] studied the fear's impact from the predator on prey in case of providing additional food for the predator. Recently, Kumar and Kumari [9] discussed the action of fear on the dynamics of the food chain model having three species. They obtained that for low fear's scale, the system stays chaotic while relatively high fear's scale leads to stability. Moreover, the rise of fear's scale further causes population extinction.
On the other hand, Aziz-Alaoui [10] investigated a dynamic system consisting of a Leslie-Gower food chain of three species. He displayed that for a realistic parameter set, chaotic dynamics could be obtained. Zhang et al. [11] studied a harvested Leslie-Gower prey-predator model. They got with the help of Pontryagin's maximal principle the optimal harvesting policy for the model. Cai et al. [12] investigated a Leslie-Gower prey-predator model with the diffusion and Allee effect on prey. They gained that the Allee effect basically rises the model spatiotemporal complexity. Later on, Abid et al. [13] were interested in the stability and optimal harvesting of a modified Leslie-Gowe preypredator model with Holling-type II functional response. Recently, Singh and Bhadauria [14] studied the dynamics of the prey-predator model with weak Allee effect II and modified Leslie-Gower. They explained that Allee effect II greatly impacts the model and can raise the risk of extinction.
Keeping the above in mind, in this paper, a tritrophic Leslie-Gower food-web system is constructed mathematically. Sokol-Howell type of function response [15] is used for describing the predation process due to the prey's group's defensive capability. The effects of fear due to the predation process are included and studied. The following section is written out as follows: the construction of the model is done. In Section 3, existence and local stability analysis are studied. In Section 4, persistence is discussed. However, in Section 5, the region of attraction for each equilibrium point of system (4) is determined by using the direct method of Lyapunov. Local bifurcation analysis is discussed in Section 6. In Section 7, the numerical simulation of the model is studied. Finally, Section 8 gives the conclusion.

The Mathematical Construction of the Model
In this section, a real-world tritrophic food-web system is constructed mathematically. It is supposed that the logistic law is applied for both the prey and the generalist-predator so that the environmental carrying capacity of the generalist-predator is not a constant but a function of the available food quantity so that there are upper limits to the rates of increase of both prey and generalist-predator. Therefore, the Leslie-Gower model is used in the construction of the proposed model, which assumes that both the prey and the generalist-predator grow according to the logistic law. Moreover, both the prey and specialist-predator have a defensive property against the predation when they are available in abundance. Therefore, a simplified Hollingtype IV functional response or Sokol-Howell type of func-tional response is used in the proposed model. This is because the per capita predation rate in the Sokol-Howell functional response rises with prey density to a maximum at a critical prey density beyond which it decreases.
Let XðTÞ represents the prey density at the time T, YðTÞ is the specialist-predator density at the time T, while ZðTÞ is the generalist-predator density at time T. Accordingly, the dynamics of such a food-web system can be described as follows: where gðXÞ = rð1 − bXÞ is the logistic growth rate of the prey species, while PðXÞ = a 0 X/ð1 + m 0 X 2 Þ and q 1 ðX, YÞ = a 1 X/ð1 are the Sokol-Howell functional responses for the specialistpredator and generalist-predator, respectively. The parameters of system (1) are supposed to be positive, and they are described as the following: the parameter r > 0 represents the intrinsic growth rate, ð1/bÞ > 0 is the carrying capacity of the environment for the prey, the positive parameters a 0 , a 1 , and a 2 are the maximum attack rates, the positive parameters m 0 , m 1 , and m 2 are the preference rates of the predator to their prey species, a 3 > 0 represents the conversion rate from individuals of the prey to specialist-predator, and c 1 > 0 is the intrinsic growth rate of a generalist-predator, while the positive parameters c 3 and c 4 are the respective food preference of generalistpredator from the prey and specialist-predator, respectively. In addition to the above, the square term c 1 Z 2 in the third equation that represents the modified Leslie-Gower term signifies the fact that mating frequency is directly proportional to the number of males as well as to that of females, see [16,17]. Finally, the parameter c 2 > 0 is a constant added in the denominator to normalize the residual reduction in the generalistpredator due to a heavy lack of the preferable food X and Y, while d > 0 is the natural death rate of the specialist-predator. Now, to add the factor of the fear for each of the prey and the specialist-predator from the predation by an upper-level predator, the growth of each of them should be reduced by a rate proportional with the abundance of their predators. Therefore, the dynamics of the food-web system given by system (1), in case of the existence of fear, can be represented as where the positive parameters n 1 , n 2 , and n 3 represent the fear 2 Journal of Applied Mathematics coefficients of the prey and specialist-predator from their direct predator. Clearly, system (2) has 16 parameters, which leads to analysis difficulty. So, in order to simplify the analysis of system (2) and minimize the number of parameters, the following nondimensional variables and parameters are utilized.
Simple computation displays that the two isoclines given by (16a) and (16b) have a unique positive intersection point, denoted by ðy * , z * Þ, provided that the following sufficient conditions hold: Here, y 1 and y 2 represent the positive root of the following two polynomials, respectively: ς 1 y 4 + ς 2 y 3 + ς 3 y 2 + ς 4 y + ς 5 = 0, where 3 12 , Now, the local stability of the above equilibrium points is studied by computing the Jacobian matrix (JM) and then determining their eigenvalues. It is simply to verify that the JM of system (4) at the point ðx, y, zÞ can be represented as Journal of Applied Mathematics with Consequently, for the TP, the JM becomes Then, the eigenvalues of J E 0 are given by λ 10 = 1 > 0, λ 20 = −α 9 < 0, and λ 30 = 0. Since there is a positive equilibrium point, then E 0 is an unstable nonhyperbolic point.
For the SPEP, the JM can be determined as where Clearly, the roots (eigenvalues) of equation (30) can be represented as It is easy to verify that all the eigenvalues of J E 3 are negative, and hence, the SPFP is locally asymptotically stable (L.A.S), presuming that the next conditions hold: The JM of system (4) at the CP can be represented as where with Consequently, the characteristic equation of J E 4 can be represented as where A 1 = −ðc 11 + c 22 Þ, with A reminder that according to the Routh-Hurwitz criterion, the characteristic equation (35) has negative real part eigenvalues, and then, the CP becomes L.A.S if and only if A 1 > 0, A 3 > 0, and Δ>0. Accordingly, the following theorem can be proved simply.

Theorem 2.
The CP is L.A.S if the below sufficient conditions hold.
Journal of Applied Mathematics

Persistence
In this part, the persistence of system (4) is investigated; it is fully recognized that the system is said to persist if and only if none of their species extinct; this means that system (4) persists if the trajectory of the system that initiates at a positive initial point does not have omega limit set on the boundary planes of its domain. System (4) has two subsystems lying in xy-plane and xz -plane, respectively, which can be represented as follows: and It can be simply confirmed that the above subsystems (39) and (40) have positive equilibrium points that coincide with those of system (4) in the interior of boundary planes xy-plane and xz-plane, respectively. Now, to discover the possibility of the existence of periodic dynamics in the boundary planes, the Dulac function approach is used.
Consider the following function B 1 ðx, yÞ = 1/xy, clearly this function B 1 ðx, yÞ > 0 and C 1 function in the interior of Therefore, Δðx, yÞ does not alter sign and not vanish under the following condition: Hence, if condition (42) holds, there are no periodic dynamics in the interior of a positive quadrant of xy-plane for subsystem (39).
Similarly, with the help of Dulac function B 2 ðx, zÞ = 1/x z 2 , which is B 2 ðx, zÞ > 0, and C 1 function in the interior of ℝ 2 + of xz-plane, it is observed that there are no periodic dynamics in the interior of a positive quadrant of xz-plane for subsystem (40) provided that the following condition holds: Theorem 3. Suppose that the boundary planes have no periodic dynamics; then, system (4) is uniformly persistent as long as that the next conditions hold: Proof. Define a function ϑðx, y, zÞ = x p 1 y p 2 z p 3 , where p 1 , p 2 , and p 3 are positive constants, and ϑðx, y, zÞ > 0 for all ðx, y, zÞ ∈ Int ℝ 3 + with ϑðx, y, zÞ = 0 if any one of x, y, or z approaches zero. Therefore, direct computation gives where the functions f i ; i = 1, 2, 3 are given in system (4). Now, according to the average Lyapunov method, the proof is satisfied provided that Ωðx, y, zÞ > 0 for all boundary equilibrium points. Therefore, where Γ i ; i = 1, 2, 3, 4, 5 are given in equation (20). Clearly, we have that Obviously, by choosing the arbitrary positive value of p 1 sufficiently larger than that of p 2 , it is obtained that ΩðE 0 Þ > 0.
Again, by choosing the arbitrary positive value of p 1 sufficiently larger than that of p 2 , it is obtained that Ωð E 1z Þ > 0 provided that condition (44c) holds.

Region of Attraction
In this section, the region of attraction for each equilibrium point of system (4) is determined using the direct method of Lyapunov. These regions of attraction for each of their equilibrium point are known as the domain or basin of attraction.
Theorem 4. The FAP is asymptotically stable (A.S) in the subregion of ℝ 3 + that satisfies the following sufficient conditions: where the symbol G is given in the proof.
Proof. Define L 1 ðx, y, zÞ = Ð x 1 ððu − 1Þ/uÞdu + y + z. Clearly, the function L 1 ðx, y, zÞ is a positive definite function that is L 1 ð1, 0, 0Þ = 0, while L 1 ðx, y, zÞ > 0, for all values in the region fðx, y, zÞ ∈ ℝ 3 + : x > 0, y ≥ 0, z ≥ 0 ; ðx, y, zÞ ≠ ð1, 0, 0Þg . Then, using some algebraic manipulation gives that Consequently, by using additional computation, the following is obtained: Now, according to the upper bound for x and y, so as t ⟶ ∞, we have x ≤ 1 and y ≤ ðα 6 /4α 9 Þ + α 6 . Therefore, it is obtained that Consequently, using conditions (52b) and (52c) with the bound of logistic term yields where G = 1/ðα 11 + α 12 + α 13 ðα 6 /4α 9 + α 6 ÞÞ − 1. Accordingly, condition (52d) guarantees that ðdL 1 /dtÞ < 0, which means it is a negative definite. Then, L 1 is a Lyapunov function, and hence, the FAP is A.S for any trajectory starting from a point that belongs to the region that satisfies the above conditions. ☐ Note that the subregion described in Theorem 4 represents the basin of attraction of the FAP. 8 Journal of Applied Mathematics Theorem 5. The SAP is A.S in the subregion of ℝ + 3 that satisfies the following sufficient conditions: x where the symbol N > 0 is given in the proof.
Proof. Define Clearly, the function L 5 ðx, y, zÞ is a positive definite function that is L 5 ðx * , y * , z * Þ = 0, while L 5 ðx, y, zÞ > 0, for all values in the region fðx, y, zÞ ∈ ℝ 3 + : x > 0, y > 0, x > 0 ; ðx, y, zÞ ≠ ðx * , y * , z * Þg. Thus, after some algebraic manipulation, it is obtained that where Journal of Applied Mathematics Consequently, using the above conditions gives Obviously, ðdL 5 /dtÞ < 0, which gives that L 5 is a Lyapunov function and hence the CP is A.S for any trajectory starting from a point that belongs to the region that satisfies the above conditions. ☐

Local Bifurcation
An investigation of the impact of varying the parameter values on the dynamics of system (4) around each equilibrium point is done using Sotomayor's theorem for local bifurcation.
Because the existence of a nonhyperbolic equilibrium point represents a necessary but not sufficient condition for bifurcation to occur, then the candidate bifurcation parameter is choosing so that the equilibrium point becomes a nonhyperbolic point. Rewrite system (4) in the following general form: Then, the second directional derivative of system (4) can become as where with Φ = ðv 1 , v 2 , v 3 Þ T be any nonzero vector, and the symbols Γ i ; i = 1, 2, 3, 4, 5 are given in equation (20). Accordingly, the local bifurcation near the equilibrium points is studied in the next theorems.
Proof. Direct computation shows that the JM of system (4) at ðE 3 , α * 9 Þ can be represented as where Γ i ; i = 1, 2, 3, 4, 5 are given in equation (28). Obviously, the matrix J 1 has two eigenvalues having negative real parts due to condition (32a) and the third one is zero, say λ 32 * = 0. Hence, E 3 is a nonhyperbolic point.
Proof. Direct computation shows that the JM of system (4) at ðE 4 , α * 12 Þ can be represented as where all the coefficients c ij ; i, j = 1, 2, 3 are given in equation (33). Clearly, the determinant of J 2 that is given by A 3 in equation (35) is zero at α * 12 . Therefore, the characteristic equation of J 2 has zero root, say λ * = 0, and then, the coexistence equilibrium point becomes a nonhyperbolic point.
Let Clearly, ρ 5 > 0, and ρ 6 < 0, due to conditions (38a)-(38c) with (38e). Now, according to Sotomayor's theorem, it is obtained that where Γ 5 is given in equation (20). Therefore, it is obtained that where Γ 5 * is given in equation (33). Now, by using equation (73), it is obtained that 13 Journal of Applied Mathematics According to Figure 3, increasing the value of α 3 causes decreasing in the predator's populations and increasing in prey population gradually; then, the solution of system (4) settled at the SPFP first and then at the FAP. However, increasing α 5 destabilizes system (4) and the solution persists at periodic dynamics. Moreover, it is observed that the effect of the parameter α 4 is similar to that of α 3 . Now, the effects of varying the parameters α 6 and α 9 on system (4) dynamical behavior are shown in Figures 5 and 6, respectively.
Obviously, Figures 5 and 6 show that the parameters α 6 and α 9 have opposite effects on the dynamic behavior of system (4) so that as increases α 6 (decreasing α 6 ), the population of the prey decreases while the populations of the 14 Journal of Applied Mathematics predators increase and hence system (4) loses their persistence and approaches to FAP as decreasing α 6 (increasing α 9 ). Furthermore, it is observed that varying the parameter α 10 has a quantitative effect on the dynamic behavior of system (4) and system (4) persists at the CP. The effects of varying the parameters α 11 and α 12 on system (4) dynamical behavior are explaining in Figures 7 and  8, respectively.
According to Figures 7 and 8, it is watched that increasing the parameter α 11 leads to losing the persistence of the system by approaching the SPFP first and then to SAP. However, decreasing (or increasing) the parameter α 12 leads to losing the persistence of the system and the solution approaches to GPFP (SPFP). It is observed that the parameter α 13 has similar effects on the dynamic behavior of system (4) as that for α 12 .
Keeping the above in mind, the second set of parameter values, which leads system (4) to chaotic dynamics, is given in the set of data (90): The dynamic behavior of system (4) is shown in Figures 9 and 10.
According to Figures 9 and 10, system (4) is rich in its dynamics including point attractors, periodic attractors, and chaos. It has many bifurcation points when varying its parameter values. It is well known that the most property that characterizes the chaotic attractor is the sensitivity to the varying initial points. So, Figure 11 shows clearly the sensitivity of the strange attractor given in Figure 9 to a small change in their initial point. Now, the impacts of varying some parameter values on the chaotic regions and the dynamic behavior of system (4) are studied using the bifurcation diagrams. It is known that the bifurcation diagram shows the birth and death of the attractors as a function of varying the parameter values. The bifurcation diagrams (BD) as a function of the parameters α 1 , α 7 , α 9 , α 12 , and α 13 are shown in  According to Figures 12 and 13, there are ranges around α 1 = 0:8 and α 7 = 0:8 in which the system is chaotic; however, it faces extinction in the specialist predator as these parameters increase and the trajectories approach periodic dynamics in xz-plane. Similar effects on the dynamics behavior of system (4) are observed for the parameters α 2 as those of α 1 and α 7 . Finally, the other BD are clearly shown the existence of chaotic regions; however, it is watched that   17 Journal of Applied Mathematics increasing these parameters leads to extinction in one predator species and the solution approaches to periodic dynamics in the plane.

Discussion and Conclusion
In this paper, a food-web model having specialist and generalist predators is proposed and studied. It is assumed that the prey has antipredator properties against predation and hence Sokol-Howell functional responses are used to satisfy this assumption. The generalist predator has an alternative food source, and hence, the Leslie-Gower type generalist predator is also included in the system. The impact of fear on the growth rate of the prey due to their fear of predation is also considered. All the properties of the solution of the model are discussed. The equilibria and their stability analysis are studied. The persistence conditions are established. The regions of asymptotic stability of these equilibria are  determined with the help of the Lyapunov method. Local bifurcation is investigated whenever it exists. Finally, extensive numerical simulations are performed to understand the effects of varying the parameter values on the dynamics of the system using two sets of parameters. It is watched that system (4) is wealthy in its dynamical behavior including all possible attractors' types such as point attractor, periodic attractor, and chaos. The fear coefficients have a clear control on the dynamics of the system so that as the fear increases, the chaotic regions decrease and the solution settled to CP or periodic dynamics in the boundary planes. Increasing the preference rates, which are equivalent to increasing the prey's group defense, leads to loss of the system persistence and the solution approaches asymptotically to FAP. Increasing the death rate of the specialist predator causes extinction in both the predators and the system approaches to FAP. Finally, system (4) is very sensitive to changes in the abundance of alternative food sources.

Data Availability
All supportable data are given in the main text.

Conflicts of Interest
The authors declare that there are no conflicts of interest.  Alfa-13 Max (z) Figure 16: The BD of system (4) as a function of α 13 in the range 1 < α 13 < 8 using the parameter values in (90).

20
Journal of Applied Mathematics