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-\tau))$, 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-\tau)))$, 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-\tau)p(x(t-\tau))$, 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 Harrison-type 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]:

$\left\{\begin{array}{l} \dfrac{\text{d}N}{\text{d}t} = u(N)-f(N)g(P), \\ \dfrac{\text{d}P}{\text{d}t} = \alpha f(N)g(P)-v(P), \end{array}\right.$

(1)

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-\frac{N}{K})$. The product $f(N)g(P)$ gives the rate at which prey is consumed, while $\alpha f(N)g(P)$ describes per capita production rate of the predator. When these functions are defined by

$f(N) = cN, \;\;\;g(P) = \frac{P}{mP+1}, \;\;\; c, m>0,$

and the function $v(P)$ is given by $v(P) = d P$ 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:

$\left\{\begin{array}{l} \dfrac{\text{d}N}{\text{d}t} = N(1-N) -\dfrac{cNP}{m P+1}, \\ \dfrac{\text{d}P}{\text{d}t} = P\big(-d+\dfrac{b N}{m P+1}\big), \end{array}\right.$

(2)

where $b = \alpha 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 $\tau$ ($>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:

$\left\{\begin{array}{l} \text{d}N(t) = N(t)\Big(1-N(t)-\dfrac{cP(t)}{mP(t)+1}\Big)\text{d}t +\sigma_1N(t)\text{d}B(t), \\ \text{d}P(t) = P(t)\Big(-d+\dfrac{b N(t-\tau)}{mP(t-\tau)+1}\Big)\text{d}t +\sigma_2P(t)\text{d}B(t), \end{array}\right.$

(3)

where $\sigma_1, \, \sigma_2\>(>0)$ represent the noise intensity. $B(t)$ is a standard one-dimensional independent Wiener process defined on a complete probability space $\big(\Omega, \mathcal{F}, \{\mathcal{F}_t\}_{t\geq0}$, $\mathbf{P}\big)$ with filtration $\{\mathcal{F}_t\}_{t\geq0}$ satisfying the usual conditions (i.e., it is right continuous and increasing while $\mathcal{F}_0$ contains all $\mathbf{P}$-null sets). We denote by $R_+^n$ the positive cone in $R^n$, that is $R_+^n = \{x\in R^n: x_i>0\>\text{for all}\> 1\leq i\leq n\}$, and denote by $\bar{R}_+^n$ the nonnegative cone in $R^n$, that is $\bar{R}_+^n = \{x\in R^n: x_i\geq 0\>\text{for all}\> 1\leq i\leq n\}$. Let $\mathcal{C}([-\tau, 0]; R_+^n)$ denote the family of continuous functions from $[-\tau, 0]$ to $R_+^n$. If $x\in R^n$, its norm is denoted by $|x| = (\sum_{i = 1}^nx_i^2)^{\frac{1}{2}}$ [17].

Since our proposed model (3) is the population model, we restrict the state space being the closed first quadrant in the $(N, P)$ plane $\bar{R}_+^2 = \big\{(N, P):\, N\geq0, \, P\geq0\big\}$.

Theorem 2.1. The positive quadrant $int(R_+^2)$ is invariant for model (3).

---

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

$\left\{\begin{array}{l} \dfrac{\text{d}N(t)}{\text{d}t} = N(t)(1-N(t))-\dfrac{cN(t)P(t)}{mP(t)+1}, \\ \dfrac{\text{d}P(t)}{\text{d}t} = P(t)\Big(-d+\dfrac{bN(t-\tau)}{mP(t-\tau)+1}\Big). \end{array}\right.$

(4)

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

$N_i = \dfrac{b(m-c)\pm\sqrt{4bcdm+b^2(m-c)^2}}{2bm}, \, P_i = \dfrac{bN_i-d}{dm}, i = 1, 2, N_1 < N_2,$

where equilibrium point $(N_1, P_1) = \big(\frac{b(m-c)-\sqrt{4bcdm+b^2(m-c)^2}}{2bm}, \frac{bN_1-d}{dm}\big)$ 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\geq 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\geq 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.

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 $\sigma_1 = \sigma_2 = 0$, we study the DDE model (4) with initial conditions $N(\vartheta) = \psi_1(\vartheta)\geq 0, \> P(\vartheta) = \psi_2(\vartheta)\geq 0$ for all $\vartheta\in[-\tau, 0]$, where $\psi_1(\vartheta), \psi_2(\vartheta)\in \mathcal{C}([-\tau, 0], R^2_+)$.

Theorem 2.3. For the DDE model (4),

(i) The nonnegative solution $(N(t), P(t))$ of the model satisfies $\limsup\limits_{t\rightarrow\infty}N(t)\leq 1$ and

$\limsup\limits_{t\rightarrow\infty}P(t)\leq\frac{(b-d)\exp(\tau)}{dm}$

for $\tau>0$ with $d < b$.

(ii) If the condition $d\geq b$ holds, then $\lim\limits_{t\rightarrow\infty}(N(t), P(t)) = (1, 0)$ holds true.

Remark 2. In the case of $\tau = 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.

Theorem 2.4. If the condition

$\begin{split} \frac{b\big(m+\exp(\tau)-\sqrt{(m+\exp(\tau))^2-4mc\exp(\tau)}\big)}{2m} < d\\ < \min\bigg\{b, \frac{b\big(m+\exp(\tau)+\sqrt{(m+\exp(\tau))^2-4mc\exp(\tau)}\big)}{2m}\bigg\} \end{split}$

(5)

with $0 < c\leq\frac{(m+\exp(\tau))^2}{4m\exp(\tau)}$ for $\tau>0$ holds, then the DDE model (4) is permanent.

Remark 3. Biologically, Theorem 2.4 means that all the species of the DDE model (4) are survivor under $0 < c\leq\frac{(m+\exp(\tau))^2}{4m\exp(\tau)}$ when the natural death rate of predator $d$ is more than $\frac{b\big(m+\exp(\tau)-\sqrt{(m+\exp(\tau))^2-4mc\exp(\tau)}\big)}{2m}$ and less than $\min\Big\{b, \frac{b\big(m+\exp(\tau)+\sqrt{(m+\exp(\tau))^2-4mc\exp(\tau)}\big)}{2m}\Big\}$.

Let $\bar{N} = N-N^*, \, \bar{P} = P-P^*$ and dropping the bars for simplicity if notations, then the DDE model (4) can be transformed into the following:

$\left\{\begin{array}{l} \dfrac{\text{d}N(t)}{\text{d}t} = (N^*+N(t))(1-(N^*+N(t))) -\dfrac{c(N^*+N(t))(P^*+P(t))}{m(P^*+P(t))+1}, \\ \dfrac{\text{d}P(t)}{\text{d}t} = (P^*+P(t))\Big(-d +\dfrac{b(N^*+N(t-\tau))}{m(P^*+P(t-\tau))+1}\Big), \end{array}\right.$

(6)

which can be rewritten as an abstract differential equation in the phase space $\mathcal{C}([-\tau, 0], R^2_+)$,

$\frac{\text{d}U(t)}{\text{d}t} = L(U(t))+G(U(t)),$

(7)

where $U(t) = (N(t), P(t))^\text{T}$,

$L(U(t)) = \left(\begin{array}{c} q_1N(t)+q_2P(t) \\ q_3N(t-\tau)+q_4P(t-\tau) \\ \end{array}\right),$

(8)

and

$G(U(t)) = \left(\begin{array}{c} -N^2(t)-\frac{cN(t)P(t)}{(mP^*+1)^2}+\frac{cmN^*P^2(t)}{(mP^*+1)^2}\\ \frac{bN(t-\tau)P(t)}{mP^*+1}-bm\frac{N^*P(t)P(t-\tau)+P^*N(t-\tau)P(t-\tau)}{(mP^*+1)^2} +\frac{bm^2N^*P^*P^2(t-\tau)}{(mP^*+1)^3} \\ \end{array}\right),$

(9)

here,

$q_1 = -N^*, \;\; q_2 = -\frac{cN^*}{(mP^*+1)^2}, \;\; q_3 = \frac{bP^*}{mP^*+1}, \;\;q_4 = -\frac{bmN^*P^*}{(mP^*+1)^2}.$

(10)

Referring to [10], we can know that given initial data $(N(\vartheta), P(\vartheta))\in\mathcal{C}([-\tau, 0];R_+^2)$, model (4) will have a unique local or global solution dependent on the parameters $q_i, \, i = 1, 2, 3, 4$.

The characteristic equation of the linearization model $\frac{\text{d}U(t)}{\text{d}t} = L(U(t))$ is

$\lambda^2-q_1\lambda+(q_1q_4-q_2q_3-q_4\lambda)\text{e}^{-\lambda\tau} = 0,$

(11)

where $q_i, \, i = 1, 2, 3, 4$ are defined by (10). After a little algebraic calculation, then we get a unique positive root:

$\omega_0 = \sqrt{\frac{-(q_1^2-q_4^2)+\sqrt{(q_1^2-q_4^2)^2+4(q_1q_4-q_2q_3)^2}}{2}},$

(12)

and the corresponding critical value of delay $\tau = \tau_k$ at which $E^* = (N^*, P^*)$ loses stability is given by

$\tau_k = \frac{1}{\omega_0}\arcsin\left(-\frac{q_4\omega_0^3+(q_1q_4-q_2q_3)q_1\omega_0} {(q_1q_4-q_2q_3)^2+q_4^2\omega_0^2}\right)+\frac{2k\pi}{\omega_0}, \;\;\; k = 0, 1, 2, \ldots$

(13)

Hence, we have the following theorem.

Theorem 2.5. For the DDE model (4), if $d < b$ holds, then there is a critical value $\tau_0$ which is the least positive value of $\tau_k$ given by (13), such that the positive equilibrium $E^* = (N^*, P^*)$ is asymptotically stable when $\tau\in(0, \tau_0)$ and unstable when $\tau>\tau_0$. Moreover, model (4) undergoes Hopf bifurcation around $E^* = (N^*, P^*)$ at $\tau = \tau_k$ $(k = 0, 1, 2, \ldots)$, that is, a family of periodic solutions bifurcate from $E^*$ as $\tau$ passes through the critical value $\tau_k$ $(k = 0, 1, 2, \ldots)$.

Remark 4. In the absence of delay $\tau = 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 $\tau\neq 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 $\tau = 0$, then there is a positive constant $\tau_0$ which is the least positive value of $\tau_k$ given by (13), such that for $\tau>\tau_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. 2 shows the time-series plots of prey $N$ and predator $P$ and phase portraits of the DDE model (4) with different time delay $\tau$, other parameters are taken as $c = 0.9, \, b = 0.7, \, d = 0.3, \, m = 0.1$ which are the same as in 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 $\tau_0 = 3.46$, Hopf bifurcation occurs and a bifurcation periodic solution exists at $\tau = \tau_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 $\tau = 2.8 < \tau_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 $\tau>\tau_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) = \bigl(N(t), P(t)\bigr)$ and $\big|Z(t)\big| = \sqrt{N^2(t)+P^2(t)}$.

Theorem 2.6. For any given initial value $Z(\vartheta) = \{(N(\vartheta), P(\vartheta)):-\tau\leq \vartheta\leq 0\}\in\mathcal{C}\big([-\tau, 0]; R_+^2\big)$, if $|Z(t-\tau)|< \frac{d-1}{b}$ on $t\geq-\tau$, then there is a unique positive solution $Z(t) = (N(t), P(t))$ to the SDDE model (3) on $t\geq-\tau$ and the solution will remain in $R_+^2$ with probability one.

Remark 5. Theorem 2.6 shows that the solution of the SDDE model (3) will remain in $R_+^2$. It can be said that model (3) can reserve the global stability of the positive equilibrium $E^* = (N^*, P^*)$ of the ODE model (2).

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.

Definition 2.7. The solution of the SDDE model (3) is said to be stochastically ultimately bounded if for any $\varepsilon\in(0, 1)$, there is a positive constant $\chi = \chi(\epsilon)$, such that for any initial value $Z(\vartheta) = \{(N(\vartheta), P(\vartheta)):-\tau\leq \vartheta\leq 0\}\in\mathcal{C}\big([-\tau, 0]; R_+^2\big)$, the solution of model (3) has the property that

$\limsup\limits_{t\rightarrow\infty}\mathbf{P}\big\{|(N(t), P(t))|>\chi\big\} < \varepsilon.$

(14)

Theorem 2.8. For any $\theta\in(0, 1)$, if $d\leq \frac{\theta b^2}{\theta b+m\exp(-\tau)}$ for $\tau>0$, there is a positive constant $C = C(\theta)$, which is independent of initial value $Z(\vartheta) = \{(N(\vartheta), P(\vartheta)):-\tau\leq \vartheta\leq 0\}\in\mathcal{C}\big([-\tau, 0]; R_+^2\big)$, such that the solution of the SDDE model (3) has the property that

$\limsup\limits_{t\rightarrow\infty}\mathbf{E}|(N(t), P(t))|^\theta\leq C.$

(15)

Theorem 2.9. For any $\theta\in(0, 1)$, if $d\leq \frac{\theta b^2}{\theta b+m\exp(-\tau)}$ for $\tau>0$, the SDDE model (3) is stochastically ultimately bounded.

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).

Theorem 2.10. Assume the same conditions of Theorem 2.6 hold. For any given initial value $Z(\vartheta) = \{(N(\vartheta), P(\vartheta)):-\tau\leq \vartheta\leq 0\}\in\mathcal{C}\big([-\tau, 0]; R_+^2\big)$, the solution $Z(t) = (N(t), P(t))$ to the SDDE model (3) has the property that

$\limsup\limits_{t\rightarrow\infty}\frac{\ln|Z(t)|}{t}\leq 1+d+\frac{\breve{C}^2}{2\hat{\sigma}^2} \;\;\; a.s.$

where $\hat{\sigma} = \min\{\sigma_1c_1, \, \sigma_2c_2\}, \, \check{C} = \max\{b, c\}$.

---

3. Numerical simulations

In this section, numerical simulation results of the SDDE model (3) for different values of $\tau$ and $\sigma_1, \, \sigma_2$ will be mentioned in subsequent two cases, firstly for $\tau = 0$ and secondly for $\tau\neq 0$ by increasing the magnitude of noise $\sigma_1, \, \sigma_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 $\sigma_1, \, \sigma_2$ and $\tau = 0$, other parametric values as (16) which are the same as in Fig. 1. In Fig. 3(a), without noise, i.e., $\sigma_1 = \sigma_2 = 0$, the interior equilibrium $E^* = (0.46, 0.64)$ of the ODE model (2) is asymptotically stable. Let $\sigma_1 = 0, \, \sigma_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 $\sigma_1 = 0.12, \, \sigma_2 = 0$, the fluctuation of prey species $N$ occurs, thus the predator species $P$ will become instability. Choosing $\sigma_1 = 1.3$ and $\sigma_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 $\sigma_1, \, \sigma_2$ and keeping $\tau = 2.8 < \tau_0 = 3.46$. In Fig. 4(a), let $\sigma_1 = 0$ and $\sigma_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 $\sigma_1 = 0.12$ and $\sigma_2 = 0$, one can see that the stochastic oscillation happens in both species with the weak fluctuation. In Fig. 4(c), increasing $\sigma_1, \, \sigma_2$ to $\sigma_1 = \sigma_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.

Fig. 5 shows the fluctuation in population densities of the prey and predator species of model (3) for different values of $\sigma_1, \, \sigma_2$ and $\tau = 3.9>\tau_0$ as mentioned in the caption. Fig. 5 illustrates the stochastic oscillation in population densities of both the species are relatively significant and they oscillate as if there is no external noises, see Figs. 5(a) and 5(b). Also, the amplitude of stochastic fluctuations increase appreciably for larger magnitudes of $\tau$ and $\sigma_2$, see Figs. 5(c).

---

4. Conclusions and remarks

For the stochastic delayed Harrison-type predator-prey 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 (3) and magnitudes of time delay and environmental noise. Our analytical results partially provide answers to the two questions proposed in the introduction: 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? What are the joint effects of time delay and noise on the Harrison-type predator-prey model? More specifically, we have follows:

1. The effect of delay induces instability: Theorem 2.5 indicates that the large time delay can induce instability and oscillations via Hopf bifurcation 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 $\tau$ is smaller than the critical value $\tau_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 $\tau$ is larger than the critical value $\tau_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 (3) of the delay model (4) will remain to have a positive stable solution, and this theorem gives a result on the robustness of the positive stable solution. Under the conditions of Theorems 2.8 and 2.9, the property of the ultimate boundedness of the SDDE model (3) will not change no matter the environmental noise is large or small. That is to say, the property of this boundedness is very robust under the noise. In short, under certain conditions, the DDE model (4) and the associated SDDE model (3) behave similarly in the sense that both the prey and 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) has almost surely asymptotic property. By fixing the values of $c, \, b, \, d$ and $m$ in model (3), we perform a large number of computer simulations with different values of the time-delay $\tau$ and the noise $\sigma_1, \, \sigma_2$. Without the time-delay $\tau = 0$, increasing the noise intensity can induce the fluctuation of the ODE model (2), see Fig. 3. With the time-delay, when $\tau < \tau_0$, add the noise intensity to the DDE model (4), it shows that the positive constant steady state $E^* = (N^*, P^*)$ of the model is unstable due the noise effects, see Fig. 4. When $\tau>\tau_0$, both the time-delay and noise factors can produce instability and increase the amplitude of the periodic oscillating orbits of two species, see Fig. 5.

---

5. Proofs

Proof of Theorem 2.1.

Proof. For all $t\in[0, t_0]$ where $t_0\in(0, \infty)$, we need to show that $N(t)>0$ and $P(t)>0$. Suppose that is not true. Then, there exists $T\in(0, t_0)$ such that for all $t\in[0, T]$, $N(t)>0$ and $P(t)>0$, and either $N(T) = 0$ or $P(T) = 0$. For all $t\in[0, T]$, we obtain

$\left\{\begin{array}{l} N(t) = N(0)\exp\Big(\int_0^t \big(1-N(s)-\frac{cP(s)}{mP(s)+1}\big)\text{d}s+\int_0^t \sigma_1\text{d}B(s)\Big), \\ P(t) = P(0)\exp\Big(\int_0^t\big(\frac{bN(s-\tau)}{mP(s-\tau)+1}-d\big)\text{d}s+\int_0^t \sigma_2\text{d}B(s)\Big). \end{array}\right.$

(17)

As $(N, P)$ is defined and continuous on $[-\tau, T]$ and $B(t)$ is a standard Brownian motion, there is a $C\geq 0$ such that for all $t\in[-\tau, T]$,

$\left\{\begin{array}{l} N(t) = N(0)\exp\Big(\int_0^t \big(1-N(s)-\frac{cP(s)}{mP(s)+1}\big)\text{d}s+\sigma_1B(t)\Big)\geq N(0)\exp(-CT), \\ P(t) = P(0)\exp\Big(\int_0^t\big(\frac{bN(s-\tau)}{mP(s-\tau)+1}-d\big)\text{d}s+\sigma_2B(t)\Big)\geq P(0)\exp(-CT). \end{array}\right.$

(18)

Taking the limit, as $t\rightarrow T$, we get

$N(T)\geq N(0)\exp(-CT)>0, \;\;\; P(T)\geq P(0)\exp(-CT)>0,$

which contradicts the fact that either $N(T) = 0$ or $P(T) = 0$. So, for all $t\in[0, t_0]$, $N(t)>0$ and $P(t)>0$.

Proof of Theorem 2.3.

Proof. For our full model (3), when $\sigma_1 = \sigma_2 = 0$, it is the following delayed model

$\left\{\begin{array}{l} \dfrac{\text{d}N(t)}{\text{d}t} = N(t)(1-N(t))-\dfrac{cN(t)P(t)}{mP(t)+1}, \\ \dfrac{\text{d}P(t)}{\text{d}t} = P(t)\Big(-d+\dfrac{bN(t-\tau)}{mP(t-\tau)+1}\Big), \end{array}\right.$

(19)

with initial conditions $N(\vartheta) = \psi_1(\vartheta)\geq 0, \> P(\vartheta) = \psi_2(\vartheta)\geq 0$ for all $\vartheta\in[-\tau, 0]$, where $\psi_1(\vartheta), \psi_2(\vartheta)\in \mathcal{C}([-\tau, 0], R^2_+)$. From the first equation of (19), we can get $\frac{\text{d}N(t)}{\text{d}t}\leq N(t)(1-N(t))$. A standard comparison theorem shows that $\limsup\limits_{t\rightarrow\infty}N(t)\leq1$. As a result, for any $\epsilon>0$, there exists a $T>0$, such that $N(t)\leq 1+\epsilon$ for $t>T$. Then from the second equation of (19), we define $P$ on all interval $[k\tau, (k+1)\tau]$ with $k = 0, 1, 2\ldots$, then one can see that $P$ is bounded on $[0, t_0]$ if $t_0 < \infty$.

From the predator equation of model (19), we get

$\frac{\text{d}P(t)}{\text{d}t} < bP(t),$

hence, for $t>\tau$,

$P(t)\leq P(t-\tau)\exp(b\tau),$

which is equivalent, for $t>\tau$, to

$P(t-\tau)\geq P(t)\exp(-b\tau).$

For any $\varepsilon>1$, there exists positive $T_\varepsilon$, such that for $t>T_\varepsilon$, $N(t) < \varepsilon$. According to the above inequality, we obtain, for $t>T_\varepsilon+\tau$,

$\frac{\text{d}P(t)}{\text{d}t} < P(t)\bigg(\frac{\varepsilon b}{(mP(t)\exp(-\varepsilon\tau)+1)}-d\bigg),$

which implies by the same arguments use for $N$ that,

$\limsup\limits_{t\rightarrow\infty}P(t)\leq\frac{(\varepsilon b-d)\exp(\varepsilon\tau)}{dm}.$

(20)

Conclusion of this Theorem holds by letting $\varepsilon\rightarrow 1$.

In addition, based on the above analysis, if $d\geq b$, for $t>T_\varepsilon+\tau$, we have $\lim\limits_{t\rightarrow\infty}P(t) = 0$.

Proof of Theorem 2.4.

Proof. From Theorem 2.3, there is a $A = \max\big\{1, \frac{(b-d)\exp(\tau)}{dm}\big\}>0$ independent of the initial condition such that

$\max\big\{\limsup\limits_{t\rightarrow\infty}N(t), \>\limsup\limits_{t\rightarrow\infty}P(t)\big\}\leq A.$

Then, we will show that there is a $B>0$, independent of the initial condition, such that

$\min\big\{\liminf\limits_{t\rightarrow\infty}N(t), \>\liminf\limits_{t\rightarrow\infty}P(t)\big\}\geq B.$

From model (4), for any $\varepsilon>1$ and for large enough $t$, we get $P(t) < \varepsilon\frac{(b-d)\exp(\tau)}{dm}$. Then, we have

$\frac{\text{d}N(t)}{\text{d}t}>N(t)\bigg(1-N(t)-\frac{c\varepsilon(b-d)\exp(\tau)}{dm}\bigg).$

Based on the standard comparison arguments, it follows that

$\liminf\limits_{t\rightarrow\infty}N(t)\geq\bigg(1-\frac{c\varepsilon(b-d)\exp(\tau)}{dm}\bigg)$

and letting $\varepsilon\rightarrow 1$, we have

$\liminf\limits_{t\rightarrow\infty}N(t)\geq\bigg(1-\frac{c(b-d)\exp(\tau)}{dm}\bigg).$

(21)

Set $C_1 = \big(1-\frac{c(b-d)\exp(\tau)}{dm}\big)$, if $c < \frac{dm\exp(-\tau)}{b-d}$ holds, then $C_1>0$. From (21) and Theorem 2.3, for any $\varepsilon>1$, there exists a positive constant $T_\varepsilon$, such that for $t>T_\varepsilon$, $N(t)>\frac{C_1}{\varepsilon}$ and $P(t) < \frac{\varepsilon(b-d)\exp(\tau)}{dm}$. Thus, for $t>T_\varepsilon+\tau$, we obtain

$\frac{\text{d}P(t)}{\text{d}t}>P(t)\bigg(\frac{bC_1}{\varepsilon(mP(t-\tau)+1)}-d\bigg).$

(22)

On the one hand, for $t>T_\varepsilon+\tau$, these inequalities lead to

$\frac{\text{d}P(t)}{\text{d}t}>P(t)\bigg(\frac{C_1bd}{\varepsilon(\varepsilon(b-d)\exp(\tau)+d)}-d\bigg),$

which involves, for $t>T_\varepsilon+\tau$,

$P(t-\tau) < P(t)\exp\bigg(-\bigg(\frac{C_1bd}{\varepsilon(\varepsilon(b-d)\exp(\tau)+d)}-d\bigg)\tau\bigg).$

(23)

On the other hand, from (22) and (23), for $t>T_\varepsilon+\tau$, we obtain

$\frac{\text{d}P(t)}{\text{d}t}>P(t)\bigg(\frac{bC_1}{\varepsilon mP(t)\exp\big(-(\frac{C_1bd}{\varepsilon(\varepsilon(b-d)\exp(\tau)+d)}-d)\tau\big)+\varepsilon}-d\bigg)$

which yields

$\liminf\limits_{t\rightarrow\infty}P(t)\geq\frac{bC_1-\varepsilon d}{\varepsilon dm}\exp\bigg(-\Big(d-\frac{C_1bd}{\varepsilon^2(b-d)\exp(\tau)+\varepsilon d}\Big)\tau\bigg).$

When $\varepsilon\rightarrow 1$, we have

$\liminf\limits_{t\rightarrow\infty}P(t)\geq\frac{bC_1-d}{dm}\exp\bigg(-\Big(d-\frac{C_1bd}{(b-d)\exp(\tau)+d}\Big)\tau\bigg) = C_2.$

Let be $B = \min\{C_1, C_2\}>0$. Therefore, model (4) is permanent.

Proof of Theorem 2.5.

Proof. In the absence of delay, i.e., $\tau = 0$, from Theorem 2.2, the positive equilibrium $E^* = (N^*, P^*)$ of model (4) is locally asymptotically stable provided that $d < b$.

Assume that for some $\tau>0$, $\text{i}\omega$ ($\omega>0$) is a root of Eq.(11), then

$-\omega^2-\text{i}q_1\omega+(q_1q_4-q_2q_3-\text{i}q_4\omega)\text{e}^{-\text{i}\omega\tau} = 0.$

(24)

Separating the real and imaginary parts, we have

$\left\{\begin{array}{l} \omega^2 = (q_1q_4-q_2q_3)\cos\omega\tau-q_4\omega\sin\omega\tau, \\ q_1\omega = -(q_1q_4-q_2q_3)\sin\omega\tau-q_4\omega\cos\omega\tau, \end{array}\right.$

(25)

thus,

$\sin\omega\tau = -\frac{q_4\omega^3+(q_1q_4-q_2q_3)q_1\omega} {(q_1q_4-q_2q_3)^2+q_4^2\omega^2}, \;\;\; \cos\omega\tau = -\frac{q_2q_3\omega^2}{(q_1q_4-q_2q_3)^2+q_4^2\omega^2},$

(26)

which leads to the following equation

$\omega^4+(q_1^2-q_4^2)\omega^2-(q_1q_4-q_2q_3)^2 = 0.$

(27)

Solving for $\omega^2$ in (27), we obtain

$\omega^2 = \frac{-(q_1^2-q_4^2)\pm\sqrt{(q_1^2-q_4^2)^2+4(q_1q_4-q_2q_3)^2}}{2}.$

(28)

Thus, Eq. (27) has a unique positive root

$\omega_0 = \sqrt{\frac{-(q_1^2-q_4^2)+\sqrt{(q_1^2-q_4^2)^2+4(q_1q_4-q_2q_3)^2}}{2}},$

(29)

and substituting this value in (25), then the corresponding critical value of delay $\tau = \tau_k$ at which $E^* = (N^*, P^*)$ loses stability is given by

$\tau_k = \frac{1}{\omega_0}\arcsin\left(-\frac{q_4\omega_0^3+(q_1q_4-q_2q_3)q_1\omega_0} {(q_1q_4-q_2q_3)^2+q_4^2\omega_0^2}\right)+\frac{2k\pi}{\omega_0}, \;\;\; k = 0, 1, 2, \ldots$

(30)

Let $\lambda(\tau) = \alpha(\tau)+\text{i}\omega(\tau)$ be the root of Eq.(11) such that $\alpha(\tau_k) = 0, \, \omega(\tau_k) = \omega_0$. Substituting $\lambda(\tau)$ into Eq.(11) and differentiating both sides of it with respect to $\tau$, we have

$\begin{split} \left(\frac{\text{d}\lambda(\tau)}{\text{d}\tau}\right)^{-1} & = \frac{(2\lambda-q_1)\text{e}^{\lambda\tau} -(q_1q_4-q_2q_3)\tau-(1+\lambda\tau)q_4} {\lambda(q_1q_4-q_2q_3+q_4\lambda)}\\ & = \frac{(2\lambda-q_1)\text{e}^{\lambda\tau}} {\lambda(q_1q_4-q_2q_3+q_4\lambda)} -\frac{q_4}{\lambda(q_1q_4-q_2q_3+q_4\lambda)} -\frac{\tau}{\lambda}, \end{split}$

(31)

which, together with (25), (26) and (29), leads to

$\begin{split} \text{sign}\Big\{\text{Re} \Big(\frac{\text{d}\lambda}{\text{d}\tau}\Big)\Big\}_{\tau = \tau_k} & = \text{sign}\Big\{\text{Re} \Big(\frac{\text{d}\lambda}{\text{d}\tau}\Big)^{-1}\Big\}_{\tau = \tau_k}\\ & = \frac{q_1^2+2\omega^2_0}{\omega^2_0(q_1^2+\omega^2_0)} +\frac{q_4^2}{(q_1q_4-q_2q_3)^2+q_4^2\omega^2_0}>0. \end{split}$

(32)

This implies that all the roots cross the imaginary axis at $\text{i}\omega$ from left to right as $\tau$ increases. Hence, the transversality condition is satisfied. Hopf bifurcation occurs at $\tau = \tau_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) = \text{e}^{u(t)}, \, P(t) = \text{e}^{v(t)}$ and using It$\hat{o}$'s formula, model (3) can be reformulated as follows:

$\left\{\begin{array}{l} \text{d}u(t) = \big(1-\text{e}^{u(t)}-\frac{c\text{e}^{v(t)}}{m\text{e}^{v(t)}+1}-\frac{\sigma_1^2}{2}\big)\text{d}t +\sigma_1\text{d}B(t), \\ \text{d}v(t) = \big(-d+\frac{b\text{e}^{u(t-\tau)}}{m\text{e}^{v(t-\tau)}+1}-\frac{\sigma_1^2}{2}\Big)\text{d}t +\sigma_2\text{d}B(t), \end{array}\right.$

(33)

the coefficient of (33) satisfy the local Lipschitz condition, then for any given initial values $u(\vartheta) = \ln N(\vartheta), \, v(\vartheta) = \ln P(\vartheta), \, \vartheta\in[-\tau, 0]$, there is a unique local solution $u(t), \, v(t)$ on $[-\tau, t_e)$, where $t_e$ is explosion time. Referring to [14,15], to show this solution is global, we need to show $t_e = \infty$ a.s. Let $l_0>0$ be sufficiently large for

$\frac{1}{l_0} < \min\limits_{-\tau\leq t\leq0}\big|Z(t)\big|\leq\max\limits_{-\tau\leq t\leq0}\big|Z(t)\big| < l_0.$

For each integer $l\geq l_0$, define the stopping time

$t_l = \inf\Big\{t\in[0, t_e): N(t)\not\in(\frac{1}{l}, l), \, P(t)\not\in(\frac{1}{l}, l)\Big\},$

where we set $\inf\emptyset = \infty$ ($\emptyset$ represents the empty set) in this paper. $t_l$ is increasing as $l\rightarrow\infty$. Let $t_\infty = \lim\limits_{l\rightarrow\infty}t_l$, then $t_\infty\leq t_e$ a.s.

If we can show that $t_\infty = \infty$ a.s., then $t_e = \infty$ a.s. and $Z(t)\in R_+^2$ a.s. for all $t\geq 0$. To show this statement, we define a $\mathcal{C}^2$-function $V:\> R_+^2\rightarrow R_+$ by

$V(Z) = (N-\log N-1)+(P-\log P-1),$

(34)

which is a nonnegative function on $N>0, \, P>0$. Let $l\geq l_0$ and $T>0$ be arbitrary. For $0\leq t\leq t_l\wedge T = \min\{t_l, T\}$, we apply the It$\hat{o}$'s formula and obtain

$\begin{array}{l} \text{d}\Big[\int_{t-\tau}^t(N(s)+P(s))\text{d}s+V(Z(t))\Big] = \Big[N-N(t-\tau)+(N-1)(1-N-\frac{cP}{mP+1})\\ +P-P(t-\tau)+(P-1)(\frac{bN(t-\tau)}{mP(t-\tau)+1}-d)+\frac{\sigma_1^2+\sigma_2^2}{2}\Big]\text{d}t\\ +\sigma_1(N-1)\text{d}B(t)+\sigma_2(P-1)\text{d}B(t)\\ \leq \Big[\frac{5}{4}+\frac{c}{m}+d-N(t-\tau)-(N-\frac{3}{2})^2+(1-d)P-P(t-\tau)+bN(t-\tau)P+\frac{\sigma_1^2+\sigma_2^2}{2}\Big]\text{d}t\\ +\sigma_1(N-1)\text{d}B(t)+\sigma_2(P-1)\text{d}B(t)\\ \leq \Big[\frac{5}{4}+\frac{c}{m}+d-N(t-\tau)-(N-\frac{3}{2})^2-P(t-\tau)-(d-1-b|Z(t-\tau)|)P+\frac{\sigma_1^2+\sigma_2^2}{2}\Big]\text{d}t\\ +\sigma_1(N-1)\text{d}B(t)+\sigma_2(P-1)\text{d}B(t)\\ \leq C\text{d}t+\sigma_1(N-1)\text{d}B(t)+\sigma_2(P-1)\text{d}B(t), \end{array}$

(35)

where $C$ is a positive number. Integrating both sides of the above inequality from $0$ to $t_l\wedge T = \min\{t_l, T\}$ and then taking the expectations give rise to,

$\begin{split} \mathbf{E}\bigg[\int_{t_l\wedge T-\tau}^{t_l\wedge T}(N(s)+P(s))\text{d}s+V(N(t_l\wedge T), P(t_l\wedge T))\bigg]\\ \leq\int_{-\tau}^0(N(s)+P(s))\text{d}s+V(N(0), P(0))+CT, \end{split}$

thus,

$\mathbf{E}V(N(t_l\wedge T), P(t_l\wedge T))\leq\int_{-\tau}^0(N(s)+P(s))\text{d}s+V(N(0), P(0))+CT.$

(36)

Set $\Omega_l = \{t_l\leq T\}$ for $l\geq l_1$. For every $\mu\in\Omega_l$, $N(t_l, \mu)$ and $P(t_l, \mu)$ equal either $l$ or $\frac{1}{l}$, hence,

$V(N(t_l, \mu), P(t_l, \mu))\geq\min\Big\{l-\log l-1, \frac{1}{l}+\log l-1\Big\}.$

(37)

It then follows from (36) that

$\begin{split} \int_{-\tau}^0(N(s)+P(s))\text{d}s&+V(N(0), P(0))+CT\geq\mathbf{E}\big[1_{\Omega_l}(\mu)V(N(t_l), P(t_l))\big]\\ &\geq\mathbf{P}\{t_l\leq T\}\min\Big\{l-\log l-1, \frac{1}{l}+\log l-1\Big\}, \end{split}$

(38)

where $1_{\Omega_l}$ is the indicator function of $\Omega_l$. Letting $l\rightarrow\infty$ leads to $\lim\limits_{l\rightarrow\infty}\mathbf{P}\{t_l\leq T\} = 0$. Since $T>0$ is arbitrary, we get $\mathbf{P}\{t_\infty < \infty\} = 0$, hence $\mathbf{P}\{t_\infty = \infty\} = 1$ as required. This completes the proof.

Proof of Theorem 2.8.

Proof. Define a function by

$V(N, P) = N^\theta+P^\theta, \;\;\; (N, P)\in R_+^2.$

(39)

For any $t>0$, using It$\hat{o}$'s formula it gives that

$\begin{split} \mathbf{L}V(N, P)& = \theta N^\theta\Big(1-N-\frac{cP}{mP+1}\Big)+\frac{\sigma_1^2}{2}\theta(\theta-1)N^\theta\\ &+\theta P^\theta\Big(\frac{bN(t-\tau)}{mP(t-\tau)+1}-d\Big)+\frac{\sigma_2^2}{2}\theta(\theta-1)P^\theta, \end{split}$

(40)

as $0 < \theta < 1$, then,

$\begin{split} \mathbf{L}V(N, P)&\leq \theta N^\theta(1-N)-\frac{\sigma_1^2}{2}\theta(1-\theta)N^\theta+\theta P^\theta\frac{bN(t-\tau)}{mP(t-\tau)+1}-\frac{\sigma_2^2}{2}\theta(1-\theta)P^\theta\\ &\leq \theta N^\theta-\frac{\sigma_1^2}{2}\theta(1-\theta)N^\theta+\frac{b\theta(b-d)\exp(\tau)}{dm}|Z(t-\tau)|^2-\frac{\sigma_2^2}{2}\theta(1-\theta)P^\theta\\ & = F(N, P)-V(N, P)-\text{e}^{\tau}|Z(t)|^2+\frac{b\theta(b-d)\exp(\tau)}{dm}|Z(t-\tau)|^2 \end{split}$

(41)

where $F(N, P) = (\theta+1) N^\theta+P^\theta-\frac{\sigma_1^2}{2}\theta(1-\theta)N^\theta-\frac{\sigma_2^2}{2}\theta(1-\theta)P^\theta+\text{e}^{\tau}|Z(t)|^2\leq C_0$ for $(N, P)\in R_+^2$. Thus,

$\mathbf{L}V(N, P)\leq C_0-V(N, P)-\text{e}^{\tau}|Z(t)|^2+\frac{b\theta(b-d)\exp(\tau)}{dm}|Z(t-\tau)|^2.$

Hence, we obtain

$\begin{split} \text{d}V(N, P)& = \mathbf{L}V(N, P)\text{d}t+\sigma_1\theta N^\theta\text{d}B(t)+\sigma_2\theta P^\theta\text{d}B(t)\\ &\leq \big(C_0-V(N, P)-\text{e}^{\tau}|Z(t)|^2+\frac{b\theta(b-d)\exp(\tau)}{dm}|Z(t-\tau)|^2\big)\text{d}t\\& +\sigma_1\theta N^\theta\text{d}B(t)+\sigma_2\theta P^\theta\text{d}B(t). \end{split}$

(42)

Applying It$\hat{o}$'s formula yields

$\begin{split} \text{d}(\text{e}^tV(N, P))& = \text{e}^tV(N, P)\text{d}t+\text{e}^t\text{d}V(N, P)\\ &\leq \text{e}^t\big(C_0-\text{e}^{\tau}|Z(t)|^2+\frac{b\theta(b-d)\exp(\tau)}{dm}|Z(t-\tau)|^2\big)\text{d}t\\ &+\text{e}^t\sigma_1\theta N^\theta\text{d}B(t)+\text{e}^t\sigma_2\theta P^\theta\text{d}B(t). \end{split}$

(43)

If $b\theta(b-d)\exp(\tau)\geq dm$, then

$\begin{split} \text{e}^t\mathbf{E}V(N, P)&\leq V(N(0), P(0))+C_0\text{e}^t-\mathbf{E}\int_0^t\text{e}^{s+\tau}|Z(s)|^2\text{d}s\\ &+\frac{b\theta(b-d)\exp(\tau)}{dm}\mathbf{E}\int_0^t\text{e}^s|Z(s-\tau)|^2\text{d}s\\ & = V(N(0), P(0))+C_0\text{e}^t-\mathbf{E}\int_0^t\text{e}^{s+\tau}|Z(s)|^2\text{d}s\\ &+\frac{b\theta(b-d)}{dm}\mathbf{E}\int_{-\tau}^{t-\tau}\text{e}^{s+\tau}|Z(s)|^2\text{d}s\\ &\leq V(N(0), P(0))+C_0\text{e}^t+\frac{b\theta(b-d)}{dm}\mathbf{E}\int_{-\tau}^{0}\text{e}^{s+\tau}|Z(s)|^2\text{d}s, \end{split}$

which implies that

$\limsup\limits_{t\rightarrow\infty}\mathbf{E}V(N(t), P(t))\leq C_0.$

(44)

Since

$|(N, P)|^{\theta}\leq\sqrt{2^{\theta}}\max\big\{N^{\theta}, P^{\theta}\big\}\leq\sqrt{2^{\theta}}V(N, P),$

(45)

we get

$\limsup\limits_{t\rightarrow\infty}\mathbf{E}|(N(t), P(t))|^\theta\leq\sqrt{2^{\theta}}\limsup\limits_{t\rightarrow\infty}\mathbf{E}V(N(t), P(t))\leq\sqrt{2^{\theta}}C_0 = C(\theta).$

(46)

The proof is complete.

Proof of Theorem 2.9.

Proof. The solutions of model (3) will remain in $R_+^2$ for all $t\geq 0$ with probability one from Theorem 2.6. From Theorem 2.8, there exists $C>0$ such that

$\limsup\limits_{t\rightarrow\infty}\mathbf{E}|(N(t), P(t))|^\theta\leq C.$

(47)

Then, we choose a constant $\epsilon>0$ sufficiently large such $\frac{C}{\epsilon} < 1$. By Chebyshev inequality,

$\mathbf{P}\big\{N^\theta+P^\theta>\epsilon\big\}\leq\frac{1}{\epsilon}\mathbf{E}\big[N^\theta+P^\theta\big]\leq\frac{C}{\epsilon}: = \varepsilon,$

which implies that

$1-\varepsilon\leq\mathbf{P}\big\{N^\theta+P^\theta < \epsilon\big\}\leq\mathbf{P}\Big\{\frac{1}{\epsilon}\leq N(t)+P(t)\leq\epsilon\Big\}.$

Noting $\big(N(t)+P(t)\big)^2\leq 2\big|(N(t), P(t))\big|^2\leq 2\big(N(t)+P(t)\big)^2$, we obtain

$\mathbf{P}\Big\{\frac{1}{\sqrt{2}\epsilon}\leq\frac{N(t)+P(t)}{\sqrt{2}}\leq\big|(N(t), P(t))\big|\leq N(t)+P(t)\leq\epsilon\Big\}\geq 1-\varepsilon.$

(48)

According to the Definition 2.7, the SDDE model (3) is stochastically ultimately bounded.

Proof of Theorem 2.10.

Proof. Define a function by $V(Z) = N+P$, using It$\hat{o}$'s formula it gives that

$\begin{split} &\ln V(Z(t)) = \ln V(Z(0))+\int_0^t\bigg(\frac{N}{V(Z(s))}\Big(1-N-\frac{cP}{mP+1}\Big)-\frac{\sigma_1^2N^2}{2V^2(Z(s))}\bigg)\text{d}s\\ &+\int_0^t\bigg(\frac{P}{V(Z(s))}\Big(\frac{bN(s-\tau)}{mP(s-\tau)+1})-d\Big)-\frac{\sigma_2^2P^2}{2V^2(Z(s))}\bigg)\text{d}s+M_1(t)+M_2(t), \end{split}$

(49)

where $M_i(t)$ is a martingale defined by $M_i(t) = \int_0^t\frac{\sigma_iZ_i}{V(Z(s))}\text{d}B(s), \>i = 1, 2$. The quadratic variation of this martingale is

$\langle M_i(t), M_i(t)\rangle = \int_0^t\frac{\sigma_i^2Z_i^2}{V(Z(s))}\text{d}s, \;\;\;i = 1, 2.$

Let $\varepsilon\in(0, 1)$ and for arbitrary positive constant $\theta$, by the exponential martingale inequality [14], for each $k\geq 1$,

$\mathbf{P}\bigg\{\sup\limits_{0\leq t\leq k}\Big(M_i(t)-\frac{\varepsilon}{2}\langle M_i(t), M_i(t)\rangle\bigg)>\frac{\theta \ln k}{\varepsilon}\Big\}\leq k^{-\theta}.$

Since the series $\Sigma_{k = 1}^\infty k^{-\theta}$ converges, applying the Borel-Cantelli lemma, there is $\Omega_0\subset\Omega$ with $\mathbf{P}(\Omega_0) = 1$ such that for $\rho\in\Omega_0$ there exists an integer $k_0 = k_0(\rho)$ and choosing $\theta = 2$, we have

$M_i(t)\leq\frac{\varepsilon}{2}\langle M_i(t), M_i(t)\rangle+\frac{2}{\varepsilon}\ln k$

for all $0\leq t\leq k$ and $k\geq k_0(\rho)$. Substituting this inequality into (49) we obtain

$\begin{split} \ln V(Z(t))&\leq\ln V(Z(0))+\int_0^t\bigg(\frac{N}{V(Z(s))}\Big(1-N-\frac{cP}{mP+1}\Big)-\frac{(1-\varepsilon)\sigma_1^2N^2}{2V^2(Z(s))}\bigg)\text{d}s\\ &+\int_0^t\bigg(\frac{P}{V(Z(s))}\Big(\frac{bN(s-\tau)}{mP(s-\tau)+1})-d\Big)-\frac{(1-\varepsilon)\sigma_2^2P^2}{2V^2(Z(s))}\bigg)\text{d}s+\frac{4}{\varepsilon}\ln k \end{split}$

(50)

for all $0\leq t\leq k, \, k\geq k_0(\rho)$ and $\rho\in\Omega_0$. By computing the following,

$\begin{split} &\frac{N}{V(Z(s))}\Big(1-N-\frac{cP}{mP+1}\Big)\leq 1+cP, \\ &\frac{P}{V(Z(s))}\Big(\frac{bN(s-\tau)}{mP(s-\tau)+1})-d\Big)\leq bN(s-\tau)+d, \\ &\frac{\sigma_1^2N^2}{2V^2(Z(s))}\geq\frac{\sigma_1^2c_1^2N}{2}, \;\;\; \frac{\sigma_2^2P^2}{2V^2(Z(s))}\geq\frac{\sigma_2^2c_2^2P}{2} \end{split}$

(51)

where $c_1, c_2\in(0, 1)$ are positive constants. Substituting (51) into (50) has

$\begin{split} \ln V(Z(t))&\leq\int_0^t\bigg(1+cP(s)+bN(s-\tau)+d-\frac{(1-\varepsilon)\sigma_1^2c_1^2N}{2} -\frac{(1-\varepsilon)\sigma_2^2c_2^2P}{2}\bigg)\text{d}s\\ &+\ln V(Z(0))+\frac{4}{\varepsilon}\ln k\\ &\leq\int_0^t\bigg(1+cP(s)+bN(s-\tau)+d-\frac{(1-\varepsilon)\hat{\sigma}^2|Z|^2}{2}\bigg)\text{d}s\\ &+\ln V(Z(0))+\frac{4}{\varepsilon}\ln k \end {split}$

(52)

for all $0\leq t\leq k, \, k\geq k_0(\rho)$ and $\rho\in\Omega_0$. From (52), we get

$\begin{split} \ln V(Z(t))&+\frac{(1-2\varepsilon)\hat{\sigma}^2}{4}\int_0^t|Z(s)|^2\text{d}s\\ &\leq\ln V(Z(0))+\frac{4}{\varepsilon}\ln k +\int_0^t\bigg(1+cP(s)+bN(s-\tau)+d-\frac{\hat{\sigma}^2|Z|^2}{2}\bigg)\text{d}s\\ &\leq C_1+\int_0^t\bigg(1+cP(s)+bN(s)+d-\frac{\hat{\sigma}^2|Z|^2}{2}\bigg)\text{d}s+\frac{4}{\varepsilon}\ln k, \end{split}$

where $C_1 = \ln V(Z(0))+b\int_{-\tau}^0N(s)\text{d}s$. Due to

$1+d+cP(t)+bN(t)-\frac{\hat{\sigma}^2|Z|^2}{2}\leq 1+d+\check{C}|Z|-\frac{\hat{\sigma}^2|Z|^2}{2}\leq 1+d+\frac{\check{C}^2}{2\hat{\sigma}^2} = C_2,$

thus, if $\rho\in\Omega_0$,

$\ln V(Z(t))+\frac{(1-2\varepsilon)\hat{\sigma}^2}{4}\int_0^t|Z(s)|^2\text{d}s\leq \check{C}+\frac{4}{\varepsilon}\ln k+C_2 t$

for all $0\leq t\leq k, \, k\geq k_0(\rho)$. Hence, for all $\rho\in\Omega_0$, if $k-1\leq t\leq k, \, k\geq k_0(\rho)$, it gives that

$\frac{1}{t}\Big(\ln V(Z(t))+\frac{(1-2\varepsilon)\hat{\sigma}^2}{4}\int_0^t|Z(s)|^2\text{d}s\Big)\leq\frac{1}{k-1}\big(\check{C}+\frac{4}{\varepsilon}\ln k\big)+C_2,$

which implies

$\limsup\limits_{t\rightarrow\infty}\frac{1}{t}\Big(\ln V(Z(t))+\frac{(1-2\varepsilon)\hat{\sigma}^2}{4}\int_0^t|Z(s)|^2\text{d}s\Big)\leq C_2.$

(53)

Letting $\varepsilon\rightarrow 0$, from (52), we have

$\limsup\limits_{t\rightarrow\infty}\frac{1}{t}\Big(\ln V(Z(t))+\frac{\hat{\sigma}^2}{4}\int_0^t|Z(s)|^2\text{d}s\Big)\leq C_2 \>\> a.s.$

(54)

and using $V(Z)\geq\frac{|Z|}{\sqrt{2}}$, we obtain

$\limsup\limits_{t\rightarrow\infty}\frac{1}{t}\Big(\ln |Z(t)|+\frac{\hat{\sigma}^2}{4}\int_0^t|Z(s)|^2\text{d}s\Big)\leq C_2 \>\> a.s.$

The desired assertion is derived.

---

Acknowledgments

This research is partially supported by the National Natural Science Foundation of China (Grant Nos. 11601226, 11426132, 71501094), the Natural Science Foundation of Jiangsu Province of China (Grant No. BK20140927, BK20150961), and the research funds from Nanjing Tech University and Jiangsu Government Scholarship for Overseas Studies. The work is also partially supported by NSF-DMS (1313312&1716802); NSF-IOS/DMS (1558127), and The James S. McDonnell Foundation 21st Century Science Initiative in Studying Complex Systems Scholar Award (UHC Scholar Award 220020472). This project has been partially supported by grants from the National Science Foundation (DMS-1263374 and DMS-1261211), and the Offices of the President and the Provost of Arizona State University.

---