1 Introduction

One important component of the predator-prey relation is predator’s functional response, i.e., the rate of prey consumption by an average predator. However, in many cases, when predators have to search for food and, therefore, have to share or compete for food, the functional response in prey-predator model should be predator-dependent. Skalski and Gilliam [1] pointed out that the predator-dependent model can provide better descriptions of predator feeding over a range of predator-prey abundances by comparing the statistical evidence from some predator-prey systems with the three predator-dependent functional responses (Hassell-Varley [2], Bedding-DeAngelis [3] and Crowley-Martin [4]); furthermore, the Bedding-DeAngelis functional response is performed even better. The classical prey-predator model with Bedding-DeAngelis functional response is

$$ \left \{ \begin{array}{@{}l} \frac{dx_{1}}{dt}=x_{1}[r_{1}-a_{11}x_{1}-\frac{a_{12}x_{2}}{1+mx_{1}+nx_{2}}],\\ \frac{dx_{2}}{dt}=x_{2}[-r_{2}+\frac{a_{21}x_{1}}{1+mx_{1}+nx_{2}}-a_{22}x_{2}]. \end{array} \right . $$
(1.1)

There is extensive literature concerned with the dynamics of this prey-predator model with Bedding-DeAngelis functional response, and we here only mention Liu and Yuan [5], Liu and Zhang [6], Zhao and Lv [7], Fan and Kuang [8], Hwang [9, 10], Guo and Wu [11] among many others.

As was pointed by Kuang [12], any model of species dynamics with delays is an approximation at best. More detailed arguments on the significance of time-delays in realistic models may also be found in the classical books of Macdonald [13] and Gopalsamy [14]. Many authors have studied the delay prey-predator system, see, e.g., [15–21]. For the study of delay population systems with Bedding-DeAngelis functional response, see [22–24]. Particularly, we consider the following population system with delays:

$$ \left \{ \begin{array}{@{}l} \frac{dx_{1}}{dt}=x_{1}[r_{1}-a_{11}x_{1}-\frac{a_{12}x_{2}(t-\tau )}{1+mx_{1}+nx_{2}}],\\ \frac{dx_{2}}{dt}=x_{2}[-r_{2}+\frac{a_{21}x_{1}(t-\tau)}{1+mx_{1}+nx_{2}}-a_{22}x_{2}]. \end{array} \right . $$
(1.2)

In system (1.2), \(x_{1}\) and \(x_{2}\) denote the population sizes of the prey and the predator, respectively. We are therefore not only interested in the positive solutions but also require the solutions not to explode at a finite time. To guarantee the positive solutions without explosion (i.e., the global positive solutions), some conditions are in general needed to impose on the coefficients of system (1.2). For example, Kuang [12] discussed the following delay prey-predator system:

$$ \left \{ \begin{array}{@{}l} \frac{dx_{1}}{dt}=x_{1}[r_{1}-a_{11}x_{1}-a_{12}x_{2}(t-\tau)],\\ \frac{dx_{2}}{dt}=x_{2}[-r_{2}+a_{21}x_{1}(t-\tau)-a_{22}x_{2}]. \end{array} \right . $$
(1.3)

He claimed that if \(\triangle_{2}>0\), then system (1.3) has a positive equilibrium \(x^{*} = (x_{1}^{*}, x_{2}^{*} ) = (\triangle_{1}/\triangle,\triangle_{2}/\triangle)\) which is globally asymptotically stable, where \(\triangle=a_{11}a_{22}+a_{12}a_{21}\), \(\triangle_{1} = r_{1}a_{22} + r_{2}a_{12}\), \(\triangle_{2} = r_{1}a_{21}-r_{2}a_{11}\).

However, population models are always affected by environmental noises. Therefore stochastic population models have recently been investigated by many authors; see, e.g., [25–33]. Mao et al. [34] have recently revealed an important fact that the environmental noise can suppress a potential population explosion. Recently, suppose that the parameter \(r_{i}\) is affected by environmental noises with \(r_{i}\rightarrow r_{i}+\sigma_{i}\,dB_{i}\), \(i=1,2\), then corresponding to system (1.3) the authors [35] obtained the following stochastic mode:

$$ \left \{ \begin{array}{@{}l} dx_{1}=x_{1}[r_{1}-a_{11}x_{1}-a_{12}x_{2}(t-\tau_{1})]\,dt+\sigma_{1}x_{1}\,dB_{1},\\ dx_{2}=x_{2}[-r_{2}+a_{21}x_{1}(t-\tau_{2})-a_{22}x_{2}]\,dt+\sigma_{1}x_{2}\,dB_{2} \end{array} \right . $$
(1.4)

and some sufficient and necessary conditions for stability in the mean and extinction of each population for the above stochastic system (1.4).

Motivated by the above work, we therefore wonder if the explosion problem for system (1.2) can be avoided by taking the environmental noise into account instead of imposing conditions on the coefficients of (1.2). To reveal this interesting fact, we stochastically perturb the delay predator-prey model (1.2) into the Itô stochastic differential delay equation

$$ \left \{ \begin{array}{@{}l} dx_{1}(t)=x_{1}(t)[r_{1}-a_{11}x_{1}-\frac{a_{12}x_{2}(t-\tau )}{1+mx_{1}+nx_{2}}]\,dt+\sigma_{1}x_{1}^{2}\,dB_{1}(t),\\ dx_{2}(t)=x_{2}(t)[-r_{2}+\frac{a_{21}x_{1}(t-\tau )}{1+mx_{1}+nx_{2}}-a_{22}x_{2}]\,dt+\sigma_{2}x_{2}^{2}\,dB_{2}(t), \end{array} \right . $$
(SM)

where \(x_{1}\) and \(x_{2}\) represent predator and prey densities at time t, respectively; \(r_{i}\), \(a_{ij}\), τ, m, n are positive constants, \(i,j=1,2\).

In addition, throughout the present paper, let \((\Omega,\mathcal{F},\{\mathcal {F}_{t}\}_{t\geq0},P)\) be a complete probability space with a filtration \(\{\mathcal{F}_{t}\}_{t\geq0}\) satisfying the usual conditions (i.e., it is right continuous and \(\mathcal{F}_{0}\) contains all P-null sets). Let \(|\cdot|\) denote the Euclidean norm in \(R^{n}\). For a given constant \(\tau>0\), let \(C([-\tau,0],R_{+}^{n})\) denote the family of all continuous \(R_{+}^{n}\)-valued functions ξ with its norm \(\|\xi\|=\sup\{|\xi(\theta)|:\theta\in[-\tau,0]\}\), where \(R_{+}=[0,+\infty)\). Also, denote by \(C_{\mathcal {F}_{0}}^{b}([-\tau,0];R_{+}^{n})\) the family of bounded, \(\mathcal {F}_{0}\)-measurable, \(C_{\mathcal{F}_{0}}^{b}([-\tau,0];R_{+}^{n})\)-valued random variables.

Remark 1.1

When \(m=n=0\), system (SM) becomes (1.4), so system (SM) is a more general stochastic system. For system (SM), so far as our knowledge is concerned, the work on a predator-prey model with stochastic perturbations seems rare. In this paper, we study system (SM) which is rather general, and some well-known systems may be viewed as its special cases, and obtain some properties of solutions to system (SM).

Remark 1.2

System (SM) is based on assuming that the noise affects parameters \(a_{11}\) and \(a_{22}\). In fact, the noise may affect other parameters in (SM), which results in other types of stochastic models, which are the future research topics; for more details, see [36].

Remark 1.3

For the analysis of population dynamical problems of (SM), two difficult issues arise: (i) how to handle the delays in the given model and (ii) how to handle the nonlinear terms in (SM). To deal with these problems, the construction of a Lyapunov functional V is quite crucial, and it is introduced in Sections 3 and 4.

2 Positive and global solutions

In order for a stochastic differential delay equation to have a unique global (i.e., no explosion in a finite time) solution for any given initial data, the coefficients of the equation are generally required to satisfy the linear growth condition and local Lipschitz condition [36]. However, the coefficients of system (SM) do not satisfy the linear growth condition, though they are locally Lipschitz continuous, so the solution of system (SM) may explode at a finite time. In this section we shall show that under simple hypothesis the solution of system (SM) is not only positive but will also not explode to infinity at any finite time.

Theorem 2.1

For any given initial data \(\{(x_{1}(\theta),x_{2}(\theta))^{\top}:-\tau\leq\theta\leq0\}=\xi\in C_{\mathcal{F}_{0}}^{b}([-\tau,0]; R_{+}^{0}\times R_{+}^{0})\), where \(R_{+}^{0}=(0,+\infty)\). If \(a_{21}\leq4\), there is a unique positive local solution \((x_{1}(t),x_{2}(t))\) to (SM) on \(t\geq-\tau\) with satisfying initial condition ξ and the solution will remain in \(R_{+}^{2}\) with probability 1.

Proof

By the biological meaning, we only focus on the positive solution to system (SM). Thus it is reasonable to make the following change of variables, \(x_{1}(t)=e^{u(t)}\), \(x_{2}(t)=e^{v(t)}\). By using Itô’s formula, system (SM) can be reformulated in the following form:

$$ \left \{ \begin{array}{@{}l} du=[r_{1}-a_{11}e^{u}-a_{12}\frac{e^{v(t-\tau)}}{1+e^{mu}+e^{nv}}-\frac {1}{2}\sigma_{1}^{2}]\,dt+\sigma_{1}e^{u}\,dB_{1}(t), \\ dv=[-r_{2}(t)+a_{21}\frac{e^{u(t-\tau )}}{1+e^{mu}+e^{nv}}-a_{22}e^{v}-\frac{1}{2}\sigma_{2}^{2}]\,dt+\sigma _{2}e^{v}\,dB_{2}(t), \quad t\geq0, \\ u(\theta)=\ln x_{1}(\theta), v(\theta)=\ln x_{2}(\theta),\quad \theta\in[-\tau,0]. \end{array} \right . $$
(2.1)

It is easy to see that the coefficients of (2.1) satisfy the local Lipschitz condition, then for any given initial values \(u(\theta)>0\), \(v(\theta)>0\), \(\theta\in[-\tau,0]\), there is a unique maximal local solution \(u(t)\), \(v(t)\) on \([-\tau,\tau_{e})\), where \(\tau_{e}\) is explosion time. To show that this solution is global, we need to show that \(\tau_{e}=\infty \) a.s. Let \(n_{0}>0\) be sufficiently large for

$$\frac{1}{n_{0}}< \min_{-\tau\leq t\leq0}\bigl|x(t)\bigr|\leq\max _{-\tau\leq t\leq 0}\bigl|x(t)\bigr|<n_{0}, $$

where \(x(t)=(x_{1}(t),x_{2}(t))^{\top}\). For each integer \(n>n_{0}\), define the stopping times:

$$\tau_{n}=\inf\biggl\{ t\in[0,\tau_{e}]:x_{1}(t) \notin\biggl(\frac{1}{n},n\biggr) \mbox{ or } x_{2}(t)\notin \biggl(\frac{1}{n},n\biggr)\biggr\} . $$

Throughout this paper, we set \(\inf\emptyset=\infty\). Obviously, \(\tau_{n}\) is increasing as \(n\rightarrow\infty\). Let \(\tau_{\infty}=\lim_{n\rightarrow}\tau_{n}\), whence \(\tau_{\infty}\leq\tau_{e}\) a.s. If we can show that \(\tau_{\infty}=\infty\) a.s., then \(\tau_{e}=\infty\) a.s. and \(x_{i}(t)\geq0\) a.s., \(i=1,2\) for all \(t\geq0\). To show this statement, let us define \(C^{2}\)-function \(V:R_{+}^{2}\rightarrow R_{+}\) by

$$V(x_{1},x_{2})=(\sqrt{x_{1}}-1-0.5\ln x_{1})+(\sqrt{x_{2}}-1-0.5\ln x_{2}). $$

The nonnegativity of this function can be obtained from \(u-1-\ln u\geq0\) on \(u>0\). Let \(n\geq n_{0}\) and \(T>0\) be arbitrary. For \(0\leq t\leq\tau_{n}\wedge T\), we can use Itô’s formula to \(\int_{t-\tau}^{t}(x_{1}^{2}(s)+x_{2}^{2}(s))\,ds+V(x_{1},x_{2})\) to obtain that

$$\begin{aligned} &d \biggl[\int_{t-\tau}^{t} \bigl(x_{1}^{2}(s)+x_{2}^{2}(s) \bigr)\,ds+V(x_{1},x_{2}) \biggr] \\ &\quad=\bigl[x_{1}^{2}-x_{1}^{2}(t- \tau) \\ &\qquad{}+0.5\bigl(x_{1}^{0.5}-1\bigr)F(x_{1},x_{2}) \bigr]\,dt +\frac{1}{8}\bigl(-x_{1}^{0.5}+2\bigr) \sigma_{1}^{2}x_{1}^{2}\,dt \\ &\qquad{}+0.5\bigl(x_{1}^{0.5}-1\bigr)\sigma_{1}x_{1}\,dB_{1}(t) +\bigl[x_{2}^{2}-x_{2}^{2}(t-\tau) +0.5 \bigl(x_{2}^{0.5}-1\bigr)G(x_{1},x_{2}) \bigr]\,dt \\ &\qquad{}+\frac{1}{8}\bigl(-x_{2}^{0.5}+2\bigr) \sigma_{2}^{2}x_{2}^{2}\,dt+0.5 \bigl(x_{2}^{0.5}-1\bigr)\sigma _{2}x_{2}\,dB_{2}(t), \end{aligned}$$
(2.2)

where

$$\begin{aligned}& F(x_{1},x_{2})=r_{1}-a_{11}x_{1}- \frac{a_{12}x_{2}(t-\tau)}{1+mx_{1}+nx_{2}}, \\& G(x_{1},x_{2})=-r_{2}-a_{22}x_{2}+ \frac{a_{21}x_{1}(t-\tau)}{1+mx_{1}+nx_{2}}. \end{aligned}$$

Compute

$$\begin{aligned} &\bigl(x_{1}^{0.5}-1\bigr)F(x_{1},x_{2}) \\ &\quad= \bigl(x_{1}^{0.5}-1\bigr)\biggl[r_{1}-a_{11}x_{1}- \frac {a_{12}x_{2}(t-\tau)}{1+mx_{1}+nx_{2}}\biggr] \\ &\quad=r_{1}\bigl(x_{1}^{0.5}-1 \bigr)-a_{11}x_{1}^{1.5}+a_{11}x_{1}+ \frac{a_{12}x_{2}(t-\tau )}{1+mx_{1}+nx_{2}}-\frac{a_{12}x_{1}^{0.5}x_{2}(t-\tau)}{1+mx_{1}+nx_{2}} \\ &\quad\leq r_{1}\bigl(x_{1}^{0.5}-1 \bigr)-a_{11}x_{1}^{1.5}+a_{11}x_{1}+a_{12}x_{2}(t- \tau) \end{aligned}$$
(2.3)

and

$$\begin{aligned} &\bigl(x_{2}^{0.5}-1\bigr)G(x_{1},x_{2}) \\ &\quad= \bigl(x_{2}^{0.5}-1\bigr)\biggl[-r_{2}-a_{22}x_{2}+ \frac {a_{21}x_{1}(t-\tau)}{1+mx_{1}+nx_{2}}\biggr] \\ &\quad=-r_{2}\bigl(x_{2}^{0.5}-1 \bigr)-a_{22}x_{2}^{1.5}+a_{22}x_{2} +\frac{a_{21}x_{2}^{0.5}x_{1}(t-\tau)}{1+mx_{1}+nx_{2}}-\frac{a_{21}x_{1}(t-\tau )}{1+mx_{1}+nx_{2}} \\ &\quad\leq -r_{2}\bigl(x_{2}^{0.5}-1 \bigr)-a_{22}x_{2}^{1.5}+a_{22}x_{2}+a_{21}x_{2}^{0.5}x_{1}(t- \tau ) \\ &\quad\leq -r_{2}\bigl(x_{2}^{0.5}-1 \bigr)-a_{22}x_{2}^{1.5}+(a_{22}+0.5a_{21})x_{2}+0.5a_{21}x_{1}^{2}(t- \tau). \end{aligned}$$
(2.4)

Substituting (2.3) and (2.4) into (2.2) yields

$$\begin{aligned} &d \biggl[\int_{t-\tau}^{t} \bigl(x_{1}^{2}(s)+x_{2}^{2}(s) \bigr)\,ds+V(x_{1},x_{2}) \biggr] \\ &\quad\leq H(x_{1},x_{2})\,dt+0.5\bigl(x_{1}^{0.5}-1 \bigr)\sigma_{1}x_{1}\,dB_{1}(t)+0.5 \bigl(x_{2}^{0.5}-1\bigr)\sigma _{2}x_{2}\,dB_{2}(t), \end{aligned}$$
(2.5)

where

$$\begin{aligned} H(x_{1},x_{2}) =&x_{1}^{2}-x_{1}^{2}(t- \tau )+0.5r_{1}\bigl(x_{1}^{0.5}-1 \bigr)-0.5a_{11}x_{1}^{1.5}+0.5a_{11}x_{1} \\ &{}+0.5a_{12}x_{2}(t-\tau)+\frac{1}{8} \bigl(-x_{1}^{0.5}+2\bigr)\sigma_{1}^{2}x_{1}^{2} \\ &{}+0.5x_{2}^{2}-0.5x_{2}^{2}(t- \tau)-0.5r_{2}\bigl(x_{2}^{0.5}-1\bigr) \\ &{}-0.5a_{22}x_{2}^{1.5}+0.5(a_{22}+0.5a_{21})x_{2} \\ &{}+0.25a_{21}x_{1}^{2}(t-\tau)+ \frac{1}{8}\bigl(-x_{2}^{0.5}+2\bigr) \sigma_{2}^{2}x_{2}^{2} \\ =&\bigl(1+0.25\sigma_{1}^{2}\bigr)x_{1}^{2}+(0.25a_{21}-1)x_{1}^{2}(t- \tau) \\ &{}+0.5r_{1}\bigl(x_{1}^{0.5}-1 \bigr)-0.5a_{11}x_{1}^{1.5}+0.5a_{11}x_{1} \\ &{}+\bigl(0.5+0.25\sigma_{2}^{2}\bigr)x_{2}^{2}-0.5x_{2}^{2}(t- \tau)-0.5r_{2}\bigl(x_{2}^{0.5}-1\bigr) \\ &{}-0.5a_{22}x_{2}^{1.5}+0.5(a_{22}+0.5a_{21})x_{2} \\ &{}-0.125x_{1}^{2.5}-0.125x_{2}^{2.5}. \end{aligned}$$

From \(a_{21}\leq4\), it is easy to see that \(H(x_{1},x_{2})\) is bounded, say by K. By (2.5), we have

$$\begin{aligned} &d \biggl[\int_{t-\tau}^{t}\bigl(x_{1}^{2}(s)+x_{2}^{2}(s) \bigr)\,ds+V(x_{1},x_{2}) \biggr] \\ &\quad\leq K\,dt+0.5\bigl(x_{1}^{0.5}-1\bigr) \sigma_{1}x_{1}\,dB_{1}(t)+0.5\bigl(x_{2}^{0.5}-1 \bigr)\sigma_{2}x_{2}\,dB_{2}(t). \end{aligned}$$

Integrating both sides from 0 to \(\tau_{n}\wedge T\), and then taking expectations, yields

$$\begin{aligned} &E \biggl[\int_{\tau_{n}\wedge T-\tau}^{\tau_{n}\wedge T}\bigl(x_{1}^{2}(s)+x_{2}^{2}(s) \bigr)\,ds+V\bigl(x_{1}(\tau_{n}\wedge T),x_{2}( \tau_{n}\wedge T)\bigr) \biggr] \\ &\quad\leq\int_{-\tau}^{0}\bigl(x_{1}^{2}(s)+x_{2}^{2}(s) \bigr)\,ds+V\bigl(x_{1}(0),x_{2}(0)\bigr)+KT. \end{aligned}$$

Consequently,

$$ EV\bigl(x_{1}(\tau_{n}\wedge T),x_{2}( \tau_{n}\wedge T)\bigr)\leq\int_{-\tau}^{0} \bigl(x_{1}^{2}(s)+x_{2}^{2}(s)\bigr)\,ds+V \bigl(x_{1}(0),x_{2}(0)\bigr)+KT. $$
(2.6)

Note that for each \(\omega\in\Omega_{n}=\{\tau_{n}\leq T\}\), there is some i such that \(x_{i}(\tau_{n},\omega)\) equals n or \(\frac{1}{n}\) for \(i=1,2\). Hence \(V(x_{1}(\tau_{n}\wedge T),x_{2}(\tau_{n}\wedge T))\) is no less than

$$\min\{\sqrt{n}-1-0.5\ln n,\sqrt{1/n}-1-0.5\ln1/n\}. $$

By (2.6) we have

$$\begin{aligned} &\int_{-\tau}^{0}\bigl(x_{1}^{2}(s)+x_{2}^{2}(s) \bigr)\,ds+V\bigl(x_{1}(0),x_{2}(0)\bigr)+KT \\ &\quad\geq E\bigl[1_{\Omega_{n}}(\omega)V\bigl(x_{1}( \tau_{n}),x_{2}(\tau_{n})\bigr)\bigr] \\ &\quad\geq P\{\tau_{n}\leq T\}\min\{\sqrt{n}-1-0.5\ln n, \sqrt{1/n}-1-0.5\ln1/n\}, \end{aligned}$$

where \(1_{\Omega_{n}}\) is the indicator function of \(\Omega_{n}\). Letting \(n\rightarrow\infty\), leads to the contraction

$$\lim_{k\rightarrow\infty}P\{\tau_{k}\leq T\}=0. $$

Since \(T>0\) is arbitrary, we must have

$$P\{\tau_{\infty}< \infty\}=0. $$

So \(P\{\tau_{\infty}= \infty\}=1\) as required. The proof is completed. □

Remark 2.1

It is well known that systems (1.2) and (1.3) may explode to infinity at a finite time for some system parameters, see [37]. However, the explosion will no longer happen as long as there is noise. In other words, Theorem 2.1 reveals the important property that the environmental noise suppresses the explosion for the delay equation.

3 Stochastically ultimate boundedness

One of the important issues in the study of population systems is the stochastically ultimate boundedness. System (SM) is said to be stochastically ultimately bounded if for any \(\varepsilon \in(0, 1)\), there is a positive constant \(H=H(\varepsilon)\) such that for any initial data \(\xi\in C_{\mathcal {F}_{0}}^{b}([-\tau,0];R_{+}^{0}\times R_{+}^{0})\), the solution \((x_{1},x_{2})^{\top}\) of system (SM) has the property that

$$\lim_{t\rightarrow\infty}\sup P\bigl\{ \bigl|x(t)\bigr|\leq H\bigr\} \geq1-\varepsilon, $$

where \(x=(x_{1},x_{2})^{\top}\).

In this section we shall investigate the stochastically ultimate boundedness of system (SM). The following theorem gives a sufficient criterion for the stochastically ultimate boundedness of population.

Lemma 3.1

Let \(\theta\in(0, 1)\) and \(\theta a_{21}\geq n\). Then there is a positive constant \(K=K(\theta)\), which is independent of the initial data \(\xi\in C_{\mathcal{F}_{0}}^{b}([-\tau,0];R_{+}^{0}\times R_{+}^{0})\), such that the solution \(x(t)\) of system (SM) has the property that

$$\lim_{t\rightarrow\infty}\sup E|x|^{\theta}\leq K. $$

Proof

Define

$$V(x_{1},x_{2})=x_{1}^{\theta}+x_{2}^{\theta}. $$

Applying Itô’s formula to system (SM), we have

$$\begin{aligned} dV(x_{1},x_{2}) =&\biggl[\theta x_{1}^{\theta}\biggl(r_{1}-a_{11}x_{1}- \frac{a_{12}x_{2}(t-\tau)}{1+mx_{1}+nx_{2}}\biggr)-\frac {\theta(1-\theta)\sigma_{1}^{2}}{2}x_{1}^{2+\theta} \biggr]\,dt \\ &{}+\theta\sigma_{1}x_{1}^{1+\theta}\,dB_{1}(t) \\ &{}+\biggl[\theta x_{2}^{\theta}\biggl(-r_{2}+ \frac{a_{21}x_{1}(t-\tau)}{1+mx_{1}+nx_{2}}-a_{22}x_{2}\biggr)-\frac {\theta(1-\theta)\sigma_{2}^{2}}{2}x_{2}^{2+\theta} \biggr]\,dt \\ &{}+\theta\sigma_{2}x_{2}^{1+\theta }\,dB_{2}(t). \end{aligned}$$
(3.1)

Denote

$$\begin{aligned} LV(x_{1},x_{2},y_{1},y_{2}) =& \theta x_{1}^{\theta}\biggl(r_{1}-a_{11}x_{1}- \frac{a_{12}y_{2}}{1+mx_{1}+nx_{2}}\biggr)-\frac{\theta (1-\theta)\sigma_{1}^{2}}{2}x_{1}^{2+\theta} \\ &{}+\theta x_{2}^{\theta}\biggl(-r_{2}+ \frac{a_{21}y_{1}}{1+mx_{1}+nx_{2}}-a_{22}x_{2}\biggr)-\frac{\theta (1-\theta)\sigma_{2}^{2}}{2}x_{2}^{2+\theta} \\ \leq&\theta x_{1}^{\theta}r_{1}- \frac{\theta(1-\theta)\sigma_{1}^{2}}{2}x_{1}^{2+\theta}-\frac{\theta (1-\theta)\sigma_{2}^{2}}{2}x_{2}^{2+\theta}+ \frac{\theta a_{21}}{n}|y|^{2} \\ =&F(x_{1},x_{2})-V(x_{1},x_{2})-e^{\tau}|x|^{2}+ \frac{\theta a_{21}}{n}|y|^{2}, \end{aligned}$$

where

$$\begin{aligned}& y=(y_{1},y_{2})^{\top}=\bigl(x_{1}(t- \tau),x_{2}(t-\tau)\bigr)^{\top},\qquad x=(x_{1},x_{2})^{\top}, \\& F(x_{1},x_{2})=\theta x_{1}^{\theta}r_{1}-\frac{\theta(1-\theta)\sigma_{1}^{2}}{2}x_{1}^{2+\theta}- \frac{\theta (1-\theta)\sigma_{2}^{2}}{2}x_{2}^{2+\theta} +x_{1}^{\theta}+x_{2}^{\theta}+e^{\tau}|x|^{2}. \end{aligned}$$

Note that \(F(x_{1},x_{2})\) is bounded in \(R_{+}^{2}\), namely

$$F(x_{1},x_{2})\leq K_{1}, \quad \forall x\in R_{+}^{2}. $$

Hence we have

$$LV(x_{1},x_{2},y_{1},y_{2})\leq K_{1}-V(x_{1},x_{2})-e^{\tau}|x|^{2}+ \frac{\theta a_{21}}{n}|y|^{2}. $$

Substituting this into (3.1) gives

$$\begin{aligned} dV(x_{1},x_{2}) \leq&\biggl[K_{1}-V(x_{1},x_{2})-e^{\tau}|x|^{2}+ \frac{\theta a_{21}}{n}|y|^{2}\biggr]\,dt \\ &{}+\theta\sigma_{1}x_{1}^{1+\theta}\,dB_{1}(t)+ \theta\sigma_{2}x_{2}^{1+\theta }\,dB_{2}(t). \end{aligned}$$
(3.2)

From (3.2) and once again by Itô’s formula, we have

$$\begin{aligned} d\bigl[e^{t}V(x_{1},x_{2})\bigr]={}&e^{t} \bigl[V(x_{1},x_{2})\,dt+dV(x_{1},x_{2}) \bigr] \\ \leq{}& e^{t}\biggl[K_{1}-e^{\tau}|x|^{2}+ \frac{\theta a_{21}}{n}\bigl|x(t-\tau)\bigr|^{2}\biggr]\,dt \\ &{}+e^{t}\theta\sigma_{1}x_{1}^{1+\theta}\,dB_{1}(t) +e^{t}\theta\sigma_{2}x_{2}^{1+\theta}\,dB_{2}(t). \end{aligned}$$

From \(\theta a_{21}\geq n\), we hence derive that

$$\begin{aligned} e^{t}EV(x_{1},x_{2}) \leq& V\bigl(x_{1}(0),x_{2}(0) \bigr)+K_{1}e^{t}-E\int_{0}^{t}e^{s+\tau}\bigl|x(s)\bigr|^{2}\,ds \\ &{}+\frac{\theta a_{21}}{n}E\int_{0}^{t}e^{s}\bigl|x(s- \tau)\bigr|^{2}\,ds \\ =&V\bigl(x_{1}(0),x_{2}(0)\bigr)+K_{1}e^{t}-E \int_{0}^{t}e^{s\tau}\bigl|x(s)\bigr|^{2}\,ds \\ &{}+\frac{\theta a_{21}}{n}E\int_{-\tau}^{t-\tau}e^{s}\bigl|x(s)\bigr|^{2}\,ds \\ \leq& V\bigl(x_{1}(0),x_{2}(0)\bigr)+K_{1}e^{t}+ \frac{\theta a_{21}}{n}\int_{-\tau}^{0}e^{s}\bigl|x(s)\bigr|^{2}\,ds. \end{aligned}$$

This implies immediately that

$$\lim_{t\rightarrow\infty}\sup EV(x_{1},x_{2})\leq K_{1}. $$

On the other hand, we have

$$|x|^{2}\leq2\max\{x_{1},x_{2}\}. $$

Thus

$$|x|^{\theta}\leq2^{\theta/2}\max\bigl\{ x_{1}^{\theta},x_{2}^{\theta}\bigr\} \leq2^{\theta /2}V(x_{1},x_{2}). $$

Hence, we have

$$\lim_{t\rightarrow\infty}\sup E|x|^{\theta}\leq2^{\theta/2}K_{1}:=K. $$

The proof is completed. □

Theorem 3.1

Let \(\theta\in(0, 1)\) and \(\theta a_{21}\geq n\). System (SM) is stochastically ultimately bounded.

Proof

From Lemma 3.1, there is \(K>0\) such that

$$\lim_{t\rightarrow\infty}\sup E|x|^{\theta}\leq K. $$

Then, for any \(\varepsilon>0\), let \(H=K^{2}/\varepsilon^{2}\). Then by Chebyshev’s inequality

$$P\bigl\{ \bigl|x(t)\bigr|>H\bigr\} \leq\frac{E(\sqrt{|x(t)|})}{\sqrt{H}}. $$

Thus

$$\lim_{t\rightarrow\infty}\sup P\bigl\{ \bigl|x(t)\bigr|>H\bigr\} \leq\frac{K}{\sqrt {H}}= \varepsilon. $$

This implies

$$\lim_{t\rightarrow\infty}\sup P\bigl\{ \bigl|x(t)\bigr|\leq H\bigr\} \geq1-\varepsilon. $$

The proof is completed. □

4 Almost surely asymptotic properties

In this section, we study the pathwise properties of system (SM).

Theorem 4.1

Let the same conditions of Theorem  2.1 hold. For any given initial data \(\{(x_{1}(\theta),x_{2}(\theta))^{\top}:-\tau\leq\theta\leq0\}=\xi\in C_{\mathcal{F}_{0}}^{b}([-\tau,0];R_{+}^{0}\times R_{+}^{0})\), the solution x(t) of system (SM) has the property that

$$\lim_{t\rightarrow\infty}\sup\frac{\ln|x(t)|}{\ln t}\leq1 \quad \textit{a.s.} $$

Proof

Define

$$V(x)=x_{1}+x_{2}. $$

It is easy to see that

$$\begin{aligned} dV\bigl(x(t)\bigr)={}&x_{1}(t)\biggl[r_{1}-a_{11}x_{1}- \frac{a_{12}x_{2}(t-\tau )}{1+mx_{1}+nx_{2}}\biggr]\,dt+\sigma_{1}x_{1}^{2}\,dB_{1}(t) \\ &{}+x_{2}(t)\biggl[-r_{2}+\frac{a_{21}x_{1}(t-\tau)}{1+mx_{1}+nx_{2}}-a_{22}x_{2} \biggr]\,dt+\sigma _{2}x_{2}^{2}\,dB_{2}(t). \end{aligned}$$

Let \(\gamma>0\) be arbitrary. By Itô’s formula, we have

$$\begin{aligned} e^{\gamma t}\ln\bigl(V\bigl(x(t)\bigr)\bigr) ={}&\ln\bigl(V\bigl(x(0)\bigr)\bigr) \\ &{}+\int_{0}^{t}e^{\gamma s} \biggl( \frac{x_{1}}{V(x(s))} \biggl[r_{1}-a_{11}x_{1}- \frac{a_{12}x_{2}(t-\tau )}{1+mx_{1}+nx_{2}} \biggr] -\frac{\sigma_{1}^{2}x_{1}^{4}}{2V^{2}(x(s))} \biggr)\,ds \\ &{}+\int_{0}^{t}e^{\gamma s} \biggl( \frac{x_{2}}{V(x(s))} \biggl[-r_{2}+\frac{a_{21}x_{1}(t-\tau )}{1+mx_{1}+nx_{2}}-a_{22}x_{2} \biggr] -\frac{\sigma_{2}^{2}x_{2}^{4}}{2V^{2}(x(s))} \biggr)\,ds \\ &{}+\gamma\int_{0}^{t}e^{\gamma s}\ln \bigl(V\bigl(x(s)\bigr)\bigr)\,ds+M_{1}(t)+M_{2}(t), \end{aligned}$$
(4.1)

where

$$M_{i}(t)=\int_{0}^{t} \frac{e^{\gamma s}\sigma_{i}x_{i}^{2}}{V(x(s))}\,dB_{i}(s),\quad i=1,2, $$

is a real-valued continuous local martingale vanishing at \(t=0\) and its quadratic form is given by

$$\bigl\langle M_{i}(t),M_{i}(t)\bigr\rangle =\int _{0}^{t}\frac{e^{2\gamma s}\sigma_{i}^{2}x_{i}^{4}}{V(x(s))}\,ds. $$

Let \(\varepsilon\in(0,1)\) and \(\theta>1\) be arbitrary. By the exponential martingale inequality (see, e.g., [36]), for each \(k\geq1\),

$$P \biggl\{ \sup_{0\leq t\leq k} \biggl[M_{i}(t)- \frac{\varepsilon}{2}e^{-\gamma k}\bigl\langle M_{i}(t),M_{i}(t) \bigr\rangle \biggr]>\frac{\theta e^{\gamma k}}{\varepsilon}\ln k \biggr\} \leq k^{-\theta}. $$

Since the series \(\sum_{k=1}^{\infty}k^{-\theta}\) converges, the well-known Borel-Cantelli lemma yields that there exists \(\Omega_{0}\subset\Omega\) with \(P(\Omega_{0})=1\) such that for any \(\omega\in\Omega_{0}\) there exists an integer \(k_{0}=k_{0}(\omega)\) such that

$$ M_{i}(t)\leq\frac{\varepsilon}{2}e^{-\gamma k}\bigl\langle M_{i}(t),M_{i}(t)\bigr\rangle +\frac{\theta e^{\gamma k}}{\varepsilon}\ln k $$
(4.2)

for all \(0\leq t\leq k\) and \(k\geq k_{0}(\omega)\). Substituting this into (4.1) and then using (4.2), we have

$$\begin{aligned} e^{\gamma t}\ln\bigl(V\bigl(x(t)\bigr)\bigr) ={}&\ln\bigl(V\bigl(x(0)\bigr)\bigr) \\ &{}+\int_{0}^{t}e^{\gamma s} \biggl( \frac{x_{1}}{V(x(s))} \biggl[r_{1}-a_{11}x_{1}- \frac{a_{12}x_{2}(t-\tau )}{1+mx_{1}+nx_{2}} \biggr] -\frac{(1-\varepsilon)\sigma_{1}^{2}x_{1}^{4}}{2V^{2}(x(s))} \biggr)\,ds \\ &{}+\int_{0}^{t}e^{\gamma s} \biggl( \frac{x_{2}}{V(x(s))} \biggl[-r_{2}+\frac{a_{21}x_{1}(t-\tau )}{1+mx_{1}+nx_{2}}-a_{22}x_{2} \biggr] -\frac{(1-\varepsilon)\sigma_{2}^{2}x_{2}^{4}}{2V^{2}(x(s))} \biggr)\,ds \\ &{}+\gamma\int_{0}^{t}e^{\gamma s}\ln \bigl(V\bigl(x(s)\bigr)\bigr)\,ds+\frac{\theta e^{\gamma k}}{\varepsilon}\ln k \end{aligned}$$
(4.3)

for all \(0\leq t\leq k\), \(k\geq k_{0}(\omega)\) and \(\omega\in\Omega_{0}\). Compute

$$\begin{aligned}& \frac{x_{1}}{V(x(s))}\biggl[r_{1}-a_{11}x_{1}- \frac{a_{12}x_{2}(t-\tau )}{1+mx_{1}+nx_{2}}\biggr] \leq r_{1}+a_{12}x_{2}(t-\tau), \end{aligned}$$
(4.4)
$$\begin{aligned}& \frac{x_{2}}{V(x(s))}\biggl[-r_{2}+\frac{a_{21}x_{1}(t-\tau )}{1+mx_{1}+nx_{2}}-a_{22}x_{2} \biggr]\leq r_{2}+a_{21}x_{1}(t- \tau), \end{aligned}$$
(4.5)
$$\begin{aligned}& \frac{\sigma_{1}^{2}x_{1}^{4}}{2V^{2}(x(s))}\geq\frac{\sigma_{1}^{2}k_{1}^{2}}{2}x_{1}^{2},\qquad \frac{\sigma_{2}^{2}x_{2}^{4}}{2V^{2}(x(s))}\geq\frac{\sigma _{2}^{2}k_{1}^{2}}{2}x_{2}^{2}, \end{aligned}$$
(4.6)

where \(k_{1},k_{2}\in(0,1)\) are positive constants. Substituting (4.4)-(4.6) into (4.3) gives

$$\begin{aligned} &e^{\gamma t}\ln\bigl(V\bigl(x(t)\bigr)\bigr) \\ &\quad=\ln\bigl(V\bigl(x(0)\bigr)\bigr)+\gamma\int_{0}^{t}e^{\gamma s} \ln\bigl(V\bigl(x(s)\bigr)\bigr)\,ds+\frac{\theta e^{\gamma k}}{\varepsilon}\ln k \\ &\qquad{}+\int_{0}^{t}e^{\gamma s} \biggl(r_{1}+a_{12}x_{2}(t-\tau)+r_{2}+a_{21}x_{1}(t- \tau)] -\frac{(1-\varepsilon)\sigma_{1}^{2}k_{1}^{2}x_{1}^{2}}{2} \\ &\qquad{}-\frac{(1-\varepsilon)\sigma_{2}^{2}k_{2}^{2}x_{2}^{2}}{2} \biggr)\,ds \end{aligned}$$
(4.7)

for all \(0\leq t\leq k\), \(k\geq k_{0}(\omega)\) and \(\omega\in\Omega_{0}\). For \(i=1,2\), note

$$\int_{0}^{t}e^{\gamma s}x_{i}(s- \tau)\,ds \leq\int_{-\tau}^{t-\tau}e^{\gamma (s+\tau)}x_{i}(s)\,ds \leq\int_{-\tau}^{0}e^{\gamma\tau}x_{i}(s)\,ds+ \int_{0}^{t}e^{\gamma (s+\tau)}x_{i}(s)\,ds. $$

We can rewrite (4.7) as

$$\begin{aligned} e^{\gamma t}\ln\bigl(V\bigl(x(t)\bigr)\bigr) ={}&\ln \bigl(V\bigl(x(0)\bigr)\bigr)+a_{12}e^{\gamma\tau}\int _{-\tau}^{0}x_{1}(s)\,ds +a_{21}e^{\gamma\tau}\int_{-\tau}^{0}x_{2}(s)\,ds+ \frac{\theta e^{\gamma k}}{\varepsilon}\ln k \\ &{}+\int_{0}^{t}e^{\gamma s} \biggl(r_{1}+a_{12}e^{\gamma \tau}x_{2}+r_{2}+a_{21}e^{\gamma\tau}x_{1} -\frac{(1-\varepsilon)\sigma_{1}^{2}k_{1}^{2}x_{1}^{2}}{2} \\ &{}-\frac{(1-\varepsilon)\sigma_{2}^{2}k_{2}^{2}x_{2}^{2}}{2} \biggr)\,ds +\gamma\int_{0}^{t}e^{\gamma s} \ln\bigl(V\bigl(x(s)\bigr)\bigr)\,ds. \end{aligned}$$
(4.8)

Obviously, the following polynomial is bounded by a positive constant, say \(K_{1}\),

$$ \gamma\ln V(x)+r_{1}+a_{12}e^{\gamma \tau}x_{2}+r_{2}+a_{21}e^{\gamma\tau}x_{1} -\frac{(1-\varepsilon)\sigma_{1}^{2}k_{1}^{2}x_{1}^{2}}{2}-\frac{(1-\varepsilon )\sigma_{2}^{2}k_{2}^{2}x_{2}^{2}}{2}\leq K_{1}. $$
(4.9)

Then by (4.8) and (4.9) we get

$$e^{\gamma t}\ln\bigl(V\bigl(x(t)\bigr)\bigr)\leq C_{1}+ \frac{\theta e^{\gamma k}}{\varepsilon}\ln k+\frac{K_{1}}{\gamma}e^{\gamma t} $$

for all \(0\leq t\leq k\), \(k\geq k_{0}(\omega)\) and \(\omega\in\Omega_{0}\), where

$$C_{1}=\ln\bigl(V\bigl(x(0)\bigr)\bigr)+a_{12}e^{\gamma\tau} \int_{-\tau}^{0}x_{1}(s)\,ds +a_{21}e^{\gamma\tau}\int_{-\tau}^{0}x_{2}(s)\,ds. $$

Consequently, for any \(\omega\in\Omega_{0}\), if \(k-1\leq t\leq k\) and \(k\geq k_{0}(\omega)\), we have

$$\frac{\ln(V(x(t)))}{\ln t}\leq\frac{1}{\ln(k-1)} \biggl[e^{-\gamma (k-1)}C_{1}+ \frac{\theta e^{\gamma}}{\varepsilon}\ln k+\frac{K_{1}}{\gamma} \biggr]. $$

This implies

$$\lim_{t\rightarrow\infty}\sup\frac{\ln(V(x(t)))}{\ln t}\leq\frac{\theta e^{\gamma}}{\varepsilon} \quad \mbox{a.s.} $$

Letting \(\varepsilon\rightarrow1\), \(\theta\rightarrow1\) and \(\gamma\rightarrow0\), we have

$$ \lim_{t\rightarrow\infty}\sup\frac{\ln(V(x(t)))}{\ln t}\leq1 \quad \mbox{a.s.} $$
(4.10)

From \(V(x)\leq\sqrt{2}|x|\) and (4.10), we get

$$\lim_{t\rightarrow\infty}\sup\frac{\ln|x(t)|}{\ln t}\leq1 \quad\mbox{a.s.} $$

 □

Remark 4.1

Similar to [37], \(\forall\varepsilon>0\), we have

$$|x|\leq\kappa+t^{1+\varepsilon} \quad\mbox{for } t\geq0, $$

where \(\kappa=\sup_{0\leq t\leq T}|x(t)|\). This means that \(|x(t)|\) will grow at most polynomially with order close to 1.

Theorem 4.2

Let the same conditions of Theorem  2.1 hold. For any given initial data \(\{(x_{1}(\theta),x_{2}(\theta))^{\top}:-\tau\leq\theta\leq0\}=\xi\in C_{\mathcal{F}_{0}}^{b}([-\tau,0];R_{+}^{0}\times R_{+}^{0})\), the solution x(t) of system (SM) has the property that

$$\lim_{t\rightarrow\infty}\sup\frac{\ln|x(t)|}{\ln t}\leq r_{1}+ \frac{\hat {\sigma}^{2}}{\check{a}^{2}} \quad \textit{a.s.}, $$

where \(\hat{\sigma}=\min\{\sigma_{1}k_{1},\sigma_{2}k_{2}\}\), \(\check{a}=\max\{ a_{12},a_{21}\}\).

Proof

Define

$$V(x)=x_{1}+x_{2}. $$

By Itô’s formula, we have

$$\begin{aligned} \ln\bigl(V\bigl(x(t)\bigr)\bigr) ={}&\ln \bigl(V\bigl(x(0)\bigr)\bigr) \\ &{}+\int_{0}^{t} \biggl(\frac{x_{1}}{V(x(s))} \biggl[r_{1}-a_{11}x_{1}-\frac {a_{12}x_{2}(s-\tau)}{1+mx_{1}+nx_{2}} \biggr] -\frac{\sigma_{1}^{2}x_{1}^{4}}{2V^{2}(x(s))} \biggr)\,ds \\ &{}+\int_{0}^{t} \biggl(\frac{x_{2}}{V(x(s))} \biggl[-r_{2}+\frac{a_{21}x_{1}(s-\tau )}{1+mx_{1}+nx_{2}}-a_{22}x_{2} \biggr] -\frac{\sigma_{2}^{2}x_{2}^{4}}{2V^{2}(x(s))} \biggr)\,ds \\ &{}+M_{1}(t)+M_{2}(t), \end{aligned}$$
(4.11)

where

$$M_{i}(t)=\int_{0}^{t} \frac{\sigma_{i}x_{i}^{2}}{V(x(s))}\,dB_{i}(s), \quad i=1,2, $$

is a real-valued continuous local martingale vanishing at \(t=0\) and its quadratic form is given by

$$\bigl\langle M_{i}(t),M_{i}(t)\bigr\rangle =\int _{0}^{t}\frac{\sigma_{i}^{2}x_{i}^{4}}{V(x(s))}\,ds. $$

Let \(\varepsilon\in(0,0.5)\) and \(\theta=2\) be arbitrary. By the exponential martingale inequality (see, e.g., [36]), for each \(k\geq1\),

$$P \biggl\{ \sup_{0\leq t\leq k} \biggl[M_{i}(t)- \frac{\varepsilon}{2} \bigl\langle M_{i}(t),M_{i}(t)\bigr\rangle \biggr]>\frac{2\ln k}{\varepsilon} \biggr\} \leq k^{-2}. $$

Since the series \(\sum_{k=1}^{\infty}k^{-\theta}\) converges, the well-known Borel-Cantelli lemma yields that there exists \(\Omega_{0}\subset\Omega\) with \(P(\Omega_{0})=1\) such that for any \(\omega\in\Omega_{0}\) there exists an integer \(k_{0}=k_{0}(\omega)\) such that

$$M_{i}(t)\leq\frac{\varepsilon}{2}\bigl\langle M_{i}(t),M_{i}(t) \bigr\rangle +\frac{2 }{\varepsilon}\ln k $$

for all \(0\leq < t\leq k\) and \(k\geq k_{0}(\omega)\). Substituting this into (4.11) we derive that

$$\begin{aligned} \ln\bigl(V\bigl(x(t)\bigr)\bigr) ={}&\ln\bigl(V\bigl(x(0)\bigr)\bigr) \\ &{}+\int_{0}^{t} \biggl(\frac{x_{1}}{V(x(s))} \biggl[r_{1}-a_{11}x_{1}-\frac{a_{12}x_{2}(s-\tau)}{1+mx_{1}+nx_{2}}\biggr] -\frac{(1-\varepsilon)\sigma_{1}^{2}x_{1}^{4}}{2V^{2}(x(s))} \biggr)\,ds \\ &{}+\int_{0}^{t} \biggl(\frac{x_{2}}{V(x(s))} \biggl[-r_{2}+\frac{a_{21}x_{1}(s-\tau )}{1+mx_{1}+nx_{2}}-a_{22}x_{2}\biggr] \\ &{}-\frac{(1-\varepsilon)\sigma_{2}^{2}x_{2}^{4}}{2V^{2}(x(s))} \biggr)\,ds +\frac{4 }{\varepsilon}\ln k \end{aligned}$$

for all \(0\leq t\leq k\), \(k\geq k_{0}(\omega)\) and \(\omega\in\Omega_{0}\). Rearranging the above inequality and using (4.4)-(4.6) give

$$\begin{aligned} \ln\bigl(V\bigl(x(t)\bigr)\bigr) =&\ln\bigl(V\bigl(x(0)\bigr) \bigr)+\frac{4 }{\varepsilon}\ln k \\ &{} +\int_{0}^{t} \biggl(r_{1}+a_{12}x_{2}(s- \tau)+r_{2}+a_{21}x_{1}(s-\tau) \\ &{}-\frac{(1-\varepsilon)\sigma_{1}^{2}k_{1}^{2}x_{1}^{2}}{2} -\frac{(1-\varepsilon)\sigma_{2}^{2}k_{2}^{2}x_{2}^{2}}{2} \biggr)\,ds \\ \leq&\ln\bigl(V\bigl(x(0)\bigr)\bigr)+\frac{4 }{\varepsilon}\ln k +\int _{0}^{t} \biggl(r_{1}+a_{12}x_{2}(s- \tau)+r_{2} \\ &{}+a_{21}x_{1}(s-\tau) -\frac{(1-\varepsilon)\hat{\sigma}^{2}|x|^{2}}{2} \biggr)\,ds \end{aligned}$$
(4.12)

for all \(0\leq t\leq k\), \(k\geq k_{0}(\omega)\) and \(\omega\in\Omega_{0}\). From (4.12), we have

$$\begin{aligned} &\ln\bigl(V\bigl(x(t)\bigr)\bigr)+\frac{(1-2\varepsilon)\hat{\sigma}^{2}}{4} \int_{0}^{t}\bigl|x(s)\bigr|^{2}\,ds \\ &\quad\leq\ln\bigl(V\bigl(x(0)\bigr)\bigr)+\frac{4 }{\varepsilon}\ln k \\ &\qquad{}+\int_{0}^{t} \biggl(r_{1}+a_{12}x_{2}(s- \tau)+r_{2}+a_{21}x_{1}(s-\tau) -\frac{\hat{\sigma}^{2}|x|^{2}}{4} \biggr)\,ds \\ &\quad\leq C_{2}+\frac{4 }{\varepsilon}\ln k +\int_{0}^{t} \biggl(r_{1}+a_{12}x_{2}(s)+r_{2}+a_{21}x_{1}(s) -\frac{\hat{\sigma}^{2}|x|^{2}}{4} \biggr)\,ds, \end{aligned}$$

where

$$C_{2}=\ln \bigl(V\bigl(x(0)\bigr)\bigr)+a_{12}\int _{-\tau}^{0}x_{1}(s)\,ds +a_{21} \int_{-\tau}^{0}x_{2}(s)\,ds. $$

It is easy to see that

$$\begin{aligned} r_{1}+a_{12}x_{2}(t)+r_{2}+a_{21}x_{1}(t) -\frac{\hat{\sigma}^{2}|x|^{2}}{4}&\leq r_{1}+\check{a}|x|- \frac{\hat{\sigma}^{2}|x|^{2}}{4} \leq r_{1}+\frac{\check{a}^{2}}{\hat{\sigma}^{2}}:=K_{2}. \end{aligned}$$

Thus, if \(\omega\in\Omega_{0}\),

$$\ln\bigl(V\bigl(x(t)\bigr)\bigr)+\frac{(1-2\varepsilon)\hat{\sigma}^{2}}{4} \int_{0}^{t}\bigl|x(s)\bigr|^{2}\,ds \leq C_{2}+\frac{4 }{\varepsilon}\ln k+K_{2}t $$

for all \(0\leq t\leq k\), \(k\geq k_{0}(\omega)\). Consequently, for any \(\omega\in\Omega_{0}\), if \(k-1\leq t\leq k\) and \(k\geq k_{0}(\omega)\), we have

$$\frac{1}{t} \biggl[\ln\bigl(V\bigl(x(t)\bigr)\bigr)+\frac{(1-2\varepsilon)\hat{\sigma }^{2}}{4} \int_{0}^{t}\bigl|x(s)\bigr|^{2}\,ds \biggr] \leq \frac{1}{k-1} \biggl[C_{2}+\frac{4 }{\varepsilon}\ln k \biggr]+K_{2}, $$

which implies

$$ \lim_{t\rightarrow\infty}\sup\frac{1}{t} \biggl[\ln\bigl(V\bigl(x(t) \bigr)\bigr) +\frac{(1-2\varepsilon)\hat{\sigma}^{2}}{4}\int_{0}^{t}\bigl|x(s)\bigr|^{2}\,ds \biggr] \leq K_{2}. $$
(4.13)

Using (4.13) and letting \(\varepsilon\rightarrow0\), we get

$$ \lim_{t\rightarrow\infty}\sup\frac{1}{t} \biggl[\ln\bigl(V\bigl(x(t) \bigr)\bigr) +\frac{\hat{\sigma}^{2}}{4}\int_{0}^{t}\bigl|x(s)\bigr|^{2}\,ds \biggr] \leq K_{2} \quad\mbox{a.s.} $$
(4.14)

Using \(V(x)\geq\frac{|x|}{\sqrt{2}}\) and (4.14), we get

$$\lim_{t\rightarrow\infty}\sup\frac{1}{t} \biggl[\ln\bigl|x(t)\bigr| + \frac{\hat{\sigma}^{2}}{4}\int_{0}^{t}\bigl|x(s)\bigr|^{2}\,ds \biggr] \leq K_{2} \quad\mbox{a.s.} $$

The proof is completed. □

5 Numerical simulations

Now let us use Milstein’s numerical method (see, e.g., [38]) to support our results. In Figure 1, we choose \(r_{1}=0.02\), \(a_{11}=0.1\), \(a_{12}=0.03\), \(r_{2}=3\), \(a_{22}=0.3\), \(\sigma_{1}=\sigma _{2}=0.02\), \(m=n=1.25\), \(\tau=1\). The difference between the conditions of Figures 1(a) and (b) is that the values of \(a_{21}\) are different. In Figure 1(a) we choose \(a_{21}=0.3<4\), then the conditions of Theorem 2.1 hold. Making Theorems 4.1 and 4.2 lead to system (SM) has almost surely asymptotic properties. In Figure 1(b) we choose \(a_{21}=500>4\), then the conditions of Theorem 2.1 are not satisfied; furthermore, the conditions of Theorems 4.1 and 4.2 do not hold. Hence, the population \(x_{1}\) goes to extinction, and \(x_{2}\) has no almost surely asymptotic properties, see Figure 1(b).

Figure 1
figure 1

Solutions of system ( SM ) for \(\pmb{r_{1}=0.02}\) , \(\pmb{a_{11}=0.1}\) , \(\pmb{a_{12} =0.03}\) , \(\pmb{r_{2}=3}\) , \(\pmb{a_{22}=0.3}\) , \(\pmb{\sigma_{1}=\sigma_{2}=0.02}\) , \(\pmb{m=n=1.25}\) , \(\pmb{\tau=1}\) . The horizontal axis represents the time t, the vertical axis represents the population sizes. (a) is with \(a_{21}=0.3\); (b) is with \(a_{21}=500\).

Full size image

6 Conclusions and future directions

A stochastic delay predator-prey system is considered and system (SM) is more general than the classical predator-prey system with the Beddington-DeAngelis functional response. We have established a sufficient condition under which system (SM) has a global positive solution. We have also discussed the asymptotic properties for the moments as well as the sample paths of the solution. In particular, we have studied a fundamental problem in population systems, namely stochastically ultimate boundedness. Our key contributions in this paper are the following:

  • This paper deals with a kind of delay stochastic population system, while most of existing results (see, e.g., [4, 5, 9–11, 22, 31]) are concerned with the non-delayed cases.

  • Our stochastic population system is therefore more complicated and the mathematics presented is more difficult.

There are still many interesting and challenging questions that need to be studied. In this paper, we only consider the growth rate \(a_{11}\), \(a_{22}\) to be stochastic; other parameters, for example, \(r_{i}\), \(i=1,2\), is stochastic, which is not studied. We wish that such questions will be investigated by some authors.