1 Introduction

Predator–prey system is very popular in the world. In order to reveal the dynamical relationship between predator and prey, a lot of predator–prey systems have been widely investigated and many good results have been obtained in the last decades, which has long been one of the hot topics in ecology [13]. Since two-species ecological models cannot describe the natural phenomena accurately and many vital behaviors can only be exhibited by systems with three or more species, for example, in the natural world, the predator often feeds on some competing prey, and hence, a three- or multi-species population system attracts more and more attention [47].

On the other hand, all species are inevitably affected by environmental noise. To better describe ecological phenomena, the white noise is introduced into a predator–prey model to reveal richer and more complex dynamics [815]. There are many kinds of stochastic perturbation. Considering the stochastic influence on the intrinsic growth rates of populations, we have \(a_{i}\rightarrow a_{i} +\xi _{i} \,d{\omega } (t)\), where \(\omega (t)\) is a standard Brownian motion defined on a complete probability space \((\Omega , \mathcal{{F} },\mathcal{{P})}\) with a filtration \(\{{\mathcal{{F}}}_{t}\}_{t\geq 0} \), \(\xi ^{2}\) is the intensity of white noise. For example, Liu [6] proposed the following three-species predator–prey model:

$$ \textstyle\begin{cases} dN_{1}(t)=N_{1}(t) (a_{1}-d_{11}N_{1}(t)-d_{12}N_{2}(t)-d_{13}N_{3}(t))\,dt+ \xi _{1}N_{1}(t)\,d\omega (t), \\ dN_{2}(t)=N_{2}(t) (a_{2}-d_{21}N_{1}(t)-d_{22}N_{2}(t)-d_{23}N_{3}(t))\,dt+ \xi _{2}N_{2}(t)\,d\omega (t), \\ dN_{3}(t)=N_{3}(t) (a_{3}-d_{31}N_{1}(t)-d_{32}N_{2}(t)-d_{33}N_{3}(t))\,dt+ \xi _{3}N_{3}(t)\,d\omega (t), \end{cases} $$

where \(N_{1}(t)\) and \(N_{2}(t)\) are the population sizes of prey species, \(N_{3}(t)\) is the population size of predator species, \(a_{i}>0\) (\(i=1,2\)) are the intrinsic rates of increase, \(a_{3}<0\) is the intrinsic rate of decrease, \(d_{12}>0\) and \(d_{21}>0\) are the parameters representing competitive effects between two prey, \(d_{13}>0\) and \(d_{23}>0\) are the coefficients of decrease of prey species due to predation, \(d_{31}<0\) and \(d_{32}<0\) are the predation rate of predator, \(d_{ii}>0\) (\(i=1,2,3\)) are the rate of competition within the same species.

As we know, predator–prey interaction is a frequently observed phenomenon. Almost all species should exhibit some delays. Considering the inevitability, more and more researchers have taken delay into an ecological model and obtained some nice results [1619]. Recently, infinite delay has been widely introduced into the ecological model since the works of Volterra to translate the cumulative effect of the past history of a system [2024]. Chen [22] et al. proposed the following model with distributed delays:

$$ \textstyle\begin{cases} dN_{1}(t)/dt= b_{1} N_{1}(t) (1- \frac{N_{1}(t)}{K} )-a_{12}N_{1}(t)N_{2}(t), \\ dN_{2}(t)/dt= -b_{2}N_{2}(t)+a_{21} \int _{-\infty }^{t}K(t-s)N_{1}(s)N_{2}(s)\,ds, \end{cases} $$

where the kernel \(K:[0,\infty )\rightarrow [0,\infty )\) is a normalized \(L^{1}\) function such that \(\int _{0}^{\infty }K(s)\,ds=1\). For distributed delay, MacDonald [25] initially proposed that it was reasonable to use gamma distribution as a kernel function, that is, \(f(t)=\frac{t^{n}\sigma ^{n+1}e^{-\sigma t}}{n!}\), where \(\sigma >0, n\) is a nonnegative integer. If \(n=0\), then the kernel \(f(t)=\sigma e^{-\sigma t}\) is called a weak kernel, otherwise it is called a strong kernel.

Motivated by the above discussion, in this paper, we consider a stochastic two-prey one-predator system with distributed delays. For convenience, we mainly consider the weak kernel case, i.e., \(f(t)=\sigma e^{-\sigma t}\). Our model is as follows:

$$ \textstyle\begin{cases} dN_{1}(t)= N_{1}(t) (r_{1}-a_{11}N_{1}(t)-a_{12} \int _{-\infty }^{t} \sigma _{2}e^{-\sigma _{2}(t-s)}N_{2}(s)\,ds \\ \hphantom{dN_{1}(t)=}{} -a_{13} \int _{-\infty }^{t}\sigma _{3}e^{- \sigma _{3}(t-s)} N_{3}(s)\,ds )\,dt+\xi _{1}N_{1}(t)\,d\omega _{1}(t), \\ dN_{2}(t)= N_{2}(t^{-}) (r_{2}-a_{21} \int _{- \infty }^{t}\sigma _{1}e^{-\sigma _{1}(t-s)} N_{1}(s)\,ds-a_{22}N_{2}(t) \\ \hphantom{dN_{2}(t)=}{} -a_{23} \int _{-\infty }^{t}\sigma _{3}e^{- \sigma _{3}(t-s)} N_{3}(s)\,ds )\,dt+\xi _{2}N_{2}(t)\,d\omega _{2}(t), \\ dN_{3}(t) = N_{3}(t) (-r_{3}+a_{31} \int _{- \infty }^{t} \sigma _{1}e^{-\sigma _{1}(t-s)}N_{1}(s)\,ds -a_{33}N_{3}(t) \\ \hphantom{dN_{3}(t) =}{} +a_{32} \int _{-\infty }^{t}\sigma _{2}e^{- \sigma _{2}(t-s)} N_{2}(s)\,ds )\,dt+\xi _{3}N_{3}(t)\,d\omega _{3}(t), \end{cases} $$
(1.1)

with the initial data

$$ N_{i}(\theta )=\varphi _{i}(\theta )\in C ((-\infty ,0],R_{+} ),\quad i=1,2,3, $$

where \(C((-\infty ,0],R_{+})\) is the set of all continuous functions from \((-\infty ,0)\) to \(R_{+}=(0,\infty )\), \(\omega _{i}(t)\) (\(i=1,2,3\)) is a standard and independent Brownian motion defined as above. All parameters are positive constants and their biological meanings refer to [6].

Define

$$ y_{i}(t)= \int _{-\infty }^{t}\sigma _{i}e^{-\sigma _{i}(t-s)}N_{i}(s) \,ds,\quad i=1,2,3. $$

Computing the derivative of \(y_{i}(t)\), then \(dy_{i}(t)=\sigma _{i}(N_{i}(t)-y_{i}(t))\,dt\), \(i=1,2,3\). Using the linear chain technique to (1.1) yields

$$ \textstyle\begin{cases} dN_{1}(t)=N_{1}(t) (r_{1}-a_{11}N_{1}(t)-a_{12}y_{2}(t) -a_{13}y_{3}(t) )\,dt+\xi _{1}N_{1}(t)\,d\omega _{1}(t), \\ dN_{2}(t)= N_{2}(t) (r_{2}-a_{21}y_{1}(t)-a_{22}N_{2}(t) -a_{23}y_{3}(t) )\,dt+\xi _{2}N_{2}(t)\,d\omega _{2}(t), \\ dN_{3}(t)= N_{3}(t) (-r_{3}+a_{31}y_{1}(t)+a_{32}y_{2}(t)-a_{33}N_{3}(t) )\,dt +\xi _{3}N_{3}(t)\,d\omega _{3}(t), \\ dy_{1}(t) = \sigma _{1}(N_{1}(t)-y_{1}(t))\,dt, \\ dy_{2}(t)= \sigma _{2}(N_{2}(t)-y_{2}(t))\,dt, \\ dy_{3}(t) = \sigma _{3}(N_{3}(t)-y_{3}(t))\,dt. \end{cases} $$
(1.2)

According to the equivalent property of (1.1) and (1.2), in what follows, we mainly consider (1.2) to reveal the dynamical properties of (1.1). Our main aims are as follows.

Firstly, we study the stability in mean and extinction of all species of (1.2), which have long been and will still be two important topics for the study of stochastic population systems.

Secondly, for a stochastic population system, instead of the positive equilibrium state of the determinate system, it is interesting and important to study the existence and uniqueness of the distribution of (1.2).

The rest work of this paper is organized as follows. Section 2 begins with some notations, definitions, and important lemmas. Section 3 focuses on the stability in mean and extinction of species of (1.2). Section 4 is devoted to the existence and uniqueness of distribution. Some numerical simulations are given in Sect. 5. Finally, we conclude the paper with a brief conclusion and discussion in Sect. 6.

2 Preliminaries

For simplicity, we give the following notations.

$$\begin{aligned}& \alpha _{1}=(a_{11}, a_{21}, -a_{31})^{T},\qquad \alpha _{2}=(a_{12}, a_{22}, -a_{32})^{T},\qquad \alpha _{3}=(a_{13}, a_{23}, a_{33})^{T}, \\& r=(r_{1},r_{2},-r_{3})^{T},\qquad \xi = \bigl(\xi _{1}^{2}/2,\xi _{2}^{2}/2, \xi _{3}^{2}/2 \bigr),\qquad A=\texttt{det}(\alpha _{1}, \alpha _{2}, \alpha _{3}), \\& b_{1}=r_{1}-\xi _{1}^{2}/2,\qquad b_{2}=r_{2}-\xi _{2}^{2}/2,\qquad b_{3}=-r_{3}- \xi _{3}^{2}/2, \\& A_{1}=\texttt{det}(r,\alpha _{2},\alpha _{3}), \qquad A_{2}= \texttt{det}(\alpha _{1},r,\alpha _{3}),\qquad A_{3}=\texttt{det}( \alpha _{1}, \alpha _{2},r), \\& \tilde{A}_{1}=\texttt{det}(\xi ,\alpha _{2},\alpha _{3}),\qquad \tilde{A}_{2}=\texttt{det}(\alpha _{1},\xi ,\alpha _{3}),\qquad \tilde{A}_{3}= \texttt{det}(\alpha _{1},\alpha _{2},\xi ), \\& \Delta _{1}=r_{2} a_{32}+r_{3} a_{22}, \qquad \Delta _{2}=r_{1}a_{31}+r_{3}a_{11}, \qquad \Delta _{3}=-r_{1}a_{21}+r_{2}a_{11}, \\& \tilde{\Delta }_{1}=\frac{\xi _{2}^{2}}{2} a_{32}+ \frac{\xi _{3}^{2}}{2} a_{22}, \qquad \tilde{\Delta }_{2}= \frac{\xi _{1}^{2}}{2}a_{31}+\frac{\xi _{3}^{2}}{2}a_{11},\qquad \tilde{ \Delta }_{3}=-\frac{\xi _{1}^{2}}{2}a_{21}+ \frac{\xi _{2}^{2}}{2}a_{11}, \\& \Delta _{1}^{*}=r_{2} a_{33}-r_{3} a_{23},\qquad \Delta _{2}^{*}=r_{1}a_{33}-r_{3}a_{13}, \qquad \Delta _{3}^{*}=r_{1}a_{22}-r_{2}a_{12}, \\& \tilde{\Delta }_{1}^{*}=\frac{\xi _{2}^{2}}{2} a_{33}- \frac{\xi _{3}^{2}}{2} a_{23}, \qquad \tilde{\Delta }_{2}^{*}= \frac{\xi _{1}^{2}}{2}a_{33}- \frac{\xi _{3}^{2}}{2}a_{13}, \qquad \tilde{\Delta }_{3}^{*}= \frac{\xi _{1}^{2}}{2}a_{22}- \frac{\xi _{2}^{2}}{2}a_{12}. \end{aligned}$$

Throughout this paper, we denote the complement minor of \(a_{ij}\) in determinant A by \(A_{ij}\) (\(i,j=1,2,3\)), and assume that \(A>0\), \(A_{i}>0\), i.e., when there is no stochastic perturbation, a positive equilibrium state exists for model (1.1). Further, for convenience, we always assume that K stands for a generic positive constant whose value may be different at different places. And for any function \(x(t)\), \(t>0\), we denote

$$ \bigl\langle x(t) \bigr\rangle =t^{-1} \int _{0}^{t}x(s)\,ds, \qquad x^{*}= \limsup_{t\rightarrow \infty } x(t),\qquad x_{*}=\liminf _{t\rightarrow \infty } x(t). $$

Now we give assumptions, definitions, and some important lemmas, which are used in our main proof.

Assumption 2.1

\(A_{13}>0\), \(A_{23}<0\), \(A_{33}>0\), \(A_{31}<0\), \(A_{32}>0 \).

Assumption 2.2

\(a_{ii}>\sum_{j=1,j\neq i}^{3}a_{ji}\), \(i,j=1,2,3\).

Remark 2.1

Assumption 2.2 means that the intra-specific competitive rates are stronger than the interaction competitive rates or predation rates among different species.

Definition 2.1

Let \(P(t)=(N_{1}(t),N_{2}(t),N_{3}(t),y_{1}(t),y_{2}(t),y_{3}(t))^{T} \in C((-\infty ,0],R_{+}^{6})\) be a solution of system (1.2), then

  1. (I)

    The population \(P(t)\) is said to be extinct if \(\lim_{t\rightarrow \infty }P(t)=0\);

  2. (II)

    The population \(P(t)\) is said to be stable in mean if \(\lim_{t\rightarrow \infty }\langle P(t)\rangle =K\), a.s., where K is a constant.

Definition 2.2

Let \(P(t)=(N_{1}(t),N_{2}(t),N_{3}(t),y_{1}(t),y_{2}(t),y_{3}(t))^{T}\in C((-\infty ,0], R_{+}^{6})\) and \(\bar{P}(t)=(\bar{N}_{1}(t),\bar{N}_{2}(t),\bar{N}_{3}(t),\bar{y}_{1}(t), \bar{y}_{2}(t),\bar{y}_{3}(t))^{T} \in C((-\infty ,0],R_{+}^{6})\) be any two positive solutions of (1.2) with the initial value \(P(0)>0\), \(\bar{P}(0)>0\), then system (1.2) is said to be globally attractive if

$$ \lim_{t\rightarrow \infty } \bigl\vert N_{i}(t)- \bar{N}_{i}(t) \bigr\vert =0, \qquad \lim_{t\rightarrow \infty } \bigl\vert y_{i}(t)-\bar{y}_{i}(t) \bigr\vert =0,\quad i=1,2,3. $$

Lemma 2.1

([26])

Suppose that \(Z(t)\in C[\Omega \times [0,+\infty ), R_{+}] \) and \(\lim_{t\rightarrow \infty }F(t)/t=0\), a.s.

  1. (a)

    If there exist two positive constants \(T>0\), \(\lambda _{0}>0\) such that, for all \(t>T\),

    $$\begin{aligned}& \ln Z(t)\leq \lambda t-\lambda _{0} \int _{0}^{t}z(s)\,ds+F(t), \quad \textit{a.s.}, \\& \textit{then}\quad \textstyle\begin{cases} \langle Z\rangle ^{*}\leq \lambda /\lambda _{0}, &\textit{a.s.}, \textit{if } \lambda \geq 0, \\ \lim_{t\rightarrow +\infty } Z(t)= 0, & \textit{a.s.}, \textit{if } \lambda < 0. \end{cases}\displaystyle \end{aligned}$$
  2. (b)

    If there exist some constants \(T>0\), \(\lambda _{0}>0,\lambda \) such that, for all \(t>T\),

    $$\begin{aligned}& \ln Z(t)\geq \lambda t-\lambda _{0} \int _{0}^{t}z(s)\,ds+F(t), \quad \textit{a.s.}, \\& \textit{then} \quad \langle Z\rangle _{*}\geq \lambda /\lambda _{0}, \quad \textit{a.s.} \end{aligned}$$

Lemma 2.2

System (1.2) has a unique solution \(P(t)=(N_{1}(t),N_{2}(t),N_{3}(t),y_{1}(t), y_{2}(t), y_{3}(t))^{T}\in C((-\infty ,0],R_{+}^{6})\) for any given initial data \(P(t_{0})\in C((-\infty ,0],R_{+}^{6})\), almost surely.

Proof

The proof is standard. For the readers’ convenience, we give the proof in Appendix A.1.

As to the expectation boundedness and asymptotical properties of the solution of (1.2), we have the following lemma. The proof is similar to that of references [20, 27, 28] and is presented in Appendix A.2. □

Lemma 2.3

Let \(P(t)=(N_{1}(t),N_{2}(t),N_{3}(t),y_{1}(t),y_{2}(t),y_{3}(t))^{T}\) be the solution of (1.2), then for any initial data \(P(t_{0})\in C((-\infty ,0],R_{+}^{6})\), there exists a positive constant \(K(p)\) such that

$$ \limsup_{t\rightarrow +\infty }\mathbb{E} \bigl(N_{i}(t)^{p} \bigr)\leq K(p), \qquad \limsup_{t\rightarrow +\infty }\mathbb{E} \bigl(y_{i}(t)^{p} \bigr) \leq K(p), $$

further,

$$ \lim_{t\rightarrow +\infty }\frac{y_{i}(t)}{ t}= 0,\qquad \limsup _{t\rightarrow +\infty } \frac{\ln N_{i}(t)}{ t} \leq 0, \quad \textit{a.s.}, i=1,2,3. $$

For the following integral equation

$$ x(t)=x(t_{0})+ \int _{t_{0}}^{t}a \bigl(s,x(s) \bigr)\,ds+\sum _{r=1}^{k} \int _{t_{0}}^{t}b_{r} \bigl(s,x(s) \bigr)\,d \sigma _{r}(s), $$
(2.1)

there is a result as follows.

Lemma 2.4

([29])

Suppose that the coefficients of (2.1) are independent of t and satisfy:

$$\begin{aligned}& \bigl\vert a(s,x)-a(s,y) \bigr\vert +\sum_{r=1}^{k} \bigl\vert b_{r}(s,x)-b_{r}(s,y) \bigr\vert \leq K \vert x-y \vert , \\& \bigl\vert a(s,x) \bigr\vert +\sum_{r=1}^{k} \bigl\vert b_{r}(s,x) \bigr\vert \leq K \bigl(1+ \vert x \vert \bigr), \end{aligned}$$

in \(U_{R}\) for any \(R>0\), and there exists a nonnegative \(C^{2}\) function \(V(x)\) in \(R^{l}\) such that

$$ LV(x)\leq -1 $$
(2.2)

outside some compact set, then system (2.1) has a solution, which is a stationary Markov process.

Lemma 2.5

([30])

Let \(f(t)\) be a nonnegative function defined on \([0,+\infty )\) such that \(f(t)\) is integrable on \([0,+\infty )\) and is uniformly continuous on \([0,+\infty )\), then \(\lim_{t\rightarrow \infty }f(t)=0\).

Lemma 2.6

Let \(P(t)=(N_{1}(t),N_{2}(t),N_{3}(t),y_{1}(t),y_{2}(t),y_{3}(t))^{T}\) be a solution of (1.2) with the initial value \(P(0)>0\), then almost every sample path of \(P(t)\) is uniformly continuous on \(t\geq 0\).

Proof

For the first equation of (1.2), it is equivalent to the following stochastic integral equation:

$$ N_{1}(s)=N_{1}(0)+ \int _{0}^{t}N_{1}(s) \bigl(r_{1}-a_{11}N_{1}(s)+a_{12}y_{2}(s) +a_{13}y_{3}(s) \bigr)\,ds+\xi _{1} \int _{0}^{t}N_{1}(s)\,d\omega _{1}(s). $$

By computation, we have

$$\begin{aligned}& \mathbb{E} \bigl\vert N_{1}(s) \bigl(r_{1}-a_{11}N_{1}(s)+a_{12}y_{2}(s)+a_{13}y_{3}(s) \bigr) \bigr\vert ^{p} \\& \quad \leq \mathbb{E} \vert N_{1} \vert ^{2p}/2+ \mathbb{E} \bigl\vert r_{1}-a_{11}N_{1}(s) +a_{12}y_{2}(s)+a_{13}y_{3}(s) \bigr\vert ^{2p}/2 \\& \quad \leq \mathbb{E} \vert N_{1} \vert ^{2p}/2+3^{2p-1} \bigl(r_{1}^{2p} +a_{11} \mathbb{E}N_{1}^{2p}+a_{12} \mathbb{E}y_{2}^{2p}+a_{13}\mathbb{E}y_{3}^{2p} \bigr)/2 \\& \quad \leq K. \end{aligned}$$

Using the moment inequality for stochastic integrals, for any \(0\leq t_{1}\leq t_{2}\), \(p>2\), we have

$$\begin{aligned} \mathbb{E} \biggl\vert \int _{t_{1}}^{t_{2}}\xi _{1}N_{1}(s)\,d \omega _{1}(s) \biggr\vert ^{p} &\leq \sigma _{1}^{p} \biggl( \frac{p(p-1)}{2} \biggr)^{ \frac{p}{2}}(t_{2}-t_{1})^{ \frac{p-2}{2}} \int _{t_{1}}^{t_{2}}E \bigl\vert N_{1}(s) \bigr\vert ^{p}\,ds \\ &\leq \sigma _{1}^{p} \biggl( \frac{p(p-1)}{2} \biggr)^{ \frac{p}{2}}(t_{2}-t_{1})^{ \frac{p-2}{2}}K. \end{aligned}$$

In the same manner, we can discuss the following five equations of (1.2) and obtain similar inequalities as above. Therefore, by Lemma 2.4 of Refs. [31], we conclude that almost every sample path of \(P(t)\) is uniformly continuous. The proof is completed. □

3 Stability in mean and extinction of species

Firstly, we give the following result on stability in mean and extinction of species of model (1.2).

Theorem 3.1

If Assumptions 2.1and 2.2hold, then for system (1.2), we have:

  1. (i)

    If \(b_{1}<0\), \(b_{2}<0\), then \(\lim_{t\rightarrow \infty }N_{i}(t)=0\), \(i=1,2,3\);

  2. (ii)

    If \(b_{1}<0\), \(b_{2}>0\), \(\Delta _{1}<\tilde{\Delta }_{1}\), then

    $$ \lim_{t\rightarrow \infty }N_{1}(t)=0, \qquad \lim_{t \rightarrow \infty } \bigl\langle N_{2}(t) \bigr\rangle =\frac{b_{2}}{a_{22}}, \qquad \lim _{t\rightarrow \infty } N_{3}(t)=0; $$

    If \(b_{1}<0\), \(b_{2}>0\), \(\Delta _{1}>\tilde{\Delta }_{1}\), then

    $$ \lim_{t\rightarrow \infty }N_{1}(t)=0,\qquad \lim_{t \rightarrow \infty } \bigl\langle N_{2}(t) \bigr\rangle = \frac{\Delta _{1}^{*}-\tilde{\Delta }_{1}^{*}}{A_{11}},\qquad \lim _{t\rightarrow \infty } \bigl\langle N_{3}(t) \bigr\rangle = \frac{\Delta _{1}-\tilde{\Delta }_{1}}{A_{11}}; $$
  3. (iii)

    If \(b_{1}>0\), \(b_{2}<0\), \(\Delta _{2}<\tilde{\Delta }_{2}\), then

    $$ \lim_{t\rightarrow \infty } \bigl\langle N_{1}(t) \bigr\rangle = \frac{b_{1}}{a_{11}},\qquad \lim_{t\rightarrow \infty } N_{2}(t)=0,\qquad \lim _{t\rightarrow \infty } \bigl\langle N_{3}(t) \bigr\rangle =0; $$

    If \(b_{1}>0\), \(b_{2}<0\), \(\Delta _{2}>\tilde{\Delta }_{2}\), then

    $$ \lim_{t\rightarrow \infty }N_{2}(t)=0, \qquad \lim_{t \rightarrow \infty } \bigl\langle N_{1}(t) \bigr\rangle = \frac{\Delta _{2}^{*}-\tilde{\Delta }_{2}^{*}}{A_{22}},\qquad \lim _{t\rightarrow \infty } \bigl\langle N_{3}(t) \bigr\rangle = \frac{\Delta _{2}-\tilde{\Delta }_{2}}{A_{22}}; $$
  4. (iv)

    If \(b_{1}>0\), \(b_{2}>0\), \(A_{i}>\tilde{A}_{i}\) (\(i=1,2,3\)), then

    $$ \lim_{t\rightarrow \infty } \bigl\langle N_{1}(t) \bigr\rangle = \frac{A_{1}-\tilde{A}_{1}}{A},\qquad \lim_{t\rightarrow \infty } \bigl\langle N_{2}(t) \bigr\rangle =\frac{A_{2}-\tilde{A}_{2}}{A},\qquad \lim_{t\rightarrow \infty } \bigl\langle N_{3}(t) \bigr\rangle = \frac{A_{3}-\tilde{A}_{3}}{A}; $$

    If \(b_{1}>0\), \(b_{2}>0\), \(A_{3}<\tilde{A}_{3}\), \(\Delta _{3}>\tilde{\Delta }_{3}\), then

    $$ \lim_{t\rightarrow \infty } \bigl\langle N_{3}(t) \bigr\rangle =0,\qquad \lim_{t\rightarrow \infty } \bigl\langle N_{1}(t) \bigr\rangle = \frac{\Delta _{3}^{*}-\tilde{\Delta }_{3}^{*}}{A_{33}},\qquad \lim_{t\rightarrow \infty } \bigl\langle N_{2}(t) \bigr\rangle = \frac{\Delta _{3}-\tilde{\Delta }_{3}}{A_{33}}; $$

Proof

For (1.2), integrating the forth to the sixth equations from 0 to t leads to

$$ \frac{y_{i}(t)-y_{i}(0)}{t}=\sigma _{i} \bigl( \bigl\langle N_{i}(t) \bigr\rangle - \bigl\langle y_{i}(t) \bigr\rangle \bigr),\quad i=1,2,3. $$

Taking the limit as \(t\rightarrow \infty \), combining with Lemma 2.3, we have

$$ \lim_{t\rightarrow \infty } \bigl\langle N_{i}(t) \bigr\rangle =\lim _{t\rightarrow \infty } \bigl\langle y_{i}(t) \bigr\rangle . $$

By utilizing Itô’s formula to \(\ln N_{i}(t)\) (\(i=1,2,3\)) and integrating both sides of the first three equations of (1.2) from 0 to t, we obtain

$$ \textstyle\begin{cases} \ln N_{1}(t)-\ln N_{1}(0) = b_{1}t-a_{11} \int _{0}^{t}N_{1}(s)\,ds -a_{12} \int _{0}^{t}N_{2}(s)\,ds \\ \hphantom{\ln N_{1}(t)-\ln N_{1}(0) =}{}-a_{13} \int _{0}^{t}N_{3}(s)\,ds +\xi _{1}\omega _{1}(t), \\ \ln N_{2}(t)-\ln N_{2}(0) = b_{2}t-a_{21} \int _{0}^{t}N_{1}(s)\,ds-a_{22} \int _{0}^{t}N_{2}(s)\,ds \\ \hphantom{\ln N_{2}(t)-\ln N_{2}(0) =}{} -a_{23} \int _{0}^{t}N_{3}(s)\,ds +\xi _{2}\omega _{2}(t), \\ \ln N_{3}(t)-\ln N_{3}(0) = b_{3}t+a_{31} \int _{0}^{t}N_{1}(s)\,ds+a_{32} \int _{0}^{t}N_{2}(s)\,ds \\ \hphantom{\ln N_{3}(t)-\ln N_{3}(0) =}{}-a_{33} \int _{0}^{t}N_{3}(s)\,ds +\xi _{3}\omega _{3}(t). \end{cases} $$
(3.1)

Denote \(\xi _{i}(t)\omega _{i}(t)=\vartheta _{i}(t)\), then

$$\begin{aligned}& \frac{\ln N_{1}(t)-\ln N_{1}(0)}{t}=b_{1}-a_{11} \bigl\langle N_{1}(t) \bigr\rangle -a_{12} \bigl\langle N_{2}(t) \bigr\rangle -a_{13} \bigl\langle N_{3}(t) \bigr\rangle + \frac{\vartheta _{1}(t)}{t}, \end{aligned}$$
(3.2)
$$\begin{aligned}& \frac{\ln N_{2}(t)-\ln N_{2}(0)}{t}=b_{2}-a_{21} \bigl\langle N_{1}(t) \bigr\rangle -a_{22} \bigl\langle N_{2}(t) \bigr\rangle -a_{23} \bigl\langle N_{3}(t) \bigr\rangle + \frac{\vartheta _{2}(t)}{t}, \end{aligned}$$
(3.3)

and

$$ \frac{\ln N_{3}(t)-\ln N_{3}(0)}{t}=b_{3}+a_{31} \bigl\langle N_{1}(t) \bigr\rangle +a_{32} \bigl\langle N_{2}(t) \bigr\rangle -a_{33} \bigl\langle N_{3}(t) \bigr\rangle \frac{\vartheta _{3}(t)}{t}. $$
(3.4)

We begin with the proof of (i).

It follows from (3.2) and (3.3) that

$$ t^{-1} \bigl[\ln N_{1}(t)-\ln N_{1}(0) \bigr]\leq b_{1}-a_{11} \bigl\langle N_{1}(t) \bigr\rangle +t^{-1}\vartheta _{1}(t) $$
(3.5)

and

$$ t^{-1} \bigl[\ln N_{2}(t)-\ln N_{2}(0) \bigr]\leq b_{2}-a_{22} \bigl\langle N_{2}(t) \bigr\rangle +t^{-1}\vartheta _{2}(t). $$
(3.6)

Using Lemma 2.1 to (3.5) and (3.6), then

$$ \lim_{t\rightarrow \infty }N_{i}(t)=0,\quad i=1,2. $$

Since \(b_{3}<0\), (3.4) implies \(\lim_{t\rightarrow \infty }N_{3}(t)=0\), and hence

$$ \lim_{t\rightarrow \infty }N_{i}(t)=0 \quad \text{for } i=1,2,3. $$

Now we prove (ii).

By the proof of (i), if \(b_{1}<0\), then \(\lim_{t\rightarrow \infty }N_{1}(t)=0\), and hence (3.3) and (3.4) imply that

$$ t^{-1} \bigl[\ln N_{2}(t)-\ln N_{2}(0) \bigr]=b_{2}-a_{22} \bigl\langle N_{2}(t) \bigr\rangle -a_{23} \bigl\langle N_{3}(t) \bigr\rangle +t^{-1}\vartheta _{2}(t) $$
(3.7)

and

$$ t^{-1} \bigl[\ln N_{3}(t)-\ln N_{3}(0) \bigr]=b_{3}+a_{32} \bigl\langle N_{2}(t) \bigr\rangle -a_{33} \bigl\langle N_{3}(t) \bigr\rangle +t^{-1}\vartheta _{3}(t). $$
(3.8)

By the elimination method, adding (3.7) applied by \(a_{33}\) and (3.8) applied by \(-a_{23}\) gives

$$ \begin{aligned} & t^{-1}\bigl(a_{33} \bigl[\ln N_{2}(t)-\ln N_{2}(0)\bigr]-a_{23}\bigl[ \ln N_{3}(t)- \ln N_{3}(0)\bigr]\bigr) \\ &\quad = (b_{2}a_{33}-b_{3}a_{23})-(a_{22}a_{33}+a_{32}a_{23}) \bigl\langle N_{2}(t) \bigr\rangle +t^{-1} \bigl(a_{33}\vartheta _{2}(t)-a_{23}\vartheta _{3}(t)\bigr). \end{aligned} $$
(3.9)

Applying Lemma 2.1 and Lemma 2.3 to (3.9), we get

$$ \bigl\langle N_{2}(t) \bigr\rangle ^{*}\leq \frac{b_{2}a_{33}-b_{3}a_{23}}{a_{22}a_{33}+a_{32}a_{23}} = \frac{\Delta _{1}^{*}-\tilde{\Delta }_{1}^{*}}{A_{11}}. $$

Substituting \(\langle N_{2}(t)\rangle ^{*}\) into (3.8), we obtain

$$ t^{-1} \bigl[\ln N_{3}(t)-\ln N_{3}(0) \bigr]\leq b_{3} +a_{32} \frac{b_{2}a_{33}-b_{3}a_{23}}{a_{22}a_{33}+a_{32}a_{23}} -a_{33} \bigl\langle N_{3}(t) \bigr\rangle +t^{-1}\vartheta _{3}(t) $$
(3.10)

and

$$ \bigl\langle N_{3}(t) \bigr\rangle ^{*}\leq \frac{b_{2}a_{32}+b_{3}a_{22}}{a_{22}a_{33}+a_{32}a_{23}}= \frac{\Delta _{1}-\tilde{\Delta }_{1}}{A_{11}}. $$

If \(b_{1}<0\), \(b_{2}>0\), \(\Delta _{1}<\tilde{\Delta }_{1}\), then (3.10) implies \(\lim_{t\rightarrow \infty } N_{3}(t) =0\). By use of (3.7) again, we have

$$ t^{-1} \bigl[\ln N_{2}(t)-\ln N_{2}(0) \bigr]\leq b_{2}+\varepsilon -a_{22} \bigl\langle N_{2}(t) \bigr\rangle +t^{-1}\vartheta _{2}(t). $$
(3.11)

Applying Lemma 2.1 to (3.11) gives

$$ \bigl\langle N_{2}(t) \bigr\rangle ^{*}\leq \frac{b_{2}}{a_{22}}. $$

Similarly, we have

$$ t^{-1} \bigl[\ln N_{2}(t)-\ln N_{2}(0) \bigr]\geq b_{2}-\varepsilon -a_{22} \bigl\langle N_{2}(t) \bigr\rangle +t^{-1}\vartheta _{2}(t) $$

and \(\langle N_{2}(t)\rangle _{*}\geq \frac{b_{2}}{a_{22}}\), and hence

$$ \lim_{t\rightarrow \infty } \bigl\langle N_{2}(t) \bigr\rangle = \frac{b_{2}}{a_{22}}. $$

If \(b_{1}<0\), \(b_{2}>0\), \(\Delta _{1}>\tilde{\Delta }_{1}\), then we can derive from (3.7) and (3.8) that

$$ t^{-1} \bigl[\ln N_{2}(t)-\ln N_{2}(0) \bigr]\geq b_{2}-a_{22} \bigl\langle N_{2}(t) \bigr\rangle -a_{23} \bigl\langle N_{3}(t) \bigr\rangle ^{*}+t^{-1}\vartheta _{2}(t) $$
(3.12)

and

$$ t^{-1} \bigl[\ln N_{3}(t)-\ln N_{3}(0) \bigr]\geq b_{3}+a_{32} \bigl\langle N_{2}(t) \bigr\rangle -a_{33} \bigl\langle N_{3}(t) \bigr\rangle ^{*}+t^{-1}\vartheta _{3}(t). $$
(3.13)

Using Lemma 2.1 to (3.12) yields

$$ \bigl\langle N_{2}(t) \bigr\rangle _{*}\geq \frac{b_{2}-a_{23}\langle N_{3}(t)\rangle ^{*}}{a_{22}} = \frac{b_{2}a_{33}-b_{3}a_{23}}{a_{22}a_{33} +a_{32}a_{23}}=\frac{\Delta _{1}^{*}-\tilde{\Delta }_{1}^{*}}{A_{11}}. $$
(3.14)

From (3.13) and (3.14), then

$$ t^{-1} \bigl[\ln N_{3}(t)-\ln N_{3}(0) \bigr]\geq b_{3}+a_{32} \bigl\langle N_{2}(t) \bigr\rangle _{*} -a_{33} \bigl\langle N_{3}(t) \bigr\rangle +t^{-1}\vartheta _{3}(t). $$
(3.15)

Applying Lemma 2.1 to (3.15) yields

$$ \bigl\langle N_{2}(t) \bigr\rangle _{*}\geq \frac{b_{3}+a_{32}\langle N_{2}(t)\rangle _{*}}{a_{33}} = \frac{b_{2}a_{32}+b_{3}a_{22}}{a_{22}a_{33}+a_{32}a_{23}} = \frac{\Delta _{1}-\tilde{\Delta }_{1}}{A_{11}}. $$

Then we have

$$ \lim_{t\rightarrow \infty }N_{1}(t)=0,\qquad \lim_{t \rightarrow \infty } \bigl\langle N_{2}(t) \bigr\rangle = \frac{\Delta _{1}^{*}-\tilde{\Delta } _{1}^{*}}{A_{11}},\qquad \lim _{t\rightarrow \infty } \bigl\langle N_{3}(t) \bigr\rangle = \frac{\Delta _{1}-\tilde{\Delta }_{1}}{A_{11}}. $$

Therefore, case (ii) is proved. The proof of case (iii) is similar to case (ii) and we omit it here.

Next we enter the proof of case (iv). We begin to eliminate \(\langle N_{1}(t)\rangle \), \(\langle N_{2}(t)\rangle \) from (3.2)–(3.4) by the elimination method. By analysis, there exist positive constants \(p=A_{13}/A_{33}>0\), \(q=-A_{23}/A_{33}>0\), multiplying both sides of (3.2)–(3.4) by p, q, and 1, respectively, adding the three inequalities yields

$$ \begin{aligned} & \frac{ \ln N_{3}(t)-\ln N_{3}(0)+p (\ln N_{1}(t) -\ln N_{1}(0) )+q (\ln N_{2}(t)-\ln N_{2}(0) )}{t} \\ &\quad = b_{1}p+b_{2}q+b_{3}-(a_{13}p+a_{23}q-a_{33}) \bigl\langle N_{3}(t) \bigr\rangle +t^{-1} \bigl(\vartheta _{1}(t)p+\vartheta _{2}(t)q+\vartheta _{3}(t) \bigr) \\ & \quad = \frac{A_{3}-\tilde{A}_{3}}{A_{33}}- \frac{A}{A_{33}}\bigl\langle N_{3}(t) \bigr\rangle +t^{-1} \bigl( \vartheta _{1}(t)p+\vartheta _{2}(t)q+\vartheta _{3}(t)\bigr). \end{aligned}$$
(3.16)

Similarly, by the elimination method, there exist constants

$$ \tilde{p}=\frac{- A_{21}}{A_{11}}< 0, \qquad \tilde{q}= \frac{A_{31}}{A_{11}}< 0,\qquad \bar{p}= \frac{- A_{12}}{A_{22}}< 0,\qquad \bar{q}=\frac{ -A_{32}}{A_{22}}< 0 $$

such that

$$ \begin{aligned} & \frac{ \ln N_{1}(t)-\ln N_{1}(0)+\tilde{p} (\ln N_{2}(t)-\ln N_{2}(0) )+\tilde{q} (\ln N_{3}(t)-\ln x_{3}(0) )}{t} \\ & \quad = \frac{A_{1}-\tilde{A}_{1}}{A_{11}} - \frac{A}{A_{11}}\bigl\langle N_{1}(t) \bigr\rangle +t^{-1} \bigl( \vartheta _{1}(t)\tilde{p}+ \vartheta _{2}(t)\tilde{q}+\vartheta _{3}(t)\bigr) \end{aligned}$$
(3.17)

and

$$ \begin{aligned} & \frac{\ln N_{2}(t)-\ln N_{2}(0)+ \bar{p} (\ln N_{1}(t)-\ln N_{1}(0) ) +\bar{q} (\ln N_{3}(t)-\ln N_{3}(0) )}{t} \\ &\quad = \frac{A_{2}-\tilde{A}_{2}}{A_{22}} - \frac{A}{A_{22}}\bigl\langle N_{2}(t) \bigr\rangle +t^{-1} \bigl( \vartheta _{1}(t)\bar{p}+\vartheta _{2}(t)\bar{q}+\vartheta _{3}(t)\bigr). \end{aligned} $$
(3.18)

Using Lemma 2.3 in equality (3.16), for arbitrarily \(\varepsilon >0\), there exists \(T>0\), for all \(t>T\), we have

$$\begin{aligned}& t^{-1} \bigl(p \bigl[\ln N_{1}(t)-\ln N_{1}(0) \bigr]+q \bigl[\ln N_{2}(t)-\ln N_{2}(0) \bigr] \bigr) \\& \quad \leq t^{-1} \bigl(p\ln N_{1}(t)+q\ln N_{2}(t) \bigr)+\varepsilon \leq \varepsilon . \end{aligned}$$
(3.19)

Substituting (3.19) into (3.16) leads to

$$\begin{aligned}& t^{-1} \bigl[\ln N_{3}(t)-\ln N_{3}(0) \bigr] \\& \quad \geq \frac{A_{3}-\tilde{A}_{3}}{A_{33}}-\varepsilon -\frac{A}{A_{33}} \bigl\langle N_{3}(t) \bigr\rangle +t^{-1} \bigl( \vartheta _{1}(t)p+\vartheta _{2}(t)q+\vartheta _{3}(t) \bigr). \end{aligned}$$
(3.20)

Since \(A_{3}>\tilde{A}_{3}\), letting \(\varepsilon >0\) be small enough such that \(A_{3}-\tilde{A}_{3}-\varepsilon >0\), then by Lemma 2.1, we have

$$ \bigl\langle N_{3}(t) \bigr\rangle _{*}\geq \frac{A_{3}-\tilde{A}_{3}}{A}. $$

Similarly, we derive from (3.2) and (3.3) that

$$\begin{aligned}& t^{-1} \bigl[\ln N_{1}(t)-\ln N_{1}(0) \bigr] \\& \quad \leq \frac{A_{1}-\tilde{A}_{1}}{A_{11}} +\varepsilon -\frac{A}{A_{11}} \bigl\langle N_{1}(t) \bigr\rangle +t^{-1} \bigl(\vartheta _{1}(t)\tilde{p}+ \vartheta _{2}(t)\tilde{q}+\vartheta _{3}(t) \bigr) \end{aligned}$$
(3.21)

and

$$\begin{aligned}& t^{-1} \bigl[\ln N_{2}(t)-\ln N_{2}(0) \bigr] \\& \quad \leq \frac{A_{2}-\tilde{A}_{2}}{A_{22}}+\varepsilon -\frac{A}{A_{22}} \bigl\langle N_{2}(t) \bigr\rangle +t^{-1} \bigl(\vartheta _{1}(t)\bar{p}+\vartheta _{2}(t) \bar{q}+\vartheta _{3}(t) \bigr). \end{aligned}$$
(3.22)

Applying Lemma 2.1 to (3.21) and (3.22) again, for sufficiently large t, we obtain

$$ \bigl\langle N_{1}(t) \bigr\rangle ^{*}\leq \frac{A_{1}-\tilde{A}_{1}}{A}, \qquad \bigl\langle N_{2}(t) \bigr\rangle ^{*} \leq \frac{A_{2}-\tilde{A}_{2}}{A}. $$

By the definition of sup limit, we deduce from (3.4) that

$$ t^{-1} \bigl[\ln N_{3}(t)-\ln N_{3}(0) \bigr] \leq b_{3}+a_{31} \bigl\langle N_{1}(t) \bigr\rangle ^{*} +a_{32} \bigl\langle N_{2}(t) \bigr\rangle ^{*}-a_{33} \bigl\langle N_{3}(t) \bigr\rangle +t^{-1}\vartheta _{3}(t). $$

Therefore, Lemma 2.1 implies

$$ \bigl\langle N_{3}(t) \bigr\rangle ^{*}\leq \frac{A_{3}-\tilde{A}_{3}}{A}. $$

By the same way, from (3.2) and (3.3), we obtain

$$\begin{aligned}& t^{-1} \bigl[\ln N_{1}(t)-\ln N_{1}(0) \bigr] \\& \quad \geq b_{1}-a_{11} \bigl\langle N_{1}(t) \bigr\rangle -a_{12} \bigl\langle N_{2}(t) \bigr\rangle ^{*}-a_{13} \bigl\langle N_{3}(t) \bigr\rangle ^{*}+t^{-1}\vartheta _{1}(t) \end{aligned}$$
(3.23)

and

$$\begin{aligned}& t^{-1} \bigl[\ln N_{2}(t)-\ln N_{2}(0) \bigr] \\& \quad \geq b_{2}-a_{21} \bigl\langle N_{1}(t) \bigr\rangle ^{*}-a_{22} \bigl\langle N_{2}(t) \bigr\rangle -a_{23} \bigl\langle N_{3}(t) \bigr\rangle ^{*}+t^{-1}\vartheta _{2}(t). \end{aligned}$$
(3.24)

Substituting \(\langle N_{1}(t)\rangle ^{*}\leq \frac{A_{1}-\tilde{A}_{1}}{A}\), \(\langle N_{2}(t)\rangle ^{*}\leq \frac{A_{2}-\tilde{A}_{2}}{A}\), \(\langle N_{3}(t)\rangle ^{*} \leq \frac{A_{3}-\tilde{A}_{3}}{A}\) into (3.23) and (3.24) and using Lemma 2.1, we have

$$ \bigl\langle N_{1}(t) \bigr\rangle _{*}\geq \frac{A_{1}-\tilde{A}_{1}}{A},\qquad \bigl\langle N_{2}(t) \bigr\rangle _{*} \geq \frac{A_{2}-\tilde{A}_{2}}{A}. $$

Therefore,

$$ \lim_{t\rightarrow \infty } \bigl\langle N_{1}(t) \bigr\rangle = \frac{A_{1}-\tilde{A}_{1}}{A}, \qquad \lim_{t\rightarrow \infty } \bigl\langle N_{2}(t) \bigr\rangle =\frac{A_{2}-\tilde{A}_{2}}{A},\qquad \lim_{t\rightarrow \infty } \bigl\langle N_{3}(t) \bigr\rangle = \frac{A_{3}-\tilde{A}_{3}}{A}, $$

which is the required assertion.

If \(b_{1}>0\), \(b_{2}>0\), \(A_{3}<\tilde{A}_{3}\), then the proof is similar to case (iii), and we omit it here. The proof is completed. □

Remark 3.1

By the process of our proof, if considering the effect of Lévy jumps, one can also establish sufficient conditions preserving the stability in mean and extinction of all species. Here we move some restricting conditions like \(R>0\) and \(b_{1}>b_{2}\), which appeared in [6].

4 Stability in distribution

Theorem 4.1

The solution of model (1.2) is a stationary Markov process, that is, there exists a stationary distribution for system (1.2) if Assumption 2.2holds.

Proof

Define

$$ \hat{V}=R_{1}\frac{N_{1}^{p}}{p}+R_{2}\frac{N_{2}^{p}}{p}+ R_{3} \biggl( \frac{N_{3}^{p}}{p}+\frac{a_{31}y_{1}^{p+1}}{\sigma _{1}(p+1)} + \frac{a_{32}y_{2}^{p+1}}{\sigma _{2}(p+1)} \biggr), $$

where \(R_{1}\), \(R_{2}\), \(R_{3}\) are positive constants defined later. By Itô’s formula, we have

$$\begin{aligned} L\hat{V}(t) = {}& R_{1}N_{1}(t)^{p} \biggl(r_{1}-a_{11}N_{1}(t)-a_{12}y_{2}(t)-a_{13}y_{3}(t)+ \frac{p-1}{2}\xi _{1}^{2} \biggr) \\ & {} +R_{2}N_{2}(t)^{p} \biggl(r_{2}-a_{21}y_{1}(t)-a_{22}N_{2}(t)-a_{23}y_{3}(t)+ \frac{p-1}{2}\xi _{2}^{2} \biggr) \\ & {} + R_{3}N_{3}(t)^{p} \biggl(-r_{3}+a_{31}y_{1}(t)+a_{32}y_{2}(t)-a_{33}N_{3}(t)+ \frac{p-1}{2}\xi _{1}^{2} \biggr) \\ & {} +R_{3} \frac{a_{31}}{p+1} \bigl(N_{1}^{p+1}(t)-y_{1}^{p+1}(t) \bigr) +R_{3} \frac{a_{32}}{p+1} \bigl(N_{2}^{p+1}(t)-y_{2}^{p+1}(t) \bigr) \\ \leq {}& R_{1}N_{1}(t)^{p} \biggl(r_{1}+ \frac{p-1}{2}\xi _{1}^{2}-a_{11}N_{1}(t) \biggr) +R_{2}N_{2}(t)^{p} \biggl(r_{2}+ \frac{p-1}{2}\xi _{2}^{2}-a_{22}N_{2}(t) \biggr) \\ & {} +R_{3}N_{3}(t)^{p} \biggl(-r_{3}+ \frac{p-1}{2} \xi _{1}^{2}-a_{33}N_{3}(t) \biggr) +R_{3}a_{31} \frac{pN_{3}^{p+1}(t)+y_{1}^{p+1}(t)}{p+1} \\ & {} +R_{3}a_{32} \frac{pN_{3}^{p+1}(t)+y_{2}^{p+1}(t)}{p+1} +R_{3} \frac{a_{31}}{p+1}\bigl(N_{1}^{p+1}(t) -y_{1}^{p+1}(t) \bigr) \\ & {} +R_{3} \frac{a_{32}}{p+1}\bigl(N_{2}^{p+1}(t)-y_{2}^{p+1}(t) \bigr) \\ ={} & R_{1}N_{1}(t)^{p} \biggl(r_{1}+ \frac{p-1}{2} \xi _{1}^{2}-a_{11}N_{1}(t) \biggr) +R_{2}N_{2}(t)^{p} \biggl(r_{2}+ \frac{p-1}{2}\xi _{2}^{2}-a_{22}N_{2}(t) \biggr) \\ & {} +R_{3}N_{3}(t)^{p} \biggl(-r_{3}+ \frac{p-1}{2} \xi _{1}^{2}-a_{33}N_{3}(t) \biggr) +R_{3} \frac{pN_{3}^{p+1}(t)}{p+1} (a_{31} +a_{32}) \\ & {} +R_{3} \frac{a_{31}}{p+1}N_{1}^{p+1}(t) +R_{3} \frac{a_{32}}{p+1}N_{2}^{p+1}(t) \\ = {}& R_{1}N_{1}(t)^{p} \biggl(r_{1}+ \frac{p-1}{2}\xi _{1}^{2} \biggr) +R_{2}N_{2}(t)^{p} \biggl(r_{2}+ \frac{p-1}{2}\xi _{2}^{2} \biggr) \\ & {} +R_{3}N_{3}(t)^{p} \biggl(-r_{3}+ \frac{p-1}{2} \xi _{1}^{2} \biggr)+ \biggl(-R_{1}a_{11}+R_{3} \frac{a_{31}}{p+1} \biggr)N_{1}^{p+1}(t) \\ & {} + \biggl(-R_{2}a_{22}+R_{3} \frac{a_{32}}{p+1} \biggr)N_{2}^{p+1}(t)+ \biggl(-R_{3}a_{33}+R_{3} \frac{p(a_{31} +a_{32})}{p+1} \biggr)N_{3}^{p+1}(t). \end{aligned}$$

By Assumption 2.2, there exist positive constants

$$\begin{aligned}& R_{1}=\frac{p+1+R_{3}a_{31}}{(p+1)a_{11}}>0,\qquad R_{2}= \frac{p+1+R_{3}a_{32}}{(p+1)a_{22}}>0, \\& R_{3}=\frac{p+1}{(p+1)a_{33}-p(a_{31}+a_{32})}>0 \end{aligned}$$

such that

$$\begin{aligned} L\hat{V} \leq{} & - \bigl(N_{1}^{p+1}(t)+N_{2}^{p+1}(t)+N_{3}^{p+1}(t) \bigr) +R_{1}N_{1}(t)^{p} \biggl(r_{1}+ \frac{p-1}{2}\xi _{1}^{2} \biggr) \\ & {} +R_{2}N_{2}(t)^{p} \biggl(r_{2}+ \frac{p-1}{2} \xi _{2}^{2} \biggr) +R_{3}N_{3}(t)^{p} \biggl(-r_{3}+ \frac{p-1}{2}\xi _{1}^{2} \biggr). \end{aligned}$$

Define \(\breve{V}=\hat{V}+\sum_{i=1}^{3}\frac{y_{i}^{p+1}}{2\sigma _{i}}\), then

$$\begin{aligned} L\breve{V} \leq {}& - \frac{1}{2} \bigl(N_{1}^{p+1}(t)+N_{2}^{p+1}(t)+N_{3}^{p+1}(t)+y_{1}^{p+1}(t)+y_{2}^{p+1}(t)+y_{3}^{p+1}(t) \bigr) \\ & {} +R_{1}N_{1}(t)^{p} \biggl(r_{1}+ \frac{p-1}{2} \xi _{1}^{2} \biggr) +R_{2}N_{2}(t)^{p} \biggl(r_{2}+ \frac{p-1}{2}\xi _{2}^{2} \biggr) \\ & {} +R_{3}N_{3}(t)^{p} \biggl(-r_{3}+ \frac{p-1}{2} \xi _{1}^{2} \biggr). \end{aligned}$$

Let

$$ \tilde{V}(t)=\sum_{i=1}^{3} \biggl( \frac{1}{N_{i}^{\iota }(t)}- \ln y_{i}(t) \biggr). $$

Applying Itô’s formula to \(\tilde{V}(t)\) yields

$$\begin{aligned} L\tilde{V} = & - \frac{\sigma _{1}(N_{1}(t)-y_{1}(t))}{y_{1}(t)} - \frac{\sigma _{2}(N_{2}(t)-y_{2}(t))}{y_{2}(t)} - \frac{\sigma _{3}(N_{3}(t)-y_{3}(t))}{y_{3}(t)} \\ & {}-\iota \bigl(N_{1}(t)\bigr)^{-\iota } \biggl(r_{1}-a_{11}N_{1}(t)-a_{12}y_{2}(t) -a_{13}y_{3}(t)- \frac{\iota +1}{2}\xi _{1}^{2} \biggr) \\ & {}-\iota \bigl(N_{2}(t)\bigr)^{-\iota } \biggl(r_{2}-a_{21}y_{1}(t)-a_{22}N_{2}(t) -a_{23}y_{3}(t)- \frac{\iota +1}{2}\xi _{2}^{2} \biggr) \\ & {}-\iota \bigl(N_{3}(t)\bigr)^{-\iota } \biggl(-r_{3}+a_{31}y_{1}(t) +a_{32}y_{2}(t)-a_{33}N_{3}(t)- \frac{\iota +1}{2}\xi _{3}^{2} \biggr) \\ \leq & -\iota \frac{r_{1}-\frac{\iota +1}{2}\xi _{1}^{2}}{N_{1}^{\iota }} +a_{11} \iota N_{1}^{1-\iota }(t) +a_{12}\iota \frac{y_{2}^{2}+N_{1}^{-2\iota }}{2} +a_{13}\iota \frac{y_{3}^{2}+N_{1}^{-2\iota }}{2} \\ & {}- \frac{\sigma _{1}(N_{1}(t)-y_{1}(t))}{y_{1}(t)} +a_{22}\iota N_{2}^{1- \iota }(t) +a_{21}\iota \frac{y_{1}^{2}+N_{2}^{-2\iota }}{2} +a_{23}\iota \frac{y_{3}^{2}+N_{2}^{-2\iota }}{2} \\ & {}- \frac{\sigma _{2}(N_{2}(t)-y_{2}(t))}{y_{2}(t)} -\iota \frac{-r_{3}- \frac{\iota +1}{2}\xi _{3}^{2}}{N_{3}^{\iota }} +a_{33}\iota N_{3}^{1-\iota }(t) - \frac{\sigma _{3}(N_{3}(t)-y_{3}(t))}{y_{3}(t)} \\ = & \sigma _{1}+\sigma _{2}+\sigma _{3}-\iota \frac{r_{1} - \frac{\iota +1}{2}\xi _{1}^{2}}{N_{1}^{\iota }} -\iota \frac{r_{2}- \frac{\iota +1}{2}\xi _{2}^{2}}{N_{2}^{\iota }} -\iota \frac{-r_{3}- \frac{\iota +1}{2}\xi _{3}^{2}}{N_{3}^{\iota }} \\ & {}+a_{21}\iota y_{1}^{2}/2+a_{12} \iota y_{2}^{2}/2 +(a_{13}+a_{23}) \iota y_{3}^{2}/2 +(a_{12}+a_{13})\iota N_{1}^{-2\iota }/{2} \\ & {}+ (a_{21}+a_{23})\iota N_{2}^{-2\iota }/{2} +a_{11}\iota N_{1}^{1- \iota }(t) \\ &{} +a_{22}\iota N_{2}^{1-\iota }(t)+a_{33} \iota N_{3}^{1-\iota }(t) - \frac{\sigma _{1}N_{1}(t)}{y_{1}(t)} - \frac{\sigma _{2}N_{2}(t) }{y_{2}(t)} - \frac{\sigma _{3}N_{3}(t)}{y_{3}(t)}. \end{aligned}$$

Define \(V(t)=\breve{V}(t)+\tilde{V}(t)\), then

$$\begin{aligned} LV = & L\breve{V}+L\tilde{V} \\ \leq & \sigma +M- \frac{1}{4} \bigl(N_{1}^{p+1}(t) +N_{2}^{p+1}(t)+N_{3}^{p+1}(t)+y_{1}^{p+1}(t) +y_{2}^{p+1}(t)+y_{3}^{p+1}(t) \bigr) \\ & {}-\iota \frac{r_{1}- \frac{\iota +1}{2}\xi _{1}^{2}}{N_{1}^{\iota }} -\iota \frac{r_{2}- \frac{\iota +1}{2}\xi _{2}^{2}}{N_{2}^{\iota }} -\iota \frac{-r_{3}- \frac{\iota +1}{2}\xi _{3}^{2}}{N_{3}^{\iota }} \\ & {}- \frac{\sigma _{1}N_{1}(t)}{y_{1}(t)} - \frac{\sigma _{2}N_{2}(t) }{y_{2}(t)} - \frac{\sigma _{3}N_{3}(t)}{y_{3}(t)}, \end{aligned}$$

where \(\sigma =\sigma _{1}+\sigma _{2}+\sigma _{3}\), and

$$\begin{aligned} M = & \max_{N_{i},y_{i}\in R_{+}^{6}} \Biggl\{ - \frac{1}{4}\sum _{i=1}^{3} \bigl(N_{i}^{p+1}(t)+y_{i}^{p+1}(t) \bigr) +R_{1}N_{1}(t)^{p} \biggl(r_{1}+ \frac{p-1}{2}\xi _{1}^{2} \biggr) \\ & {}+R_{2}N_{2}(t)^{p} \biggl(r_{2}+ \frac{p-1}{2} \xi _{2}^{2} \biggr) +R_{3}N_{3}(t)^{p} \biggl(-r_{3}+ \frac{p-1}{2}\xi _{1}^{2} \biggr) \\ &{} + \frac{a_{21}\iota y_{1}^{2}}{2}+ \frac{a_{12}\iota y_{2}^{2}}{2} + \frac{(a_{13}+a_{23})\iota y_{3}^{2}}{2} + \frac{(a_{12}+a_{13})\iota N_{1}^{-2\iota }}{2} \\ &{} + \frac{(a_{21}+a_{23})\iota N_{2}^{-2\iota }}{2} +a_{11}\iota N_{1}^{1- \iota }(t) +a_{22}\iota N_{2}^{1-\iota }(t) +a_{33}\iota N_{3}^{1-\iota }(t) \Biggr\} . \end{aligned}$$

Choose \(\varepsilon >0\) small enough such that

$$ 0< \varepsilon < \min \biggl\{ \biggl( \frac{\iota (r_{i}-\frac{\xi _{i}^{2}}{2}(\iota +1)}{ \sigma +M+1} \biggr)^{\frac{1}{\iota }}, \biggl( \frac{1}{4(\sigma +M+1)} \biggr)^{ \frac{1}{p+1}}, \frac{\sigma _{i}}{\sigma +M+1},i=1,2,3 \biggr\} . $$

Define the following bounded closed set:

$$ D_{\varepsilon }= \biggl\{ (N_{1},N_{2},N_{3},y_{1},y_{2},y_{3}) \in R_{+}^{6} \Big|\varepsilon < N_{i}< \frac{1}{\varepsilon }, {\varepsilon }^{2}< y_{i}< \frac{1}{\varepsilon },i=1,2,3 \biggr\} , $$

and for \(i=1,2,3\), denote

$$\begin{aligned}& D_{\varepsilon }^{i}= \biggl\{ (N_{1},N_{2},N_{3},y_{1},y_{2},y_{3}) \in R_{+}^{6}\Big|N_{i}> \frac{1}{\varepsilon } \biggr\} , \\& D_{\varepsilon }^{i+3}= \biggl\{ (N_{1},N_{2},N_{3},y_{1},y_{2},y_{3}) \in R_{+}^{6}\Big|y_{i}> \frac{1}{\varepsilon } \biggr\} , \\& D_{\varepsilon }^{i+6}= \bigl\{ (N_{1},N_{2},N_{3},y_{1},y_{2},y_{3}) \in R_{+}^{6}\Big|0< N_{i}< \varepsilon \bigr\} , \\& D_{\varepsilon }^{i+9}= \biggl\{ (N_{1},N_{2},N_{3},y_{1},y_{2},y_{3}) \in R_{+}^{6}\Big|\varepsilon < N_{i}< \frac{1}{\varepsilon }, 0< y_{i}< { \varepsilon }^{2} \biggr\} . \end{aligned}$$

Denote the complement of \(D_{\varepsilon }\) by \(D_{\varepsilon }^{C}\), then it is easy to get \(D_{\varepsilon }^{C}=\bigcup_{j=1}^{12} D_{\varepsilon }^{j}\). For all \((N_{1},N_{2},N_{3},y_{1},y_{2},y_{3})\in D_{\varepsilon }^{C}\), we discuss as follows.

  1. (i)

    If \((N_{1},N_{2},N_{3},y_{1},y_{2},y_{3})\in D_{\varepsilon }^{i}\), \(i=1,2,3\), then

    $$ LV\leq \sigma +M-\frac{1}{4}N_{i}^{p}\leq \sigma +M- \frac{1}{4\varepsilon ^{p+1}}< -1; $$
  2. (ii)

    If \((N_{1},N_{2},N_{3},y_{1},y_{2},y_{3})\in D_{\varepsilon }^{i+3}\), \(i=1,2,3\), then

    $$ LV\leq \sigma +M-\frac{1}{4}y_{i}^{p}\leq \sigma +M- \frac{1}{4\varepsilon ^{p+1}}< -1; $$
  3. (iii)

    If \((N_{1},N_{2},N_{3},y_{1},y_{2},y_{3})\in D_{\varepsilon }^{i+6}\), \(i=1,2,3\), then

    $$ LV\leq \sigma +M-\iota \frac{r_{i}-\frac{\xi _{i}^{2}}{2}(\iota +1)}{N_{i}^{\iota }}\leq \sigma +M-\iota \frac{r_{i}-\frac{\xi _{i}^{2}}{2}(\iota +1)}{{\varepsilon }^{\iota }} < -1; $$
  4. (iv)

    If \((N_{1},N_{2},N_{3},y_{1},y_{2},y_{3})\in D_{\varepsilon }^{i+9}\), \(i=1,2,3\), then

    $$ LV\leq \sigma +M-\sigma _{i} N_{i}/y_{i}\leq \sigma +M-\sigma _{i} \varepsilon /\varepsilon ^{2} < -1. $$

Consequently, for any \((N_{1},N_{2},N_{3},y_{1},y_{2},y_{3})\in D_{\varepsilon }^{C}\), we have

$$ \sup_{(N_{1},N_{2},N_{3},y_{1},y_{2},y_{3})\in R_{+}^{6}} LV(N_{1},N_{2},N_{3},y_{1},y_{2},y_{3}) \leq -1. $$

Therefore, it follows from Lemma 2.4 that there exists a stationary distribution for system (1.2). The proof is completed. □

Theorem 4.2

Under Assumption 2.2, solutions of model (1.2) are globally attractive.

Proof

Firstly, let \(N(t)=N(t,N(\phi ))\) and \(\bar{N}(t)=\bar{N}(t,\bar{N}(\phi ))\) be any two solutions of model (1.1) with the initial data \(N(\phi ), \bar{N}(\phi )\in C([-\tau ,0],R_{+}^{3})\). We only need to prove \(\lim_{t\rightarrow \infty }\mathbb{ E}|N_{i}(t)-\bar{N}_{i}(t)|=0\) for \(i=1,2,3\).

Define

$$ V(t)=\sum_{i=1}^{3}D_{i} \bigl\vert \ln N_{i}(t)-\ln \bar{N}_{i}(t) \bigr\vert , $$

where \(D_{i}\) (\(i=1,2,3\)) is defined later in the proof. By computing the upper right derivative of \(V(t)\), then

$$ \begin{aligned} D^{+}V(t)\leq {}& D_{1} \mathrm{sign} \bigl(N_{1}(t)-\bar{N}_{1}(t) \bigr)\bigl[-a_{11} \bigl(N_{1}(t)-\bar{N}_{1}(t) \bigr)-a_{12}\bigl(y_{2}(t)-\bar{y}_{2}(t)\bigr) \\ & {}-a_{13}\bigl(y_{3}(t)-\bar{y}_{3}(t)\bigr) \bigr]\,dt +D_{2} \mathrm{sign} \bigl(N_{2}(t)- \bar{N}_{2}(t)\bigr) \\ & {}\times \bigl[-a_{21} \bigl(y_{1}(t)- \bar{y}_{1}(t)\bigr)-a_{22}\bigl(N_{2}(t)- \bar{N}_{2}(t)\bigr) -a_{23}\bigl(y_{3}(t)- \bar{y}_{3}(t)\bigr)\bigr]\,dt \\ & {}+D_{3} \mathrm{sign} \bigl(N_{3}(t)- \bar{N}_{3}(t)\bigr)\bigl[a_{31} \bigl(y_{1}(t)- \bar{y}_{1}(t)\bigr)+a_{32}\bigl(y_{2}(t)- \bar{y}_{2}(t)\bigr) \\ & {}-a_{33}\bigl(Ny_{3}(t)-\bar{N}_{3}(t)\bigr) \bigr]\,dt. \end{aligned} $$
(4.1)

On the other hand, by (1.2), we have

$$ \frac{d(y_{i}-\bar{y}_{i})}{dt}=\sigma _{i}(N_{i}-\bar{N}_{i})- \sigma _{i}(y_{i}-\bar{y}_{i}),\quad i=1,2,3, $$

that is,

$$ y_{i}(t)-\bar{y}_{i}(t)=e^{-\sigma _{1}t} \bigl(y_{i}(0)-\bar{y}_{i}(0) \bigr)+ \sigma _{1}e^{-\sigma _{1}t} \int _{0}^{t} e^{\sigma _{1}s} \bigl(N_{i}(s)- \bar{N}_{i}(s) \bigr)\,ds,\quad i=1,2,3. $$

Therefore,

$$ \bigl\vert y_{i}(t)-\bar{y}_{i}(t) \bigr\vert \leq e^{-\sigma _{1}t} \bigl\vert y_{i}(0)-\bar{y}_{i}(0) \bigr\vert + \sigma _{1}e^{-\sigma _{1}t} \int _{0}^{t} e^{\sigma _{1}s} \bigl\vert N_{i}(s)- \bar{N}_{i}(s) \bigr\vert \,ds,\quad i=1,2,3. $$

Integrating two sides of the above inequality from 0 to t, we have

$$\begin{aligned} \int _{0}^{t} \bigl\vert y_{i}(t)- \bar{y}_{i}(t) \bigr\vert \leq & \frac{1-e^{-\sigma _{1}t}}{\sigma _{i}} \bigl\vert y_{i}(0)- \bar{y}_{i}(0) \bigr\vert +\sigma _{i} \int _{0}^{t}\,dv \int _{0}^{v} e^{\sigma _{1}(s-v)} \bigl\vert N_{i}(s)-\bar{N}_{i}(s) \bigr\vert \,ds \\ = & \frac{1-e^{-\sigma _{1}t}}{\sigma _{i}} \bigl\vert y_{i}(0)- \bar{y}_{i}(0) \bigr\vert + \int _{0}^{t}\bigl(1- e^{\sigma _{1}(s-t)}\bigr) \bigl\vert N_{i}(s)- \bar{N}_{i}(s) \bigr\vert \,ds \\ \leq & \frac{1}{\sigma _{i}} \bigl\vert y_{i}(0)- \bar{y}_{i}(0) \bigr\vert + \int _{0}^{t} \bigl\vert N_{i}(s)- \bar{N}_{i}(s) \bigr\vert \,ds,\quad i=1,2,3. \end{aligned}$$

Integrating both sides of (4.1) from 0 to t and taking expectations give

$$\begin{aligned} V(t) \leq & V(0)+D_{1} \biggl(-a_{11} \int _{0}^{t} \bigl\vert N_{1}(s)- \bar{N}_{1}(s) \bigr\vert \,ds+a_{12} \int _{0}^{t} \bigl\vert y_{2}(s)- \bar{y}_{2}(s) \bigr\vert \,ds \\ & {}+a_{13} \int _{0}^{t} \bigl\vert y_{3}(s)- \bar{y}_{3}(s) \bigr\vert \,ds \biggr) +D_{2} \biggl(-a_{22} \int _{0}^{t} \bigl\vert N_{2}(s)- \bar{N}_{2}(s) \bigr\vert \,ds \\ & {} +a_{21} \int _{0}^{t} \bigl\vert y_{1}(s) - \bar{y}_{1}(s) \bigr\vert \,ds+a_{23} \int _{0}^{t} \bigl\vert y_{3}(s)- \bar{y}_{3}(s) \bigr\vert \,ds \biggr) \\ & {}+D_{3} \biggl(-a_{33} \int _{0}^{t} \bigl\vert N_{3}(s)- \bar{N}_{3}(s) \bigr\vert \,ds+a_{31} \int _{0}^{t} \bigl\vert y_{1}(s) - \bar{y}_{1}(s) \bigr\vert \,ds\\ & {} +a_{32} \int _{0}^{t} \bigl\vert y_{2}(s)- \bar{y}_{2}(s) \bigr\vert \,ds \biggr) \\ \leq & V(0)+(D_{2}a_{31}+D_{3}a_{21}-D_{1}a_{11}) \int _{0}^{t} \bigl\vert N_{1}(s)- \bar{N}_{1}(s) \bigr\vert \,ds \\ &{} +(D_{1}a_{12}+D_{3}a_{32}-D_{2}a_{22}) \int _{0}^{t} \bigl\vert N_{2}(s)- \bar{N}_{2}(s) \bigr\vert \,ds \\ & {}+(D_{1}a_{13}+D_{2}a_{23}-D_{3}a_{33}) \int _{0}^{t} \bigl\vert N_{3}(s)- \bar{N}_{3}(s) \bigr\vert \,ds \\ & {}+ \frac{D_{2}a_{21}+D_{3}a_{31}}{\sigma _{1}} \bigl\vert y_{1}(0)- \bar{y}_{1}(0) \bigr\vert + \frac{D_{1}a_{12}+D_{3}a_{32}}{\sigma _{2}} \bigl\vert y_{2}(0)-\bar{y}_{2}(0) \bigr\vert \\ & {}+ \frac{D_{1}a_{13}+D_{2}a_{23}}{\sigma _{3}} \bigl\vert y_{3}(0)- \bar{y}_{3}(0) \bigr\vert . \end{aligned}$$

For the following equations,

$$ \textstyle\begin{cases} D_{1}a_{11}-D_{2}a_{21}-D_{3}a_{31}=1, \\ -D_{1}a_{12}+D_{2}a_{22}-D_{3}a_{32}=1, \\ -D_{1}a_{13}-D_{2}a_{23}+D_{3}a_{33}=1, \end{cases} $$

under Assumption 2.2, the coefficient matrix of \(D_{1}\), \(D_{2}\), and \(D_{3}\) is a nonsingular M-matrix, then by M-matrix theory, there exists \(D_{i}>0\) (\(i=1,2,3\)) satisfying the equation. Therefore,

$$\begin{aligned}& V(t)+ \int _{0}^{t} \bigl\vert N_{1}(s)- \bar{N}_{1}(s) \bigr\vert \,ds+ \int _{0}^{t} \bigl\vert N_{2}(s) - \bar{N}_{2}(s) \bigr\vert \,ds+ \int _{0}^{t} \bigl\vert N_{3}(s)- \bar{N}_{3}(s) \bigr\vert \,ds \\& \quad \leq V(0)+ \frac{D_{2}a_{21}+D_{3}a_{31}}{\sigma _{1}} \bigl\vert y_{1}(0)- \bar{y}_{1}(0) \bigr\vert + \frac{D_{1}a_{12}+D_{3}a_{32}}{\sigma _{2}} \bigl\vert y_{2}(0)- \bar{y}_{2}(0) \bigr\vert \\& \qquad {}+ \frac{D_{1}a_{13}+D_{2}a_{23}}{\sigma _{3}} \bigl\vert y_{3}(0)- \bar{y}_{3}(0) \bigr\vert \\& \quad < +\infty , \end{aligned}$$

which means \(|N_{i}(t)-\bar{N}_{i}(t)|\in \mathtt{L}_{1}[0,+\infty )\). Consequently, we can derive from Lemma 2.5 and Lemma 2.6 that

$$ \lim_{t\rightarrow \infty } \bigl\vert N_{i}(t)- \bar{N}_{i}(t) \bigr\vert =0, \quad i=1,2,3. $$

The proof is completed. □

Remark 4.1

Combining the existence of distribution and the global attractivity of solutions of (1.2), we conclude that system (1.2) has a unique distribution, which is stable.

5 Numerical simulations

In this section, we give some numerical simulations to validate our theoretical results. By the Milstein higher order method proposed by Higham [32], we numerically simulate the solutions of system (1.2). Using discretization Brownian path over \([0,T]\) and writing efficient Matlab codes, we can obtain the corresponding simulation figures one by one.

Let

$$\begin{aligned}& A=\begin{vmatrix} 0.4 & 0.1 & 0.15 \\ 0.2 & 0.4 & 0.1 \\ -0.15 & -0.1 & 0.4 \end{vmatrix} =0.0645,\qquad \sigma _{1}=\sigma _{2}=\sigma _{3}=1, \\& r_{1}=0.31,\qquad r_{2}=0.31,\qquad r_{3}=0.01. \end{aligned}$$

In the following, without special mention, we only change the parameter of white noise and keep the rest of parameters unchanged so as to clearly see the dynamical effect of white noise.

Case (i) \(b_{1}<0\), \(b_{2}<0\).

Let \(\xi _{1}=0.8\), \(\xi _{2}=0.8\), \(\xi _{3}=0.5292\), then an easy computation yields \(b1=-0.01\), \(b_{2}=-0.01\), \(b_{3}=-0.15\). It follows from Theorem 3.1 that all species are extinct, illustrated in Fig. 1.

Figure 1
figure 1

Extinction of all species of system (1.2) with \(\xi _{1}=0.8\), \(\xi _{2}=0.8\), \(\xi _{3}=0.5292\). (a) Extinction of \(N_{1}\), \(N_{2}\), \(N_{3}\). (b) Extinction of \(y_{1}\), \(y_{2}\), \(y_{3}\)

Case (ii) \(b_{1}<0\), \(b_{2}>0\).

Let \(\xi _{1}=0.8\), \(\xi _{2}=0.4\), \(\xi _{3}=0.01\), then \(b1=-0.01\), \(b_{2}=0.23\), \(b_{3}=-0.01\), and \(\Delta _{1}^{*}=0.019\), \(\tilde{\Delta }_{1}^{*}=0.093\), \(A_{11}=0.17\), \(\frac{\Delta _{1}-\tilde{\Delta }_{1}}{A_{11}}=0.019/0.17=0.1118\), \(\frac{\Delta _{1}^{*}-\tilde{\Delta }_{1}^{*}}{A_{11}}=0.093/0.17=0.5471\). By Theorem 3.1, then

$$\begin{aligned}& \lim_{t\rightarrow \infty }N_{1}(t)=0, \qquad \lim_{t \rightarrow \infty } \bigl\langle N_{2}(t) \bigr\rangle =0.5471, \\& \lim_{t\rightarrow \infty } \bigl\langle N_{3}(t) \bigr\rangle =0.1118, \end{aligned}$$

which is illustrated in Fig. 2(a).

Figure 2
figure 2

Dynamical behavior of (1.2). (a) is for (1.2) with \(\xi _{1}=0.8\), \(\xi _{2}=0.4\), \(\xi _{3}=0.01\), then \(\lim_{t\rightarrow \infty }N_{1}(t)=0\), \(\lim_{t \rightarrow \infty }\langle N_{2}(t)\rangle =0.5471\), \(\lim_{t \rightarrow \infty }\langle N_{3}(t)\rangle =0.1118\), (b) is for (1.2) with \(r_{3}=0.1\), \(\xi _{1}=0.8\), \(\xi _{2}=0.469\), \(\xi _{3}=0.01\), then \(N_{1}(t)\), \(N_{3}(t) \) are both extinct, \(\lim_{t\rightarrow \infty }\langle N_{2}(t)\rangle =0.5\)

If \(r_{3}=0.1\), \(\xi _{1}=0.8\), \(\xi _{2}=0.469\), \(\xi _{3}=0.01\), then \(b1=-0.01\), \(b_{2}=0.2\), \(b_{3}=-0.1\), \(b_{2}a_{32}+b_{3}a_{22} =-0.02<0\). Hence, Theorem 3.1 implies \(N_{1}(t)\), \(N_{3}(t) \) are both extinct, and species \(N_{2}(t)\) is stable in mean, and \(\lim_{t\rightarrow \infty }\langle N_{2}(t)\rangle =0.5\), see Fig. 2(b).

Case (iii) \(b_{1}>0\), \(b_{2}<0\).

Let \(\xi _{1}=0.4\), \(\xi _{2}=0.8\), \(\xi _{3}=0.01\), then \(b1=0.23\), \(b_{2}=-0.01\), \(b_{3}=-0.01\). By computation, \(\Delta _{2}-\tilde{\Delta }_{2} =0.0305\), \(\Delta _{2}^{*}-\tilde{\Delta }_{2}^{*}=0.0935\), \(A_{22}=0.1825\). Hence, it follows from Theorem 3.1 that \(N_{2}(t)\) is extinct, and

$$ \lim_{t\rightarrow \infty } \bigl\langle N_{1}(t) \bigr\rangle = \frac{\Delta _{2}^{*}-\tilde{\Delta }_{2}^{*}}{A_{22}}=0.1671, \qquad \lim_{t\rightarrow \infty } \bigl\langle N_{3}(t) \bigr\rangle = \frac{\Delta _{2}-\tilde{\Delta }_{2}}{A_{22}}=0.5123. $$

Figure 3(a) verifies it correctly.

Figure 3
figure 3

Dynamical behavior of (1.2). (a) is for system (1.2) with \(\xi _{1}=0.4\), \(\xi _{2}=0.8\), \(\xi _{3}=0.01\), then \(N_{2}(t)\) is extinct, and \(\lim_{t\rightarrow \infty }\langle N_{1}(t)\rangle =0.1671\), \(\lim_{t\rightarrow \infty }\langle N_{3}(t)\rangle =0.5123\). (b) is for system (1.2) with \(r_{3}=0.1\), \(\xi _{1}=0.469\), \(\xi _{2}=0.8\), \(\xi _{3}=0.01\), then \(N_{2}(t)\), \(N_{3}(t)\) are both extinct and \(\lim_{t\rightarrow \infty }\langle N_{1}(t)\rangle =0.5\)

If \(r_{3}=0.1\), \(\xi _{1}=0.469\), \(\xi _{2}=0.8\), \(\xi _{3}=0.01\), then \(b1=0.2\), \(b_{2}=-0.01\), \(b_{3}=-0.1\), and \(b_{1}a_{31}+b_{3}a_{11}-0.01\). Theorem 3.1 indicates that \(N_{2}(t)\), \(N_{3}(t)\) are both extinct and \(\lim_{t\rightarrow \infty }\langle N_{1}(t)\rangle =0.2/0.4=0.5\), see Fig. 3(b).

Case (iv) \(b_{1}>0\), \(b_{2}>0\).

We choose \(\xi _{1}=0.1414\), \(\xi _{2}=0.1414\), \(\xi _{3}=0.2\) such that \(b1=0.3\), \(b_{2}=0.3\), \(b_{3}=-0.03\). By computation, then \(A_{1}-\tilde{A}_{1} =0.036\), \(A_{2}-\tilde{A}_{2} =0.0266\), \(A_{3}- \tilde{A}_{3}=0.0153 \). Therefore, by Theorem 3.1, we have

$$\begin{aligned}& \lim_{t\rightarrow \infty } \bigl\langle N_{1}(t) \bigr\rangle = \frac{A_{1}-\tilde{A}_{1}}{A}=0.5581,\qquad \lim_{t \rightarrow \infty } \bigl\langle N_{2}(t) \bigr\rangle = \frac{A_{2}-\tilde{A}_{2}}{A}=0.4124, \\& \lim_{t\rightarrow \infty } \bigl\langle N_{3}(t) \bigr\rangle = \frac{A_{3}-\tilde{A}_{3}}{A}=0.2372, \end{aligned}$$

which is illustrated in Fig. 4.

Figure 4
figure 4

Stable case of (1.2) with \(\xi _{1}=0.1414\), \(\xi _{2}=0.1414\), \(\xi _{3}=0.2\). Theorem 3.1 shows that \(\lim_{t\rightarrow \infty }\langle N_{1}(t)\rangle = \frac{A_{1}-\tilde{A}_{1}}{A}=0.5581\), \(\lim_{t\rightarrow \infty }\langle N_{2}(t)\rangle =\frac{A_{2}-\tilde{A}_{2}}{A}=0.4124\), \(\lim_{t\rightarrow \infty }\langle N_{3}(t)\rangle = \frac{A_{3}-\tilde{A}_{3}}{A}=0.2372\)

By use of Theorems 4.1 and 4.2, we know that system (1.2) has a unique distribution, which is revealed in Figs. 5 and 6. Figure 5 is the probability density function of preys \(N_{1}(t)\), \(N_{2}(t)\) and predator \(N_{3}(t)\), respectively. Figure 6 shows the attractivity of the solutions. They both indicate the existence and stability of stationary distribution function. The simulation results verify that when the condition is satisfied, that is, the white noise is relatively small, system (1.2) is stable.

Figure 5
figure 5

Probability density function of preys \(N_{1}(t)\), \(N_{2}(t)\) and predator \(N_{3}(t)\) of (1.2) with \(\xi _{1}= 0.1414\), \(\xi _{2}=0.1414\), \(\xi _{3}=0.2\)

Figure 6
figure 6

Attractivity of solutions of system (1.2)

If \(\xi _{1}=0.1414\), \(\xi _{2}=0.1414\), \(\xi _{3}=0.5292\), then \(b1=0.3\), \(b_{2}=0.3\), \(b_{3}=-0.15\). By computation, we have \(A_{3}-\tilde{A}_{3}=-0.0015<0 \), \(\Delta _{3}-\tilde{\Delta }_{3}=0.06\), \(\Delta _{3}^{*}-\tilde{\Delta }_{3}^{*}=0.09\), \(A_{33}=0.14 \), which guarantees that the condition of case (iv) holds, and hence Theorem 3.1 implies the two-prey are stable in mean and the predator \(N_{3}\) is extinct, further,

$$ \lim_{t\rightarrow \infty } \bigl\langle N_{1}(t) \bigr\rangle = \frac{\Delta _{3}-\tilde{\Delta }_{3}}{A_{33}}=0.4286,\qquad \lim_{t\rightarrow \infty } \bigl\langle N_{2}(t) \bigr\rangle = \frac{\Delta _{3}^{*}-\tilde{\Delta }_{3}^{*}}{A_{33}}=0.6429 $$

and

$$ \lim_{t\rightarrow \infty } \bigl\langle N_{3}(t) \bigr\rangle =0. $$

Figure 7 indicates the result is true.

Figure 7
figure 7

Dynamical behavior of (1.2) with \(\xi _{1}=0.1414\), \(\xi _{2}=0.1414\), \(\xi _{3}=0.5292\). Theorem 3.1 implies the two prey are stable in mean and the predator \(N_{3}\) is extinct, and \(\lim_{t\rightarrow \infty }\langle N_{1}(t)\rangle = \frac{\Delta _{3}-\tilde{\Delta }_{3}}{A_{33}}=0.4286\), \(\lim_{t \rightarrow \infty }\langle N_{2}(t)\rangle = \frac{\Delta _{3}^{*}-\tilde{\Delta }_{3}^{*}}{A_{33}}=0.6429\), \(\lim_{t\rightarrow \infty }\langle N_{3}(t)\rangle =0\)

6 Conclusion and discussion

In this paper, we consider a three-species stochastic predator–prey system with distributed delays. Theorem 3.1 gives sufficient conditions of the stability in mean and extinction of each species. Theorems 4.1 and 4.2 give the existence and uniqueness of distribution of each species. Finally, by numerical simulations, we illustrate the validity of our theoretical results.

Theorem 3.1 implies that stochastic parameter \(\xi _{i}\) (\(i=1,2,3\)) has some important influences to the extinction, stability in mean of all species of (1.2), which is illustrated by our simulations clearly. Simulations reveal that small intensity of white noise strengthens the stability of (1.2), while large intensity of white noise will bring serious influence to the dynamical behavior.

Recently, regime switching appears in a biological system frequently, and many nice results have been obtained by many researchers. How about the white noise affecting the dynamical behavior of a predator–prey system with regime switching? We believe it is very interesting and leave it for our future work.