DYNAMICS OF A STOCHASTIC DELAYED HARRISON-TYPE PREDATION MODEL: EFFECTS OF DELAY AND STOCHASTIC COMPONENTS

. This paper investigates the complex dynamics of a Harrison-type predator-prey model that incorporating: (1) A constant time delay in the func- tional response term of the predator growth equation; and (2) environmental noise in both prey and predator equations. We provide the rigorous results of our model including the dynamical behaviors of a positive solution and Hopf bifurcation. We also perform numerical simulations on the eﬀects of delay or/and noise when the corresponding ODE model has an interior solution. Our theoretical and numerical results show that delay can either remain sta- bility or destabilize the model; large noise could destabilize the model; and the combination of delay and noise could intensify the periodic instability of the model. Our results may provide us useful biological insights into population managements for prey-predator interaction models.


1.
Introduction. Dynamics of prey-predator systems has been an important research theme in ecology. It is well known that prey-predator interaction is not only one of the basic interspecies relations for ecological and social systems, but also the basic block of more complicated food chains, food webs, and biochemical network structures [21]. In the biological world, there are many examples of precesses that involve significant delays that cannot be ignored [13]. Examples include the time between fertilization and birth in the case of sexual reproduction, the time between initiation of cellular division and effective division in the case of mitosis, and the time required for digestion in the case of consumption of nutrient and its conversion to viable biomass [1]. Time delay plays a vital role in population dynamics of prey-predator, which has been recognized to contribute critically to the stable or unstable outcomes of prey population due to predation [10]. Population systems are exposed to the impact of a large number of random unpredictable factors. Then, environmental noise component is another important factor that affects prey-predator dynamics and, it related to climate, geographical distribution, geological features, disaster, human intervention, and other environmental factors [14]. Thus, a realistic prey-predator model should include both time delays and environmental noise.
For example, it is reasonable to assume that the death of prey is instantaneous when attacked by their predator but its contribution to the growth of predator population must be delayed by some time delay [25]. Kuang [10] mentioned that animals must take time to digest their food before further activities and responses taking place. Samanta [26] argued that the effect of time delay due to the time required in going from the egg stage to the adult stage, gestation period etc, has to be taken into account. Three different ways of incorporating a constant time delay into the prey-predator models were presented by Martin and Ruan [18]. In the first case, a time delay appears in the prey specific growth term, which is based on the assumption that in the absence of predators prey satisfies the Hutchinson's equation x (t) = x(t)(1−x(t−τ )), that is also called the delayed logistic equation [19]. In the second case, there is a delay in the predator response term of the predator equation, that is, y (t) = y(t)(−d + p(x(t − τ ))), which can be regarded as a gestation period or reaction time of predators. The third one is a time delay in the interaction term of the predator equation, that is, y (t) = −dy(t)+y(t−τ )p(x(t−τ )), which assumes that the change rate of predators depends on the number of preys and of predators present at some previous time. In this work, we incorporate the second type of delay in the functional response term of the predator growth equation.
Since the nature is full of uncertainty and random phenomena, the natural growth of species often does not follow strictly deterministic laws but rather oscillate randomly about some average, so that the population density never attains a fixed value with the advancement of time rather exhibits continuous oscillation around some average values [14]. The basic mechanism and factors of population growth like resources and vital rates-birth, death, immigration and emigration, change non-deterministically due to continuous fluctuations in the environment (e.g. variation in intensity of sunlight, temperature, water level, etc.) [20]. It is necessary and important to consider the corresponding stochastic population model, i.e., the effects of environmental noise, which undeniably arise from either environmental variability or internal species. For example, in [3], the authors extended a classical epidemic model from a deterministic framework to a stochastic differential equation one through introducing random fluctuations. Among the various ways of constructing a stochastic model systems for a given deterministic system, Cai et al. [4] propose a stochastic version of the epidemic model with nonlinear incident rate. Mao et al. [16] revealed that given population systems are subject to environmental noise, it can suppress a potential population explosion. That is to say, different structures of environmental noise may have different effects on the population systems. In our work, we would include noise processes in both prey and predator.
As we mentioned earlier, a realistic prey-predator model should include both time delays and environmental noise. There is a considerate amount of work done by many scholars, for instance [17,2,22,11,23,27,6,12,28,24,5]. Of them, Mao et al. considered the effects of environmental noise on the delay Lotka-Volterra in [17]. Vasilova [27] studied a stochastic Gilpin-Ayala predator-prey model with time-dependent delay. He established sufficient conditions for the existence of a global positive solution of the considered system, and proved stochastically ultimate boundedness, the long-time behavior of trajectories and extinction of species. In [6], Han et al. investigated two-species Lotka-Volterra delayed stochastic predator-prey system, with and without pollution, respectively. They revealed that there exists a unique nonnegative solution in each system that is permanent in time average under certain conditions. The convergence of the distributions of the solutions of a stochastic two-predator one-prey model with time delay was considered by Liu, Bai and Jin et al. [12]. Wang et al. [28] investigated the dynamics of a stochastic FIV model with seasonality analytically and numerically and, established sufficient criteria for extinction and weak persistence of the disease in the mean.
Motivated by the works of [17,27,12], we propose to study a delayed Harrisontype predator-prey model with noise terms that is lack in literature. We aim to explore and address the following questions through our analytic and numerical results via comparisons to the corresponding ODE model: 1. What are the dynamical effects of time delay in the functional response term of the predator growth equation on the Harrison-type predator-prey model? 2. What are the joint effects of time delay and noise on the Harrison-type predator-prey model? The remaining of this paper is organized as follows. In Section 2, we derive our model and provide a summary on the dynamics of the corresponding ODE model. And we analyze the stability of a positive equilibrium under the effect of time delay or noise or both of them, and provide sufficient conditions that guarantee the existence of Hopf bifurcation at a positive equilibrium. In Section 3, numerical simulations are given to verify the theoretic analysis. In Section 4, the paper provides a brief conclusion and remarks. All proofs of our theoretical results are presented in the last section.

2.
Model derivations and the related dynamics. The general predator-prey model can be expressed by the model of nonlinear ordinary differential equations [7]: where N = N (t) and P = P (t) denote the prey and predator population at time t, respectively. The function u(N ) represents the growth rate of the prey in the absence of predation and is given by the traditional logistic form u(N ) = rN (1− N K ). The product f (N )g(P ) gives the rate at which prey is consumed, while αf (N )g(P ) describes per capita production rate of the predator. When these functions are defined by f (N ) = cN, g(P ) = P mP + 1 , c, m > 0, and the function v(P ) is given by v(P ) = dP where d > 0 is the natural death rate of the predator population. Then, we make a change of variables and the according Harrison-type functional response predator-prey model is described by the following form: where b = αc stands for conversion rate from prey to predator. Model (2) and its many kinds of extended forms have been widely studied in ecological literature [8,31,29,30,9]. In [29], Wang and Cai presented a theoretical analysis of processes of pattern formation that involves organisms distribution and their interaction of spatially distributed population with cross-diffusion in a Harrison-type predator-prey model. Wang et al. [30] investigated the complex dynamics induced by Allee effect in the reaction-diffusion predator-prey model. Jin [9] studied the solution moment stability for a Harrison-type predator-prey model with parametric dichotomous noises.
In this paper, we will study the effects of fluctuating environment on the dynamical behaviors of the Harrison-type predation model (2) with time delay τ (> 0) which is introduced into the functional response term of the predator growth equation, and incorporate white noise in each equations of the model. Moreover, we assume that fluctuations in the environment show themselves as fluctuations in the growth rate of prey N and in the death rate of predator P , then the corresponding stochastic delayed Harrison-type predator-prey model can be described as: where σ 1 , σ 2 (> 0) represent the noise intensity. B(t) is a standard one-dimensional independent Wiener process defined on a complete probability space Ω, F, {F t } t≥0 , P with filtration {F t } t≥0 satisfying the usual conditions (i.e., it is right continuous and increasing while F 0 contains all P-null sets). We denote by R n + the positive cone in R n , that is R n + = {x ∈ R n : x i > 0 for all 1 ≤ i ≤ n}, and denote byR n + the nonnegative cone in R n , that isR n + = {x ∈ R n : x i ≥ 0 for all 1 ≤ i ≤ n}. Let C([−τ, 0]; R n + ) denote the family of continuous functions from [−τ, 0] to R n + . If x ∈ R n , its norm is denoted by |x| = (

2.1.
With only delay in the predator growth term. Before further studying the dynamics of our full model (3), we provide the following results regarding its corresponding ODE model (2) and the following delayed model Model (3) has the trivial constant steady state E 0 = (0, 0) and the boundary equilibrium E 1 = (1, 0), and the other nonnegative constant steady states of (3) are the nonnegative constant steady states of the corresponding ODE model (2) in the absence of delay and noise. Let is always negative.
Theorem 2.2 (Dynamics of the ODE model). The ODE model (2) always has two boundary equilibria: the extinction equilibrium E 0 = (0, 0) which is always saddle; and the prey-only-equilibrium E 1 = (1, 0) which is locally stable if d ≥ b holds while it is a saddle otherwise. As a summary, the conditions of the number of equilibria of the ODE model (2) are provided as follows: 1. If d ≥ b holds, then the ODE model (2) has no interior equilibrium with two boundary equilibria: E 0 = (0, 0) being a saddle and E 1 = (1, 0) being locally stable. 2. If d < b holds, then the ODE model (2) has a unique interior equilibrium E * = (N * , P * ) which is locally stable. Moreover, it is globally asymptotically stable in domain R 2 + . The existence and stability conditions of the equilibrium of the ODE model (2) is listed in Table 1.

Equilibrium Existence Condition Stability Condition
(0, 0) Always exists Always saddle Always sink Table 1. The existence and stability of equilibria for model (2) where The proof of Theorem 2.2 is rather standard, and detailed investigations relating to the ODE model may be found in Ref. [29], hence is omitted. Remark 1. We would like to point out that Wang and Cai [29] provided similar results of the ODE model (2). For the ODE model (2), some preliminary results are presented, including dissipativeness, boundedness, permanence of the solutions, and the equilibria stability analysis of the model. Model (2) has two boundary equilibria E 0 = (0, 0), E 1 = (1, 0) and a unique positive equilibrium E * = (N * , P * ). From Theorem 2.2, for the equilibrium E 1 = (1, 0), ecologically it means that if the natural death rate of predator d is more than the conversion rate from prey to predator b, then it leads to the extinction of predator population. For the equilibrium E * = (N * , P * ), ecologically it implies that if the natural death rate of predator d is less than the conversion rate from prey to predator b, then model (2) will be locally asymptotically stable around E * . We also know that E * = (N * , P * ) of the ODE model (2) is globally asymptotically stable if it attracts all positive solutions of the model. That is, the boundedness of the solution together with the saddle nature of boundary equilibrium point E 1 = (1, 0) and asymptotical stability of E * = (N * , P * ) leads to the conclusion that all the trajectories will approach E * with increasing time t. Hence E * is a global attractor and the model is globally asymptotically stable. When we take the parameters as c = 0.9, b = 0.7, d = 0.3, m = 0.1 for model (2), the positive equilibrium E * = (0.46, 0.64) exists which is spiral sink and shown in Fig.1. For model (3), when σ 1 = σ 2 = 0, we study the DDE model (4) Remark 2. In the case of τ = 0, i.e., model (4) without time-delay, or model (2), the dissipation and boundedness of the solution are similar to the results in Theorem 2.3. That's to say, the solutions of the model with or without time-delay effect are always dissipative and uniformly bounded. Moreover, when the natural death rate of predator d is more than the conversion rate from prey to predator b, the boundary equilibrium E 1 = (1, 0) of model (4) is stable.
and less than . LetN = N − N * ,P = P − P * and dropping the bars for simplicity if notations, then the DDE model (4) can be transformed into the following: which can be rewritten as an abstract differential equation in the phase space where and here, Referring to [10], we can know that given initial data (N (ϑ), P (ϑ)) ∈ C([−τ, 0]; R 2 + ), model (4) will have a unique local or global solution dependent on the parameters where q i , i = 1, 2, 3, 4 are defined by (10). After a little algebraic calculation, then we get a unique positive root: and the corresponding critical value of delay τ = τ k at which E * = (N * , P * ) loses stability is given by Hence, we have the following theorem.

Remark 4.
In the absence of delay τ = 0, when the natural death rate of predator d is less than the conversion rate from prey to predator b, then the ODE model (2) remain locally asymptotically stable around E * = (N * , P * ). In the presence of delay τ = 0, when the conditions of Theorem 2.5 hold for the DDE model (4) implies that the interior equilibrium E * is conditionally stable. Based on the above analysis, we can see that if E * = (N * , P * ) of model (4) is stable for τ = 0, then there is a positive constant τ 0 which is the least positive value of τ k given by (13), such that for τ > τ 0 , E * becomes unstable. It is worthy to mention here that the amplitudes of both prey and predator population increase with the increase of time delay (see our simulations below for details).  Fig. 1, and the initial value is set to (N (0), P (0)) = (1.2, 1.0). According to Theorem 2.5, there is a critical value τ 0 = 3.46, Hopf bifurcation occurs and a bifurcation periodic solution exists at τ = τ 0 , which displays that both prey and predator population reach periodic oscillations around the equilibrium E * = (0.47, 0.64) in finite time. The equilibrium E * is asymptotically stable when τ = 2.8 < τ 0 (see Fig. 2(a)), that is to say, the delay model (4) is conditionally stable at the unique positive equilibrium E * . The equilibrium E * loses its stability and becomes unstable when τ > τ 0 (see Figs. 2(b) and 2(c)). Comparing Figs. 2(a), 2(b) with 2(c), one can see that a stable limit cycle arises, which is caused by the Hopf bifurcation. It implies that time delay can induce instability and oscillations via Hopf bifurcation in the model. We also realize that a large enough time delay can increase the amplitude of the oscillating orbit of species.

2.2.
With both delay and noise effect. In the following, we turn our focus on the dynamical behaviors of the SDDE model (3) with both time delay and noise effect. Denote Z(t) = N (t), P (t) and Z(t) = N 2 (t) + P 2 (t).
Theorem 2.6. For any given initial value The property of Theorem 2.6 makes us continue to discuss how the solution varies in R 2 + in more detail. And we present that the definition of stochastic ultimate boundedness [14] is one of the important topics in population dynamics and defined as follows.  Remark 6. The above analyses show that under certain conditions, the original ODE model (2) and the associated SDDE model (3) behave similarly in the sense that both have positive solutions which will not explode to infinity in a finite time, and actually, will be ultimately bounded. That is to say, under certain conditions the noise will not spoil these nice properties.
The following result gives the almost surely asymptotic property of model (3). 3. Numerical simulations. In this section, numerical simulation results of the SDDE model (3) for different values of τ and σ 1 , σ 2 will be mentioned in subsequent two cases, firstly for τ = 0 and secondly for τ = 0 by increasing the magnitude of noise σ 1 , σ 2 step by step. Other parameters are taken as c = 0.9, b = 0.7, d = 0.3, m = 0.1 (16) and the initial value is set to (N 0 , P 0 ) = (1.2, 1.0). Fig. 3 shows the time-series plots of model (3) only with different noise intensities σ 1 , σ 2 and τ = 0, other parametric values as (16) which are the same as in Fig. 1. In Fig. 3(a), without noise, i.e., σ 1 = σ 2 = 0, the interior equilibrium E * = (0.46, 0.64) of the ODE model (2) is asymptotically stable. Let σ 1 = 0, σ 2 = 0.12 in Fig. 3(b), we can see that the positive solution of model (3) is stochastic fluctuation. In Fig. 3(c), take σ 1 = 0.12, σ 2 = 0, the fluctuation of prey species N occurs, thus the predator species P will become instability. Choosing σ 1 = 1.3 and σ 2 = 1.5 in Fig. 3(d), it shows the stochastically fluctuating population distribution for both prey and predator species when the magnitude of environmental noises increase. Fig. 4 shows the evolution of both the prey and predator species of the SDDE model (3) with different noise intensities σ 1 , σ 2 and keeping τ = 2.8 < τ 0 = 3.46. In Fig. 4(a), let σ 1 = 0 and σ 2 = 0.12, it indicates that the prey and predator species fluctuate around the steady state value (0.46, 0.64) with small oscillation. In Fig. 4(b), take σ 1 = 0.12 and σ 2 = 0, one can see that the stochastic oscillation happens in both species with the weak fluctuation. In Fig. 4(c), increasing σ 1 , σ 2 to σ 1 = σ 2 = 0.5, it illustrates the prey and predator species fluctuation with strong fluctuation over time t and the amplitude of oscillation is amplified due to environmental noise increases.

Conclusions and remarks.
For the stochastic delayed Harrison-type predatorprey model (3), one of our interesting findings is that our proposed delayed predation model can exhibit Hopf bifurcation. And we prove that the positive equilibrium is asymptotically stable when the delay is less than a certain critical value and unstable when the delay is greater than the critical value. These findings have also been reported in the delayed Harrison-type predator-prey system with diffusive effects under Neumann boundary conditions in [31]. The work of [31] shows that the system has a Hopf bifurcation and, the authors analyze the direction of Hopf bifurcation and the stability of bifurcating periodic solution. In addition, we provide the conditions of dissipative, uniformly bounded, and permanent of our proposed delayed model.
We also provide additional results on the case when the Harrison-type predation model is added by both time delay and environmental noise under certain conditions. Our analysis and numerical simulations suggest that large delay could destabilize the model; and the combination of delay and noise could intensify the periodic instability of the SDDE model (3). Complicated analytical conditions for the above results depend upon certain parametric restrictions involving the parameters of model  whereas the corresponding model without delay is asymptotically stable at the interior equilibrium. Delay induced instability is shown in Figs. 2(b) and 2(c) and cause the populations to fluctuate. If the time delay τ is smaller than the critical value τ 0 , the positive equilibrium E * = (N * , P * ) of the DDE model (4) maintains its stability when E * is stable in the corresponding ODE model (2), see Fig. 2(a). If the time delay τ is larger than the critical value τ 0 , the time delay can induce a stable limit cycle generated through the Hopf bifurcation and the larger time delay can increase the amplitude of the oscillating orbits of two species. 2. The synergistic effects of delay and noise: We illustrate some relevant properties of the corresponding stochastic delayed predator-prey model (3) and reveal the effect of environmental noise on the model. The increase of the noise intensity has a drastic impact on the dynamical behavior of both species with or without the delay effect. Theorem 2.6 shows that under some conditions, the stochastically perturbed SDDE model (  predator species have positive solutions which will not explode to infinity in a finite time and, in fact, will be ultimately bounded. That is, under certain conditions environmental noise will not spoil these nice properties. Moreover, Theorem 2.10 shows model (3)  Proof. For all t ∈ [0, t 0 ] where t 0 ∈ (0, ∞), we need to show that N (t) > 0 and P (t) > 0. Suppose that is not true. Then, there exists T ∈ (0, t 0 ) such that for all t ∈ [0, T ], N (t) > 0 and P (t) > 0, and either N (T ) = 0 or P (T ) = 0. For all As (N, P ) is defined and continuous on [−τ, T ] and B(t) is a standard Brownian motion, there is a C ≥ 0 such that for all t ∈ [−τ, T ], , Proof. For our full model (3), when σ 1 = σ 2 = 0, it is the following delayed model with initial conditions N (ϑ) = ψ 1 (ϑ) ≥ 0, P (ϑ) = ψ 2 (ϑ) ≥ 0 for all ϑ ∈ [−τ, 0], where ψ 1 (ϑ), ψ 2 (ϑ) ∈ C([−τ, 0], R 2 + ). From the first equation of (19), we can get . A standard comparison theorem shows that lim sup t→∞ N (t) ≤ 1. As a result, for any > 0, there exists a T > 0, such that N (t) ≤ 1 + for t > T . Then from the second equation of (19), we define P on all interval [kτ, (k + 1)τ ] with k = 0, 1, 2 . . ., then one can see that P is bounded on [0, t 0 ] if t 0 < ∞. From the predator equation of model (19), we get hence, for t > τ , P (t) ≤ P (t − τ ) exp(bτ ), which is equivalent, for t > τ , to For any ε > 1, there exists positive T ε , such that for t > T ε , N (t) < ε. According to the above inequality, we obtain, for t > T ε + τ , which implies by the same arguments use for N that, Conclusion of this Theorem holds by letting ε → 1.
In addition, based on the above analysis, if d ≥ b, for t > T ε + τ , we have lim t→∞ P (t) = 0. Then, we will show that there is a B > 0, independent of the initial condition, such that min lim inf From model (4), for any ε > 1 and for large enough t, we get Based on the standard comparison arguments, it follows that Set holds, then C 1 > 0. From (21) and Theorem 2.3, for any ε > 1, there exists a positive constant T ε , such that for t > T ε , Thus, for t > T ε + τ , we obtain On the one hand, for t > T ε + τ , these inequalities lead to dP (t) dt > P (t) which involves, for t > T ε + τ , On the other hand, from (22) and (23), for t > T ε + τ , we obtain When ε → 1, we have Let be B = min{C 1 , C 2 } > 0. Therefore, model (4) is permanent.
Proof of Theorem 2.5.
Assume that for some τ > 0, iω (ω > 0) is a root of Eq.(11), then Separating the real and imaginary parts, we have thus, which leads to the following equation Solving for ω 2 in (27), we obtain Thus, Eq. (27) has a unique positive root and substituting this value in (25), then the corresponding critical value of delay τ = τ k at which E * = (N * , P * ) loses stability is given by Let λ(τ ) = α(τ ) + iω(τ ) be the root of Eq.(11) such that α(τ k ) = 0, ω(τ k ) = ω 0 . Substituting λ(τ ) into Eq.(11) and differentiating both sides of it with respect to τ , we have which, together with (25), (26) and (29), leads to sign Re dλ dτ This implies that all the roots cross the imaginary axis at iω from left to right as τ increases. Hence, the transversality condition is satisfied. Hopf bifurcation occurs at τ = τ k given by (30), which are Hopf bifurcation values of model (4). This completes the proof.
Proof of Theorem 2.6.
Proof. From the biological meaning, we only consider the positive solution to the SDDE model (3). Taking the change of variables as N (t) = e u(t) , P (t) = e v(t) and using Itô's formula, model (3) can be reformulated as follows: the coefficient of (33) satisfy the local Lipschitz condition, then for any given initial values u(ϑ) = ln N (ϑ), v(ϑ) = ln P (ϑ), ϑ ∈ [−τ, 0], there is a unique local solution where t e is explosion time. Referring to [14,15], to show this solution is global, we need to show t e = ∞ a.s. Let l 0 > 0 be sufficiently large for For each integer l ≥ l 0 , define the stopping time where we set inf ∅ = ∞ (∅ represents the empty set) in this paper. t l is increasing as l → ∞. Let t ∞ = lim l→∞ t l , then t ∞ ≤ t e a.s.
Proof of Theorem 2.9.
Proof. The solutions of model (3) will remain in R 2 + for all t ≥ 0 with probability one from Theorem 2.6. From Theorem 2.8, there exists C > 0 such that lim sup t→∞ E|(N (t), P (t))| θ ≤ C.
According to the Definition 2.7, the SDDE model (3) is stochastically ultimately bounded.
Proof of Theorem 2.10.