1. Introduction

The dynamical relationship between prey and predator has long been and will continue to be a dominant theme in ecology due to its universal importance and existence. One important component of the predator-prey is functional response, i.e. the rate of prey consumption by an average predator. The functional response can be classified into two types: predator-dependent and prey-dependent. The classical Holling types I-III [1] [2] are strictly prey-dependent functional response; The main predator-dependent functional response has Crowley-Martin type [3], Hassell-Varley type [4], as well as Beddington-DeAngelis type by Beddington [5] and DeAngelis et al. [6]. There is much significant evidence to suggest that Beddington-DeAngelis functional response occurs quite frequently in natural systems and laboratory (see e.g. [7] [8] ). The classical predator-prey model with Beddington-DeAngelis functional response can be expressed as follows

$\{\begin{array}{l}\text{d}x\left(t\right)=x\left(t\right)\left[a_1-b_{1}x\left(t_1\right)-\frac{c_{1}y}{1+m_{1}x\left(t\right)+m_{2}y\left(t\right)}\right]\text{d}t,\hfill \\ \text{d}y\left(t\right)=y\left(t\right)\left[a_2-b_{2}y\left(t\right)+\frac{c_{2}x\left(t\right)}{1+m_{1}x\left(t\right)+m_{2}y\left(t\right)}\right]\text{d}t.\hfill \end{array}$ (1.1)

where $x\left(t\right)$ and $y\left(t\right)$ represent the size of the prey and predator populations at time t, respectively. The parameter $a_1$ denotes the intrinsic growth rate of the prey population and $a_2$ denotes the death date of the predator population. The parameter $b_1$ and $b_2$ are the density-dependent coefficients of the prey and predator populations, respectively. The parameter $c_1$ and $c_2$ represent the capturing rate of the predator and the rate of conversion of nutrients into the reproduction for the predator, respectively.

However, the model is deterministic, and does not incorporate the effect of environmental noise, which is always present. In the real world, population models are always affected by the environmental noise, which is an important component in an ecosystem [9] [10]. Thus, it is interesting to study how the environmental noise affects the population models. To fit the reality better, many authors have introduced white noise into the population dynamics to reveal the effects of the white noise [11] [12]. Inspired by the above facts, in this paper, we assume that fluctuations in the environment mainly affect the intrinsic growth rate $a_1$ and the death rate $a_2$, that is

$a_1\to a_1+{\alpha }_1{\stackrel{\dot{}}{W}}_1\left(t\right),a_2\to a_2+{\alpha }_2{\stackrel{\dot{}}{W}}_2\left(t\right).$

Then we obtain the following stochastic system

$\{\begin{array}{l}\text{d}x\left(t\right)=x\left(t\right)\left[a_1-b_{1}x\left(t\right)-\frac{c_{1}y}{1+m_{1}x\left(t\right)+m_{2}y\left(t\right)}\right]\text{d}t+{\alpha }_{1}x\left(t\right)\text{d}{W}_1\left(t\right),\hfill \\ \text{d}y\left(t\right)=y\left(t\right)\left[a_2-b_{2}y\left(t\right)+\frac{c_{2}x\left(t\right)}{1+m_{1}x\left(t\right)+m_{2}y\left(t\right)}\right]\text{d}t+{\alpha }_{2}y\left(t\right)\text{d}{W}_2\left(t\right).\hfill \end{array}$ (1.2)

On the other hand, more realistic and interesting models of population interactions should take the effects of time delay into account [13] [14] [15] [16]. In general, delay differential equations can exhibit much more complicated dynamics than differential equations without delay. Liu [17] has investigated global asymptotic stability of the positive equilibrium about stochastic predator-prey system with Beddingtons-DeAngelis and time delay. However, so far as we know a very little amount of work has been done with the stochastic predator-prey system with Beddingtons-DeAngelis and time delay. Therefore it is interesting and important to study the following stochastic delayed predator-prey model with Beddington-DeAngelis functional response.

$\{\begin{array}{l}\text{d}x\left(t\right)=x\left(t\right)\left[a_1-b_{1}x\left(t-{\tau }_1\right)-\frac{c_{1}y}{1+m_{1}x\left(t\right)+m_{2}y\left(t\right)}\right]\text{d}t+{\alpha }_1\left(t\right)\text{d}{W}_1\left(t\right),\hfill \\ \text{d}y\left(t\right)=y\left(t\right)\left[a_2-b_{2}y\left(t-{\tau }_2\right)+\frac{c_{2}x\left(t-{\tau }_3\right)}{1+m_{1}x\left(t-{\tau }_3\right)+m_{2}y\left(t-{\tau }_3\right)}\right]\text{d}t+{\alpha }_{2}y\left(t\right)\text{d}{W}_2\left(t\right).\hfill \end{array}$ (1.3)

with the initial conditions

$x_0\left(\theta \right)={\varphi }_1\left(\theta \right)>0,y_0\left(\theta \right)={\varphi }_2\left(\theta \right)>0,\theta \in \left[-\tau ,0\right],\tau =\max \left\{{\tau }_1,{\tau }_2,{\tau }_3\right\}.$

where $\tau >0$ denotes the delay;

$\varphi =\left({\varphi }_1,{\varphi }_2\right)\in C\left(\left[-\tau ,0\right],R_{+}^2\right)$, $R_{+}^2=\left\{\left(x,y\right):x\ge 0,y\ge 0\right\}$,

$\Vert \varphi \Vert =\max \left\{\left|\varphi \left(\theta \right)\right|:\theta \in \left[-\tau ,0\right]\right\}.$

and $\left|\varphi \right|$ is any norm in $R_{+}^2$. As usual, we use the notation $x_t\left(\theta \right)=x\left(t+\theta \right)$ for $\theta \in \left[-\tau ,0\right]$.

The rest of the paper is organized as follows. In Section 2, we show that system (1.3) has a global positive solution. In Section 3, stochastic ultimate boundedness is studied. In Section 4, we investigate the asymptotic moment estimation. In Section 5, we present numerical simulations to illustrate our mathematical findings. We close the paper with conclusions and discussions in Section 6.

2. Global Positive Solutions

Throughout this paper, unless otherwise specified, let $\left(\text{Ω},{\left\{{\mathcal{F}}_t\right\}}_{t\ge 0},\mathcal{P}\right)$ be a complete probability space with a filtration ${\left\{{\mathcal{F}}_t\right\}}_{t\ge 0}$ satisfying the usual conditions (i.e. it is right continuous and ${\mathcal{F}}_0$ contains all $\mathcal{P}$ -null sets). Moreover, let $W_i\left(t\right),\left(i=1,2\right)$ be standard Brownian motions defined on this probability space. Also let $R_{+}^n=\left\{x\in R^n:x_i>0\text{for}\text{all}1\le i\le n\right\}$.

In order for a stochastic differential equation to have a global solution for any given initial condition, it is generally necessary to data the coefficients of the equation are generally required to satisfy the liner growth condition and local Lipschitz condition (see e.g. [18] ). However, the coefficients of (1.3) neither obey the linear growth condition nor local Lipschitz condition. The existence of local positive solutions is given by variable substitution and Itô’s formula.

Lemma 2.1. For any initial value $\left\{\left(x\left(t\right),y\left(t\right)\right):-\tau \le t\le 0\right\}\in C\left(\left[-\tau ,0\right];R_{+}^2\right)$, there is a unique positive local solution $\left(x\left(t\right),y\left(t\right)\right),t\in \left[-\tau ,{\tau }_e\right)$ of system (1.3), where $\tau =\max \left\{{\tau }_1,{\tau }_2,{\tau }_3\right\}$ and ${\tau }_e$ is the explosion time.

Proof. Consider the following system

$\{\begin{array}{l}\text{d}f\left(t\right)=\left[a_1-b_1{\text{e}}^{f\left(t-{\tau }_1\right)}-\frac{1}{2}{\alpha }_1^2{\text{e}}^{2f\left(t\right)}-\frac{c_1{\text{e}}^{g\left(t\right)}}{1+m_1{\text{e}}^{f\left(t\right)}+m_2{\text{e}}^{g\left(t\right)}}\right]\text{d}t+{\alpha }_1\text{d}{W}_1\left(t\right),\hfill \\ \text{d}g\left(t\right)=\left[a_2-b_2{\text{e}}^{g\left(t-{\tau }_2\right)}-\frac{1}{2}{\alpha }_2^2{\text{e}}^{2g\left(t\right)}+\frac{c_2{\text{e}}^{f\left(t-{\tau }_3\right)}}{1+m_1{\text{e}}^{f\left(t-{\tau }_3\right)}+m_2{\text{e}}^{g\left(t-{\tau }_3\right)}}\right]\text{d}t+{\alpha }_2\text{d}{W}_2\left(t\right).\hfill \end{array}$ (2.1)

with initial value $f\left(0\right)=\log x_0,g\left(0\right)=\log y_0$. It is clear that the coefficient of system (2.1) satisfy local Lipschitz condition, then there is an unique local solution $\left(f\left(t\right),g\left(t\right)\right),t\in \left[0,{\tau }_e\right)$ of system (2.1). Therefore, by Itô’s formula, it is easy to find that $\left(x\left(t\right)={\text{e}}^{f\left(t\right)},y\left(t\right)={\text{e}}^{g\left(t\right)}\right)$ is the unique positive local solution of the system (1.3) with the initial value $x_0>0,y_0>0$.

Next, give the existence of the positive solution.

Theorem 2.1. For any given initial value

$\left\{\left(x\left(t\right),y\left(t\right)\right):-\tau \le t\le 0\right\}\in C\left(\left[-\tau ,0\right];R_{+}^2\right)$,

there is a unique solution $\left(x\left(t\right),y\left(t\right)\right)$ of system (1.3) on $t\ge 0$, and the solution will remain in $R_{+}^2$ with probability 1.

Proof. Since, Lemma 2.1 shows that there is a positive local solution $\left(x\left(t\right),y\left(t\right)\right),t\in \left[0,{\tau }_e\right)$ of system (1.3), then to show this solution is global, we only need to show that ${\tau }_e=\infty ,a.s.$, Let $m_0\ge 0$ be sufficiently large so that both $x_0$ and $y_0$ lie within the interval $\left[1/u_0,u_0\right]$. For each integer $u>u_0$, define the stopping time

${\tau }_u=\inf \left\{t\in \left[0,{\tau }_e\right):x\left(t\right)\notin \left(1/u,u\right),\text{or}y\left(t\right)\notin \left(1/u,u\right)\right\}.$

where throughout this paper, we set $\inf \Phi =\infty$ (as usual $\Phi$ denotes the empty set). Clearly, ${\tau }_u$ is increasing as $u\to \infty$. Set ${\tau }_{\infty }={\lim }_{u\to \infty }{\tau }_u$, Whence ${\tau }_{\infty }<{\tau }_e,a.s.$. If we can show that ${\tau }_{\infty }=\infty ,a.s.$, Then ${\tau }_e=\infty$ and $\left(x\left(t\right),y\left(t\right)\right)\in R_{+}^2,a.s.$. For if this statement is false, then there are a pair of constants $T>0$ and $\epsilon \in \left(0,1\right)$, such that

$P\left\{{\tau }_u\le T\right\}>\epsilon .$

Hence there is an integer $U_1\ge U_0$ such that

$P\left\{{\tau }_u\le T\right\}\ge \epsilon ,u\ge U_1.$ (2.2)

Define a C²-function $V:R_{+}^2\to {\bar{R}}_{+}$ by

$V_1\left(x\left(t\right),y\left(t\right)\right)=\left(\sqrt{x}-1-0.5\log x\right)+\left(\sqrt{y}-1-0.5\log y\right).$

$V_2\left(x\left(t\right),y\left(t\right)\right)=V_1\left(x\left(t\right),y\left(t\right)\right)+{\displaystyle {\int }_{t-{\tau }_1}^t{\left|x\left(s\right)\right|}^2\text{d}s}+{\displaystyle {\int }_{t-{\tau }_2}^t{\left|y\left(s\right)\right|}^2\text{d}s}.$

The non-negativity of $V_1\left(x\left(t\right),y\left(t\right)\right)$ can be seen from

$\sqrt{u}-1-0.5\log u\ge 0,\forall u>0$.

Using Itô’s formula, we get

$\begin{array}{l}\text{d}{V}_2=\text{d}\left[V_1\left(x\left(t\right),y\left(t\right)\right)+{\displaystyle {\int }_{t-{\tau }_1}^t{\left|x\left(s\right)\right|}^2\text{d}s}+{\displaystyle {\int }_{t-{\tau }_2}^t{\left|y\left(s\right)\right|}^2\text{d}s}\right]\\ =0.5\left(x^{-0.5}\left(t\right)-x^{-1}\left(t\right)\right)x\left(t\right)[(a_1-b_{1}x\left(t-{\tau }_1\right)\begin{array}{c}\\ \end{array}\\ -\frac{c_{1}y\left(t\right)}{1+m_{1}x\left(t\right)+m_{2}y\left(t\right)})\text{d}t+{\alpha }_1\text{d}{W}_1\left(t\right)]\\ +0.5{\alpha }_1^2{x}^2\left(t\right)\left(-0.25x^{-1.5}\left(t\right)+0.5x^{-2}\left(t\right)\right)\text{d}t\end{array}$

$\begin{array}{l}+0.5\left(y^{-0.5}\left(t\right)-y^{-1}\left(t\right)\right)y\left(t\right)[(a_2-b_{2}y\left(t-{\tau }_1\right)\begin{array}{c}\\ \end{array}\\ -\frac{c_{2}x\left(t-{\tau }_3\right)}{1+m_{1}x\left(t-{\tau }_3\right)+m_{2}y\left(t-{\tau }_3\right)})\text{d}t+{\alpha }_2\text{d}{W}_2\left(t\right)]\\ +0.5{\alpha }_2^2{y}^2\left(t\right)\left(-0.25y^{-1.5}\left(t\right)+0.5y^{-2}\left(t\right)\right)\text{d}t\\ +[{\left|x\left(t\right)\right|}^2-{\left|x\left(t-{\tau }_1\right)\right|}^2+{\left|y\left(t\right)\right|}^2-{\left|y\left(t-{\tau }_2\right)\right|}^2]\text{d}t\end{array}$

$\begin{array}{l}=0.5\left(x^{0.5}\left(t\right)-1\right)\left[\left(a_1-b_{1}x\left(t-{\tau }_1\right)-\frac{c_{1}y\left(t\right)}{1+m_{1}x\left(t\right)+m_{2}y\left(t\right)}\right)\text{d}t+{\alpha }_1\text{d}{W}_1\left(t\right)\right]\\ +0.5{\alpha }_1^2\left(-0.25x^{0.5}\left(t\right)+0.5\right)\text{d}t+0.5\left(y^{0.5}\left(t\right)-1\right)[(a_2-b_{2}y\left(t-{\tau }_1\right)\begin{array}{c}\\ \end{array}\\ -\frac{c_{2}x\left(t-{\tau }_3\right)}{1+m_{1}x\left(t-{\tau }_3\right)+m_{2}y\left(t-{\tau }_3\right)})\text{d}t+{\alpha }_2\text{d}{W}_2\left(t\right)]\\ +0.5{\alpha }_2^2\left(-0.25y^{0.5}\left(t\right)+0.5\right)\text{d}t\\ +\left[{\left|x\left(t\right)\right|}^2-{\left|x\left(t-{\tau }_1\right)\right|}^2+{\left|y\left(t\right)\right|}^2-{\left|y\left(t-{\tau }_2\right)\right|}^2\right]\text{d}t\end{array}$

$\begin{array}{l}=\{{\left|x\left(t\right)\right|}^2-{\left|x\left(t-{\tau }_1\right)\right|}^2+0.5[a_1{x}^{0.5}\left(t\right)-b_{1}x\left(t-{\tau }_1\right)x^{0.5}\left(t\right)\begin{array}{c}\\ \end{array}\\ -\frac{c_{1}y\left(t\right)x^{0.5}\left(t\right)}{1+m_{1}x\left(t\right)+m_{2}y\left(t\right)}-a_1+\frac{c_{1}y\left(t\right)}{1+m_{1}x\left(t\right)+m_{2}y\left(t\right)}-0.25{\alpha }_1^2{x}^{0.5}\left(t\right)+0.5{\alpha }_1^2]\\ \begin{array}{c}\\ \end{array}+0.5b_{1}x\left(t-{\tau }_1\right)\}\text{d}t+\{{\left|y\left(t\right)\right|}^2-{\left|y\left(t-{\tau }_2\right)\right|}^2\begin{array}{c}\\ \end{array}\end{array}$

$\begin{array}{l}+0.5[a_2{y}^{0.5}\left(t\right)-b_{2}y\left(t-{\tau }_2\right)y^{0.5}\left(t\right)+\frac{c_{2}x\left(t-{\tau }_3\right)y^{0.5}\left(t\right)}{1+m_{1}x\left(t-{\tau }_3\right)+m_{2}y\left(t-{\tau }_3\right)}-a_2\\ -\frac{c_{2}x\left(t-{\tau }_3\right)y^{0.5}\left(t\right)}{1+m_{1}x\left(t-{\tau }_3\right)+m_{2}y\left(t-{\tau }_3\right)}-0.25{\alpha }_2^2{y}^{0.5}\left(t\right)+0.5{\alpha }_2^2]+0.5b_{2}y\left(t-{\tau }_2\right)\}\text{d}t\\ +0.5\left(x^{0.5}\left(t\right)-1\right)\text{d}{W}_1\left(t\right)+0.5\left(y^{0.5}\left(t\right)-1\right)\text{d}{W}_2(t)\end{array}$

$\begin{array}{l}\le \{{\left|x\left(t\right)\right|}^2-{\left|x\left(t-{\tau }_1\right)\right|}^2+0.5[a_1{x}^{0.5}\left(t\right)+\frac{c_{1}y\left(t\right)}{1+m_{1}x\left(t\right)+m_{2}y\left(t\right)}\\ \begin{array}{c}\\ \end{array}-0.25{\alpha }_1^2{x}^{0.5}\left(t\right)+0.5{\alpha }_1^2]+{\left(0.25b_1\right)}^2+{\left|x\left(t-{\tau }_1\right)\right|}^2\}\text{d}t\\ +\{{\left|y\left(t\right)\right|}^2-{\left|y\left(t-{\tau }_2\right)\right|}^2+0.5[a_2{y}^{0.5}\left(t\right)+\frac{c_{2}x\left(t-{\tau }_3\right)y^{0.5}\left(t\right)}{1+m_{1}x\left(t-{\tau }_3\right)+m_{2}y\left(t-{\tau }_3\right)}\\ \begin{array}{c}\\ \end{array}-0.25{\alpha }_2^2{y}^{0.5}\left(t\right)+0.5{\alpha }_2^2]+{\left(0.25b_2\right)}^2+{\left|y\left(t-{\tau }_2\right)\right|}^2\}\text{d}t\\ +0.5\left(x^{0.5}\left(t\right)-1\right){\alpha }_1\text{d}{W}_1\left(t\right)+0.5\left(y^{0.5}\left(t\right)-1\right){\alpha }_2\text{d}{W}_2(t)\end{array}$

$\begin{array}{l}\le \left\{{\left|x\left(t\right)\right|}^2+0.5\left[a_1{x}^{0.5}\left(t\right)+c_1/m_2+0.5{\alpha }_1^2+0.125b_1^2-0.25{\alpha }_1^2{x}^{0.5}\left(t\right)\right]\right\}\text{d}t\\ +\left\{{\left|y\left(t\right)\right|}^2+0.5\left[a_2{y}^{0.5}\left(t\right)+c_2{y}^{0.5}\left(t\right)/m_1+0.5{\alpha }_2^2+0.125b_2^2-0.25{\alpha }_2^2{y}^{0.5}\left(t\right)\right]\right\}\text{d}t\\ +0.5\left(x^{0.5}\left(t\right)-1\right){\alpha }_1\text{d}{W}_1\left(t\right)+0.5\left(y^{0.5}\left(t\right)-1\right){\alpha }_2\text{d}{W}_2\left(t\right)\\ =M\left(x\left(t\right),y\left(t\right)\right)\text{d}t+0.5\left(x^{0.5}\left(t\right)-1\right){\alpha }_1\text{d}{W}_1\left(t\right)+0.5\left(y^{0.5}\left(t\right)-1\right){\alpha }_2\text{d}{W}_2\left(t\right).\end{array}$ (2.3)

where

$\begin{array}{l}M\left(x\left(t\right),y\left(t\right)\right)\\ =\left\{{\left|x\left(t\right)\right|}^2+0.5\left[a_1{x}^{0.5}\left(t\right)+c_1/m_2+0.5{\alpha }_1^2+0.125b_1^2-0.25{\alpha }_1^2{x}^{0.5}\left(t\right)\right]\right\}\\ +\left\{{\left|y\left(t\right)\right|}^2+0.5\left[a_2{y}^{0.5}\left(t\right)+c_2{y}^{0.5}\left(t\right)/m_1+0.5{\alpha }_2^2+0.125b_2^2-0.25{\alpha }_2^2{y}^{0.5}\left(t\right)\right]\right\}.\end{array}$

which implies that

$M\left(x\left(t\right),y\left(t\right)\right)\le M^{*},$

Because next inequality exists, we can get (2.3)

$\begin{array}{l}\frac{1}{2}{b}_{1}x\left(t-{\tau }_1\right)\le {\left(\frac{1}{4}{b}_1\right)}^2+{\left|x\left(t-{\tau }_1\right)\right|}^2,\\ \frac{1}{2}{b}_{2}y\left(t-{\tau }_2\right)\le {\left(\frac{1}{4}{b}_2\right)}^2+{\left|y\left(t-{\tau }_2\right)\right|}^2.\end{array}$

To sum up, we can get

$\begin{array}{l}\text{d}{V}_2\left(x\left(t\right),y\left(t\right)\right)\\ =\text{d}\left[V_1\left(x\left(t\right),y\left(t\right)\right)+{\displaystyle {\int }_{t-{\tau }_1}^t{\left|x\left(s\right)\right|}^2\text{d}s}+{\displaystyle {\int }_{t-{\tau }_2}^t{\left|y\left(s\right)\right|}^2\text{d}s}\right]\\ \le M\left(x\left(t\right),y\left(t\right)\right)+0.5\left(x^{0.5}\left(t\right)-1\right){\alpha }_1\text{d}{W}_1\left(t\right)+0.5\left(y^{0.5}\left(t\right)-1\right){\alpha }_2\text{d}{W}_2\left(t\right).\end{array}$ (2.4)

Integrating both sides of the above inequality from 0 to ${\tau }_m\wedge T$ and then taking the expectations leads to

$E\left\{{\displaystyle {\int }_0^{{\tau }_u\wedge T}\text{d}\left[V_1\left(x\left(t\right),y\left(t\right)\right)+{\displaystyle {\int }_{t-{\tau }_1}^t{\left|x\left(s\right)\right|}^2\text{d}s}+{\displaystyle {\int }_{t-{\tau }_2}^t{\left|y\left(s\right)\right|}^2\text{d}s}\right]}\right\}\le M^{*}E\left({\tau }_u\wedge T\right),$

So

$\begin{array}{l}E\left[{\displaystyle {\int }_{{\tau }_u\wedge T-{\tau }_1}^{{\tau }_u\wedge T}{\left|x\left(s\right)\right|}^2\text{d}s}+{\displaystyle {\int }_{{\tau }_u\wedge T-{\tau }_2}^{{\tau }_u\wedge T}{\left|y\left(s\right)\right|}^2\text{d}s}\right]+E\left[x\left({\tau }_u\wedge T\right),y\left({\tau }_u\wedge T\right)\right]\\ \le E{\displaystyle {\int }_{-{\tau }_1}^0{\left|x\left(s\right)\right|}^2\text{d}s}+E{\displaystyle {\int }_{-{\tau }_2}^0{\left|y\left(s\right)\right|}^2\text{d}s}+V_1\left(x\left(0\right),y\left(0\right)\right)+M^{*}E\left({\tau }_u\wedge T\right),\end{array}$

Hence

$\begin{array}{l}E\left[x\left({\tau }_u\wedge T\right),y\left({\tau }_u\wedge T\right)\right]\\ \le E{\displaystyle {\int }_{-{\tau }_1}^0{\left|x\left(s\right)\right|}^2\text{d}s}+E{\displaystyle {\int }_{-{\tau }_2}^0{\left|y\left(s\right)\right|}^2\text{d}s}+V_1\left(x\left(0\right),y\left(0\right)\right)+M^{*}T<+\infty .\end{array}$ (2.5)

Set ${\Omega }_u=\left\{{\tau }_u\le T\right\}$ for $u\ge U_1$, then by (2.2), we know $P\left({\Omega }_u\right)\ge \epsilon$. Note that for every $\omega \in {\Omega }_u$, there is at last one of $x\left({\tau }_u,\omega \right),y\left({\tau }_u,\omega \right)$ equal either u or 1/u, then $V_1\left(x\left({\tau }_u\right),y\left({\tau }_u\right)\right)$ is no less then $\min \left\{\left(u-1-\log u\right),\left(\frac{1}{u}-1+\log u\right)\right\}$.

It then follows from (2.2) and (2.5) that

$\begin{array}{l}E{\displaystyle {\int }_{-{\tau }_1}^0{\left|x\left(s\right)\right|}^2\text{d}s}+E{\displaystyle {\int }_{-{\tau }_2}^0{\left|y\left(s\right)\right|}^2\text{d}s}+V_1\left(x\left(0\right),y\left(0\right)\right)+M^{*}E\left({\tau }_u\wedge T\right)\\ \ge E\left[1_{{\Omega }_u\left(\omega \right)}\left(x\left({\tau }_u,\omega \right),y\left({\tau }_u,\omega \right)\right)\right]\ge \epsilon \min \left\{\left(u-1-\log u\right),\left(\frac{1}{u}-1+\log u\right)\right\}.\end{array}$

where $1_{{\Omega }_u\left(\omega \right)}$ is the indicator function of ${\Omega }_u$. Letting $u\to \infty$ leads to the contradiction that

$+\infty >E{\displaystyle {\int }_{-{\tau }_1}^0{\left|x\left(s\right)\right|}^2\text{d}s}+E{\displaystyle {\int }_{-{\tau }_2}^0{\left|y\left(s\right)\right|}^2\text{d}s}+V_1\left(x\left(0\right),y\left(0\right)\right)+M^{*}E\left({\tau }_u\wedge T\right)=+\infty .$

So we must have ${\tau }_{\infty }=\infty ,a.s.$.

3. Stochastic Ultimate Boundedness

Define 3.1. The solution of system (1.3) is random and ultimately bounded, if there exists an any positive constant $H=H\left(\epsilon \right)$ so that for any initial value $\left\{\left(x\left(t\right),y\left(t\right)\right):-\tau \le t\le 0\right\}\in C\left(\left[-\tau ,0\right];R_{+}^2\right)$, it satisfies

$\underset{t\to \infty }{\lim \sup }P\left\{\left|x\left(t\right),y\left(t\right)\right|>H\right\}<\epsilon .$

Lemma 3.1. For any initial value $\left\{\left(x\left(t\right),y\left(t\right)\right):-\tau \le t\le 0\right\}\in C\left(\left[-\tau ,0\right];R_{+}^2\right)$, $\left(x\left(t\right),y\left(t\right)\right)$ is a solution of the system (1.3), there exists positive constants $H\left(\rho \right),0<\rho <1$ satisfies

$\underset{t\to \infty }{\lim \sup }E{\left|\left(x\left(t\right),y\left(t\right)\right)\right|}^{\rho }=H\left(\rho \right).$

Proof. Define $V_3\left(x\right)=x^{\rho }\left(t\right)+y^{\rho }\left(t\right)$, If $\left(x\left(t\right),y\left(t\right)\right)\in R_{+}^2$, we have

$\text{d}{V}_3\left(x\left(t\right),y\left(t\right)\right)=LV\left(x\left(t\right),y\left(t\right)\right)\text{d}t+\rho {\alpha }_1{x}^{\rho }\left(t\right)\text{d}{W}_1\left(t\right)+\rho {\alpha }_2{y}^{\rho }\left(t\right)\text{d}{W}_2\left(t\right).$ (3.1)

where

$\begin{array}{l}L{V}_3\left(x\left(t\right),y\left(t\right)\right)\\ =\rho x^{\rho }\left[a_1-b_{1}x\left(t-{\tau }_1\right)-\frac{c_{1}y\left(t\right)}{1+m_{1}x\left(t\right)+m_{2}y\left(t\right)}\right]+\frac{\rho \left(\rho -1\right)}{2}{\alpha }_1^2{x}^{\rho }\left(t\right)\\ +\rho y^{\rho }\left[a_2-b_{2}y\left(t-{\tau }_2\right)+\frac{c_{2}y\left(t-{\tau }_3\right)}{1+m_{1}x\left(t-{\tau }_3\right)+m_{2}y\left(t-{\tau }_3\right)}\right]+\frac{\rho \left(\rho -1\right)}{2}{\alpha }_2^2{y}^{\rho }\left(t\right)\\ \le a_1\rho x^{\rho }\left(t\right)-\frac{\rho \left(1-\rho \right){\alpha }_1^2}{2}{x}^{\rho }\left(t\right)+a_2\rho y^{\rho }\left(t\right)\\ +\frac{c_2\rho x\left(t-{\tau }_3\right)y^{\rho }\left(t\right)}{1+m_{1}x\left(t-{\tau }_3\right)+m_{2}y\left(t-{\tau }_3\right)}-\frac{\rho \left(1-\rho \right){\alpha }_2^2}{2}{y}^{\rho }\left(t\right).\end{array}$

Because of $0<\rho <1$, $\rho \left(1-\rho \right)>0$, so

$\begin{array}{l}L{V}_3\left(x\left(t\right),y\left(t\right)\right)\\ \le a_1\rho x^{\rho }\left(t\right)-\frac{\rho \left(1-\rho \right){\alpha }_1^2}{2}{x}^{\rho }\left(t\right)+x^{\rho }\left(t\right)+a_2\rho y^{\rho }\left(t\right)-\frac{\rho \left(1-\rho \right){\alpha }_2^2}{2}{y}^{\rho }\left(t\right)\\ +\frac{c_2\rho x\left(t-{\tau }_3\right)y^{\rho }\left(t\right)}{1+m_{1}x\left(t-{\tau }_3\right)+m_{2}y\left(t-{\tau }_3\right)}+y^{\rho }\left(t\right)-V_3\left(x\left(t\right),y\left(t\right)\right)\\ \le H-V_3\left(x\left(t\right),y\left(t\right)\right).\end{array}$

where H is a positive number, substitute it into Equation (3.2) to get

$\begin{array}{l}\text{d}{V}_3\left(x\left(t\right),y\left(t\right)\right)\\ \le \left[H-V_3\left(x\left(t\right),y\left(t\right)\right)\right]\text{d}t+\rho {\alpha }_1{x}^{\rho }\left(t\right)\text{d}{W}_1\left(t\right)+\rho {\alpha }_2{y}^{\rho }\left(t\right)\text{d}{W}_2\left(t\right).\end{array}$

Applying Itô’s formula again, get

$\begin{array}{l}\text{d}\left[{\text{e}}^t{V}_3\left(x\left(t\right),y\left(t\right)\right)\right]\\ ={\text{e}}^t\left[V_3\left(x\left(t\right),y\left(t\right)\right)+\text{d}{V}_3\left(x\left(t\right),y\left(t\right)\right)\right]\\ \le {\text{e}}^{t}H\text{d}t+\rho {\alpha }_1{\text{e}}^t{x}^{\rho }\left(t\right)\text{d}{W}_1\left(t\right)+\rho {\alpha }_2{\text{e}}^t{y}^{\rho }\left(t\right)\text{d}{W}_2\left(t\right).\end{array}$ (3.2)

Taking the expectation of both sides of above inequality (3.2)

${\text{e}}^{t}E{V}_3\left(x\left(t\right),y\left(t\right)\right)\le V_3\left(x\left(0\right),y\left(0\right)\right)+H\left({\text{e}}^t-1\right).$

Namely

$\underset{t\to \infty }{\lim \sup }\text{E}{V}_3\left(x\left(t\right),y\left(t\right)\right)\le H.$

And because

${\left|\left(x\left(t\right),y\left(t\right)\right)\right|}^{\theta }\le 2^{\rho /2}\max \left\{x^{\rho }\left(t\right),y^{\rho }\left(t\right)\right\}\le 2^{\rho /2}{V}_3\left(x\left(t\right),y\left(t\right)\right),$

So

$\underset{t\to \infty }{\lim \sup }\text{E}{\left|\left(x\left(t\right),y\left(t\right)\right)\right|}^{\rho }\le 2^{\rho /2}\underset{t\to \infty }{\lim \sup }\text{E}{V}_3\left(x\left(t\right),y\left(t\right)\right)\le 2^{\rho /2}H.$

Therefore

$\underset{t\to \infty }{\lim \sup }\text{E}{\left|\left(x\left(t\right),y\left(t\right)\right)\right|}^{\rho }\le H_1,$

Among them

$H_1=2^{\rho /2}H.$

Further considering the stochastic ultimate boundedness of the solution, the following propreties hold true.

Theorem 3.1. The solution of system (1.3) is finally bounded by randomness.

Proof. Applying Lemma 3.1, set $\rho =1/2$, then there exists $K>0$, so that

$\underset{t\to \infty }{\lim \sup }\text{E}{\left|\left(x\left(t\right),y\left(t\right)\right)\right|}^{\frac{1}{2}}\le K.$

For any $\epsilon >0$, setting $H_1={\left[K/\epsilon \right]}^{\frac{1}{\rho }}$, an application of Chebyshev’s inequality, there is

$P\left\{\left|\left(x\left(t\right),y\left(t\right)\right)\right|>H_1\right\}<\frac{E\left[{\left|\left(x\left(t\right),y\left(t\right)\right)\right|}^{\rho }\right]}{H_1^{\rho }}\le \epsilon ,$

So

$\underset{t\to \infty }{\lim \sup }P\left\{\left|\left(x\left(t\right),y\left(t\right)\right)\right|>H_1\right\}\le K\cdot \frac{\epsilon }{K}=\epsilon .$

Namely

$\underset{t\to \infty }{\lim \sup }P\left\{\left|\left(x\left(t\right),y\left(t\right)\right)\right|\le H_1\right\}>1-\epsilon .$

which is the desired assertion.

4. Asymptotic Moment Estimation

Theorem 2.1 and Theorem 3.1 show that, for any given initial condition, system (1.3) has a unique global positive solution and the solution is random and finally has upper bounded. The asymptotic moment of the solution is estimated below.

Theorem 4.1. For any given $\theta \in \left(0,1\right)$, there is positive constant $K=K\left(\theta \right)$ such that the solutions of system (1.3) with the initial condition $\left\{\left(x\left(t\right),y\left(t\right)\right):-\tau \le t\le 0\right\}\in C\left(\left[-\tau ,0\right];R_{+}^2\right)$ , have the following property

$\underset{t\to \infty }{\lim \sup }\frac{1}{T}{\displaystyle {\int }_0^{T}E\left[x^{\theta }\left(s\right)+y^{\theta }\left(s\right)\right]\text{d}s}\le K.$ (4.1)

where $K\left(\theta \right)$ is in dependent of the initial value.

Proof. Define a C²-function

$V_4\left(x\left(t\right),y\left(t\right)\right)=x^{\theta }\left(t\right)+y^{\theta }\left(t\right),\left(x\left(t\right),y\left(t\right)\right)\in R_{+}^2.$

By Itô’s formula, one can see that

$\text{d}{V}_4\left(x\left(t\right),y\left(t\right)\right)=L{V}_4\left(x\left(t\right),y\left(t\right)\right)\text{d}t+\theta {\alpha }_1{x}^{\theta }\left(t\right)\text{d}{W}_1\left(t\right)+\theta {\alpha }_2{y}^{\theta }\left(t\right)\text{d}{W}_2\left(t\right).$ (4.2)

where

$\begin{array}{l}L{V}_4\left(x\left(t\right),y\left(t\right)\right)\\ =\theta x^{\theta }\left(t\right)\left[a_1-b_1\left(t-{\tau }_1\right)-\frac{c_{1}y\left(t\right)}{1+m_{1}x\left(t\right)+m_{2}y\left(t\right)}\right]+\frac{\theta \left(\theta -1\right){\alpha }_1^2}{2}{x}^{\theta }\left(t\right)\\ +\theta y^{\theta }\left(t\right)\left[a_2-b_2\left(t-{\tau }_2\right)+\frac{c_{2}x\left(t-{\tau }_3\right)}{1+m_{1}x\left(t-{\tau }_3\right)+m_{2}y\left(t-{\tau }_3\right)}\right]+\frac{\theta \left(\theta -1\right){\alpha }_2^2}{2}{y}^{\theta }\left(t\right)\\ \le a_1\theta x^{\theta }\left(t\right)-\frac{\theta \left(1-\theta \right){\alpha }_1^2}{2}{x}^{\theta }\left(t\right)+a_2\theta y^{\theta }\left(t\right)\\ +\frac{c_2\theta x\left(t-{\tau }_3\right)y^{\theta }\left(t\right)}{1+m_{1}x\left(t-{\tau }_3\right)+m_{2}y\left(t-{\tau }_3\right)}-\frac{\theta \left(1-\theta \right){\alpha }_2^2}{2}{y}^{\theta }\left(t\right).\end{array}$

Therefore

$\begin{array}{l}L{V}_4\left(x\left(t\right),y\left(t\right)\right)+\frac{\theta \left(1-\theta \right){\alpha }_1^2}{4}{x}^{\theta }\left(t\right)+\frac{\theta \left(1-\theta \right){\alpha }_2^2}{4}{y}^{\theta }\left(t\right)\\ \le a_1\theta x^{\theta }\left(t\right)-\frac{\theta \left(1-\theta \right){\alpha }_1^2}{4}{x}^{\theta }\left(t\right)+a_2\theta y^{\theta }\left(t\right)\\ +\frac{c_2\theta x\left(t-{\tau }_3\right)y^{\theta }\left(t\right)}{1+m_{1}x\left(t-{\tau }_3\right)+m_{2}y\left(t-{\tau }_3\right)}-\frac{\theta \left(1-\theta \right){\alpha }_2^2}{4}{y}^{\theta }\left(t\right)\\ \le M.\end{array}$ (4.3)

where M is a positive number, if we take ${\alpha }^2=\min \left\{{\alpha }_1^2,{\alpha }_2^2\right\}$, from (4.1), we can get

$\begin{array}{l}L{V}_4\left(x\left(t\right),y\left(t\right)\right)+\frac{\theta \left(1-\theta \right){\alpha }^2}{4}\left(x^{\theta }\left(t\right)+y^{\theta }\left(t\right)\right)\\ \le L{V}_4\left(x\left(t\right),y\left(t\right)\right)+\frac{\theta \left(1-\theta \right){\alpha }_1^2}{4}{x}^{\theta }\left(t\right)+\frac{\theta \left(1-\theta \right){\alpha }_2^2}{4}{y}^{\theta }\left(t\right)\le M.\end{array}$ (4.4)