Abstract
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.
1. Introduction
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 [4], habitat alterations, vigilance [5] 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 [7]. Snowshoe hares shift to relatively safe but less profitable microhabitats [8]. Birds react to the sound of predator species with antipredator defenses, and they flee from their nests at the first sight of danger [1]. Other pieces of evidence suggest that fear can also affect populations like mule deer [9] and dugongs [10].
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 [11]. 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 [12]. 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 [13]. 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 [14]. 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 [15]. Sarkar and Khajanchi [16] 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 [17], 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 [18]) 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 .
Lemma 2. There exists an such that for any solution of system (1) with an initial condition (2), there must be a such that and for .
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 [19], 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.
Theorem 1 implies that loses its stability to the occurrence of the positive equilibrium when (11) holds. The local stability of is addressed in the following theorem.
When (11) holds, linearizing system (1) at the positive equilibrium , we get the characteristic equation as follows:where
When , (21) yields
It is easy to see that and . To ensure the local stability of , we need .
Let , is a root of (21); then, (21) implies
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. [20].
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.
Suppose is the eigenvector of corresponding to ; then . Then, from the definition of , (34), (37), and (38), we have
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 [20] 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 [20], 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 .
5. Permanence
In order to prove the uniform persistence of system (1), we need the uniform persistence theorem for infinite dimensional systems from [24]. 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).
Theorem 4. Let (11) hold; then there exists an such that every solutionof system (1) with initial conditions ,and for somesatisfies
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.
In Lemma 2, we have shown that the solutions of system (1) are ultimately bounded; together with Theorem 4, we now give the main result in this section.
Theorem 5. System1 is permanent provided11 holds and the initial conditions are shown inTheorem 4.
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 [12].
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.
(a)
(b)
(a)
(b)
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 .
(a)
(b)
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 .
(a)
(b)
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.
(a)
(b)
(c)
(d)
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.
(a)
(b)
(a)
(b)
In Figure 8, we explore the effects of delay on the dynamics of system (64). Figure 8 shows that as increases, at first the positive equilibrium is stable; then when is greater than a certain critical value, periodic solution occurs; this is coincident with Theorem 3, which implies that delay destabilizes the equilibrium and induces oscillatory dynamics. By comparing Figure 3 and Figure 8, we see that the difference is big in the dynamics when the delay-dependent factor is incorporated or not.
(a)
(b)
In Figure 9, taking as a bifurcation parameter, we found that the fear effects can stabilize the system. In Figure 10, we draw the bifurcation diagram with respect to ; the shapes of the bifurcation diagrams for system (64) in Figure 10 and the original system (1) in Figure 5 are very similar.
(a)
(b)
(a)
(b)
7. Discussion
In this paper, taking into account the time delay corresponding to the gestation period and the effect that the fear of predators has on prey, we have proposed a predator-prey model with Holling-type II functional response. We have shown that the boundary equilibrium is globally asymptotically stable if (11) is reversed and pointed out that the instability of implies the permanence of the system and existence of the positive equilibrium . Furthermore, for the system without delay-dependent factor, by choosing the delay as a bifurcation parameter, sufficient conditions for the existence of Hopf bifurcation have been presented. By using the theory of normal form and center manifold theorem, the explicit formulae determining the direction, stability, and other properties of bifurcating periodic solutions have been derived.
The novelty of this paper is that we have incorporated both the fear effect and a delay-dependent factor that ensures the predators that do not survive the delay interval do not contribute to the growth of the predator population. Our investigation shows that the dynamics of system (1) are affected significantly by not only the biomass transfer delay, but also the level of the fear effect. In [25], the delay-dependent factor was neglected.
Numerical simulations are performed to show the roles that time delay and fear effects can play in predator-prey interactions. For the original system (1), we have shown that both large time delay and high levels of fear (or equivalently strong antipredator responses) can stabilize the predator-prey system by excluding the existence of periodic solutions. And also, we observe that if the time delay as well as the fear level increases, the predator population goes to extinction. While for system (64) which neglects the delay-dependent factor, numerical simulations show that the delay has destabilizing effect while the large cost of fear has a stabilizing effect on the predator-prey interactions. Figure 8 indicates that ignoring the delay-dependent factor results in model in (63) that does not predict the extinction of the predators, even if the delay is very large. Therefore, time delay and fear effect significantly influence the dynamics of the system and play important roles in ecosystems stability.
In the current paper, we assumed that the fear effects only reduce the birth rate of the prey, but it may also increase the death rate of adult prey due to long-term physiological impacts [11]. To make the model more realistic, some other important mechanisms can also be incorporated; for example, spatial effect is an important factor due to foraging behaviors of both prey and predator. The functional response plays a vital role in determining the different complex dynamical behaviors; different functional response functions can be chosen to reflect the group defense. It would be interesting and challenging to consider the aforementioned possible extensions; we leave these problems for further investigation.
Data Availability
All data generated or analyzed during this study are included in this article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (11801431), the International Cooperation on Science and Technology of the Key Research and Development Program of Shaanxi Province (2018KW-021), and the Natural Science Basic Research Plan in Shaanxi Province of China (2021JM-445).