Effect of fear on two predator-one prey model in deterministic and fluctuating environment

Recent ecological studies on predator-prey interactions has concentrated on determining the impacts of antipredator behavior due to fear of predators. These studies are mainly confined into one predator-one prey system. But in case of multiple predator attack on single prey species, fear mechanism is still unknown. The combined impact of multiple predator often cannot be anticipated from their independent effects. So coexistence of multiple predators and prey’s fitness becomes an important issue from an ecological point of view. Based on the above observations, we proposed and analyzed a model consisting of two competing predator sharing a common prey where prey’s reproduction rate is affected due to fear generated by the predators. We first study the boundedness, uniform persistence, stability and Hopf bifurcation of the deterministic model. Thereafter, we have investigated the existence and uniqueness of the global positive solution, boundedness, asymptotic stability of the stochastic model.  Numerical examples are provided to support our obtained  results.


Introduction
Predator-prey interactions have been studying over a long period of time through a different perspective. Mathematical models play a crucial role to explore the dynamics of such interactions. Recently, the term ecology of fear creates a major interest to the researchers to include this fact in mathematical model describing the predator-prey interactions. Pioneer work of Zanette et al. [36] inspired many scientists to study the fear effect in predator-prey model [18,19,20,21,22,27,32,33,34,35,36]. Qi and Meng [24] investigated the threshold behaviour of a stochastic predator-prey system with prey refuge, fear effect and non-constant mortality rate. Sarkar and Khajanchi [26] studied the impact of fear on the growth of prey in a predator-prey model without environmental noise. They discussed permanence, Hopf bifurcation and stability of the coexistence equilibrium point. Xia and Yuan [34] analysed a stochastic predator-prey model with prey refuge and fear effect. They obtained stationary distribution and performed survival analysis of the model. Liu and Jiang [12] addressed the influence of fear factor on the dynamics of stochastic predator-prey model. They derived sufficient criteria for the existence and uniqueness of an ergodic stationary distribution of positive solution to the system. These studies mainly demonstrated the two dimensional predator-prey models which cannot give real picture and lots of complex behaviour can be observed in three or more interacting population system [17,30]. Smith [29] showed the existence of Hopf bifurcation in a two-predator-one prey competition model. Savitri et al. [28] studied the dynamical behaviour of modified Leslie-Gower one prey-two predator with competition and found permanence, global stability and bifurcation. Farkas and Freedman [5] derived stability condition in two-predator and one-prey system. More studies on two predators and one prey system can be found in [1,4,10]. The above studies did not consider the fear effect for obtaining coexistence results.
In view of above, it is reasonable to develop a model that predicts how prey respond to multiple predators when the fear factor comes into the scenario. By formulating deterministic models in ecology, it is often difficult to anticipate the future behaviour of the system properly. As deterministic models ignore the effect of random environmental fluctuation, the actual behaviour of a natural system cannot be well understood, even small noise can suppresses explosions in population dynamics [15]. It is well known that the natural growth of population is always affected by random perturbations and hence this should be incorporated into deterministic model. In [11], the authors demonstrated the impacts of white noises on the persistence and extinction of a one prey-two predator system and also obtained the condition for the existence and uniqueness of an ergodic stationary distribution of the positive solutions. The effect of the white noise on the prey-predator system are described in [9,14,23]. There are several approaches to include random effects in the model, both from ecological and mathematical point of view [8]. In this work, we follow the approach discussed in [2].
The main target of this paper is to investigate the impact of fear on prey population when two predators are competing for the same prey species in the deterministic and stochastic environment.
The paper is structured as follows. In the next section, we propose our model. The biological validity of the model is examined in Section 3. We demonstrate some dynamical properties of the model (2.2) in Section 4. The mathematical model with environmental noise is presented and analysed in Section 5. Computer simulations are given in Section 6. A brief discussion concludes the paper in Section 7.

Deterministic model
In [36], the authors investigated the following dynamical behaviour of the following predator-prey model employing the cost of fear into prey reproduction rate in the presence of predator.
where the variables x and y denote the densities of prey and predator population at time t respectively. All the parameters are positive. r is the birth rate, α is the natural death rate of prey, β denotes the decay rate due to intra-species competition, p represents the per capita predator consumption rate, m is the constant handling time for each prey captured, b is the conversion rate of preys biomass to predators biomass. k is the level of fear and d is the mortality rate of predator. In above study of prey antipredator strategies concern responses to single predator. But in real ecosystem, complex indirect interactions are common and a pair of species interaction cannot describe the patterns observed [17]. In prey-predator interactions, antipredator behavior strategy of prey can change the fitness and coexistence of predators [16]. So we have incorporated another predator z in system (2.1) and then system (2.1) becomes (2.2) Here k and l denote the level of fear due to the predator y and z respectively. ρ and σ are interspecific competition coefficient between the two predator species. n is the constant handling time for each prey captured. δ and µ are the intraspecies competition among the predator species. h represents the death rate of predator z. Specific example illustrates the above situation. The fresh water snail Physella gyrina faces two types of predators that exhibits different behavioral responses. The snail take shelter in a covered area when pumpkinseed sunfish comes in the water column. But another predator crayfish waits under the rocks in search of food. Then snail escapes from this place and moved to the water surface [31]. Furthermore, in grassland system, grasshopper Melanoplus femurrubrum are attacked by the different types of spider namely Pisaurina mira, Rabidosa rabida and Phidippus rimator. The grasshopper change their habitat from grasses to herbs in the presence of first two types of predators.
But it faces problems when third type of predator is present as it is available everywhere.

Positivity and boundedness of solutions
In this section, we first examine positivity and boundedness of solutions of system (2.2) which concern the validity and well behaved nature of the model. We first check positivity.
Proof. The positivity of x(t), y(t), z(t) can be verified by the equations with x 0 , y 0 , z 0 > 0. As x(t) > 0 for all t > 0. The same argument is valid for component y(t) and z(t).
Hence the interior of R 3 + , is an invariant set of system (2.2).
Proof. Let us consider the function The time derivative along a solution of (2.2) is Thus the the following inequality is satisfied. (3.1) Note that µ = min{d, h} is no larger than d and h. Thus (3.1) can further lead to By using the comparision theorem [3], we get Taking limit when t → ∞, we have Hence system (2.2) is bounded.

Existence of equilibria and stability analysis
Evidently, system (2.2) has five non-negative equilibrium points. The population free equilibrium point E 0 = (0, 0, 0). The predator free equilibrium point E 1 = ( r−α β , 0, 0). Here E 0 always exists and E 1 exists if r > α. The predator free equilibrium point E 2 = (x,ȳ, 0). The equilibrium point E 2 can be found in xy-plane provided it satisfies the following equations From Eqn. (4.2), we find the value of y as Now using the value of y in Eqn. (4.3), we get the following equation in x, These conditions imply that Eqn. (4.4) has a positive root. Similarly, we can show that The local stability properties of the equilibrium points can be examined through the Jacobian matrix around the each equilibrium point. .
It is very difficult to determine the coordinates of the coexistence equilibrium point. The existence of the coexistence equilibrium point can be obtained by showing uniform persistence of the system which in turn implies the existence [6,7].
Persistence. If there exists a compact set E ⊂ B = {(x, y, z) : x > 0, y > 0, z > 0} where all solutions of system (2.2) ultimately lie, the system is said to be persistent. The following theorem establishes uniform persistence of system (2.2).
Theorem 4.1. Suppose E 1 , E 2 and E 3 exist. If there are no limit cycles in the x − y and x − z plane respectively. If , h + σȳ < cqx 1 + nx and d + ρẑ < bpx 1 + mx , then system (2.2) is uniformly persistent.
Proof. Proceeding along the lines in [19], we can prove the theorem and is omitted here.
2) is uniformly persistent under the assumptions of the Theorem 4.1. Now we conclude that there exist a time T such that x(t), y(t), z(t) > L, 0 < L < 1, for t > T.

Remark 2.
If there are finite number of limit cycles , then persistence condition becomes for each limit cycle (φ(t),ψ(t)) in the xy-plane and each limit cycle (φ(t),ψ(t)) in the xz-plane where τ is the appropriate period.
The above theorem implies uniform persistence of system (2.2) under certain restrictions. Further, it is shown that in [7], uniform persistence ensures the existence of a positive equilibrium point. Hence E 4 = (x * , y * , z * ) exists provided the conditions of Theorem 4.1 are fulfilled.
To discuss the stability of E 4 , we find out the characteristic equation around E 4 which is given by where From Routh Hurwitz criterion, we can say that E 4 = (x * , y * , z * ) is locally asymptotically stable if the following conditions are satisfied. a 1 > 0, a 2 > 0 and a 1 a 2 > a 3 .
and detD > 0 where D is defined in the proof. Then E 4 is globally asymptotically stable.
Bifurcation study. Next, we present a bifurcation result. For convenience of presentation, we define f (k) = a 1 (k)a 2 (k) − a 3 (k).
then the positive equilibrium point E 4 is unstable when k < k * but is stable when k > k * and a Hopf bifurcation of periodic solutions emerges at k = k * .
Proof. Using the method in [25], we observe that f (k) is monotonic increasing in the neighborhood of k = k * . As a i (k) > 0, i = 1, 2, 3, f (k) < 0 when k < k * then E 4 becomes unstable. Again,f (k) > 0 when k > k * then E 4 becomes stable. Applying a result in [13], we can establish Hopf bifurcation. Similarly, one can show bifurcation taking l as a bifurcation parameter.

Stochastic model
In our system (2.2), we consider the random effects to address some biological issues in future. This leads to a stochastic generalization of our system. The random perturbations are dependent on the white noises, which are directly proportional to the distances of x(t), y(t), z(t) from the equilibrium values of x * , y * , z * respectively. So our system (2.2) becomes where σ i > 0, i = 1, 2, 3 and represent the intensity of environmental fluctuations and B i (t), i = 1, 2, 3 are independent Brownian motions.
Throughout this analysis, we take (Ω, F, P ) as a complete probability space with a filtration {F t } t∈R+ satisfying the conventional condition, namely right continuity and increasing, whereas F 0 consists of all P -void sets.
Proof. Using local Lipschitz property, one can show existence of unique solution on [0, τ e ] where τ e represents the explosion time. We now prove that this solution is global i.e.,τ e = ∞ almost surely.
Let a 0 > 0 be large enough for x(0), y(0) and z(0) lying within the interval [ 1 a0 , a 0 ]. Given any integer a ≥ a 0 , define the stopping time by where we take inf Ø = ∞ where Ø denotes the empty set. Since τ a is nondecreasing as a → ∞, then we have τ ∞ = lim a→∞ τ a . Then τ ∞ ≤ τ e a.s. We have to establish that τ ∞ = ∞ a.s. If this is not true, ∃T > 0 and ∈ (0, 1) such that P {τ a ≤ T } > . Consequently, there is an integer a 1 ≥ a 0 such that Consider a function V :

FEAR EFFECT IN A PREDATOR-PREY MODEL IN DETERMINISTIC AND FLUCTUATING ENVIRONMENT 9
Clearly, the above function is non-negative. If U (t) ∈ R 3 + , by applying Itô formula, we obtain where M is a positive constant. Let τ a ∧ T = min{τ a , T }. Then Thus it follows that Take Ω a = {τ a ≤ T } for a ≥ a 1 and from (5.2), we get P (Ω a ) ≥ . Observe that for each ξ ∈ Ω a , there is x(τ a , ξ), y(τ a , ξ), z(τ a , ξ) equals either a or 1 a , hence V (x(τ a , ξ), y(τ a , ξ), z(τ a , ξ)) is no less than So we get where I Ωa(ξ) is the indicator function of Ω a . Taking a → ∞, so ∞ = V (x(0), y(0), z(0)) + M T < ∞, a contradiction. Hence we have, τ e = ∞ almost surely.
Now we show stochastic ultimate boundedness.

This implies that lim
So we can find a positive constant π 1 such that lim t→∞ supE| U (t)| < π 1 .
We now establish global stability of the positive steady state E 4 in stochastic environment. Lyapunov functional will be constructed to prove global stability. In the remark after Theorem 4.1, we note that x(t), y(t), z(t) > L.
where the matrix H is given in the proof. Then E 4 is globally asymptotically stable with probability one.
Proof. Consider the same Lyapunov function used in Theorem 4.2.
Using Itô formula, we have where X T = {|x − x * |, |y − y * |, |z − z * |} and the entries of the matrix H = [h ij ] 3×3 are given by: Hence, one can verify that under the conditions of the theorem, H is positive definite, leading to dV (x, y, z) < 0. This conclude globally asymptotical stability of E 4 , completing the proof.

Verification of analytical results
In this section, we present some numerical computations to justify our analytical results. We mainly focus here to show the effect of fear and intraspecific competition rate among the prey species on the  3). Comparing Figs. 1 and 3, we observe that the value of x * , y * , z * of E 4 decreases as long as the level of fear increases. The corresponding stochastic solutions with σ 1 = 12, σ 2 = 9, σ 3 = 10 are shown in Fig. 4. Now we change the parameter β in (6.1) and assuming the value as 5, we note that the coexistence equilibrium point E 4 (0.2939, 0.7908, 0.6600) is globally stable (see Fig. 5). The stochastic solutions with σ 1 = 12, σ = 9, σ 3 = 10 are shown in Fig. 6. In this case, the conditions of Theorem 6 are not satisfied, noises make the solutions fluctuating in random manner. In Fig. 7, we choose σ 1 = 0.1, σ 2 = 0.2, σ 3 = 0.3, the other parameters are same as in Fig. 6, the conditions of Theorem 6 are satisfied and hence the coexistence equilibrium point E 4 (0.2939, 0.7908, 0.6600) is globally asymptotically stable with probability one.

Discussion
The impact of fear in predator-prey model demands a lot of attention. In case multiple predators attacking a single prey species, antipredator strategy is still unknown. With this fact in mind, we have proposed a one prey-two predator species incorporating fear factor in the preys reproduction rate. Here the two predators are competing for the same resource. For the biological validity of the model, we have examined positivity and boundedness of the system. The model analysis shows that the level of fear has a major influences on the dynamical behavior of the system. It is observed that increase amount of fear can stabilize the system. The local stability, persistence, bifurcation and global stability are investigated in deterministic system.   We also explore the impact of environmental noises on the system. By applying Itô formula, we are able to determine the existence of a unique global solution for any positive initial point. Boundedness is demonstrated. Global stability condition is derived by constructing a Lyapunov function. It is noted that when intensities are small, the system remains stable. These results are biologically important as it concerns the coexistence of all the species.
But the major difference between our work and the other recent works is the inclusion of one more prey in predator-prey system. In this case, predator may also modify prey population dynamics by changing the preys interactions with the other species. Thus the impacts on species coexistence depend not only by the predation but also by the interspecific competition between the predator species.  To our understanding, the dynamical study of deterministic and stochastic version of two predator competing for a single resource with fear has not done yet. Permanence study remains for future work.
Conflict of interest. The author declares that there is no conflict of interest in publishing this paper.
Funding. This research received no specific grant from any funding agency in the public, commercial or not-for-profit sectors.