Abstract
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 Leslie-Gower 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.
1. 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–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 prey-predator 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.
2. 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 Holling-type IV functional response or Sokol-Howell type of functional 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 represents the prey density at the time , is the specialist-predator density at the time , while is the generalist-predator density at time . Accordingly, the dynamics of such a food-web system can be described as follows: where is the logistic growth rate of the prey species, while and , with , are the Sokol-Howell functional responses for the specialist-predator and generalist-predator, respectively. The parameters of system (1) are supposed to be positive, and they are described as the following: the parameter represents the intrinsic growth rate, is the carrying capacity of the environment for the prey, the positive parameters , and are the maximum attack rates, the positive parameters , and are the preference rates of the predator to their prey species, represents the conversion rate from individuals of the prey to specialist-predator, and is the intrinsic growth rate of a generalist-predator, while the positive parameters are the respective food preference of generalist-predator from the prey and specialist-predator, respectively. In addition to the above, the square term 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 is a constant added in the denominator to normalize the residual reduction in the generalist-predator due to a heavy lack of the preferable food and , while 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 , and represent the fear 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.
Therefore, the nondimensional system corresponding to system (2) is given by
Here, the interaction functions are defined on . Moreover, since the right-hand side functions of system (4) are continuous and have continuous partial derivatives, hence they are Lipschitz functions. Thus, the solution of system (4) exists and is unique.
Theorem 1. The solutions of system (4), which start in are uniformly bounded.
Proof. Let be a solution of system (4) with nonnegative initial condition . From the first equation in system (4), it is obtained that
hence, by solving this differential inequality, it is obtained that as . Define the function ; then, the following is resulting
or .
Then, solving this differential inequality gives for :
Assume that , then by differentiating it with respect to time gives
where .
Therefore, rearranging and substituting the upper bound of the logistic term give
where .
Then, solving the last differential inequality as gives
Hence, all solutions of system (4) that initiate in are confined in the region , that is, all the solution are uniformly bounded.☐
3. Equilibria and Their Local Stability
System (4) has at most six equilibrium points, which can be described as follows: (1)The trivial equilibrium point (TP) denoted by always exist(2)The first axial equilibrium point (FAP) denoted by always exist(3)The second axial equilibrium point (SAP) denoted by , where is any positive real number, exists if and only if(4)The generalist-predator free equilibrium point (GPFP) that denoted by , wherewith , and , exists uniquely in the positive quadrant of -plane if and only if the following conditions hold: (5)The specialist-predator free equilibrium point (SPFP) denoted by can be represented aswhere , and . Obviously, SPFP exists uniquely in the -plane if and only if the following condition holds: (6)The coexistence equilibrium point (CP) denoted by can be represented as follows:The point represents the intersection point of the following two isoclines in the positive quadrant of the -plane.
Simple computation displays that the two isoclines given by (16a) and (16b) have a unique positive intersection point, denoted by , provided that the following sufficient conditions hold:
Here, and represent the positive root of the following two polynomials, respectively: where
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 can be represented as where with , , , , and .
Consequently, for the TP, the JM becomes
Then, the eigenvalues of are given by . Since there is a positive equilibrium point, then is an unstable nonhyperbolic point.
For the FAP, the JM can be represented as
Then, the eigenvalues of are clearly given by . Hence, the FAP is a nonhyperbolic point due to the existence of zero eigenvalue.
For the SAP, the JM can be represented as
Then, the eigenvalues of are clearly given by , , and . Hence, the SAP is a nonhyperbolic point due to the existence of zero eigenvalue.
For the GPFP, the JM can be represented as where with , , and . Clearly, the characteristic equation of can be represented as
Because the JM has a zero eigenvalue (), then is a nonhyperbolic equilibrium point.
For the SPEP, the JM can be determined as where with , , , , and . Therefore, the characteristic equation of can be represented as
Clearly, the roots (eigenvalues) of equation (30) can be represented as
It is easy to verify that all the eigenvalues of 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 , , , , and . Consequently, the characteristic equation of can be represented as where , 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 , , and . Accordingly, the following theorem can be proved simply.
Theorem 2. The CP is L.A.S if the below sufficient conditions hold.
4. 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 -plane and -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 -plane and -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 , clearly this function and function in the interior of of -plane. Moreover, it is obvious that
Therefore, 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 -plane for subsystem (39).
Similarly, with the help of Dulac function , which is , and function in the interior of of -plane, it is observed that there are no periodic dynamics in the interior of a positive quadrant of -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 , where are positive constants, and for all with if any one of , or approaches zero. Therefore, direct computation gives where the functions are given in system (4). Now, according to the average Lyapunov method, the proof is satisfied provided that for all boundary equilibrium points. Therefore, where are given in equation (20). Clearly, we have that Obviously, by choosing the arbitrary positive value of sufficiently larger than that of , it is obtained that . Note that conditions (44a) and (44b) guarantee that . Again, by choosing the arbitrary positive value of sufficiently larger than that of , it is obtained that provided that condition (44c) holds. Clearly, condition (44d) guarantees that . Finally, under condition (44e). Hence, system (4) is uniformly persistent, and the proof is done.☐
5. 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 that satisfies the following sufficient conditions: where the symbol is given in the proof.
Proof. Define . Clearly, the function is a positive definite function that is , while , for all values in the region . Then, using some algebraic manipulation gives that
Consequently, by using additional computation, the following is obtained:
Now, according to the upper bound for and , so as , we have and . Therefore, it is obtained that
Consequently, using conditions (52b) and (52c) with the bound of logistic term yields
where .
Accordingly, condition (52d) guarantees that , which means it is a negative definite. Then, 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.
Theorem 5. The SAP is A.S in the subregion of that satisfies the following sufficient conditions: where the symbol is given in the proof.
Proof. Define . Clearly, the function is a positive definite function that is , while , for all values in the region . Then, using some algebraic manipulation gives that
where ; clearly, is positive at the SAP.
Consequently, by using conditions (57a) and (57b), it is obtained that . Hence, the derivative of the function is a negative definite, and then, is a Lyapunov function. Hence, the SAP becomes A.S point for any point that belongs to the subregion that satisfies the given conditions and the proof is done.☐
Theorem 6. The GPFP is A.S in the subregion of that satisfies the following sufficient conditions: where all the symbols are given in the proof.
Proof. Define . Clearly, the function is a positive definite function that is , while , for all values in the region . Then, using some algebraic manipulation gives that
Furthermore, computation gives
Therefore, using conditions (59a), (59b), (59c), and (59d) gives that
where , , , , .
Accordingly, using condition (59f), it is observed that , which means that is a Lyapunov function and hence the GPFP is an A.S for any trajectory starting from a point that belongs to the region that satisfies the above conditions.☐
Theorem 7. Assume that the SPFP is L.A.S; then, it has a basin of attraction that belongs to the that satisfies the following conditions: where the symbols , and are given in the proof.
Proof. Define . Clearly, the function is a positive definite function that is , while , for all values in the region . Then, using some algebraic manipulation gives that Consequently, using conditions (63a), (63b), (63d), (63e), and (63f) gives Moreover, using condition (63c) leads to Accordingly, , which means that is a Lyapunov function and hence the SPFP is A.S for any trajectory starting from a point that belongs to the region that satisfies the above conditions.☐
Theorem 8. Assume that the CP is L.A.S; then, it has a basin of attraction that belongs to the that satisfies the following conditions: where the symbols are given in the proof.
Proof. Define Clearly, the function is a positive definite function that is , while , for all values in the region . Thus, after some algebraic manipulation, it is obtained that where Consequently, using the above conditions gives Obviously, , which gives that 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.☐
6. 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 be any nonzero vector, and the symbols are given in equation (20). Accordingly, the local bifurcation near the equilibrium points is studied in the next theorems.
Theorem 9. Suppose that condition (32a) is satisfied; then, system (4) at the SPFP undergoes a transcritical bifurcation when the parameter passes through the value , provided that the following condition holds: where the symbols and are given in the proof.
Proof. Direct computation shows that the JM of system (4) at can be represented as
where are given in equation (28). Obviously, the matrix has two eigenvalues having negative real parts due to condition (32a) and the third one is zero, say . Hence, is a nonhyperbolic point.
Let be the eigenvector corresponding to the eigenvalue . Thus, gives that , where be any real number.
Now, let represent the eigenvector corresponding to the eigenvalue of the matrix . Thus, gives that , where be any real number. Now, according to Sotomayor’s theorem, it is obtained that
Therefore, ; hence, system (4) has no saddle-node bifurcation. Moreover, since
Then, .
Also, by using equation (73), it is obtained that
where
Accordingly, it is easy to verify that
Clearly, , under condition (75), and then, a transcritical bifurcation take place in the sense of Sotomayor.☐
Theorem 10. Suppose that conditions (38a)–(38c) and (38e) are satisfied; then, system (4) at CP undergoes a saddle-node bifurcation when the parameter passes through the value , provided that the following condition holds: where are given in the proof.
Proof. Direct computation shows that the JM of system (4) at can be represented as
where all the coefficients are given in equation (33). Clearly, the determinant of that is given by in equation (35) is zero at . Therefore, the characteristic equation of has zero root, say , and then, the coexistence equilibrium point becomes a nonhyperbolic point.
Let be the eigenvector corresponding to the eigenvalue . Thus, gives that , where be any real number, and , while . Clearly, , and , due to conditions (38a)–(38c) with (38e).
Now, let represent the eigenvector corresponding to the eigenvalue of the matrix . Thus, gives that , where be any real number, and , while . Clearly, , and , due to conditions (38a)–(38c) with (38e).
Now, according to Sotomayor’s theorem, it is obtained that
where is given in equation (20). Therefore, it is obtained that
where is given in equation (33). Now, by using equation (73), it is obtained that
where
Accordingly, it is simple to verify that
Clearly, , under condition (82), and then, a saddle-node bifurcation takes place in the sense of Sotomayor.☐
7. Numerical Simulation
In this part, the dynamics of system (4) are studied numerically. The objective is to complete the vision of the global dynamics of the system and detect the effect of varying the parameter values on the system’s dynamical behavior. Two sets of hypothetical parameter values are used to satisfy our objective. For the following set of parameter values, other sets can be used to study the system:
It is watched that system (4) has a globally asymptotically stable CP as shown in Figure 1.
(a)
(b)
(c)
(d)
According to Figure 1, system (4) is globally asymptotically stable given by which means system (4) is persistent at the CP. However, for the data (89) with increasing , system (4) loses its persistence and approaches asymptotically to GPFP as shown in Figure 2 for the typical value of .
(a)
(b)
(c)
(d)
On the other hand, varying the other fear’s parameters leads to a quantitative change in the position of the CP and system (4) still persistent. Now, the effects of varying the parameters and are presented in Figures 3 and 4, respectively.
(a)
(b)
(c)
(d)
(a)
(b)
According to Figure 3, increasing the value of 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 destabilizes system (4) and the solution persists at periodic dynamics. Moreover, it is observed that the effect of the parameter is similar to that of .
Now, the effects of varying the parameters and on system (4) dynamical behavior are shown in Figures 5 and 6, respectively.
Obviously, Figures 5 and 6 show that the parameters and have opposite effects on the dynamic behavior of system (4) so that as increases (decreasing ), the population of the prey decreases while the populations of the predators increase and hence system (4) loses their persistence and approaches to FAP as decreasing (increasing ). Furthermore, it is observed that varying the parameter has a quantitative effect on the dynamic behavior of system (4) and system (4) persists at the CP.
The effects of varying the parameters and on system (4) dynamical behavior are explaining in Figures 7 and 8, respectively.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
According to Figures 7 and 8, it is watched that increasing the parameter leads to losing the persistence of the system by approaching the SPFP first and then to SAP. However, decreasing (or increasing) the parameter leads to losing the persistence of the system and the solution approaches to GPFP (SPFP). It is observed that the parameter has similar effects on the dynamic behavior of system (4) as that for .
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.
(a)
(b)
(c)
(d)
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.
(a)
(b)
(c)
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 , and are shown in Figures 12–16.
According to Figures 12 and 13, there are ranges around and 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 -plane. Similar effects on the dynamics behavior of system (4) are observed for the parameters as those of and . Finally, the other BD are clearly shown the existence of chaotic regions; however, it is watched that increasing these parameters leads to extinction in one predator species and the solution approaches to periodic dynamics in the plane.
8. 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.
Acknowledgments
This article is funded by the authors themself only.