The impact of fear factor and self-defence on the dynamics of predator-prey model with digestion delay

Abstract: In this paper, we propose both deterministic and stochastic predator-prey models with digestion delay, incorporating fear factor and self-defence. For the deterministic model, the existence and stability of the equilibrium are investigated and the occurrence of Hopf bifurcation is studied. For the stochastic model, we investigate the existence of a unique global positive solution of the model and analyze the asymptotic behavior of the global solution around the equilibriums of the deterministic model. Finally, numerical simulations are carried out to verify our analytical results, which indicate that the intensity of white noise, fear factor and self-defence have a significant relationship with the dynamics of the predator-prey model and expand the theoretical analyses.


Introduction
Predator-prey dynamics model is one of the most important topics in mathematical ecology research, which is universal and significant. Since the early model proposed by Lotka [1] and Volterra [2], the studies of predator-prey interactions have attracted rapidly increasing attention. Generally, a predator-prey model can be described by a system of ordinary differential equations as follows: , y(t))y(t), dy(t) dt = f 2 (y(t)) + cp(x(t), y(t))y(t), where x(t) and y(t) are the population size of the prey and predator at time t, respectively. f 1 (x(t)) is the growth rate of prey population; f 2 (y(t)) is the growth rate of predator population. c is a positive constant (c < 1 due to the biological significance) which represents the coefficient in converting prey into a new predator. p(x(t), y(t)) represents the functional response, i.e., how predator consumes the prey species. Vast and rich studies [3][4][5][6][7] have been done by taking different kinds of functional responses into consideration to investigate the effect of direct predation. However, in the real world, the presence of predators also has indirect effects on the prey population. The outcomes of recent theoretical studies and experimental shreds of evidence [8][9][10][11] have unveiled that the indirect impact of predator on prey species sometimes becomes more crucial in altering the dynamics of ecosystems. The threat of predation alone is enough to reduce the growth, survival and fecundity of prey, which is arbitrated through an alternation of prey's behavioural, morphological or psychological trait [12,13]. For instance, the sounds of predator will induce continuous fear on birds which compels them to leave their primary nest and move to a safer zone. This kind of behavior leads to serious damage to the birth rate of the new offspring inevitably. Frequent changes in habitat will not only affect the basic reproduction of the individuals but also cause a serious impact on the survival of adults [8]. In 2009, Sheriff et al. [14] reported that when pregnant snowshoe hares were exposed to a trained dog (predator), the birth of cubs was significantly reduced. In addition, Zanette et al. [15] performed an experiment on song sparrows and the outcome exhibited that the quantity of posterity of the song sparrows was decreased by 40% due to fear induced by predators. Shreds of evidence suggest that the indirect impact of predator should not be ignored when analyzing the dynamics of prey-predator interaction. In this context, Wang et al. [16] were the first to establish an ordinary differential equation model including the fear effect. Through their work, we know that the level of fear would not affect the dynamic of the system if the predator consumes the prey according to linear functional response, while for the model with Holling type II functional response, the level of fear affected predator-prey interactions in several ways. Subsequently, they proposed a predator-prey model with the cost of fear and adaptive avoidance of predators in [17]. Inspired by their research, Sasmal [18] investigated a model with fear effect and Allee effect where the prey is taken as infected. Sha et al. [19] and Hossain [20] focused on eco-epidemiological models with fear effect. Few other works have been found in this direction, such as discrete [21], delayed [22] and diffusive [23,24] models. Recently, Wang and Zou [25] investigated a model with the benefits of the anti-predation response and digestion time which is necessary for biomass transfer from prey to predator after predation. The chance of the prey being vulnerable by predators was thought to be decreased for the presence of anti-predation response, which makes the model more in line with the actual situation. Obviously, such anti-predation responses are beneficial to the survival of the prey. In addition, it is worth noting that in reality, predation comes with risks. The prey will resist when it encounters a predator instead of just escaping. For example, leopards run the risk of being stabbed by the hard spines on the back and tail of the porcupine when preying on the porcupine. Once stabbed, the leopard may die due to wound infection or unable to prey. The death of the predator caused by self-defense of the prey can be regarded as the predation of the predator by the prey, which will reduce the growth of the predator population. The reduction is the cost of predation and it may exist whether the predation is successful or not. Noting that this kind of passive predation does not need time to digest and will not affect the next predation, the linear function response is suitable. What needs to be pointed out is that self-defence is different from conservative and protective behaviour and group defence, as it is the fight for survival after the failure of avoidance strategy. By fighting against the predator, the chance of being predated will obviously be decreased, which will benefit the prey population. Moreover, the prey's self-defence will definitely consume its energy and this will clearly affect the production rate and death rate of the prey. To avoid making the system too complicated, we only consider the cost of the prey's self-defence for the predator. Hence, we get the following model with fear factor, digestion delay and self-defence: where x(t) and y(t) have the same meaning as the system (1.1). τ ≥ 0 is the average time needed for biomass transfer after predation from prey to predator. All the parameters are positive constants: r represents the intrinsic growth rate of the prey; d i (i = 1, 2) represent the death rate; a and h donate the intra-specific competition of the prey and predator, respectively; k represents the level of fear; c i (i = 1, 2) represent the decreasing rate of reproduction and predation, respectively; c is the biomass transform efficiency constant and θ represents the level of self-defence. From the realistic point of view, system (1.2) extend model (3.24) in [25] by considering the self-defence of prey and intra-specific of predator. Therefore, for the dynamics of the deterministic model, we will focus on both of these parameters in addition to the fear factor. On the other hand, parameters involved in the system always fluctuate around some average values since the environmental fluctuations, such as temperature, humidity, rainfall and so on. In view of this, many researchers [26][27][28][29] have introduced stochastic environmental perturbations into deterministic models. Since the death rate d i (i = 1, 2) are most vulnerable by the environmental fluctuations and most influential in determining the dynamics of the system, we will focus on the effects of perturbations on these parameters. In this paper, we perturb −d 1 and −d 2 by where B 1 (t) and B 2 (t) are standard Brownian motion with B 1 (0) = B 2 (0) = 0, σ 2 1 and σ 2 2 are the intensity of white noise. Whereupon, we obtain the following stochastic model: System (1.3) satisfies the initial conditions where φ 1 (ς) and φ 2 (ς) are both nonnegative continuous functions on [−τ, 0]. In addition, we suppose that φ i (0) > 0 (i = 1, 2) for biological meaning. Throughout this paper, unless otherwise specified, let (Ω, F , {F t } t≥0 , P) denote a complete probability space with a filtration {F t } t≥0 satisfying the usual conditions (i.e. it is right continuous and F 0 contains all P-null sets). The framework of the paper is arranged as follows. In Section 2, we analyze the existence and stability of the equilibrium points and the existence of Hopf bifurcations of system (1.2). In Section 3, we investigate the dynamic behavior of system (1.3), including the existence and uniqueness of the global positive solution and asymptotic property of system (1.3) around the equilibrium points of the corresponding deterministic system (1.2). Several numerical simulations are introduced to support our theoretical results in Section 4. Finally, we give a brief discussion and summarize the main results in Section 5.

Dynamics of the deterministic model
Obviously, system (1.2) has the following three possible non-negative equilibrium points: (i) The trivial equilibrium E 0 = (0, 0), which always exists; (ii) The predator-free equilibrium E 1 = ( r−d 1 a , 0), which exists when r − d 1 > 0, matching biological meaning; (iii) The positive (coexistence) h and x * is the positive solution of equation One can see that A > 0, but the sign of B, C, D, E and F is hard to confirm. By Descartes' Rule of Signs [30], we have the following existence criterion of the positive root of the equation in view of the sign of the coefficient. (I)The equation has exactly one positive root if The equation has an even number of positive roots, including none , for the rest of conditions. Next, we study the local stability for each of the equilibrium points. The linearized form of system (1.2) at an equilibrium (x,ȳ) is as follows: (2.1) where , We can derive the characteristic equation corresponding to system (2.1) as Theorem 2.1. The trivial equilibrium E 0 is locally asymptotically stable if r < d 1 and unstable if r > d 1 . When r > d 1 , the predator-free equilibrium E 1 is locally asymptotically stable if cp Hence, E 0 is locally asymptotically stable if r < d 1 and unstable if r > d 1 .
When r > d 1 , the predator-free equilibrium E 1 exists, at which the characteristic equation becomes Noting that at E 1 , J 11 = d 1 − r < 0, the stability of E 1 is determined by equation Calculation shows that Eq (2.3) is equivalent to The proof is complete.
Remark 2.1. When θ = 0, (2.4) reduces to d 2 (1 + c 1 k) > (cp − d 2 q(1 + c 1 k)) r−d 1 a , which is consistent with the conclusion obtained in [25]. Equation (2.4) indicates that the stability of the predator-free equilibrium E 1 is related to the value of θ. Unfortunately, as previously analyzed, the presence of h and θ makes the conditions for the existence of the positive equilibrium less intuitive. Numerical simulations in Section 4 shows that when E 1 loss its stability, the positive equilibrium will exist.
(2) If (H2) is satisfied, Eq (2.5) only has roots with negative real parts, which indicates that the positive equilibrium E * is locally asymptotically stable.
(3) If (H3) is satisfied, Eq (2.5) possesses roots with positive real parts, which indicates that the positive equilibrium E * is unstable.
Now, we check whether the roots of Eq (2.2) will cross the pure imaginary axis as τ increases. When τ 0, plugging λ = ωi (ω > 0) into Eq (2.2) to obtain the critical value where Hopf bifurcation occurs. Separating the real and imaginary parts, we obtain By trigonometric identity, we obtain Recalling that has no positive root when Hence no root of Eq (2.2) will cross the pure imaginary axis for all τ ≥ 0, that is, the stability of the positive equilibrium E * will not change as τ increases.
Lemma 2.2. Assume that the positive equilibrium E * of system (1.2) exists.
If Eq (2.8) is reversed, Eq (2.7) has a unique positive root where Then we obtain a sequence of τ n corresponding to ω 0 : (2.9) Hence we obtain possible value for τ at which Hopf bifurcation may occur. When τ = τ 0 , ±iω 0 is a pair of pure imaginary roots of Eq (2.2). For the sake of confirming whether Hopf bifurcation occurs at the first critical value τ = τ 0 , the transversal condition needs to be verified. Differentiating Eq (2.2) with respect to τ gives Therefore, at τ = τ 0 , we have . where The proof about the transversal condition is complete and a Hopf bifurcation exists when τ = τ 0 .
It is worth noting that Eq (2.7) has a unique positive root and By the insightful work of Cooke and Grossman [31], one can see that (i) The unstable positive equilibrium of system (1.2) never becomes stable; (ii) If the positive equilibrium is stable for τ = 0, it becomes unstable at τ 0 and remains so as τ is increased; (iii) Only crossing of the imaginary axis from left to right is possible as τ increases. Therefore, any roots on the right half plane will remain on the right plane and the stability of the positive equilibrium can only be lost but not regained. Then, we obtain the lemma as follows.
Lemma 2.3. Assume that the positive equilibrium E * of system (1.2) exists.
then hypothesis (H1) will never hold, which will simplify Theorem 2.2. Furthermore, Eq (2.8) reduces to 2d 2 J 11 J 12 K 21 − J 2 12 K 2 21 > 0, which happens to be Eq (3.40) in [25]. Hence, Theorem 3.4 in [25] is a special case of Theorem 2.2 under conditions h = 0 and θ = 0. The formula of hypotheses indicates that both the value of h and θ will affect the stability of the positive equilibrium E * . We shall numerically investigate the impact of them in Section 4.

Dynamics of stochastic model
In this section, we mainly study the dynamics of the stochastic system (1.3). We first prove that there exists a unique global positive solution of model system (1.3) with any given positive initial value. Then, we investigate the asymptotic property of the stochastic system (1.3) around the equilibriums of its corresponding deterministic system (1.2).

Existence and uniqueness of the global positive solution
The following result is related to the existence and uniqueness of the global positive solution, which is a prerequisite for researching the long term behavior of system (1.3).
Theorem 3.1. For any given initial value Eq (1.4), system (1.3) has a unique global positive solution (x(t), y(t)) for t ≥ −τ. Furthermore, the solution will remain in R 2 + with probability one, namely (x(t), y(t)) ∈ R 2 + for all t ≥ −τ almost surely. Proof. Note that the coefficients of system (1.3) are locally Lipschitz continuous and for any given initial value Eq (1.4), there exists a unique local solution (x(t), y(t)) on t ∈ [−τ, τ e ), where τ e is the explosion time [32]. In order to demonstrate that this solution is global, it's sufficient to prove that τ e = ∞ a.s.. Let k 0 ≥ 0 be a sufficiently large constant for every component of (x(0), y(0)) all lying within the interval [ 1 For each integer k ≥ k 0 , we define the stopping time as follows: Throughout this paper, we set inf ∅ = ∞(as usual ∅ denotes the empty set). Clearly, τ k is increasing as k → ∞. Set τ ∞ = lim k→∞ τ k , hence τ ∞ ≤ τ e a.s.. To complete the proof, we only need to show that τ ∞ = ∞ a.s.. If this statement is false, there is a pair of constants T > 0 and ε ∈ (0, 1) such that P{τ k ≤ T } ≥ ε for any k ≥ k 0 .
We define a C 2 −function V : R 2 + → R + as follows: Applying Itô formula yields where K is a positive number, the remainder of the proof follows that in Li et al. [33], here, we omit it. The proof is complete.

Asymptotic Property
Due to the interference of stochastic noise, system (1.3) does not exist any equilibriums. However, we are interested in whether system (1.3) has stability and what is the influence of white noise. In this subsection, we discuss the asymptotic behaviors of system (1.3) around the predator-free equilibrium and the positive equilibrium of the corresponding deterministic system (1.2). Next, we study the asymptotic property of system (1.3) around the predator-free equilibrium E 1 = ( r−d 1 a , 0) when r − d 1 > 0 and the positive equilibrium, respectvely. Theorem 3.2. Let (x(t), y(t)) be a solution of model (1.3) with initial value (1.4). If the condition r − d 1 > σ 2 1 holds, then lim sup Proof. Define the function where ω i (i=1,2,3) are positive constants to be determined later.
Applying Itô formula yields: For where Using the inequality (a + b) 2 ≤ 2a 2 + 2b 2 for all a, b ∈ R, we have For where For where For where Choosing ω 1 = c 2x 0 , ω 2 = c 2 , ω 3 = cpx 0 d 2 (1+c 1 k) , one can see that Integrating both sides of Eq (3.1) from 0 to t and taking the expectation yield Divide both sides by t and take the limit superior, and then we obtain lim sup This completes the proof.

Proof.
Recalling that E * = (x * , y * ) is the positive equilibrium of system (1.2), we have An application of Itô formula implies where Using the inequality (a + b) 2 ≤ 2a 2 + 2b 2 for all a, b ∈ R, we have where where the inequality a(x − r a ) 2 ≥ 0 is applied. Define Define V 5 = y − y * − y * ln y y * , Moreover, define Additionally, define Construct the function where ω i (i = 1, · · · , 4) are positive constants to be determined as follows.
Therefore, we have Choosing Integrating both sides from 0 to t and taking the expectation yield Divide both sides by t and take the limit superior, and then we obtain lim sup where This completes the proof.
Remark 3.2. From Theorem 3.3, one can see that the solution of system (1.3) will be swinging near the positive equilibrium under certain conditions. In addiction, there is a positive correlation between the swing intensity and the intensity of environmental noise.
In order to show the role of τ on results more intuitively, we choose θ = 0.01 and k = 10, and consider the set of values 0. 1,2,8,15, and 30 (as demonstrated in Figure 5). The results indicate that when τ > τ 0 = 5.7998, the solution of system (1.2) shows more oscillating behaviors and the positive equilibrium is unstable. However, once the parameter is reduced to less than the threshold, the system solution will quickly converge to the equilibrium point, which verifies Theorem 2.2-(II)-(ii)) very well.
To explore the impact of self-defence on the dynamics of system (1.2) and to compare with the results obtained in [25], we choose the same parameter values as in [25]: h = 0, k = 1, d 2 = 0.05, c = 0.4, c 1 = 1, c 2 = 1 and τ = 120, and consider the set of values 0, 0.01, 0.03, 0.036, and 0.04 (as demonstrated in Figure 6). The curve of θ = 0 is in line with that in [25] and a periodic solution occurs. While as θ increases, the oscillation amplitude of both the prey and the predator gradually decreases and the equilibrium tends to stable eventually. Hence the level of self-defence acts as a stabilizing factor in system (1.2) (Remark 2.2).   Furthermore, Theorem 2.2 demonstrates that when Φ 1 > 0 and J 11 > hy * , the positive equilibrium E * of system (1.2) is unstable for all τ > 0. We set a = 0.05, q = 0.8, θ = 0.01 and k = 2 to ensure that the prerequisites are established and the numeric solutions are demonstrated in Figure 7. In such a case, the solutions are periodic and as τ increases, the amplitude of the oscillation (periodic solution) increases.  For the stochastic model, we shall adopt Milstein's Higher Order Method mentioned in [29] to verify our theoretical results. The corresponding discretization equations are where the time increment ∆t > 0, η is the integer part of τ ∆t − 1, ξ k and k (k = 1, 2, ..., n) are independent Gaussian random variables that follow the standard normal distribution N(0, 1).
In Figure 8, we consider the initial value (x(0), y(0)) = (0.2, 0.1) with related parameter values: θ = 0.01, k = 30, σ i = 0.01(i = 1, 2) and other parameters are given in Eq (4.1). In this case, r − d 1 = 0.02 > σ 2 1 = 1 × 10 −4 , which means the condition in Theorem 3.2 is satisfied. Hence, all positive solutions of system (1.3) will fluctuate around the predator-free equilibrium. As clearly shown in Figure 8, the curve of the stochastic model goes around E 1 = (2, 0), which supports Theorem 3.2. Besides, we take the parameters in Figure 1(b) and σ i = 0.01(i = 1, 2). In this situation, the conditions in Theorem 3.3 are satisfied:   Figure 11. Trajectories of stochastic system (1.3) and its corresponding deterministic system (1.2) with σ 2 = 0.01, 0.05, 0.2 and σ 1 = 0.01. keep increasing, the population affected will go to extinction (such as σ i = 0.2 (i = 1, 2)), which can be observed clearly from Figures 10(a) and Figure 11(b). Notice that Theorem 2.2 claims the coexistence of predator and prey. This discrepancy highlights the impact of stochastic fluctuations on the predator-prey model. What's more, the numerical simulations also verify a basic fact that the survival of predator is based on the survival of the prey while the prey can persist in the absence of the predator. Ultimately, we focus on the impacts of the response level k and the self-defense level θ on the prey, respectively. As stated at the beginning of the paper, we firmly believe that the scale of prey must have a significant relationship with anti-predation response and self-defence. Here, under the condition that σ i = 0.01 (i = 1, 2) and other factors are given in Eq (4.1), we choose different response level k = 1, 10, 20 (with θ = 0.01 fixed) and self-defense level θ = 0.001, 0.01, 0.05 (with k = 20 fixed). As is clearly demonstrated in Figure 12, both the response level and the self-defense level exert farreaching influence on the scale of prey. With the increase of the self-defense ability and anti-predator consciousness, the population of prey will gradually increase and trend to the predator-free equilibrium, which has two implications: (i) High-level self-defense will prevent the predator from being preyed on by the predator, otherwise the predator will eventually become extinct; (ii) Compared with the negative effects, the benefits of high-level response are more significant. One of the possible reasons is that the increase in escaped prey compensates for the decline in newborns. These are consistent with what we know: predators prefer safer and more efficient predation.

Conclusions
This paper is mainly related to a deterministic predator-prey model with fear factor and self-defence, considering the existence of digestion delay, and its corresponding stochastic model. Rather than simply considering the negative impact of the fear effect, we take the benefits of anti-predation behavior and resistance into consideration. Both avoidance and resistance of prey along with digestion delay all affect the stability of the equilibrium point of the deterministic system and may lead to Hopf bifurcation. In addition, as researchers have not paid much attention to the corresponding stochastic model, we mainly investigate the asymptotic behavior around the solution of the deterministic model. In contrast with existing results, our research shows that the level of fear alone is not enough to completely determine the dynamics of the predator-prey model. The level of self-defence acts as a stabilizing factor, as well. The results obtained can be applied to keep the coexistence of biological populations and maintain ecological balance. Furthermore, numerical results support and expand our analyses, indicating that: (i) The introduction of disturbance of noise may drift the coexistence equilibrium to the predator-free equilibrium; (ii) Small noise will not change the stability of the equilibrium of system (1.3) and the amplitude of the population is positively correlated with the intensity of withe noise; (iii) Large environmental fluctuations will lead to the extinction of predator and prey; (iv) The scale of prey has a significant relationship with both anti-predation behavior and self-defence, so is predator.
In addition to the conclusions that have been reached, there are some possible extensions of our work. As mentioned before, we are of great interest in the impact of self-defence on the prey population. A replacement of 1 1+c 1 k (respectively, r 1+c 2 ky(t) ) by g(k, θ) (respectively, f (k, θ, y(t))) which is decreasing in both k and θ may be a meaningful attempt. Furthermore, as pointed out in Xiao and Chen [34], stage structure has a significant influence in the stability of periodic orbits. The probability of self-defense is different when the prey encounters immature or mature predators. Moreover, Lan et al. [35] proposed a model with psychological effects and impulsive toxicants in polluted environments. These models are more complex but realistic. We can consider models with factors mentioned above for further investigations.