Dynamics of a Predator-Prey Model with Fear Effect and Time Delay
In this paper, we propose a time-delayed predator-prey model with Holling-type II functional response, which incorporates the gestation period and the cost of fear into prey reproduction. The dynamical behavior of this system is both analytically and numerically investigated from the viewpoint of stability, permanence, and bifurcation. We found that there are stability switches, and Hopf bifurcations occur when the delay passes through a sequence of critical values. The explicit formulae which determine the direction, stability, and other properties of the bifurcating periodic solutions are given by using the normal form theory and center manifold theorem. We perform extensive numerical simulations to explore the impact of some important parameters on the dynamics of the system. Numerical simulations show that high levels of fear have a stabilizing effect while relatively low levels of fear have a destabilizing effect on the predator-prey interactions which lead to limit-cycle oscillations. We also found that the model with or without a delay-dependent factor can have a significantly different dynamics. Thus, ignoring the delay or not including the delay-dependent factor might result in inaccurate modelling predictions.
Predator-prey interaction is a central topic in ecology and evolutionary biology which has been extensively studied from different aspects by many researchers over the last few decades. Mathematical models have played a vital role in understanding these complex scenarios. Most of the existing predator-prey models are based upon classical Lotka–Volterra formalism, and it is usually assumed that the predators can impact the prey population only through direct killing as direct predation is relatively easy to observe in nature. However, the presence of predator may significantly change the behaviour and physiology of prey to such an extent that it could affect the prey population more effective than direct predation [1–3].
In the presence of predators, many prey species change their behaviour due to predation risk and show a variety of antipredator responses, which incorporates foraging activity , habitat alterations, vigilance  and some physiological changes, etc. [3, 6]. Many researchers observed that due to fear of predation risk, the reproduction of the scared prey decreases. For example, in the Greater Yellowstone Ecosystem, wolves influence the reproductive physiology of Elks . Snowshoe hares shift to relatively safe but less profitable microhabitats . Birds react to the sound of predator species with antipredator defenses, and they flee from their nests at the first sight of danger . Other pieces of evidence suggest that fear can also affect populations like mule deer  and dugongs .
In 2011, Zanette et al. conducted a manipulation on song sparrows during a whole breeding season and showed that there is a reduction in the offspring production of the song sparrows due to the fear from predator . Very few researchers have investigated the impact of fear effect in a predator-prey model with the help of mathematical modelling. Wang et al. first proposed a predator-prey model by incorporating the cost of fear into prey production. Their analysis showed that fear effect can stabilize the oscillation of the system . Wang and Zou formulated a predator-prey model with the cost of fear and adaptive avoidance of predators. The prey population is divided into juvenile and mature stages; this model is represented by a system of delay differential equations with maturation delay . They theoretically analyzed the model dynamics for constant defense level and adaptive defense level. Numerical simulations showed that both strong adaptation of adult prey and the large cost of fear have destabilizing effect, but large population of predators has a stabilizing effect on the system. Sasmal and Takeuchi studied a predator-prey system with fear effect, group defense, and Holling-type IV functional response; they observed that the model can exhibit multistability . Panday et al. studied the impact of fear in a tritrophic food chain model with Holling-type II functional response and observed that fear can stabilize the system from chaos to stable focus through the period-halving phenomenon, and chaotic dynamics can be controlled by the fear factors . Sarkar and Khajanchi  proposed and analyzed a prey-predator system introducing the cost of fear into prey reproduction with Holling-type II functional response. Their mathematical analysis suggests that strong antipredator responses can stabilize the prey-predator interactions by ignoring the existence of periodic behaviors.
Inspired by the aforementioned works, the aim of the current work is to develop a predator-prey model, which incorporates the fear effects of predators on prey as well as the time delay. Further, we investigate the impacts of important parameters on the dynamics of system. The rest of the paper is organized as follows. In the next section, we will give the underlying assumptions and formulate the model. In the following section, we present some results on the stability of equilibria and existence of Hopf bifurcation. In Section 4, we investigate the properties of the bifurcating periodic solutions by using the normal form theory and center manifold theorem. We show the permanence of the system in Section 5. In Section 6, we provide some numerical simulation results which reveal some potential roles that the fear effect and the time delay may play in predator-prey interactions. A brief discussion completes the paper.
2. Model Formulation
Let and be the densities of prey and predator populations at time , respectively. We assume that the prey population follows a logistic growth without any predator population, which can be split into three different parts, namely, the birth rate, the natural death rate, and the density dependent death due to intraspecific competition among preys. Let be the birth rate of prey, be the natural death rate of prey, and be the death rate due to intraprey competition. Since the fear effect will reduce the production, we will multiply the production term of prey by a factor which accounts for the cost of antipredator defense due to predation fear. Also, it is reasonable to assume that the production of predators after predation is not instantaneous; for the sake of simplicity, we suppose that the time lag required for gestation of predator is . Therefore, our predator-prey model with fear effect is given byHere, the parameter refers to the level of fear, which is due to the antipredator response of prey. In this study, we adopt the widely used Holling-type II functional response of predators , where is the predation rate and is the half-saturation constant. is the conversion efficiency from prey biomass to predator biomass and is the death rate of predator. Since predator goes through a gestation stage, and it endures the death rate, then a factor appears in the first term of the second equation of system (1). The parameters , and are all positive constants, and and are nonnegative constants.
The initial conditions take the form ofwhere denotes the space of continuous functions mapping the interval into and denotes the norm of an element in by .
It is straightforward to show (see ) the following.
Lemma 1. The solution of system (1) with an initial condition (2) is unique and nonnegative for all .
Note that from the first equation in (1), we obtain ; by a comparison argument, we see that if , then , and applying the theory of asymptotically autonomous systems to the second equation in (1), we also have . That is, when , both prey and predator species go to extinction. So, throughout this paper, we assume that .
Proof. We first consider the systemwhose solution satisfiesSince , then by the comparison principle, we haveHence, for any small enough, there exists a sufficiently large such thatFor and , we haveIt follows thatThis completes the proof. □
Lemma 3 (see , Lemma 2.1). Consider the following equation:where;for. We have(i)if , then (ii)if , then It is easy to check that (1) has a trivial extinction equilibrium . This system also has a semitrivial equilibrium under the condition . To find the positive (coexistence) equilibrium , letting the right-hand side be zero, we then get , which exists if , and is a positive root of the equationwhere ,,. Clearly, (10) has one positive root if and only if ; this inequality is equivalent to
Then, we conclude that (1) has a unique positive equilibrium if (11) holds, and the expression for is as follows:
3. Stability Analysis and Hopf Bifurcation
Theorem 1. The trivial extinction equilibrium is always an unstable saddle point. The semitrivial equilibriumis globally asymptotically stable inif (11) is reversed and is unstable if (11) holds; meanwhile, exists.
Proof. For equilibrium , the characteristic equation isObviously, (13) has one positive root and one negative root . Hence, is a saddle point.
The characteristic equation at is as follows:Let ; then it is clear that (14) has one root and the other roots are determined by the equation .
If (11) is reversed, then . Let with be a root of . Then, we haveIf , thenA Contradiction. This shows that all roots of must have negative real parts. Therefore, is locally asymptotically stable.
Next, we will show that is globally attractive if (11) is reversed. Since , we can choose small enough such thatFrom (6), we know thatfor . But for the equationby (17) and Lemma 3, it follows that . Therefore, an application of the standard comparison argument yields . Then, the limiting equation for isTherefore, , and is globally asymptotically stable in .
Note also that if (11) holds, then , and . Hence, has at least one positive root and is unstable.
When , (21) yields
It is easy to see that and . To ensure the local stability of , we need .
Squaring and adding both equations of (24), it follows thatwhere . A straightforward computation shows that
Now, we can summarize the above result on the local stability of the positive equilibrium .
Theorem 2. Assume that (11) holds; thenexists. Moreover, (i) if and, then (21) has no purely imaginary root, and is locally asymptotically stable for all . (ii) If , then (21) has a unique pair of purely imaginary roots , and may lose stability as increases from zero, where and is the unique positive root of (25).
To further investigate the stability of under condition (ii) of Theorem 2, we need to check the transversality condition at . Since and are both time delay-dependent, it is hard to discuss the transversality condition, if we suppose that the coefficient in conversing prey into predator is a constant, independent of delay , then we can extend the result in Theorem 2 as follows.
Theorem 3. Ifis a constant, and , then there exists a sequence of values of :such that is locally asymptotically stable whenand unstable when. Furthermore, system (1) undergoes a Hopf bifurcation at when, , where is defined by (28).
Proof. From (24), we know that corresponding to isTo complete the proof of Theorem 3, we only need to show thatDifferentiating (21) with respect to , we haveIt follows thatand hence,Then, we have . Therefore, the transversality condition holds, and the Hopf bifurcation occurs at .
4. Properties of Hopf Bifurcations
In this section, we will study the properties of the Hopf bifurcations by using the normal theory and the center manifold theorem due to Hassard et al. .
Let ,, and dropping the bars for simplification of notations, system (1) is transformed into a functional differential equation in aswhere , and , are given, respectively, byandwhere
By the Riesz representation theorem, there exists a function of bounded variation for , such that
In fact, we can choosewhere denotes the Dirac delta function. For , definex
Then, system (1) is equivalent towhere for . For , defineand a bilinear inner productwhere . Then, and are adjoint operators. By the discussion in Section 3, we know that are eigenvalues of . Hence, they are also eigenvalues of . We first need to compute the eigenvectors of and corresponding to and , respectively.
Similarly, we can obtain the eigenvector of corresponding to , where
Choosing as , then by (42), we see .
Now, we will use the algorithms given in  and using a computation process similar to that in [21–23], we get the coefficients used in determining the qualities of bifurcating periodic solutions:where
Thus, we can determine and . Furthermore, we can compute by (45). Then, we can compute the following values:
The properties of bifurcating periodic solution in the center manifold at the critical value are determined by the values in (47). From , we know that determines the directions of the Hopf bifurcation: if , then the Hopf bifurcation is supercritical (subcritical); determines the stability of the bifurcating periodic solutions: the bifurcating periodic solutions are stable (unstable) if ; and determines the period of the bifurcating periodic solutions: the period increase (decrease) if .
In order to prove the uniform persistence of system (1), we need the uniform persistence theorem for infinite dimensional systems from . Let be a complete metric space. Suppose that is open and dense in , with , . Assume that is a semigroup on satisfying
Let and let be the global attractor for . Assume further that(i)There is a such that is compact for ;(ii) is point dissipative in ;(iii) is isolated and has an acyclic covering , where ;(iv) for .
Then, is a uniform repeller with respect to , i.e., there is an such that for any ,, where is the distance of from .
Now, we consider the persistence of system (1).
Proof. LetIt suffices to show that there exists an such that for any solution of system (1) initiating from ,. To this end, we verify below that the conditions of the above uniform persistence theorem are satisfied. Firstly, we show is positively invariant. By (1), we haveandThen, and for all if. This implies that is positively invariant.
Similarly, we can show that is positively invariant, and condition (48) is satisfied. We have verified the point dissipativeness of the semiflow of system (1) in 2.
Denote the limit set of the solution of system (1) starting in by . LetFor any given , we have for all . If for some , then for all , and as . If for some , we have the following two cases.
Case 1. . Then from (52), we have for all , and , it follows that or as .
Case 2. . Then we must have ; from the first equation of (1), we have for all ; again from the second equation of (1), we have as .
Finally, we see that for any given , as , or as . Thus, . If invariant sets and are isolated, then is isolated and is an acyclic covering. The isolated invariance of and will follow from the following proof.
We now show thatwhere () denotes the stable manifold of (). Suppose ; then there exists a solution ,, of (1) with initial conditions in , such thatSince , we can choose sufficiently small so thatLet sufficiently large such thatThen, for , we haveConsider the following equation:By the comparison principle, we have for all . It is easy to see that (59) has a positive equilibrium , which is globally attractive. Note that for all and ; this is a contradiction to (55).
Suppose ; then there exists a solution ,, of (1) with initial conditions in , such thatThen, for any , by (60), there exists a such thatHence, by (1),Since (11) holds, we can choose sufficiently small so that . Then, by Lemma 3, the solution of the linear delay differential equation with positive initial condition must converge to infinity as . Therefore, by the comparison principle, we have as , which is a contradiction to (60). Thus, (54) holds. This completes the proof.
6. Numerical Simulations
In order to investigate the role of time delay and fear in our model (1), we perform some numerical simulations to illustrate the analytical results observed in the previous sections. All the simulations are carried out in MATLAB.
First, we choose the parameter values aswhere most of the parameters are taken from .
For the above set of parameter values, in Figure 1, we observe that the system shows a limit cycle around the interior equilibrium . If we further increase the value of , for , we see that the solutions approach the coexistence equilibrium; see Figure 2.
To observe the long-term behavior for a range of values of , in Figure 3, we draw the bifurcation diagram with respect to . From Figure 3, we can see that as we increase the value of , the prey and predator populations show oscillatory coexistence for the small values of ; then, both populations show stable coexistence. This indicates that increasing delay may induce a transition from the state where the populations of the prey and predator oscillate periodically to a stability situation. And also, we note that the prey population increases, while the predator population decreases as we gradually increase , and finally the predator population goes to extinction; the system becomes stable around the semitrivial equilibrium .
To investigate the impact of fear when the system is unstable around the positive equilibrium, we fix . We draw the bifurcation diagram considering as a bifurcation parameter. Figure 4 shows that as we increase the level of the fear effect in prey, the coexistence of prey-predator changes from stable periodic oscillations to stable equilibrium. Figure 4 indicates that the prey population tends to a steady state. The fear effect will not affect the prey population over the long term. Biologically, after a certain level of fear with prey population, the fear has no effect due to physiological impact when the prey populations are habituated with fear from predator species, while the predator population eventually will decrease with increasing .
If we choose parameter values as those in Figure 1 and let vary in the interval , Figure 5 gives the bifurcation diagram. From this figure, we can see that the predator population can survive only for the large values of . If the value of the predation rate is low, then the predator population cannot survive. As we increase the value of , the coexistence of prey-predator changes at stable equilibrium to stable oscillatory coexistence.
To investigate the properties of the bifurcating periodic solution, we neglect the delay-dependent factor ; that is, we consider the following system:
We use the following parameter values:
Then, (64) has a positive equilibrium and satisfies the conditions indicated in Theorem 3; then, is asymptotically stable when . Take ; we then have , , , , and . Thus, is stable for as illustrated in Figure 6. When passes through the critical value , will lose its stability and a Hopf bifurcation will occurs; that is, a family of periodic solutions bifurcate from . Since , , then the Hopf bifurcation is supercritical and its direction is . These bifurcating periodic solutions from at are stable; see Figure 7.