Dynamic analysis of a stochastic four species food-chain model with harvesting and distributed delay

Abstract

A stochastic four species food-chain model is proposed in this paper. Here, artificial harvest in each species and the effect of time delay for interaction between species are considered, which makes the model more applicable in real situations. Specifically, we address the stochastic global dynamics behavior, including the existence of global positive solutions, stochastic ultimate boundedness, extinction with probability one, persistence in mean and global stability. The asymptotic stability in the probability distribution is obtained, and the criterion for the existence and non-existence of the optimal harvesting strategy is also derived. Furthermore, this paper can provide reference for the research of general n-species stochastic food-chain models.

1 Introduction

Ecosystem of one species is very rare in nature. In natural ecosystems, the coexistence of a large number of species is almost universal (see [1]). Over the last few decades, two or three species systems such as predator–prey and food-chain systems have long been the main topic of mathematical ecology and ecology (see [2–5]). However, we have realized that many phenomena in nature cannot be described by an ecosystem with two species or three interacting species. It is extremely important to develop theoretical methods with four or more species (see [6, 7]). El-Owaidy et al. in [7] studied a four-level generalized food-chain model, they analyzed the existence of a bounded solution and investigated the stability of various equilibrium points. However, the complex dynamics of the model was not explored.

Nowadays, the harvesting policies and regulations of wildlife in various countries have been gradually established. The most important thing of such regulations is to formulate an optimal harvesting plan that integrates the three aspects of ecology, environment and economy. So, it is extremely important to develop theoretical methods to get optimal harvesting result (see [8–17]). Tuerxun et al. in [17] studied a stochastic two-predators one-prey system with distributed delays, harvesting and Lévy jumps. They mainly discussed the global dynamics and the optimal harvesting strategy, also obtained the optimal harvesting result was affected by environmental fluctuations.

In the presence of such a variety of environmental randomness, which can lead to crucial impact (see [18–24]). Lande et al. in [19] carried out that maybe largely because of the ignorance of randomness, the extinction of numerous species caused by over-harvesting. The next question is: if all species are affected by harvesting and environmental randomness, what role does randomness play? To answer the above question and inspired by the above literature, this article discusses the influence of environmental randomness.

Apart from randomness, time delay is another factor that is easily overlooked by scholars. Xu (see [25]) and Ma (see [26]) pointed out that systems with distributed delay are divided into two categories: discrete delay and continuous distributed delay. Whether distributed delay or discrete delay has a crucial impact on the result, because it is inevitable in nature world. Generally, when changes occur, it takes a certain amount of time for species in nature to show this effect. No species will react immediately in this situation (see [25, 26]).

In [27], the authors studied a tri-trophic stochastic food-chain model with harvesting. For each species, the threshold of persistence in mean and extinction, and the criterion for the stability in distribution of the system are obtained. Furthermore, the necessary and sufficient criterion for existence of the optimal harvesting strategy are established. The sustainable maximum yield and optimal harvesting effort are also given. In [28], taking harvesting and distributed delays into consideration, the authors investigated a class of stochastic three species food-chain models. The global dynamics of the model, including global asymptotic stability, extinction, random boundedness, and the probability distribution are obtained. Furthermore, the maximum of expectation of sustainable yield (MESY for short) and the optimal harvesting strategy are acquired.

However, we see that the similar research work for the general n-species (when $n\geq 4$) stochastic food-chain models is not found up to now. After a preliminary attempt, we find that there exists the larger difficulty to straightway investigate the dynamical behavior of the general n-species stochastic food-chain model at present. We cannot yet give a universal method or formula to establish the ideal results for the moment. For this reason, we focus on the following stochastic four species food-chain model with harvesting and distributed delay:

$$ \textstyle\begin{cases} \mathrm{d}x_{1}(t)= x_{1}(t) [r_{1}-h_{1}-a_{11}x_{1}(t)-a_{12} \int _{-\tau _{12}} ^{0}x_{2}(t+\theta )\,\mathrm{d}\mu _{12}(\theta ) ]\,\mathrm{d}t+ \sigma _{1}x_{1}(t) \,\mathrm{d}B_{1}(t), \\ \mathrm{d}x_{2}(t)= x_{2}(t) [r_{2}-h_{2}+a_{21} \int _{- \tau _{21}}^{0}x_{1}(t+\theta ) \,\mathrm{d}\mu _{21}(\theta )-a_{22}x_{2}(t) \\ \hphantom{\mathrm{d}x_{2}(t)={}}{} -a_{23} \int _{-\tau _{23}}^{0}x_{3}(t+\theta ) \,\mathrm{d} \mu _{23}(\theta ) ]\,\mathrm{d}t+\sigma _{2}x_{2}(t)\,\mathrm{d}B_{2}(t), \\ \mathrm{d}x_{3}(t)= x_{3}(t) [r_{3}-h_{3}+a_{32} \int _{- \tau _{32}}^{0}x_{2}(t+\theta ) \,\mathrm{d}\mu _{32}(\theta )-a_{33}x_{3}(t) \\ \hphantom{\mathrm{d}x_{3}(t)={}}{} -a_{34} \int _{-\tau _{34}}^{0}x_{4}(t+\theta ) \,\mathrm{d}\mu _{34}( \theta ) ]\,\mathrm{d}t+\sigma _{3}x_{3}(t)\,\mathrm{d}B_{3}(t), \\ \mathrm{d}x_{4}(t)= x_{4}(t) [r_{4}-h_{4}+a_{43} \int _{- \tau _{43}} ^{0}x_{3}(t+\theta )\,\mathrm{d}\mu _{43}(\theta )-a_{44}x_{4}(t) ] \,\mathrm{d}t+\sigma _{4}x_{4}(t)\,\mathrm{d}B_{4}(t). \end{cases} $$

(1)

We expect that the method and results introduced in this paper can help for us to investigate the general n-species stochastic food-chain models.

In model (1), the parameter $r_{1}>0$ is intrinsic growth rate of species $x_{1}$, $r_{i}\leq 0$ ($i=2,3,4$) represents death rates of species $x_{i}$, $a_{ii}>0$ ($i=1,2,3,4$) is density dependent coefficient of species $x_{i}$, $a_{12}\geq 0$, $a_{23}\geq 0$ and $a_{34}\geq 0$ are capture rates, $a_{21}\geq 0$, $a_{32}\geq 0$ and $a_{43}\geq 0$ stand for efficiency of food conversion, $h_{i}\geq 0$ ($i=1,2,3,4$) measures for the harvesting effort of species $x_{i}$, $\mu _{ij}(\theta )$ ($i,j=1,2,3,4$) is nonnegative variation function defined on $[-\tau _{ij},0]$ satisfying $\int _{-\tau _{ij}}^{0}\,\mathrm{d}\mu _{ij}(\theta )=1$, $B_{i}(t)$ ($i=1,2,3,4$) is the independent standard Brownian motion defined on the complete probability space $(\Omega ,\{\mathcal{F}_{t}\}_{t\geq 0},P)$ with a filtration $\{\mathcal{F}_{t}\}_{t\geq 0}$ satisfying the usual conditions, and $\sigma _{i}^{2}$ ($i=1,2,3,4$) is the intensity of $B_{i}(t)$.

We conduct research from two major aspects of model (1) in this article. One is the global dynamics, we mainly use the stochastic inequalities, the inequality estimation technique and Lyapunov function method to obtain. Another is about the harvesting, we consider the relation between the extinction, persistence of species and the influence of harvesting, we also obtain the optimal harvesting strategy $H^{*}=(h_{1}^{*},h_{2}^{*},h_{3}^{*},h_{4}^{*})$ and the maximal expectation of sustained yield $Y(H^{*})=\lim_{t\to \infty }\sum_{i=1}^{4}E(h^{*}_{i}x_{i}(t))$ under the premise that all species are not extinct.

The specific content of this paper is listed as follows. We first provide few necessary lemmas to prove the main results. In Sect. 2, for any positive initial value, the existence of the global unique positive solution is obtained, meanwhile, the random boundedness is also acquired. In Sect. 3, we not only establish the overall criterion of extinction and persistence in mean, but also establish the condition of the global asymptotic stability in distribution. We address the discussion that the impact of harvesting on extinction and persistence, and provide the sufficient and prerequisite criterion for the existence and non-existence of optimal harvesting strategy in Sect. 4. In order to clarify the main conclusions of the paper, we provide numerical simulations in Sect. 5. Finally, in Sect. 6, we not only give a concise conclusion, but also put forward some interesting relevant questions based on the thinking of this research.

2 Preliminaries

Firstly, introduce the following notations:

$$\begin{aligned}& b_{1}=r_{1}-h_{1}- \frac{1}{2}\sigma _{1}^{2},\qquad b_{2}=r_{2}-h_{2}- \frac{1}{2} \sigma _{2}^{2},\qquad b_{3}=r_{3}-h_{3}- \frac{1}{2}\sigma _{3}^{2}, \\& b_{4}=r_{4}-h_{4}- \frac{1}{2} \sigma _{4}^{2},\qquad \Delta _{11}=b_{1},\qquad \Delta _{21}=b_{1}a_{22}-b_{2}a_{12},\qquad \Delta _{22}=b_{1}a_{21}+b_{2}a_{11}, \\& \Delta _{31}={b_{1}(a_{22}a_{33}+a_{32}a_{23})-b_{2}a_{33}a_{12}+b_{3}a_{12}a_{23}}, \\& \Delta _{32}=a_{33}(b_{1}a_{21}+b_{2}a_{11})-b_{3}a_{11}a_{23}, \\& \Delta _{33}=(b_{1}a_{21}+b_{2}a_{11})a_{32}+b_{3}(a_{11}a_{22}+a_{12}a_{21}), \\& \Delta _{41}=b_{1}(a_{22}a_{34}a_{43}+a_{22}a_{33}a_{44}+a_{23}a_{32}a_{44})-b_{2}(a_{34}a_{43}a_{12}+a_{33}a_{44}a_{12}) \\& \hphantom{\Delta _{41}={}}{} +b_{3}a_{12}a_{23}a_{44}-b_{4}a_{12}a_{23}a_{44}, \\& \Delta _{42}=b_{1}(a_{21}a_{34}a_{43}+a_{21}a_{33}a_{44})+b_{2}(a_{34}a_{43}a_{11}+a_{33}a_{44}a_{11}) \\& \hphantom{\Delta _{42}={}}{}-b_{3}a_{11}a_{23}a_{44}+b_{4}a_{11}a_{23}a_{34}, \\& \Delta _{43}=b_{1}a_{21}a_{32}a_{44}+b_{2}a_{11}a_{32}a_{44}+b_{3}a_{44}(a_{11}a_{22}+a_{12}a_{21})-b_{4}a_{34}(a_{11}a_{22}+a_{12}a_{21}), \\& \Delta _{44}=b_{1}a_{21}a_{32}a_{43}+b_{2}a_{11}a_{32}a_{43}+b_{3}a_{43}(a_{11}a_{22}+a_{12}a_{21}) \\& \hphantom{\Delta _{44}={}}{}+b_{4} \bigl(a_{33}(a_{11}a_{22}+a_{12}a_{21})+a_{11}a_{23}a_{32} \bigr), \\& H_{1}=a_{11}, \qquad H_{2}=a_{11}a_{22}+a_{12}a_{21},\qquad H_{3}=a_{11}a_{22}a_{33}+a_{33}a_{12}a_{21}+a_{11}a_{32}a_{23}, \\& H_{4}=a_{11}a_{22}a_{33}a_{44}+a_{12}a_{21}a_{33}a_{44}+a_{11}a_{32}a_{23}a_{44} +a_{34}a_{43}a_{11}a_{22}+a_{34}a_{43}a_{12}a_{21}. \end{aligned}$$

It is clear that $b_{i}\leq 0$ for $i=2,3,4$ and when $b_{1}\geq 0$ we have $\Delta _{21}\geq 0$. Furthermore, we have the following lemma.

Lemma 1

If $\Delta _{44}>0$ (≥0), then $\Delta _{33}>0$, $\Delta _{41}>0$ (≥0), $\Delta _{42}>0$ (≥0) and $\Delta _{43}>0$ (≥0). If $\Delta _{33}>0$ (≥0), then $\Delta _{22}>0$, $\Delta _{31}>0$ (≥0) and $\Delta _{32}>0$ (≥0).

Proof

Let $\Delta _{44}>0$. Obviously, we have $\Delta _{33}>0$. Since $\Delta _{44}a_{44}-\Delta _{43}a_{43}=-H_{2}[-b_{4}(a_{34}a_{43} +a_{33}a_{44})+a_{11}a_{23}a_{32}]$, we obtain $\Delta _{44}a_{44}\leq \Delta _{43}a_{43}$, which implies $\Delta _{43}>0$.

By calculating we furthermore have

$$ \begin{aligned} &\Delta _{43}a_{33}+\Delta _{44}a_{34}-b_{3} \bigl[(a_{34}a_{43}+a_{33}a_{44}) (a_{11}a_{22}+a_{12}a_{21}) \\ &\quad {}+a_{11}a_{23}a_{32}a_{44} \bigr]-b_{4}a_{33}a_{34}(a_{11}a_{22}+a_{12}a_{21})= \Delta _{42}a_{32}. \end{aligned} $$

Hence, we obtain $\Delta _{42}>0$.

Furthermore, by calculating we also have

$$ \begin{aligned} &a_{22}\Delta _{42}+a_{23} \Delta _{43}-b_{2} \bigl[a_{11}a_{44}(a_{22}a_{33}+a_{23}a_{32}) \\ &\quad {}+a_{34}a_{43}(a_{11}a_{22}+a_{12}a_{21}) +a_{12}a_{21}a_{33}a_{44} \bigr]=a_{21} \Delta _{41}. \end{aligned} $$

Hence, we obtain $\Delta _{41}>0$.

Let $\omega _{1}^{*}=\frac{\Delta _{31}}{H_{3}}$, $\omega _{2}^{*}=\frac{\Delta _{32}}{H_{3}}$, $\omega _{3}^{*}=\frac{\Delta _{33}}{H_{3}}$. Then $\omega _{3}^{*}>0$. By calculating, we can obtain

$$ a_{32}\omega _{2}^{*}=-b_{3}+a_{33} \omega _{3}^{*}>0,\qquad a_{21}\omega _{1}^{*}=-b_{2}+a_{22} \omega _{2}^{*}+a_{23}\omega _{3}^{*}>0. $$

Therefore, we have $\Delta _{31}>0$ and $\Delta _{32}>0$. Obviously, we have $\Delta _{22}>0$. Similarly, we can prove the case of “≥0”. □

Lemma 2

For all real numbers $P\geq 0$, $Q\geq 0$, $P_{j}\geq 0$, and $a>0$, $b>0$ with $\frac{1}{a}+\frac{1}{b}=1$, where $1\leq j\leq n$, one has

$$ \Biggl(\sum_{j=1}^{n}P_{j} \Biggr)^{a}\leq n^{a}\sum_{j=1}^{n}P_{j}^{a},\qquad PQ \leq \frac{P^{a}}{a}+\frac{Q^{b}}{b}. $$

Lemma 3

Assume that positive constants $\alpha _{i}$, $i=1,2,3,4$ and an integer $n>0$ such that

$$ \alpha _{i} \biggl(-a_{ii}+ \frac{a_{ii-1}}{2n} \biggr)+\alpha _{i-1} \frac{a_{i-1i}}{2n^{2}}+ \alpha _{i+1}\frac{na_{i+1i}}{2}< 0,\quad i=1,2,3,4, $$

(2)

where we stipulate $\alpha _{0}=\alpha _{5}=0$ and $a_{10}=a_{01}=a_{45}=a_{54}=0$.

Proof

Define the matrix as follows:

$$ P= \begin{pmatrix} a_{11}& -\frac{n}{2}a_{21}& 0& 0 \\ -\frac{1}{2n^{2}}a_{12}& a_{22}- \frac{1}{2n}a_{21}& -\frac{n}{2}a_{32}&0 \\ 0& -\frac{1}{2n^{2}}a_{23}& a_{33}- \frac{1}{2n}a_{32}& -\frac{n}{2}a_{43} \\ 0&0& -\frac{1}{2n^{2}}a_{34}& a_{44}- \frac{1}{2n}a_{43} \end{pmatrix}. $$

Calculating the principal minors of P, we obtain

$$\begin{aligned}& P_{2}=\det \begin{pmatrix} a_{11}& -\frac{n}{2}a_{21} \\ -\frac{1}{2n^{2}}a_{12}& a_{22}- \frac{1}{2n}a_{21} \end{pmatrix} =a_{11} \biggl(a_{22}-\frac{1}{2n}a_{21} \biggr)- \frac{1}{4n}a_{12}a_{21}, \\& \begin{aligned} P_{3}&=\det \begin{pmatrix} a_{11}& -\frac{n}{2}a_{21}& 0 \\ -\frac{1}{2n^{2}}a_{12}& a_{22}- \frac{1}{2n}a_{21}& -\frac{n}{2}a_{32} \\ 0& -\frac{1}{2n^{2}}a_{23}& a_{33}- \frac{1}{2n}a_{32} \end{pmatrix} \\ &=a_{11} \biggl(a_{22}-\frac{1}{2n}a_{21} \biggr) \biggl(a_{33}-\frac{1}{2n}a_{32} \biggr)- \biggl(a_{33}- \frac{1}{2n}a_{32} \biggr) \frac{1}{4n}a_{12}a_{21}- \frac{1}{4n}a_{11}a_{23}a_{32}, \end{aligned} \end{aligned}$$

and

$$ P_{4}=\det P= \biggl(a_{44}-\frac{1}{2n}a_{43} \biggr)P_{3}-\frac{1}{4n}a_{34}a_{43}P_{2}. $$

From the expressions of $P_{i}$ ($i=2,3,4$), we easily see that there are enough large integers n such that $P_{i}>0$ ($i=2,3,4$). Therefore, P is a M-matrix. By the properties of M-matrix, there are the positive constant vector $\alpha =(\alpha _{1},\alpha _{2},\alpha _{3},\alpha _{4})^{T}$ such that $P\alpha >0$. Therefore, $(-P)\alpha <0$ which is equivalent to the inequality (2). This completes the proof. □

Let $r=\max \{\tau _{12},\tau _{21},\tau _{23},\tau _{32},\tau _{34}, \tau _{43}\}$. From the biological background, the initial data of any solution $x(t)=(x_{1}(t),x_{2}(t),x_{3}(t),x_{4}(t))$ for model (1) is defined as follows:

$$ x(\theta )= \bigl(\varsigma (\theta ),\xi (\theta ),\kappa ( \theta ),\eta ( \theta ) \bigr) ,\quad -r\leq \theta \leq 0. $$

(3)

For model (1), in regard to ultimate boundedness and existence of the positive global solution, we obtain the conclusions shown below.

Lemma 4

For all initial data $x(\theta )=(\varsigma (\theta ),\xi (\theta ),\kappa (\theta ),\eta ( \theta ))\in C([-\gamma ,0],R_{+}^{4})$, model (1) with condition (3) has a unique global solution $x(t)=(x_{1}(t),x_{2}(t),x_{3}(t),x_{4}(t))\in R^{4}_{+}$ a.s. for all $t\geq 0$. Moreover, for any $p>0$ there exist constants $K_{i}(p)>0$, $i=1,2,3,4$ such that

$$ \limsup_{t\to \infty }E \bigl[x_{i}^{p}(t) \bigr]\leq K_{i}(p),\quad i=1,2,3,4. $$

Proof

The proof of Lemma 4 is similar to Lemma 3 given in [28]. But, here we will give an improvement. Obviously, the model investigated here has local Lipschitz continuous coefficients. Then, for any $(\varsigma (\theta ),\xi (\theta ),\kappa (\theta ),\eta (\theta )) \in C([-r,0],R_{+}^{4})$, there is a unique solution $x(t)=(x_{1}(t),x_{2}(t),x_{3}(t),x_{4}(t))\in R^{4}_{+}$ on $t\in [-r,\tau _{e})$, here $\tau _{e}$ represents the explosion time. In order to obtain that the solution is global, $\tau _{e}=\infty$ a.s. should be proved. We first assume a large enough $k_{0}>0$ to let $\varsigma (0),\xi (0),\kappa (0),\eta (0)\in (\frac{1}{k_{0}},k_{0})$. Then, for any integer $k>k_{0}$, we define the following stopping time:

$$ \tau _{k}=\inf \biggl\{ t\in [0,\tau _{e}): \min_{1\leq i\leq 4} \bigl\{ x_{i}(t) \bigr\} \leq \frac{1}{k} \mbox{ or } \max_{1\leq i\leq 4} \bigl\{ x_{i}(t) \bigr\} \geq k \biggr\} . $$

(4)

$\tau _{k}$ is increasing as $k\rightarrow \infty$. Let $\tau _{\infty }=\lim_{k\to \infty }\tau _{k}$. $\tau _{\infty }\leq \tau _{e}$ a.s. is obtained. Therefore, we just have to clarify $\tau _{\infty }=\infty$ a.s.

Suppose that the assertion is wrong, then constants $\varepsilon \in (0,1)$ and $T>0$ such that $P(\tau _{\infty }\leq T)>\varepsilon$. Thus, an integer $k_{1}>k_{0}$ satisfying

$$ P(\tau _{k}\leq T)>\varepsilon $$

(5)

for any $k>k_{1}$. Define $V_{i}(x_{i})=x_{i}-1-\ln x_{i}$ ($i=1,2,3,4$). Using the Itô formula, we obtain

$$ \mathrm{d}V_{i}(x_{i})= \mathcal{L}\bigl[V_{i}(x_{i})\bigr]\,\mathrm{d}t+\sigma _{i}(x_{i}-1) \,\mathrm{d}B_{i}(t),\quad i=1,2,3,4, $$

(6)

here

$$\begin{aligned}& \mathcal{L} \bigl[V_{1}(x_{1}) \bigr]= (x_{1}-1) \biggl(r_{1}-h_{1}-a_{11}x_{1}(t)-a_{12} \int _{-\tau _{12}}^{0}x_{2}(t+\theta ) \, \mathrm{d}\mu _{12}(\theta ) \biggr) + \frac{1}{2}\sigma _{1}^{2}, \\& \begin{aligned} \mathcal{L} \bigl[V_{i}(x_{i}) \bigr]={}& (x_{i}-1) \biggl(r_{i}-h_{i}+a_{ii-1} \int _{- \tau _{ii-1}}^{0}x_{i-1}(t+\theta ) \, \mathrm{d}\mu _{i-1i}(\theta )-a_{ii}x_{i}(t) \\ &{} -a_{ii+1} \int _{-\tau _{ii+1}}^{0}x_{i+1}(t+\theta ) \, \mathrm{d}\mu _{ii+1}( \theta ) \biggr)+\frac{1}{2}\sigma _{i}^{2},\quad i=2,3, \end{aligned} \\& \mathcal{L} \bigl[V_{4}(x_{4}) \bigr]= (x_{4}-1) \biggl(r_{4}-h_{4}-a_{44}x_{4}(t)+a_{43} \int _{-\tau _{43}}^{0}x_{3}(t+\theta ) \, \mathrm{d}\mu _{43}(\theta ) \biggr) + \frac{1}{2}\sigma _{4}^{2}. \end{aligned}$$

By Lemma 2 and for an integer $n>0$, we get

$$ \begin{aligned} &\mathcal{L} \bigl[V_{1}(x_{1}) \bigr]\leq \frac{\sigma _{1}^{2}}{2}-(r_{1}-h_{1})+ \frac{n^{2}}{2}a_{12}+(r_{1}-h_{1})x_{1}+a_{11}x_{1}-a_{11}x_{1}^{2} \\ &\hphantom{\mathcal{L} \bigl[V_{1}(x_{1}) \bigr]\leq{}}{} +\frac{1}{2n^{2}}a_{12} \int _{-\tau _{12}}^{0}x_{2}^{2}(t+ \theta ) \,\mathrm{d}\mu _{12}(\theta ), \\ &\mathcal{L} \bigl[V_{i}(x_{i}) \bigr]\leq \frac{\sigma _{i}^{2}}{2}-(r_{i}-h_{i}) + \frac{n}{2}a_{ii-1} \int _{-\tau _{ii-1}}^{0}x_{i-1}^{2}(t+ \theta ) \,\mathrm{d}\mu _{ii-1}(\theta ) \\ &\hphantom{\mathcal{L} \bigl[V_{i}(x_{i}) \bigr]\leq{}}{} +(r_{i}-h_{i})x_{i}+a_{ii}x_{i}-a_{ii}x_{i}^{2}+ \frac{x_{i}^{2}}{2n}a_{ii-1}+ \frac{n^{2}}{2}a_{ii+1} \\ &\hphantom{\mathcal{L} \bigl[V_{i}(x_{i}) \bigr]\leq{}}{} +\frac{1}{2n^{2}}a_{ii+1} \int _{-\tau _{ii+1}}^{0}x_{i+1}^{2}(t+ \theta )\,\mathrm{d}\mu _{ii+1}(\theta ),\quad i=2,3, \\ &\mathcal{L} \bigl[V_{4}(x_{4}) \bigr]\leq \frac{\sigma _{4}^{2}}{2}-(r_{4}-h_{4})+ \frac{x_{4}^{2}}{2n}a_{43}+(r_{4}-h_{4})x_{4}+a_{44}x_{4} \\ &\hphantom{\mathcal{L} \bigl[V_{4}(x_{4}) \bigr]\leq{}}{} -a_{44}x_{4}^{2}+\frac{n}{2}a_{43} \int _{-\tau _{43}}^{0}x_{3}^{2}(t+ \theta )\,\mathrm{d}\mu _{43}(\theta ). \end{aligned} $$

(7)

Define $V_{0}(x)=\sum_{i=1}^{4}\alpha _{i}V_{i}(x_{i})+ V_{5}(t)$, here $x=(x_{1},x_{2},x_{3},x_{4})$ and

$$ \begin{aligned} V_{5}(t)={}&\sum _{i=1}^{3}\alpha _{i} \frac{1}{2n^{2}}a_{ii+1} \int _{- \tau _{ii+1}}^{0} \int _{t+\theta }^{t}x_{i+1}^{2}(s)\, \mathrm{d}s\, \mathrm{d} \mu _{ii+1}( \theta ) \\ &{}+\sum_{i=1}^{3}\alpha _{i+1}\frac{n}{2}a_{i+1i} \int _{-\tau _{i+1i}}^{0} \int _{t+\theta }^{t}x_{i}^{2}(s)\, \mathrm{d}s\, \mathrm{d} \mu _{i+1i}(\theta ). \end{aligned} $$

(8)

From the Itô formula

$$ \mathrm{d} \bigl[V_{0}(x) \bigr]= \mathcal{L}V_{0}(x)\,\mathrm{d}t+\sum_{i=1}^{4} \alpha _{i}\sigma _{i}(x_{i}-1) \, \mathrm{d}B_{i}(t). $$

From (7) and (8), we obtain

$$ \begin{aligned} \mathcal{L} \bigl[V_{0}(x) \bigr] ={}&\sum _{i=1}^{4}\alpha _{i} \mathcal{L}V_{i}(x_{i})+ \frac{\mathrm{d}}{\mathrm{d}t}V_{5}(t) \\ \leq {}&\sum_{i=1}^{4}\alpha _{i} \biggl\{ \frac{\sigma _{i}^{2}}{2}-(r_{i}-h_{i})+a_{ii+1} \frac{n^{2}}{2}+(r_{i}-h_{i})x_{i}+a_{ii}x_{i}-a_{ii}x_{i}^{2}+a_{ii-1} \frac{1}{2n}x_{i}^{2} \biggr\} \\ &{}+\sum_{i=1}^{4}\alpha _{i}a_{ii+1}\frac{1}{2n^{2}}x_{i+1}^{2}+ \sum_{i=1}^{4} \alpha _{i+1}a_{i+1i}\frac{n}{2}x_{i}^{2}, \end{aligned} $$

where we stipulate $\alpha _{5}=0$, $a_{10}=a_{01}=0$ and $a_{45}=a_{54}=0$. From (2) in Lemma 3, it is easy to find that there is a constant $K>0$ so that

$$ \mathrm{d} \bigl[V_{0}(x) \bigr] \leq K\, \mathrm{d}t+\sum_{i=1}^{4}\alpha _{i}\sigma _{i}(x_{i}-1) \, \mathrm{d}B_{i}(t). $$

(9)

Then, from (5) and (9), then we can get the following contradiction:

$$ \infty >V_{0} \bigl(x(0) \bigr)+KT\geq \infty . $$

Hence, we derive $\tau _{\infty }=\infty$ a.s., as a result, $\tau _{e}=\infty$ a.s.

For a constant $p>0$, assume $R_{1}(t)=e^{t}x_{1}^{p}(t)$. Using the Itô formula again,

$$ \mathrm{d}R_{1}(t)=\mathcal{L}R_{1}(t)\, \mathrm{d}t+pe^{t}x_{1}^{p} \sigma _{1}\,\mathrm{d}B_{1}(t), $$

(10)

where

$$ \begin{aligned} \mathcal{L}R_{1}(t)&=e^{t}x_{1}^{p} \biggl\{ 1+ \frac{p(p-1)\sigma _{1}^{2}}{2}+p \biggl[r_{1}-h_{1}-a_{11}x_{1} -a_{12} \int _{-\tau _{12}}^{0}x_{2}(t+\theta ) \, \mathrm{d}\mu _{12}(\theta ) \biggr] \biggr\} \\ &\leq \biggl\{ \biggl[p(r_{1}-h_{1})+1+ \frac{p(p-1)\sigma _{1}^{2}}{2} \biggr]x_{1}^{p}-pa_{11}x_{1}^{p+1} \biggr\} e^{t}. \end{aligned} $$

(11)

Assume that an integer $n>0$ and a constant $p>0$ satisfy $a_{22}-a_{21}\frac{p}{p+1}n^{-\frac{p+1}{p}}>0$, we define $R_{2}(t)$ as follows:

$$ R_{2}(t)=C_{1}^{*}R_{1}(t)+e^{t}x_{2}^{p}(t) +e^{\tau _{21}} \frac{pn^{p+1}}{p+1}a_{21} \int _{-\tau _{21}}^{0} \int _{t+\theta }^{t}e^{s}x_{1}^{p+1}(s) \,\mathrm{d}s\,\mathrm{d}\mu _{21}(\theta ), $$

(12)

where $C_{1}^{*}=a_{11}^{-1}e^{\tau _{21}}n^{p+1}a_{21}$. By the Itô formula, we get

$$ \mathrm{d}R_{2}(t)=\mathcal{L}R_{2}(t)\, \mathrm{d}t+C_{1}^{*}pe^{t}x_{1}^{p} \sigma _{1} \,\mathrm{d}B_{1}(t)+pe^{t}x_{2}^{p} \sigma _{2}\,\mathrm{d}B_{2}(t), $$

(13)

where

$$ \begin{aligned} \mathcal{L}R_{2}(t)={}&C^{*}_{1}e^{t}x_{1}^{p} \biggl\{ 1+ \frac{p(p-1)\sigma _{1}^{2}}{2}+p \biggl[r_{1}-h_{1}-a_{11}x_{1} -a_{12} \int _{-\tau _{12}}^{0}x_{2}(t+\theta ) \, \mathrm{d}\mu _{12}(\theta ) \biggr] \biggr\} \\ &{}+e^{t}x_{2}^{p} \biggl\{ 1+ \frac{p(p-1)\sigma _{2}^{2}}{2}+p \biggl[r_{2}-h_{2}+a_{21} \int _{-\tau _{21}}^{0}x_{1}(t+\theta ) \, \mathrm{d}\mu _{21}(\theta ) \\ &{}-a_{22}x_{2}(t)-a_{23} \int _{-\tau _{23}}^{0}x_{3}(t+\theta ) \, \mathrm{d}\mu _{23}(\theta ) \biggr] \biggr\} \\ &{}+e^{\tau _{21}}\frac{pn^{p+1}}{p+1}a_{21} \biggl(e^{t}x_{1}^{p+1}(t)- \int _{- \tau _{21}}^{0}e^{t+\theta }x_{1}^{p+1}(t+ \theta )\,\mathrm{d}\mu _{21}( \theta ) \biggr) \\ \leq {}&C_{1}^{*}e^{t} \biggl\{ \biggl[1+ \frac{p(p-1)\sigma _{1}^{2}}{2}+p(r_{1}-h_{1}) \biggr]x_{1}^{p}-pa_{11}x_{1}^{p+1} \biggr\} \\ &{}+e^{t} \biggl\{ \biggl[1+\frac{p(p-1)\sigma _{2}^{2}}{2}+p(r_{2}-h_{2}) \biggr]x_{2}^{p}-p \biggl[a_{22}-a_{21} \frac{p}{p+1}n^{-\frac{{p+1}}{p}} \biggr]x_{2}^{p+1} \\ &{}+\frac{p}{p+1}n^{p+1}a_{21} \int _{-\tau _{21}}^{0}x_{1}^{p+1}(t+ \theta )\,\mathrm{d}\mu _{21}(\theta ) \biggr\} \\ &{}+e^{\tau _{21}}\frac{pn^{p+1}}{p+1}a_{21} \biggl(e^{t}x_{1}^{p+1}(t) -e^{- \tau _{21}} \int _{-\tau _{21}}^{0}e^{t}x_{1}^{p+1}(t+ \theta )\,\mathrm{d} \mu _{21}(\theta ) \biggr) \\ \leq {}&e^{t} \biggl\{ \biggl[1+\frac{p(p-1)\sigma _{2}^{2}}{2}+p(r_{2}-h_{2}) \biggr]x_{2}^{p}-p \biggl[a_{22}-a_{21} \frac{p}{p+1}n^{-\frac{{p+1}}{p}} \biggr]x_{2}^{p+1} \\ &{}+C_{1}^{*} \biggl[1+\frac{p(p-1)\sigma _{1}^{2}}{2}+p(r_{1}-h_{1}) \biggr]x_{1}^{p}-e^{ \tau _{21}} \frac{p^{2}}{p+1}n^{p+1}a_{21}x_{1}^{p+1} \biggr\} . \end{aligned} $$

(14)

Assume that an integer $n>0$ and a constant $p>0$ satisfy $a_{33}-a_{32}\frac{p}{p+1}n^{-\frac{p+1}{p}}>0$, we define $R_{3}(t)$ as follows:

$$ R_{3}(t)=C_{2}^{*}R_{2}(t)+e^{t}x_{3}^{p} +e^{\tau _{32}} \frac{pn^{p+1}}{p+1}a_{32} \int _{-\tau _{32}}^{0} \int _{t+\theta }^{t}e^{s}x_{2}^{p+1}(s) \,\mathrm{d}s\,\mathrm{d}\mu _{32}(\theta ), $$

(15)

where $C_{2}^{*}=a_{22}^{-1}e^{\tau _{32}}n^{p+1}a_{32}$. By the Itô formula, we obtain

$$ \mathrm{d}R_{3}(t)=\mathcal{L}R_{3}(t)\, \mathrm{d}t+C_{2}^{*} \bigl(C_{1}^{*}pe^{t}x_{1}^{p} \sigma _{1}\,\mathrm{d}B_{1}(t)+pe^{t}x_{2}^{p} \sigma _{2}\,\mathrm{d}B_{2}(t) \bigr)+pe^{t}x_{3}^{p} \sigma _{3}\,\mathrm{d}B_{3}(t), $$

(16)

where similarly to (14) we can obtain

$$ \begin{aligned} \mathcal{L}R_{3}(t)\leq{} & e^{t} \biggl\{ \biggl[1+p(r_{3}-h_{3})+ \frac{p(p-1)\sigma _{3}^{2}}{2} \biggr]x_{3}^{p}-p \biggl[a_{33} -a_{32} \frac{p}{p+1}n^{-\frac{{p+1}}{p}} \biggr]x_{3}^{p+1} \\ &{}+C_{2}^{*} \biggl[1+p(r_{2}-h_{2})+ \frac{p(p-1)\sigma _{2}^{2}}{2} \biggr]x_{2}^{p}- \frac{p^{2}}{p+1} \bigl(n^{p+1}a_{32}e^{\tau _{32}}+n^{-\frac{p+1}{p}}a_{21}C_{2}^{*} \bigr)x_{2}^{p+1} \\ &{}+C_{1}^{*}C_{2}^{*} \biggl[1+p(r_{1}-h_{1})+\frac{p(p-1)\sigma _{1}^{2}}{2} \biggr]x_{1}^{p}-C_{2}^{*}e^{ \tau _{21}} \frac{p^{2}}{p+1}n^{p+1}a_{21}x_{1}^{p+1} \biggr\} . \end{aligned} $$

Finally, assume that an integer $n>0$ and a constant $p>0$ satisfy $a_{44}-a_{43}\frac{p}{p+1}n^{-\frac{p+1}{p}}>0$, we define $R_{4}(t)$ as follows:

$$ R_{4}(t)=C_{3}^{*}R_{3}(t)+e^{t}x_{4}^{p} +e^{\tau _{43}} \frac{pn^{p+1}}{p+1}a_{43} \int _{-\tau _{43}}^{0} \int _{t+\theta }^{t}e^{s}x_{3}^{p+1}(s) \,\mathrm{d}s\,\mathrm{d}\mu _{43}(\theta ), $$

(17)

where $C_{3}^{*}=a_{33}^{-1}e^{\tau _{43}}n^{p+1}a_{43}$. From the Itô formula, we derive

$$ \begin{aligned} \mathrm{d}R_{4}(t)={}& \mathcal{L}R_{4}(t)\, \mathrm{d}t+C_{3}^{*} \bigl(C_{2}^{*} \bigl(C_{1}^{*}pe^{t}x_{1}^{p} \sigma _{1}\,\mathrm{d}B_{1}(t)+pe^{t}x_{2}^{p} \sigma _{2}\,\mathrm{d}B_{2}(t) \bigr) \\ &{}+pe^{t}x_{3}^{p}\sigma _{3} \,\mathrm{d}B_{3}(t) \bigr)+pe^{t}x_{4}^{p} \sigma _{4} \,\mathrm{d}B_{4}(t), \end{aligned} $$

(18)

where similarly to (14) we can obtain

$$ \begin{aligned} \mathcal{L}R_{4}(t)\leq {}& e^{t} \biggl\{ \biggl[1+p(r_{4}-h_{4})+ \frac{p(p-1)\sigma _{4}^{2}}{2} \biggr]x_{4}^{p}-p \biggl[a_{44}-a_{43} \frac{p}{p+1}n^{-\frac{{p+1}}{p}} \biggr]x_{4}^{p+1} \\ &{}+C_{3}^{*} \biggl[1+p(r_{3}-h_{3})+ \frac{p(p-1)\sigma _{3}^{2}}{2} \biggr]x_{3}^{p}- \frac{p^{2}}{p+1} \bigl(n^{p+1}a_{43}e^{\tau _{43}}+n^{-\frac{p+1}{p}}a_{32}C_{3}^{*} \bigr)x_{3}^{p+1} \\ &{}+C_{2}^{*}C_{3}^{*} \biggl[1+p(r_{2}-h_{2})+\frac{p(p-1)\sigma _{2}^{2}}{2} \biggr]x_{2}^{p} \\ &{}-C_{3}^{*} \frac{p^{2}}{p+1} \bigl(n^{p+1}a_{32}e^{\tau _{32}} +n^{-\frac{p+1}{p}}a_{21}C_{2}^{*} \bigr)x_{2}^{p+1} \\ &{}+C_{1}^{*}C_{2}^{*}C_{3}^{*} \biggl[1+p(r_{1}-h_{1})+ \frac{p(p-1)\sigma _{1}^{2}}{2} \biggr]x_{1}^{p}-C_{2}^{*}C_{3}^{*}e^{\tau _{21}} \frac{p^{2}}{p+1}n^{p+1}a_{21}x_{1}^{p+1} \biggr\} . \end{aligned} $$

For any $t\geq 0$, we derive here exists a constant $K_{4}(p)>0$ so that $\mathcal{L}R_{4}(t)\leq K_{4}(p)e^{t}$ by the above inequality, Hence, $E[R_{4}(t)]\leq E[R_{4}(0)]+K_{4}(p)(e^{t}-1)$ for all $t\geq 0$ is obtained. Consequently, from the definitions of $Q_{i}(t)$ ($i=1,2,3,4$) we furthermore have

$$ \begin{aligned} & C_{3}^{*}E \bigl[R_{3}(t) \bigr]\leq E \bigl[R_{4}(0) \bigr]+K_{4}(p) \bigl(e^{t}-1 \bigr), \\ & C_{2}^{*}C_{3}^{*}E \bigl[R_{2}(t) \bigr]\leq E \bigl[R_{4}(0) \bigr]+K_{4}(p) \bigl(e^{t}-1 \bigr), \\ & C_{1}^{*}C_{2}^{*}C_{3}^{*}E \bigl[R_{1}(t) \bigr]\leq E \bigl[R_{4}(0) \bigr]+K_{4}(p) \bigl(e^{t}-1 \bigr), \end{aligned} $$

for any $t\geq 0$. Since $E[e^{t}x_{i}^{p}(t)]\leq E[R_{i}(t)]$ ($i=1,2,3,4$) is also acquired for all $t\geq 0$, this shows that there are constants $K_{i}(p)>0$, $i=1,2,3,4$, satisfying $\limsup_{t\to \infty }E[x_{i}^{p}(t)]\leq K_{i}(p)$ ($i=1,2,3,4$). □

Remark 1

Observing the proof process from Lemma 4, we easily find that Lemma 4 seemingly can be extended to the general n-species stochastic food-chain system with distributed delay and harvesting.

Lemma 5

Suppose that the functions $P\in C(R_{+}\times \Omega , R_{+})$ and $Q\in C(R_{+}\times \Omega , R)$ satisfy $\lim_{t\rightarrow \infty }\frac{Q(t)}{t}=0$ a.s.

• (1) Assume that there exist a few constants $\beta >0$, $T>0$ and $\beta _{0}>0$ such that for $t\geq T$

  $$ \ln P(t)=\beta t-\beta _{0} \int _{0}^{t}P(s)\,\mathrm{d}s+Q(t)\quad \textit{a.s.}, $$

  then $\lim_{t\to \infty }\langle P(t)\rangle =\frac{\beta }{\beta _{0}}$ a.s., and $\lim_{t\to \infty }\frac{\ln P(t)}{t}=0$ a.s.

• (2) Assume that there are constants $T>0$, $\beta _{0}>0$ and $\beta \in R$ such that for $t\geq T$

  $$ \ln P(t)\leq \beta t-\beta _{0} \int _{0}^{t}P(s)\,\mathrm{d}s+Q(t)\quad \textit{a.s.}, $$

  then $\limsup_{t\to \infty }\langle P(t)\rangle \leq \frac{\beta }{\beta _{0}}$ a.s. as $\beta \geq 0$, and $\lim_{t\to \infty }P(t)=0$ a.s. as $\beta < 0$.

• (3) Assume that there exist constants $\beta >0$, $\beta _{0}>0$ and $T>0$ such that for $t\geq T$

  $$ \ln P(t)\geq \beta t-\beta _{0} \int _{0}^{t}P(s)\,\mathrm{d}s+Q(t)\quad \textit{a.s.}, $$

  then $\liminf_{t\to \infty }\langle P(t)\rangle \geq \frac{\beta }{\beta _{0}}$ a.s.

We consider an auxiliary system as follows:

$$ \textstyle\begin{cases} \mathrm{d}y_{1}(t)= y_{1}(t) [r_{1}-h_{1}-a_{11}y_{1}(t) ]\,\mathrm{d}t+\sigma _{1}y_{1}(t) \,\mathrm{d}B_{1}(t), \\ \mathrm{d}y_{i}(t)= y_{i}(t) [r_{i}-h_{i}+a_{ii-1} \int _{-\tau _{ii-1}}^{0}y_{i-1}(t+ \theta ) \,\mathrm{d}\mu _{ii-1}(\theta ) -a_{ii}y_{i}(t) ]\,\mathrm{d}t \\ \hphantom{\mathrm{d}y_{i}(t)={}}{} +\sigma _{i}y_{i}(t)\,\mathrm{d}B_{i}(t),\quad i=2,3,4, \end{cases} $$

(19)

and the initial value is given by

$$ \bigl(y_{1}(\theta ),y_{2}(\theta ),y_{3}(\theta ),y_{4}(\theta ) \bigr)= \bigl( \varsigma ( \theta ),\xi (\theta ),\kappa (\theta ),\eta (\theta ) \bigr),\quad -r \leq \theta \leq 0. $$

(20)

Here, we use the same argument as in the proof of Lemma 3, with the condition (20) we can easily derive model (19) has a unique global solution $(y_{1}(t),y_{2}(t),y_{3}(t),y_{4}(t))\in R^{4}_{+}$ a.s. for all $t\geq 0$. The following conclusions are derived.

Here, for convenience, we denote $\Lambda _{11}=\Delta _{11}$, $\Lambda _{22}=\Delta _{22}$, $\Lambda _{33}=\Delta _{33}-b_{3}a_{12}a_{21}$ and $\Lambda _{44}=\Delta _{44}-b_{4}(a_{33}a_{12}a_{21}+a_{11}a_{23}a_{32})$.

Lemma 6

Assume that $(y_{1}(t),y_{2}(t),y_{3}(t),y_{4}(t))$ is any positive global solution of model (19). We derive:

• (1) Suppose that $\Lambda _{11}<0$, then $\lim_{t\to \infty }y_{i}(t)=0$ a.s., $i=1,2,3,4$.

• (2) Suppose that $\Lambda _{11}=0$, then $\lim_{t\to \infty }\langle Z_{1}(t)\rangle =0$, and $\lim_{t\to \infty }y_{i}(t)=0$ a.s., $i=2,3,4$.

• (3) Suppose that $\Lambda _{11}>0$ and $\Lambda _{22}<0$, then $\lim_{t\to \infty }\langle y_{1}(t)\rangle = \frac{\Lambda _{11}}{a_{11}}$, and $\lim_{t\to \infty }y_{i}(t)=0$ a.s., $i=2,3,4$.

• (4) Suppose that $\Lambda _{22}=0$, then $\lim_{t\to \infty }\langle y_{1}(t)\rangle = \frac{\Lambda _{11}}{a_{11}}$, and $\lim_{t\to \infty }\langle y_{2}(t)\rangle =0$, $\lim_{t\to \infty }y_{i}(t)=0$ a.s., $i=3,4$.

• (5) Suppose that $\Lambda _{22}>0$ and $\Lambda _{33}<0$, then $\lim_{t\to \infty }\langle y_{i}(t)\rangle = \frac{\Lambda _{ii}}{\prod_{j=1}^{i}a_{jj}}$ a.s., $i=1,2$ and $\lim_{t\to \infty }y_{i}(t)=0$ a.s., $i=3,4$.

• (6) Suppose that $\Lambda _{33}=0$, then $\lim_{t\to \infty }\langle y_{i}(t)\rangle = \frac{\Lambda _{ii}}{\prod_{j=1}^{i}a_{jj}}$ a.s., $i=1,2$, $\lim_{t\to \infty }\langle y_{3}(t)\rangle =0$ a.s. and $\lim_{t\to \infty }y_{4}(t)=0$ a.s.

• (7) Suppose that $\Lambda _{33}>0$ and $\Lambda _{44}<0$, then $\lim_{t\to \infty }\langle y_{i}(t)\rangle = \frac{\Lambda _{ii}}{\prod_{j=1}^{i}a_{jj}}$ a.s., $i=1, 2,3$, $\lim_{t \to \infty }y_{4}(t)=0$ a.s.

• (8) Suppose that $\Lambda _{44}=0$, then $\lim_{t\to \infty }\langle y_{i}(t)\rangle = \frac{\Lambda _{ii}}{\prod_{j=1}^{i}a_{jj}}$ a.s., $i=1, 2,3$, $\lim_{t \to \infty }\langle y_{4}(t)\rangle =0$ a.s.

• (9) Suppose that $\Lambda _{44}>0$, then $\lim_{t\to \infty }\langle y_{i}(t)\rangle = \frac{\Lambda _{ii}}{\prod_{j=1}^{i}a_{jj}}$ a.s., $i=1,2,3,4$.

• (10) $\limsup_{t\to \infty }\frac{\ln y_{i}(t)}{t}\leq 0$ a.s., for $i=1,2,3,4$.

The proving process of Lemma 6 is similar to Lemma 5 given in [28]. We hence omit it here. It is clear that Lemma 6 also seemingly can be extended to the general n-species stochastic food-chain system with distributed delay and harvesting.

Lemma 7

Assume that $(x_{1}(t),x_{2}(t),x_{3}(t),x_{4}(t))$ and $(y_{1}(t),y_{2}(t),y_{3}(t),y_{4}(t))$ are the solutions of model (1) and model (19), respectively. Then, for any $-r\leq \theta \leq 0$ and $i=1,2,3,4$, we obtain:

• (1) If the initial conditions such that $x_{i}(\theta )\leq y_{i}(\theta )$, then $x_{i}(t)\leq y_{i}(t)$ for $t\geq 0$,

• (2) $\limsup_{t\to \infty }\frac{\ln x_{i}(t)}{t}\leq 0$ a.s.,

• (3) $\lim_{t\to \infty }\frac{1}{t}\int _{t-\tau }^{t}x_{i}(s)\,\mathrm{d}s=0$ a.s. when the constant $\tau >0$.

Proof

From model (1) we get

$$\begin{aligned}& \mathrm{d}x_{1}(t)\leq x_{1}(t) \bigl[r_{1}-h_{1}-a_{11}x_{1}(t) \bigr] \,\mathrm{d}t+\sigma _{1}x_{1}(t) \, \mathrm{d}B_{1}(t), \\& \begin{aligned} \mathrm{d}x_{i}(t)\leq {}& x_{i}(t) \biggl[r_{i}-h_{i}+a_{ii-1} \int _{- \tau _{ii-1}}^{0}x_{i-1}(t+\theta ) \, \mathrm{d}\mu _{ii-1}(\theta )-a_{ii}x_{i}(t) \biggr] \,\mathrm{d}t \\ &{}+\sigma _{i}x_{i}(t)\,\mathrm{d}B_{i}(t),\quad i=2,3,4. \end{aligned} \end{aligned}$$

From the comparison theorem, we obtain $x_{i}(t)\leq y_{i}(t)$ ($i=1,2,3,4$) on $t\geq 0$. Thus, for a constant $\tau >0$, we find that $\limsup_{t\to \infty }\frac{\ln x_{i}(t)}{t}\leq 0$ a.s. and $\lim_{t\to \infty }\frac{1}{t}\int _{t-\tau }^{t}x_{i}(s)\,\mathrm{d}s=0$ a.s. ($i=1,2,3,4$) hold from Lemma 6. □

Remark 2

It is clear that Lemma 7 also is satisfied for the general n-species stochastic food-chain system with distributed delay and harvesting.

3 Global dynamics

Here, we firstly introduce the following useful lemma.

Lemma 8

Assume that model (1) has the solution $(x_{1}(t),x_{2}(t),x_{3}(t),x_{4}(t))$. If there is an $i\in \{1,2,3\}$ to satisfy $\lim_{t\to \infty }\langle x_{i}(t)\rangle =0$ a.s., then, for all $j>i$, $\lim_{t\to \infty }x_{j}(t)=0$ a.s. holds.

Proof

We first use the Itô formula, then

$$\begin{aligned}& \ln x_{1}(t) = b_{1}t-a_{11} \int _{0}^{t}x_{1}(s) \, \mathrm{d}s-a_{12} \int _{0}^{t}x_{2}(s)\,\mathrm{d}s + \phi _{1}(t), \end{aligned}$$

(21)

$$\begin{aligned}& \begin{aligned} \ln x_{i}(t) ={}&b_{i}t+a_{ii-1} \int _{0}^{t}x_{i-1}(s)\,\mathrm{d}s -a_{ii} \int _{0}^{t}x_{i}(s)\,\mathrm{d}s \\ &{}-a_{ii+1} \int _{0}^{t}x_{i+1}(s)\,\mathrm{d}s+ \phi _{i}(t),\quad i=2,3, \end{aligned} \end{aligned}$$

(22)

and

$$ \ln x_{4}(t) =b_{4}t+a_{43} \int _{0}^{t}x_{3}(s) \, \mathrm{d}s-a_{44} \int _{0}^{t}x_{4}(s) \,\mathrm{d}s+ \phi _{4}(t), $$

(23)

where

$$\begin{aligned}& \begin{aligned} \phi _{1}(t)={}&\sigma _{1}B_{1}(t)+\ln x_{1}(0)+a_{12} \int _{-\tau _{12}}^{0} \int _{t+\theta }^{t}x_{2}(s)\,\mathrm{d}s \,\mathrm{d}\mu _{12}(\theta ) \\ &{}-a_{12} \int _{-\tau _{12}}^{0} \int _{\theta }^{0}x_{2}(s)\,\mathrm{d}s \,\mathrm{d}\mu _{12}( \theta ), \end{aligned} \\& \begin{aligned} \phi _{i}(t)={}&\sigma _{i}B_{i}(t)+\ln x_{i}(0) +a_{ii-1} \int _{-\tau _{ii-1}}^{0} \int _{\theta }^{0}x_{i-1}(s)\,\mathrm{d}s \,\mathrm{d}\mu _{ii-1}(\theta ) \\ &{}-a_{ii-1} \int _{-\tau _{ii-1}}^{0} \int _{t+\theta }^{t}x_{i-1}(s) \, \mathrm{d}s \,\mathrm{d}\mu _{ii-1}(\theta ) +a_{ii+1} \int _{-\tau _{ii+1}}^{0} \int _{t+\theta }^{t}x_{i+1}(s)\,\mathrm{d}s \,\mathrm{d}\mu _{ii+1}(\theta ) \\ &{}-a_{ii+1} \int _{-\tau _{ii+1}}^{0} \int _{\theta }^{0}x_{i+1}(s) \,\mathrm{d}s \,\mathrm{d}\mu _{ii+1}(\theta ),\quad i=2,3, \end{aligned} \\& \begin{aligned} \phi _{4}(t)={}&\sigma _{4}B_{4}(t)+\ln x_{4}(0)+a_{43} \int _{-\tau _{43}}^{0} \int _{\theta }^{0}x_{3}(s)\,\mathrm{d}s \,\mathrm{d}\mu _{43}(\theta ) \\ &{}-a_{43} \int _{-\tau _{43}}^{0} \int _{t+\theta }^{t}x_{3}(s)\,\mathrm{d}s \,\mathrm{d} \mu _{43}(\theta ). \end{aligned} \end{aligned}$$

Obviously, $\lim_{t\to \infty }\frac{\phi _{i}(t)}{t}=0$ a.s. is obtained for $i=1,2,3,4$ by Lemma 7. Assume $\lim_{t\to \infty }\langle x_{i}(t)\rangle =0$ a.s. Then for any constant $\varepsilon >0$ with $b_{i+1}+a_{i+1i}\varepsilon <0$ there exists a $T>0$ to satisfy $\int _{0}^{t}x_{i}(s)\,\mathrm{d}s<\varepsilon t$ for any $t\geq T$. Therefore, for $t\geq T$, by (22) and (23), the following inequality is found:

$$ \ln x_{i+1}(t)\leq b_{i+1}t+a_{i+1i}\varepsilon t -a_{i+1i+1} \int _{0}^{t}x_{i+1}(s) \, \mathrm{d}s+ \phi _{i+1}(t). $$

Thus, by Lemma 5 we derive $\lim_{t\to \infty }x_{i+1}(t)=0$ a.s. Consequently, $\lim_{t\to \infty }x_{j}(t)=0$ a.s. for any $j>i$. □

Remark 3

It is easy for us to find that Lemma 8 also seemingly can be extended to the general n-species stochastic food-chain system with distributed delay and harvesting.

In the following theorem, we state and prove a screening criterion as a main result in this paper on the extinction and persistence in mean of global positive solutions for model (1).

Theorem 1

Suppose that $(x_{1}(t),x_{2}(t),x_{3}(t),x_{4}(t))$ is any positive global solution of model (1). Then we derive:

• (1) If $\Delta _{11}<0$, then $\lim_{t\to \infty }x_{j}(t)=0$ a.s. for $j=1,2,3,4$.

• (2) If $\Delta _{11}=0$, then $\lim_{t\to \infty }\langle x_{1}(t)\rangle =0$ and $\lim_{t\to \infty }x_{j}(t)=0$ a.s. for $j=2,3,4$.

• (3) If $\Delta _{11}>0$ and $\Delta _{22}<0$, then $\lim_{t\to \infty }\langle x_{1}(t)\rangle = \frac{\Delta _{11}}{H_{1}}$ and $\lim_{t\to \infty }x_{j}(t)=0$ a.s. for $j=2,3,4$.

• (4) If $\Delta _{22}=0$, then $\lim_{t\to \infty }\langle x_{1}(t)\rangle = \frac{\Delta _{11}}{H_{1}}$, $\lim_{t\to \infty }\langle x_{2}(t)\rangle =0$ and $\lim_{t\to \infty }x_{j}(t)=0$ a.s. for $j=3,4$.

• (5) If $\Delta _{22}>0$ and $\Delta _{33}<0$, then $\lim_{t\to \infty }\langle x_{j}(t)\rangle = \frac{\Delta _{2j}}{H_{2}}$, $j=1,2$, and $\lim_{t\to \infty }x_{j}(t)=0$ a.s. for $j=3,4$.

• (6) If $\Delta _{33}=0$ and the condition

  $$ a_{33}a_{22}H_{2}-a_{12}a_{21}a_{23}a_{32}>0 $$

  (24)

  holds, then $\lim_{t\to \infty }\langle x_{j}(t)\rangle = \frac{\Delta _{2j}}{H_{2}}$, $j=1,2$, $\lim_{t\to \infty }\langle x_{3}(t)\rangle =0$ and $\lim_{t\to \infty }x_{4}(t)=0$ a.s.

• (7) If $\Delta _{33}>0$, $\Delta _{44}<0$ and the condition (24) holds, then $\lim_{t\to \infty }\langle x_{j}(t)\rangle = \frac{\Delta _{3j}}{H_{3}}$, $j=1, 2,3$, $\lim_{t\to \infty }x_{4}(t)=0$ a.s.

• (8) If $\Delta _{44}=0$ and the condition

  $$ (a_{22}a_{33}H_{2}-a_{12}a_{21}a_{23}a_{32})a_{44}H_{3} -a_{23}a_{32}a_{34}a_{43}H_{2}^{2}>0 $$

  (25)

  holds, then $\lim_{t\to \infty }\langle x_{j}(t)\rangle = \frac{\Delta _{3j}}{H_{3}}$, $j=1,2,3$ and $\lim_{t\to \infty }\langle x_{4}(t)\rangle =0$ a.s.

• (9) If $\Delta _{44}>0$ and the condition (25) holds, then $\lim_{t\to \infty }\langle x_{j}(t)\rangle = \frac{\Delta _{4j}}{H_{4}}$, $j=1,2,3,4$, a.s.

Proof

For model (1), $(x_{1}(t),x_{2}(t),x_{3}(t),x_{4}(t))$ can be the positive global solution. Let $V_{2}(t)=a_{21}\ln x_{1}(t)+a_{11}\ln x_{2}(t)$, $V_{3}(t)=a_{32}V_{2}(t)+H_{2}\ln x_{3}(t)$ and $V_{4}(t)=a_{43}V_{3}(t)+H_{3}\ln x_{4}(t)$. From (21)–(23), we obtain

$$ V_{4}(t)=\Delta _{44}t-H_{4} \int _{0}^{t}x_{4}(s)\,\mathrm{d}s+ \phi _{5}(t), $$

(26)

where $\phi _{5}(t)=a_{21}a_{32}a_{43}\phi _{1}(t)+a_{43}a_{11}a_{32}\phi _{2}(t)+a_{43}H_{2} \phi _{3}(t)+H_{3}\phi _{4}(t)$. we apply the similar method that used for $\phi _{1}(t)$, $\lim_{t\to \infty }\frac{\phi _{5}(t)}{t}=0$ a.s. is obtained.

If $\Delta _{44}>0$, then, by Lemma 7, and for any $\varepsilon >0$ with $\Delta _{44}-3\varepsilon >0$, there exists a constant $T>0$ satisfies $\ln x_{1}(t)<\frac{\varepsilon }{a_{43}a_{32}a_{21}+1}t$, $\ln x_{2}(t)<\frac{\varepsilon }{a_{43}a_{32}a_{11}+1}t$ and $\ln x_{3}(t)<\frac{\varepsilon }{a_{43}H_{2}+1}t$ for all $t\geq T$. Then, from (26) we obtain

$$ H_{3}\ln x_{4}(t)>(\Delta _{44}-3 \varepsilon )t-H_{4} \int _{0}^{t}x_{4}(s) \,\mathrm{d}s+ \phi _{5}(t) $$

for all $t\geq T$. Thus, from the arbitrary ε and Lemma 4

$$ \liminf_{t\to \infty } \bigl\langle x_{4}(t) \bigr\rangle \geq \frac{\Delta _{44}}{H_{4}} $$

(27)

is obtained.

If $\Delta _{44}\leq 0$, then since $\liminf_{t\to \infty }\langle x_{4}(t)\rangle \geq 0$, we also have $\liminf_{t\to \infty }\langle x_{4}(t)\rangle \geq \frac{\Delta _{44}}{H_{4}}$. Let $U_{2}(t)=a_{22}\ln x_{1}(t)-a_{12}\ln x_{2}(t)$ and $U_{4}(t)=a_{43}U_{2}(t)-a_{12}a_{23}\ln x_{4}(t)$. By (21), (22) and (23), we compute

$$ U_{4}(t)=(a_{43} \Delta _{21}-b_{4}a_{12}a_{23})t+a_{12}a_{23}a_{44} \int _{0}^{t}x_{4}(s) \, \mathrm{d}s-H_{2}a_{43} \int _{0}^{t}x_{1}(s) \, \mathrm{d}s+ \phi _{6}(t), $$

(28)

where $\phi _{6}(t)=a_{22}a_{43}\phi _{1}(t)-a_{12}a_{43}\phi _{2}(t)-a_{12}a_{23} \phi _{4}(t)$. we apply the similar method that used for $\phi _{1}(t)$, $\lim_{t\to \infty }\frac{\phi _{6}(t)}{t}=0$ a.s. is derived. For all $\varepsilon >0$, there exists a constant $T>0$ such that $\ln x_{2}(t)<\frac{\varepsilon }{a_{43}a_{12}+1}t$, $\ln x_{3}(t)<\frac{\varepsilon }{a_{12}a_{23}+1}t$ and

$$ \int _{0}^{t}x_{4}(s)\,\mathrm{d}s \leq \Bigl( \limsup_{t\to \infty } \bigl\langle x_{4}(t) \bigr\rangle + \varepsilon \Bigr)t $$

for any $t\geq T$.

Thus, from (28), we obtain

$$ \begin{aligned} a_{43}a_{22} \ln x_{1}(t)\leq {}& \Bigl(a_{43}\Delta _{21}-b_{4}a_{12}a_{23}+2 \varepsilon +a_{12}a_{23}a_{44} \Bigl(\limsup _{t\to \infty } \bigl\langle x_{4}(t) \bigr\rangle + \varepsilon \Bigr) \Bigr)t \\ &{}-H_{2}a_{43} \int _{0}^{t}x_{1}(s)\,\mathrm{d}s+ \phi _{6}(t) \end{aligned} $$

(29)

for any $t\geq T$. Thus, from the arbitrary ε and Lemma 4

$$ \limsup_{t\to \infty } \bigl\langle x_{1}(t) \bigr\rangle \leq \frac{(a_{43}\Delta _{21}-b_{4}a_{12}a_{23} +a_{12}a_{23}a_{44}\limsup_{t\to \infty }\langle x_{4}(t)\rangle )}{H_{2}a_{43}} $$

(30)

is furthermore obtained.

From (21) and (22), we have

$$ V_{3}(t)=\Delta _{33}t-H_{3} \int _{0}^{t}x_{3}(s)\,\mathrm{d}s+ \phi _{7}(t)-a_{34}H_{2} \int _{0}^{t}x_{4}(s)\,\mathrm{d}s, $$

(31)

where $\phi _{7}(t)=a_{21}a_{32}\phi _{1}(t)+a_{11}a_{32}\phi _{2}(t)+H_{2} \phi _{3}(t)$. we apply the similar method that used for $\phi _{1}(t)$, $\lim_{t\to \infty }\frac{\phi _{7}(t)}{t}=0$ a.s. is obtained. By Lemma 7, and for all $\varepsilon >0$, there exists a constant $T>0$ for any $t\geq T$ such that $\ln x_{1}(t)<\frac{\varepsilon }{a_{32}a_{21}+1}t$, $\ln x_{2}(t)<\frac{\varepsilon }{a_{32}a_{11}+1}t$ and $\int _{0}^{t}x_{4}(s)\,\mathrm{d}s\leq (\limsup_{t\to \infty }\langle x_{4}(t) \rangle +\varepsilon )t$. Thus, for all $t\geq T$ and by (31)

$$ H_{2}\ln x_{3}(t)>(\Delta _{33}-2 \varepsilon )t-a_{34}H_{2} \Bigl(\limsup _{t \to \infty } \bigl\langle x_{4}(t) \bigr\rangle + \varepsilon \Bigr)t-H_{3} \int _{0}^{t}x_{3}(s) \, \mathrm{d}s+ \phi _{7}(t) $$

is furthermore obtained.

If $\Delta _{33}-a_{34}H_{2}\limsup_{t\to \infty }\langle x_{4}(t) \rangle >0$, then by Lemma 5 and the arbitrary ε we furthermore have

$$ \liminf_{t\to \infty } \bigl\langle x_{3}(t) \bigr\rangle \geq \frac{\Delta _{33}}{H_{3}}- \frac{a_{34}H_{2}}{H_{3}} \Bigl(\limsup_{t\to \infty } \bigl\langle x_{4}(t) \bigr\rangle \Bigr). $$

(32)

If $\Delta _{33}-a_{34}H_{2}\limsup_{t\to \infty }\langle x_{4}(t) \rangle \leq 0$, then since $\liminf_{t\to \infty }\langle x_{3}(t)\rangle \geq 0$, we also have

$$ \liminf_{t\to \infty } \bigl\langle x_{3}(t) \bigr\rangle \geq \frac{\Delta _{33}}{H_{3}}-\frac{a_{34}H_{2}}{H_{3}} \Bigl(\limsup _{t\to \infty } \bigl\langle x_{4}(t) \bigr\rangle \Bigr). $$

Provided that, for any $\varepsilon >0$, there is a constant $T>0$ satisfying for $t\geq T$

$$ \int _{0}^{t}x_{1}(s)\,\mathrm{d}s \leq \frac{a_{43}\Delta _{21}-b_{4}a_{12}a_{23}+a_{12}a_{23}a_{44}(\limsup_{t\to \infty }\langle x_{4}(t)\rangle +\varepsilon )}{H_{2}a_{43}} $$

and

$$ \int _{0}^{t}x_{3}(s)\,\mathrm{d}s \geq \frac{\Delta _{33}}{H_{3}}- \frac{a_{34}H_{2}}{H_{3}} \Bigl(\limsup _{t\to \infty } \bigl\langle x_{4}(t) \bigr\rangle - \varepsilon \Bigr). $$

Combining with (22), for all $t\geq T$,

$$ \begin{aligned} \ln x_{2}(t) \leq {}& \biggl[b_{2}+a_{21} \frac{a_{43}\Delta _{21}-b_{4}a_{12}a_{23}+a_{12}a_{23}a_{44} (\limsup_{t\to \infty }\langle x_{4}(t)\rangle +\varepsilon )}{H_{2}a_{43}} \\ &{}-a_{23} \biggl(\frac{\Delta _{33}}{H_{3}}-\frac{a_{34}H_{2}}{H_{3}} \Bigl( \limsup_{t\to \infty } \bigl\langle x_{4}(t) \bigr\rangle - \varepsilon \Bigr) \biggr) \biggr]t+\phi _{2}(t)-a_{22} \int _{0}^{t}x_{2}(s)\,\mathrm{d}s \end{aligned} $$

(33)

is furthermore obtained.

We have $\lim_{t\to \infty }\frac{\phi _{2}(t)}{t}=0$ a.s. by Lemma 7. We denote

$$ \begin{aligned} M_{1}={}& b_{2}+a_{21} \frac{(a_{43}\Delta _{21}-b_{4}a_{12}a_{23}+a_{12}a_{23}a_{44}\limsup_{t\to \infty } \langle x_{4}(t)\rangle )}{H_{2}a_{43}} \\ &{}-a_{23} \biggl(\frac{\Delta _{33}}{H_{3}}-\frac{a_{34}H_{2}}{H_{3}}\limsup _{t \to \infty } \bigl\langle x_{4}(t) \bigr\rangle \biggr). \end{aligned} $$

If $M_{1}\geq 0$, then we can obtain

$$ \begin{aligned} \limsup_{t\to \infty } \bigl\langle x_{2}(t) \bigr\rangle \leq {}&\frac{1}{a_{22}} \biggl[b_{2}+a_{21} \frac{a_{43}\Delta _{21}-b_{4}a_{12}a_{23}}{H_{2}a_{43}}-a_{23} \frac{\Delta _{33}}{H_{3}} \\ &{}+ \biggl(\frac{a_{12}a_{21}a_{23}a_{44}}{a_{43}H_{2}}+ \frac{a_{23}a_{34}H_{2}}{H_{3}} \biggr)\limsup _{t\to \infty } \bigl\langle x_{4}(t) \bigr\rangle \biggr] = \frac{M_{1}}{a_{22}}. \end{aligned} $$

(34)

If $M_{1}<0$, then $\lim_{t\to \infty }x_{2}(t)=0$ is directly obtained. From this and Lemma 8, $\lim_{t\to \infty }x_{j}(t)=0$, $j=3,4$, is furthermore derived.

Let $M_{1}\geq 0$, for all $\varepsilon >0$, there is a constant $T>0$ such that

$$ \int _{0}^{t}x_{4}(s)\,\mathrm{d}s \geq \biggl( \frac{\Delta _{44}}{H_{4}}-\varepsilon \biggr)t,\qquad \int _{0}^{t}x_{2}(s)\,\mathrm{d}s \leq \biggl( \frac{M_{1}}{a_{22}}+\varepsilon \biggr)t $$

for any $t\geq T$. From (22), (27) and (34), we derive for any $t\geq T$

$$ \ln x_{3}(t)\leq \biggl(b_{3}+a_{32} \biggl(\frac{M_{1}}{a_{22}}+ \varepsilon \biggr)-a_{34} \biggl( \frac{\Delta _{44}}{H_{4}}-\varepsilon \biggr) \biggr)t-a_{33} \int _{0}^{t}x_{3}(s) \, \mathrm{d}s+ \phi _{4}(t). $$

(35)

We have $\lim_{t\to \infty }\frac{\phi _{3}(t)}{t}=0$ a.s. by Lemma 7. We denote

$$ M_{2}=b_{3}+\frac{a_{32}}{a_{22}}M_{1}-a_{34} \frac{\Delta _{44}}{H_{4}}. $$

If $M_{2}\geq 0$, then from the arbitrary ε and Lemma 4 we furthermore have

$$ \limsup_{t\to \infty } \bigl\langle x_{3}(t) \bigr\rangle \leq \frac{1}{a_{33}} \biggl[b_{3}+ \frac{a_{32}M_{1}}{a_{22}}-\frac{a_{34}\Delta _{44}}{H_{4}} \biggr]. $$

(36)

If $M_{2}<0$, then we obtain $\lim_{t\to \infty }x_{3}(t)=0$. From this and Lemma 8, $\lim_{t\to \infty }x_{4}(t)=0$ is furthermore obtained.

Let $M_{2}\geq 0$. From (36) and for any $\varepsilon >0$, there exists a constant $T>0$, we obtain

$$ \int _{0}^{t}x_{3}(s)\,\mathrm{d}s \leq \frac{1}{a_{33}} \biggl[b_{3}+ \frac{a_{32}M_{1}}{a_{22}}- \frac{a_{34}\Delta _{44}}{H_{4}}+ \varepsilon \biggr]t $$

for $t\geq T$. From (23), we derive

$$ \ln x_{4}(t)\leq \biggl[b_{4}+\frac{a_{43}}{a_{33}} \biggl(b_{3}+ \frac{a_{32}M_{1}}{a_{22}}-\frac{a_{34}\Delta _{44}}{H_{4}}+ \varepsilon \biggr) \biggr]t -a_{44} \int _{0}^{t}x_{4}(s)\,\mathrm{d}s+ \phi _{4}(t) $$

(37)

for any $t\geq T$. We have $\lim_{t\to \infty }\frac{\phi _{4}(t)}{t}=0$ a.s. by Lemma 7. We denote

$$ M_{3}=b_{4}+\frac{a_{43}}{a_{33}} \biggl(b_{3}+ \frac{a_{32}M_{1}}{a_{22}}- \frac{a_{34}\Delta _{44}}{H_{4}} \biggr). $$

If $M_{3}\geq 0$, then from the arbitrary ε and Lemma 4

$$ \begin{aligned} \limsup_{t\to \infty } \bigl\langle x_{4}(t) \bigr\rangle \leq {}&\frac{1}{a_{44}} \biggl[b_{4}+ \frac{a_{43}}{a_{33}} \biggl[b_{3}+ \frac{a_{32}}{a_{22}}M_{1}-a_{34} \frac{\Delta _{44}}{H_{4}} \biggr] \biggr] \\ ={}&\frac{1}{a_{44}} \biggl[b_{4}+\frac{a_{43}}{a_{33}} \biggl[b_{3}+ \frac{a_{32}}{a_{22}} \biggl[b_{2}+a_{21} \frac{a_{43}\Delta _{21}-b_{4}a_{12}a_{23}}{H_{2}a_{43}} \\ &{}-a_{23}\frac{\Delta _{33}}{H_{3}} \biggr]-a_{34} \frac{\Delta _{44}}{H_{4}} \biggr] \biggr] \\ &{}+ \frac{a_{12}a_{21}a_{23}a_{32}a_{44}H_{3}+a_{23}a_{32}a_{34}a_{43}H_{2}^{2}}{a_{22}a_{33}a_{44}H_{2}H_{3}} \limsup_{t\to \infty } \bigl\langle x_{4}(t) \bigr\rangle \end{aligned} $$

(38)

is furthermore obtained. By a detailed calculation we can obtain

$$ \begin{aligned} &\frac{1}{a_{44}} \biggl[b_{4}+ \frac{a_{43}}{a_{33}} \biggl[b_{3}+ \frac{a_{32}}{a_{22}} \biggl[b_{2}+a_{21} \frac{a_{43}\Delta _{21}-b_{4}a_{12}a_{23}}{H_{2}a_{43}}-a_{23} \frac{\Delta _{33}}{H_{3}} \biggr]-a_{34}\frac{\Delta _{44}}{H_{4}} \biggr] \biggr] \\ &\quad =\frac{1}{a_{22}a_{33}a_{44}} \bigl[a_{22}a_{33}a_{44}H_{2}H_{3}-a_{12}a_{21}a_{23}a_{32}a_{44}H_{3}-a_{23}a_{32}a_{34}a_{43}H_{2}^{2} \bigr] \frac{\Delta _{44}}{H_{4}}. \end{aligned} $$

Thus, we furthermore find that (38) is equal with the inequality as follows:

$$ \begin{aligned} & \bigl[a_{22}a_{33}a_{44}H_{2}H_{3}-a_{12}a_{21}a_{23}a_{32}a_{44}H_{3}-a_{23}a_{32}a_{34}a_{43}H_{2}^{2} \bigr] \limsup_{t\to \infty } \bigl\langle x_{4}(t) \bigr\rangle \\ &\quad \leq \bigl[a_{22}a_{33}a_{44}H_{2}H_{3}-a_{12}a_{21}a_{23}a_{32}a_{44}H_{3}-a_{23}a_{32}a_{34}a_{43}H_{2}^{2} \bigr] \frac{\Delta _{44}}{H_{4}}. \end{aligned} $$

(39)

If $M_{3}<0$, then from (37) and Lemma 5 we directly have $\lim_{t\to \infty }x_{4}(t)=0$.

Assume $\Delta _{44}>0$, then we can obtain

$$\begin{aligned}& \begin{aligned} M_{1}\geq {}& b_{2}+a_{21} \frac{a_{43}\Delta _{21}-b_{4}a_{12}a_{23}+a_{12}a_{23}a_{44}\frac{\Delta _{44}}{H_{4}}}{H_{2}a_{43}} \\ &{}-a_{23}\frac{\Delta _{33}}{H_{3}}+ \frac{a_{23}a_{34}H_{2}\Delta _{44}}{H_{3}H_{4}}=a_{22} \frac{\Delta _{41}}{H_{4}}>0, \end{aligned} \\& M_{2}\geq b_{3}+\frac{a_{32}\Delta _{41}}{H_{4}}- \frac{a_{34}\Delta _{44}}{H_{4}}=a_{33}\frac{\Delta _{43}}{H_{4}}>0, \\& M_{3}\geq b_{4}+b_{3} \frac{a_{43}}{a_{33}}- \frac{a_{32}a_{43}\Delta _{41}}{a_{33}H_{4}}- \frac{a_{43}\Delta _{44}}{a_{33}a_{34}H_{4}} =a_{44} \frac{\Delta _{44}}{H_{4}}> 0. \end{aligned}$$

Hence, from (39) and condition (25), $\limsup_{t\to \infty }\langle x_{4}(t)\rangle \leq \frac{\Delta _{44}}{H_{4}}$ is obtained. Hence, $\lim_{t\to \infty }\langle x_{4}(t)\rangle = \frac{\Delta _{44}}{H_{4}}$ is directly derived.

From the above conclusion and (30), we have

$$ \begin{aligned} \limsup_{t\to \infty } \bigl\langle x_{1}(t) \bigr\rangle \leq{} &\frac{(a_{43}\Delta _{21}-b_{4}a_{12}a_{23}+a_{12}a_{23}a_{44}\frac{\Delta _{44}}{H_{4}})}{H_{2}a_{43}} \\ ={}&\frac{b_{1}(a_{22}a_{34}a_{43}+a_{22}a_{33}a_{44}+a_{23}a_{32}a_{44})-b_{2}a_{12}(a_{33}a_{44}+a_{34}a_{43})}{H_{4}} \\ &{}+\frac{-b_{4}a_{12}a_{23}a_{34}+b_{3}a_{12}a_{23}a_{44}}{H_{4}}= \frac{\Delta _{41}}{H_{4}}. \end{aligned} $$

(40)

Then, from (32) we furthermore obtain

$$ \liminf_{t\to \infty } \bigl\langle x_{3}(t) \bigr\rangle \geq \frac{\Delta _{33}}{H_{3}}- \frac{a_{34}H_{2}}{H_{3}} \frac{\Delta _{44}}{H_{4}}=\frac{\Delta _{43}}{H_{4}}. $$

(41)

Similarly, by (33)

$$ \limsup_{t\to \infty } \bigl\langle x_{2}(t) \bigr\rangle \leq \frac{b_{2}H_{4}+a_{21}\Delta _{41}-a_{23}\Delta _{43}}{a_{22}H_{4}}= \frac{\Delta _{42}}{H_{4}} $$

(42)

is also obtained. For all $\varepsilon >0$, there is a $T>0$ for all $t\geq T$ such that $\int _{0}^{t}x_{2}(s)\,\mathrm{d}s<(\frac{\Delta _{42}}{H_{4}}+\varepsilon )t$ and $\int _{0}^{t}x_{4}(s)\,\mathrm{d}s>(\frac{\Delta _{44}}{H_{4}}-\varepsilon )t$. From (22), we compute

$$ \ln x_{3}(t) \leq b_{3}t+a_{32} \biggl(\frac{\Delta _{42}}{H_{4}}+ \varepsilon \biggr)-a_{33} \int _{0}^{t}x_{3}(s)\,\mathrm{d}s -a_{34} \biggl( \frac{\Delta _{44}}{H_{4}}-\varepsilon \biggr)+\phi _{3}(t). $$

(43)

We have $\lim_{t\to \infty }\phi _{3}(t)=0$. Thus, from the arbitrariness of ε and Lemma 5

$$ \limsup_{t\to \infty } \bigl\langle x_{3}(t) \bigr\rangle \leq \frac{b_{3}H_{4}+a_{32}\Delta _{42}-a_{34}\Delta _{44}}{a_{33}H_{4}}= \frac{\Delta _{43}}{H_{4}} $$

(44)

is furthermore derived. Hence, $\lim_{t\to \infty }\langle x_{3}(t)\rangle = \frac{\Delta _{43}}{H_{4}}$ is obtained.

From (21) and (22), we compute

$$ V_{2}(t)=\Delta _{22}t-H_{2} \int _{0}^{t}x_{2}(s) \, \mathrm{d}s-a_{11}a_{23} \int _{0}^{t}x_{3}(s)\,\mathrm{d}s+ \phi _{8}(t), $$

(45)

where $\phi _{8}(t)=a_{21}\phi _{1}(t)+a_{11}\phi _{2}(t)$. By Lemma 7, $\lim_{t\to \infty }\frac{\phi _{8}(t)}{t}=0$ a.s. is obtained. From Lemma 7 and for $\varepsilon >0$, there exists a $T>0$ satisfying $\int _{0}^{t}x_{3}(s)\,\mathrm{d}s<(\frac{\Delta _{43}}{H_{4}}+\varepsilon )t$ and $\ln x_{1}(t)<\frac{\varepsilon }{a_{21}+1}t$ for $t>T$. Thus

$$ a_{11}\ln x_{2}(t)\geq \biggl(\Delta _{22}-a_{11}a_{23} \biggl( \frac{\Delta _{43}}{H_{4}}+\varepsilon \biggr)-\varepsilon \biggr) t-H_{2} \int _{0}^{t}x_{2}(s) \,\mathrm{d}s+ \phi _{8}(t) $$

(46)

is obtained. Therefore, from the arbitrariness of ε and Lemma 5

$$ \liminf_{t\to \infty } \bigl\langle x_{2}(t) \bigr\rangle \geq \frac{H_{4}\Delta _{22}-a_{11}a_{23}\Delta _{43}}{H_{2}H_{4}}= \frac{\Delta _{42}}{H_{4}} $$

(47)

is furthermore obtained. Then, we obtain $\lim_{t\to \infty }\langle x_{2}(t)\rangle = \frac{\Delta _{42}}{H_{4}}$.

For any $\varepsilon >0$, there is a $T>0$ such that $\int _{0}^{t} x_{2}(s)\,\mathrm{d}s<(\frac{\Delta _{42}}{H_{4}}+ \varepsilon )$ for any $t>T$. Hence, from (21), we have

$$ \ln x_{1}(t)\geq \biggl(b_{1}-a_{12} \biggl(\frac{\Delta _{42}}{H_{4}}+ \varepsilon \biggr) \biggr)t-a_{11} \int _{0}^{t}x_{1}(s)\,\mathrm{d}s+ \phi _{1}(t). $$

(48)

Hence, by Lemma 5 and the arbitrariness of ε we furthermore have

$$ \liminf_{t\to \infty } \bigl\langle x_{1}(t) \bigr\rangle \geq \frac{b_{1}H_{4}-a_{12}\Delta _{42}}{a_{11}H_{4}}= \frac{\Delta _{41}}{H_{4}}. $$

(49)

Then, we also have $\lim_{t\to \infty }\langle x_{1}(t)\rangle = \frac{\Delta _{41}}{H_{4}}$. Therefore, conclusion (9) in Theorem 1 is proved.

Assume $\Delta _{44}=0$. If there an $i\in \{1,2,3\}$ such that $M_{i}<0$, then we furthermore have $\lim_{t\to \infty }x_{i+1}(t)=0$ from the above discussions. Hence, by Lemma 8, $\lim_{t\to \infty }x_{4}(t)=0$. Otherwise, we have $M_{i}\geq 0$, $i=1,2,3$. Then, from the above discussions we also have

$$ \begin{aligned} & \bigl[a_{22}a_{33}a_{44}H_{2}H_{3}-a_{12}a_{21}a_{23}a_{32}a_{44}H_{3}-a_{23}a_{32}a_{34}a_{43}H_{2}^{2} \bigr] \limsup_{t\to \infty } \bigl\langle x_{4}(t) \bigr\rangle \\ &\quad \leq \bigl[a_{22}a_{33}a_{44}H_{2}H_{3}-a_{12}a_{21}a_{23}a_{32}a_{44}H_{3}-a_{23}a_{32}a_{34}a_{43}H_{2}^{2} \bigr] \frac{\Delta _{44}}{H_{4}}=0. \end{aligned} $$

Therefore, from condition (25) we have $\lim_{t\to \infty }\langle {x_{4}}\rangle =0$ a.s.

Assume $\Delta _{44}<0$. Then from (26) we obtain

$$ V_{4}(t)\leq \Delta _{44}t+\phi _{5}(t). $$

Hence,

$$ \limsup_{t\to \infty }\frac{1}{t} \bigl[ \bigl[a_{32} \bigl(a_{21}\ln x_{1}(t)+a_{11} \ln x_{2}(t) \bigr)+H_{2}\ln x_{3}(t) \bigr]a_{43}+H_{3}\ln x_{4}(t) \bigr]\leq \Delta _{44}< 0. $$

Thus, we have

$$ \lim_{t\to \infty } \bigl[ \bigl(x_{1}(t) \bigr)^{a_{43}a_{32}a_{21}} \bigl(x_{2}(t) \bigr)^{a_{43}a_{32}a_{11}} \bigl(x_{3}(t) \bigr)^{a_{43}H_{2}} \bigl(x_{4}(t) \bigr)^{H_{3}} \bigr]=0. $$

This shows that there exists a $j\in \{1,2,3,4\}$ that satisfies $\lim_{t\to \infty }x_{j}(t)=0$. Consequently, by Lemma 8, $\lim_{t\to \infty }x_{4}(t)=0$.

In conclusion, when $\Delta _{44}\leq 0$ we always obtain $\lim_{t\to \infty }\langle x_{4}(t)\rangle =0$ or $\lim_{t\to \infty }x_{4}(t)=0$. Thus, applying the similar arguments used in the proving process of Theorem 1 listed in [28], the remaining conclusions in Theorem 1 can be proved. □

Remark 4

Observe the proving process of the above, the criterion (25) only used to obtain $\limsup_{t\to \infty }\langle x_{4}(t)\rangle \leq \frac{\Delta _{44}}{H_{4}}$ from the inequality (39). This shows that conditions (24) and (25) appear to be the supererogatory and pure mathematical conditions.

Remark 5

An important and interesting open problem is how to extend Theorem 1 to the general n-species stochastic food-chain system with distributed delay and harvesting.

In the following theorem, we mainly investigate that, for all positive global solutions of model (1), the conclusion about global attractiveness in the expectation.

Theorem 2

For initial conditions $\phi ,\phi ^{*}\in C([-r,0],R_{+}^{4})$, assume that model (1) has two solutions $(x_{1}(t;\phi ),x_{2}(t;\phi ),x_{3}(t;\phi ),x_{4}(t;\phi ))$ and $(y_{1}(t;\phi ^{*}),y_{2}(t;\phi ^{*}),y_{3}(t;\phi ^{*}), y_{4}(t;\phi ^{*}))$. If there are positive constants $w_{i}$ ($i=1,2,3,4$) such that

$$\begin{aligned}& w_{1}a_{11}-w_{2}a_{21}>0,\qquad w_{i}a_{ii}-w_{i-1}a_{i-1i}-w_{i+1}a_{i+1i}>0\quad (i=2,3), \\& w_{4}a_{44}-w_{3}a_{34}>0. \end{aligned}$$

Then

$$ \lim_{t\to \infty }E \Biggl(\sum_{i=1}^{4} \bigl\vert x_{i}(t,\phi )-y_{i} \bigl(t,\phi ^{*} \bigr) \bigr\vert ^{2} \Biggr)^{ \frac{1}{2}}=0. $$

The proof of Theorem 2 is similar to Theorem 2 from [28]. Hence it is omitted here. Now let $\mathcal{P}([-r,0],R_{+}^{4})$ represent the whole measurable probability space on $C([-r,0],R_{+}^{4})$. For $\mathcal{P}_{1},\mathcal{P}_{2}\in \mathcal{P}([-r,0],R_{+}^{4})$, set the metric as follows:

$$ d_{L}(\mathcal{P}_{1},\mathcal{P}_{2})= \sup_{f\in {L}} \biggl\vert \int _{R_{+}^{4}}f(u) \mathcal{P}_{1}(\mathrm{d}u)- \int _{R_{+}^{4}}f(u)\mathcal{P}_{2}(\mathrm{d}u) \biggr\vert , $$

where

$$ L= \bigl\{ f:\mathcal{C} \bigl([-r,0],R_{+}^{4} \bigr) \rightarrow R: \bigl\vert f(u_{1})-f(u_{2}) \bigr\vert \leq \Vert u_{1}-u_{2} \Vert , \bigl\vert f(\cdot ) \bigr\vert \leq 1 \bigr\} . $$

Let $p(t,\phi ,\mathrm{d}x)$ represents the transition probability of process $x(t)=(x_{1}(t),x_{2}(t),x_{3}(t),x_{4}(t))$. In the following theorem, we consider the condition of asymptotically stability. The results as follows are obtained.

Theorem 3

Suppose that positive constants $q_{i}$ ($i=1,2,3,4$) satisfies

$$\begin{aligned}& q_{1}a_{11}-q_{2}a_{21}>0,\qquad q_{i}a_{ii}-q_{i-1}a_{i-1i}-q_{i+1}a_{i+1i}>0\quad (i=2,3),\\& q_{4}a_{44}-q_{3}a_{34}>0. \end{aligned}$$

Then model (1) is asymptotically stable in distribution, i.e., for all initial value $\phi \in C([-\gamma ,0],R_{+}^{4})$, a unique probability measure $v(\cdot )$ satisfies the transition probability $p(t,\phi ,\cdot )$ of solution $(x_{1}(t,\phi ),x_{2}(t,\phi ),x_{3}(t,\phi ),x_{4}(t,\phi ))$ such that

$$ \lim_{t\to \infty }d_{BL} \bigl(p(t,\phi ,\cdot ),v( \cdot ) \bigr)=0. $$

Remark 6

Obviously, Theorems 2 and 3 also seemingly can be extended to the general n-species stochastic food-chain system with distributed delay and harvesting.

4 Effect of harvesting

We firstly introduce the following lemma.

Lemma 9

Assume that there exist positive constants $q_{i}$ ($i=1,2,3,4$) satisfying

$$\begin{aligned}& q_{1}a_{11}-q_{2}a_{21}>0,\qquad q_{i}a_{ii}-q_{i-1}a_{i-1i}-q_{i+1}a_{i+1i}>0\quad (i=2,3), \\& q_{4}a_{44}-q_{3}a_{34}>0. \end{aligned}$$

Then we have

$$ a_{22}a_{33}a_{44}H_{2}H_{3}-a_{12}a_{21}a_{23}a_{32}a_{44}H_{3}-a_{23}a_{32}a_{34}a_{43}H_{2}^{2}>0. $$

Proof

In fact, we obtain

$$ a_{11}>\frac{q_{2}}{q_{1}}a_{21},\qquad a_{ii}> \frac{1}{q_{i}}(q_{i-1}a_{i-1i}+q_{i+1}a_{i+1i}),\quad i=2,3, \qquad a_{44}>\frac{q_{3}}{q_{4}}a_{34}. $$

By calculating we obtain

$$ \begin{aligned} &(a_{22}a_{33}H_{2}-a_{12}a_{21}a_{23}a_{32})a_{44}H_{3} \\ &\quad > \biggl(\frac{1}{q_{2}}[q_{1}a_{12}+q_{3}a_{32}] \frac{1}{q_{3}}[q_{2}a_{23}+q_{4}a_{43}] [a_{11}a_{22}+a_{12}a_{21}]-a_{12}a_{21}a_{23}a_{32} \biggr)a_{44}H_{3} \\ &\quad \geq \biggl( \biggl(\frac{1}{q_{2}}q_{1}a_{12} \frac{1}{q_{3}}[q_{2}a_{23}+q_{4}a_{43}]+ \frac{q_{4}}{q_{2}}a_{32}a_{43} \biggr) (a_{11}a_{22}+a_{12}a_{21}) \biggr)a_{44}H_{3} \\ &\quad \geq \frac{q_{4}}{q_{2}}a_{32}a_{43}(a_{11}a_{22}+a_{12}a_{21}) \frac{q_{3}}{q_{4}}a_{34}H_{3}. \end{aligned} $$

Since

$$ H_{3}>(a_{11}a_{22}+a_{12}a_{21})a_{33}>(a_{11}a_{22}+a_{12}a_{21}) \frac{q_{2}}{q_{3}}a_{23}, $$

we furthermore obtain

$$ \begin{aligned} &(a_{22}a_{33}H_{2}-a_{12}a_{21}a_{23}a_{32})a_{44}H_{3} \\ &\quad > \frac{q_{4}}{q_{2}}a_{32}a_{43}(a_{11}a_{22}+a_{12}a_{21}) \frac{q_{3}}{q_{4}}a_{34}(a_{11}a_{22}+a_{12}a_{21}) \frac{q_{2}}{q_{3}}a_{23} \\ &\quad =a_{23}a_{32}a_{34}a_{43}H_{2}^{2}. \end{aligned} $$

This completes the proof. □

For the convenience, we define the following matrix:

$$ B= \begin{pmatrix} a_{11}&a_{12}&0&0 \\ -a_{21}&a_{22}&a_{23}&0 \\ 0&-a_{32}&a_{33}&a_{34} \\ 0&0&-a_{43}&a_{44} \end{pmatrix}. $$

It is clear that the determinant $\vert B \vert =H_{4}>0$. Hence, there exists $B^{-1}$. Let $H=(h_{1},h_{2},h_{3},h_{4})^{T}$ and $R=(r_{1}-\frac{1}{2}\sigma _{1}^{2}, r_{2}-\frac{1}{2}\sigma _{2}^{2}, r_{3}-\frac{1}{2}\sigma _{3}^{2}, r_{4}-\frac{1}{2}\sigma _{4}^{2})^{T}$. Furthermore, let $H^{*}=(h_{1}^{*},h_{2}^{*},h_{3}^{*},h_{4}^{*})^{T}=(B(B^{-1})^{T}+E)^{-1}R$, where E is the unit matrix.

Theorem 4

Suppose that the positive constants $q_{i}$ ($i=1,2,3,4$) satisfy

$$\begin{aligned}& q_{1}a_{11}-q_{2}a_{21}>0,\qquad q_{i}a_{ii}-q_{i-1}a_{i-1i}-q_{i+1}a_{i+1i}>0\quad (i=2,3), \\& q_{4}a_{44}-q_{3}a_{34}>0. \end{aligned}$$

Then the following conclusions hold.

($\mathcal{A}_{1}$):

If $h_{i}^{*}\geq 0$ for $i=1,2,3,4$, $\Delta _{44}|_{h_{i}=h_{i}^{*},i=1,2,3,4}>0$ and $B^{-1}+(B^{-1})^{T}$ is positive semi-definite, then model (1) has the optimal harvesting strategy $H=H^{*}$ and

$$ MESY\triangleq Y \bigl(H^{*} \bigr)= \bigl(H^{*} \bigr)^{T}B^{-1} \bigl(R-H^{*} \bigr). $$

(50)

($\mathcal{A}_{2}$):

If any of the following conditions holds:

($\mathcal{B}_{1}$):

$\Delta _{44}|_{h_{i}=h_{i}^{*},i=1,2,3,4}\leq 0$;

($\mathcal{B}_{2}$):

$h_{1}^{*}<0$ or $h_{2}^{*}<0$ or $h_{3}^{*}<0$ or $h_{4}^{*}<0$;

($\mathcal{B}_{3}$):

$B^{-1}+(B^{-1})^{T}$ is not positive semi-definite,

then the optimal harvesting strategy of model (1) does not exist.

Proof

Let $\mathcal{U}=\{H=(h_{1},h_{2},h_{3},h_{4})^{T}\in R^{4}:\Delta _{44}>0,h_{i} \geq 0,i=1,2,3,4\}$. Obviously, for $H\in \mathcal{U}$, the conclusion (9) of Theorem 1 stands. Meanwhile, supposing $H^{*}$ exists, then $H^{*}\in \mathcal{U}$.

At the beginning, let us prove $(\mathcal{A}_{1})$. The set $\mathcal{U}$ is not empty for $H^{*}\in \mathcal{U}$. From Theorem 3, a unique invariant measure $v(\cdot )$ for model (1) exists. And thus it yields by Corollary 3.4.3 in Prato and Zbczyk [29] that $v(\cdot )$ is strongly mixing. The measure $v(\cdot )$ is ergodic from Theorem 3.2.6 in [29]. For initial condition $(\varsigma (\theta ),\xi (\theta ),\kappa (\theta ),\eta (\theta )) \in C([-r,0],R_{+}^{4})$, model (1) has a global positive solution $x(t)=(x_{1}(t),x_{2}(t),x_{3}(t),x_{4}(t))$. In view of Theorem 3.3.1 in [29], we obtain

$$ \lim_{t\rightarrow \infty }\frac{1}{t} \int _{0}^{t}H^{T}x(s)\,\mathrm{d}s= \int _{R_{+}^{4}}H^{T}xv(\mathrm{d}x), $$

(51)

where $H=(h_{1},h_{2},h_{3},h_{4})^{T}\in \mathcal{U}$. Let $\varrho (z)$ be the stationary probability density of model (1), then we have

$$ Y(H)=\lim_{t\rightarrow \infty }E \Biggl[\sum _{i=1}^{4}h_{i}x_{i}(t) \Biggr]=\lim_{t \rightarrow \infty }E \bigl[H^{T}x(t) \bigr]= \int _{R_{+}^{4}}H^{T}x\varrho (x)\,\mathrm{d}x. $$

(52)

For model (1), in view of the invariant measure is sole, one also has a one-to-one correspondence among $\varrho (z)$ and its corresponding invariant measure;

$$ \int _{R_{+}^{4}}H^{T}x\varrho (x)\,\mathrm{d}x= \int _{R_{+}^{4}}H^{T}xv(\mathrm{d}x) $$

(53)

is deduced. Therefore, from the conclusion (9) of Theorem 1, Lemma 9, and (51)–(53)

$$\begin{aligned} Y(H)={}&\lim_{t\to +\infty } \frac{1}{t} \int _{0}^{t}H^{T}x(s)\,\mathrm{d}s \\ ={}&h_{1}\lim_{t\to +\infty }\frac{1}{t} \int _{0}^{t}x_{1}(s) \, \mathrm{d}s+h_{2} \lim_{t\to +\infty }\frac{1}{t} \int _{0}^{t}x_{2}(s)\,\mathrm{d}s +h_{3} \lim_{t\to +\infty }\frac{1}{t} \int _{0}^{t}x_{3}(s)\,\mathrm{d}s \\ &{}+h_{4}\lim_{t\to +\infty }\frac{1}{t} \int _{0}^{t}x_{4}(s)\,\mathrm{d}s \\ ={}&h_{1}\Delta _{41}+h_{2}\Delta _{42}+h_{3}\Delta _{43}+h_{4} \Delta _{44} \end{aligned}$$

is obtained. It can be carefully calculated that $Y(H)=H^{T}B^{-1}(R-H)$. Calculating the gradient of $Y(H)$, we have

$$ \frac{\partial Y(H)}{\partial H}=\frac{\partial H^{T}}{\partial H}B^{-1}(R-H)+ \frac{\partial (R-H)^{T}}{\partial H} \bigl(B^{-1} \bigr)^{T}H. $$

Since $\frac{\partial H^{T}}{\partial H}=E$ is unit matrix, we furthermore have

$$ \frac{\partial Y(H)}{\partial H}=B^{-1}(R-H)- \bigl(B^{-1} \bigr)^{T}H=B^{-1}R- \bigl(B^{-1}+ \bigl(B^{-1} \bigr)^{T} \bigr)H. $$

Solving the equation $\frac{\partial Y(H)}{\partial H}=0$, we obtain the critical value $H=(B^{-1}+(B^{-1})^{T})^{-1}B^{-1}R$. We have

$$ H= \bigl[B^{-1} \bigl(B \bigl(B^{-1} \bigr)^{T}+E \bigr) \bigr]^{-1}B^{-1}R= \bigl(B \bigl(B^{-1} \bigr)^{T}+E \bigr)^{-1}R=H^{*}. $$

Furthermore calculating the Hessian matrix of $Y(H)$, we obtain

$$ \frac{\partial }{\partial H} \biggl(\frac{\partial Y(H)}{\partial H} \biggr)=- \frac{\partial (B^{-1}+(B^{-1})^{T})H}{\partial H}=- \bigl(B^{-1}+ \bigl(B^{-1} \bigr)^{T} \bigr). $$

(54)

Since $B^{-1}+{(B^{-1})}^{T}$ is positive semi-definite, from the existence principle of extremum value for multivariable functions, we find that $Y(H)$ has the maximum global value $H=H^{*}$. Clearly, $H^{*}$ is unique, hence if $H^{*}\in \mathcal{U}$, i.e., $h^{*}_{i}\geq 0$ ($i=1,2,3,4$) and $\Delta _{44}|_{h_{i}=h^{*}_{i},i=1,2,3,4}>0$, thus we finally obtain the result that $H^{*}$ is an optimal harvesting strategy, and MESY shown in (50).

Now we need to prove $(\mathcal{A}_{2})$. We first assume that $(\mathcal{B}_{1})$ or $(\mathcal{B}_{2})$ stands. Suppose that $\widetilde{\Gamma }=(\gamma _{1},\gamma _{2},\gamma _{3},\gamma _{4})$ is the optimal harvesting strategy, thus $\Gamma \in \mathcal{U}$. That is,

$$ \Delta _{44}|_{h_{i}=\gamma _{i},i=1,2,3,4}>0,\quad \gamma _{i}\geq 0, i=1,2,3,4. $$

(55)

Then again, if $\Gamma \in \mathcal{U}$ is the optimal harvesting strategy, we find that Γ is the unique solution of the equation $\frac{\partial Y(H)}{\partial H}=0$. Therefore, we have $( h_{1}^{*},h_{2}^{*},h_{3}^{*},h_{4}^{*})=(\gamma _{1},\gamma _{2}, \gamma _{3},\gamma _{4})$. Thus, condition (55) becomes

$$ \Delta _{44}|_{h_{i}= h_{i}^{*},i=1,2,3,4}>0, \quad h_{i}^{*} \geq 0, i=1,2,3,4, $$

which is impossible.

Lastly, let us consider $(\mathcal{B}_{3})$. We first assume that $(\mathcal{B}_{1})$ and $(\mathcal{B}_{2})$ fail to stand. Thus, $h_{i}^{*}\geq 0$, $i=1,2,3,4$, and $\Delta _{44}|_{h_{i}= h_{i}^{*},i=1,2,3,4}>0$. Thus, $\mathcal{U}$ is not empty. That is to say, (51)–(54) hold. Let $B^{-1}+(B^{-1})^{T}=(b_{ij})_{4\times 4}$. Then, by calculating we have

$$ b_{11}= \frac{2(a_{22}a_{33}a_{44}+a_{22}a_{34}a_{43}+a_{23}a_{32}a_{44})}{H_{4}}. $$

Obviously, $b_{11}>0$. In other words, $B^{-1}+(B^{-1})^{T}$ is not negative semi-definite. From $\mathcal{B}_{3}$, we see that $B^{-1}+(B^{-1})^{T}$ is indefinite. Therefore, there is no optimal harvesting strategy if $\mathcal{B}_{3}$ holds. □

Remark 7

We easily observe from the above proving process of Theorem 4 that, for the general n-species stochastic food-chain system with distributed delay and harvesting, similar results can be established.

5 Numerical examples

Next, we give three examples and a few figures to illustrate our main results. The numerical methods are proposed in the numerical examples section of [28]. In model (1), we indicate the initial conditions $x_{1}(\theta )=0.3e^{\theta }$, $x_{2}(\theta )=0.2e^{\theta }$, $x_{3}(\theta )=0.3e^{\theta }$ and $x_{4}(\theta )=0.2e^{\theta }$ for all $\theta \in [-\ln 2,0]$, and $\tau _{12}=\tau _{21}=\tau _{23}=\tau _{32}=\tau _{34}=\tau _{43}= \ln 2$ in the numerical simulations as follows.

Example 1

The parameters $r_{1}=2.0$, $r_{2}=-1.0$, $r_{3}=-0.5$, $r_{4}=-0.1$ and $h_{1}=h_{2}=h_{3}=h_{4}=0$ are set. It is assumed that the parameters set for model (1) are as shown below.

Case 1.1. $a_{11}=1$, $a_{22}=1$, $a_{33}=2$, $a_{44}=0.5$, $a_{12}=2$, $a_{21}=2$, $a_{23}=1$, $a_{32}=2$, $a_{34}=1$, $a_{43}=1$, $\sigma _{1}=0.5$, $\sigma _{2}=0.3$, $\sigma _{3}=0.9$ and $\sigma _{4}=0.9$. We have $\Delta _{33}=0.8850>0$, $\Delta _{44}=-5.1750<0$ and $a_{22}a_{33}H_{2}-a_{12}a_{21}a_{23}a_{32}=2.0000>0$. Thus, based on the conclusion (7) in Theorem 1, one shows that $x_{i}(t)$ feature persistence in mean for $i=1,2,3$ and $x_{4}(t)$ is extinct. Figure 1 shows the dynamic responses of $x_{i}(t)$, for $i=1,2,3,4$.

Figure 1

The time series diagram shows that for ${i=1,2,3}$, $x_{i}(t)$ is persistent in mean, $x_{4}(t)$ is extinct

Case 1.2. $a_{11}=0.8$, $a_{22}=1$, $a_{33}=2.5$, $a_{44}=1.8$, $a_{12}=1$, $a_{21}=2$, $a_{23}=1$, $a_{32}=2$, $a_{34}=0.3$, $a_{43}=1$, $\sigma _{1}=0.1$, $\sigma _{2}=0.1$, $\sigma _{3}=0.2$ and $\sigma _{4}=0.9711$. We have $\Delta _{44}=0$ and $(a_{22}a_{33}H_{2}-a_{12}a_{21}a_{23}a_{32})a_{44}H_{3}-a_{23}a_{32}a_{34}a_{43}H_{2}^{2}=3.0360>0$. Thus, based on the conclusion (8) in Theorem 1, one shows that $x_{i}(t)$ feature persistence in mean for $i=1,2,3$ and $x_{4}(t)$ is extinct in mean. Figure 2 shows the dynamic responses of $x_{i}(t)$, for $i=1,2,3,4$.

Figure 2

The time series diagram shows that for ${i=1,2,3}$, $x_{i}(t)$ is persistent in mean, $x_{4}(t)$ is extinct in mean

Case 1.3. $a_{11}=0.5$, $a_{22}=2$, $a_{33}=2.5$, $a_{44}=1.2$, $a_{12}=1$, $a_{21}=2.5$, $a_{23}=2$, $a_{32}=2$, $a_{34}=0.6$, $a_{43}=2$, $\sigma _{1}=0.1$, $\sigma _{2}=0.2$, $\sigma _{3}=0.5$ and $\sigma _{4}=0.5$. We have $\Delta _{44}=11.1163>0$ and $(a_{22}a_{33}H_{2}-a_{12}a_{21}a_{23}a_{32})a_{44}H_{3}-a_{23}a_{32}a_{34}a_{43}H_{2}^{2}=5.7000>0$. Thus, based on the conclusion (9) in Theorem 1, one shows that $x_{i}(t)$ are persistent in mean. Figure 3 shows the dynamic responses of $x_{i}(t)$. Here, $i=1,2,3,4$.

Figure 3

The time series diagram shows that for ${i=1,2,3,4}$, $x_{i}(t)$ is persistent in mean

Example 2

In model (1), parameters $r_{1}=2.0$, $r_{2}=-1.0$, $r_{3}=-0.5$, $r_{4}=-0.1$ and $h_{1}=h_{2}=h_{3}=h_{4}=0$ are fixed. It is assumed that the parameters set for model (1) are shown below.

Case 2.1. $a_{11}=1$, $a_{22}=1$, $a_{33}=2$, $a_{44}=0.5$, $a_{12}=2$, $a_{21}=2$, $a_{23}=1$, $a_{32}=2$, $a_{34}=1$, $a_{43}=1$, $\sigma _{1}=0.5$, $\sigma _{2}=0.3$, $\sigma _{3}=0.9$ and $\sigma _{4}=0.9$. We have $\Delta _{33}=0.8850>0$, $\Delta _{44}=-2.6500<0$ and $a_{22}a_{33}H_{2}-a_{12}a_{21}a_{23}a_{32}=-3<0$. Clearly, the criterion of conclusion (7) in Theorem 1 is not met. However, the dynamic responses of $x_{i}(t)$ ($i=1,2,3,4$), which are given in Fig. 4, show that $x_{i}(t)$ feature persistence in mean and $x_{4}(t)$ is extinct, here $i=1,2,3$.

Figure 4

The time series diagram shows that for ${i=1,2,3}$, $x_{i}(t)$ is persistent in mean, $x_{4}(t)$ is extinct

Case 2.2. $a_{11}=0.8$, $a_{22}=1$, $a_{33}=2.5$, $a_{44}=1.8$, $a_{12}=1$, $a_{21}=2$, $a_{23}=1$, $a_{32}=2$, $a_{34}=1$, $a_{43}=1$, $\sigma _{1}=0.1$, $\sigma _{2}=0.1$, $\sigma _{3}=0.2$ and $\sigma _{4}=0.9711$. We have $\Delta _{44}=0$ and $(a_{22}a_{33}H_{2}-a_{12}a_{21}a_{23}a_{32})a_{44}H_{3}-a_{23}a_{32}a_{34}a_{43}H_{2}^{2}=-7.9400<0$. Clearly, the criterion of conclusion (8) in Theorem 1 is not met. However, the dynamic responses of $x_{i}(t)/$ ($i=1,2,3,4$), which are given in Fig. 5, show that $x_{i}(t)$ feature persistence in mean and $x_{4}(t)$ is extinct in mean, here $i=1,2,3$.

Figure 5

The time series diagram shows that for ${i=1,2,3}$, $x_{i}(t)$ is persistent in mean, $x_{4}(t)$ is extinct in mean

Case 2.3. $a_{11}=0.5$, $a_{22}=2$, $a_{33}=2.5$, $a_{44}=1$, $a_{12}=1$, $a_{21}=2.5$, $a_{23}=2$, $a_{32}=2$, $a_{34}=0.6$, $a_{43}=2$, $\sigma _{1}=0.1$, $\sigma _{2}=0.2$, $\sigma _{3}=0.5$ and $\sigma _{4}=0.5$. We have $\Delta _{44}=11.1163>0$ and $(a_{22}a_{33}H_{2}-a_{12}a_{21}a_{23}a_{32})a_{44}H_{3}-a_{23}a_{32}a_{34}a_{43}H_{2}^{2}=-5.0500<0$. Clearly, the criterion of conclusion (9) in Theorem 1 is not met. However, the dynamic responses of $x_{i}(t)$, which are given in Fig. 6, show that $x_{i}(t)$ for ($i=1,2,3,4$) feature persistence in mean.

Figure 6

The time series diagram shows that for ${i=1,2,3,4}$, $x_{i}(t)$ is persistent in mean

Example 3

In model (1), take parameters $r_{1}=1.5$, $r_{2}=-0.5$, $r_{3}=-0.03$ and $r_{4}=-0.01$, $m_{1}=2.5$, $m_{2}=1.3$, $m_{3}=0.8$, $m_{4}=1.4$, $a_{11}=1.6$, $a_{12}=0.2$, $a_{22}=2$, $a_{21}=2.5$, $a_{23}=1$, $a_{32}=2.5$, $a_{33}=2$, $a_{34}=0.2$, $a_{43}=0.1$, $a_{44}=2$, $\sigma _{1}=0.1$, $\sigma _{2}=0.2$, $\sigma _{3}=0.1$ and $\sigma _{4}=0.1$. Then we have $m_{1}a_{11}-m_{2}a_{21}=0.7500>0$, $m_{2}a_{22}-m_{1}a_{12}-m_{3}a_{32}=0.1000>0$, $m_{3}a_{33}-m_{2}a_{23}-m_{4}a_{43}=0.1600>0$ and $m_{4}a_{44}-m_{3}a_{34}=2.6400>0$. Hence, a condition of Theorem 4 is satisfied, and by calculating we have $h_{1}^{*}=0.3377>0$, $h_{2}^{*}=0.4811$, $h_{3}^{*}=0.5147$, $h_{4}^{*}=0.0083$, and $\Delta _{44}=0.5424>0$. Thus, the criterion of conclusion ($\mathcal{A}_{1}$) in Theorem 4 is met. Then, the optimal harvesting strategy $H^{*}=(0.3377,0.4811,0.5147,0.0083)^{T}$ is obtained, we also have $Y(H^{*})=0.4316$. The dynamic response is shown in Fig. 7.

Figure 7

The time series diagram shows the optimal harvesting strategy

6 Conclusion

By investigating the effect of harvesting and distributed delays on the stochastic model and taking four species into accounts, we extend the main investigation in [28]. Using the inequality estimation technique, the Lyapunov function method and the stochastic integrals inequalities, in Theorem 1, the critical values between persistence in mean and extinction are investigated; as the result shows, environmental randomness can affect the extinction and persistence of a species in terms of the demographics of species and lower tropical species. However, the environmental randomness affects the average abundance of a species at all trophic levels. Global attractiveness and global asymptotic stability in distribution of model (1) are discussed in Theorem 2 and Theorem 3, respectively. The existence of the maximum of sustainable yield, the optimal harvesting strategy and the optimal harvesting effort and are obtained in Theorem 4, the result shows that the optimal harvesting strategy has a close relation with environmental fluctuations. Finally, numerical simulations are provided to support theoretical findings.

There are some issues that may need to be followed up to continue the discussion. First of all, similar research work (see [28]) for the general n-species random food-chain systems is not found up to now. From Remarks 1–7, we easily find that, for the general n-species stochastic food-chain system with distributed delay and harvesting, Theorems 2–4 can be established. However, the remaining problem is how to extend Theorem 1 to the general n-species random food-chain model. Secondly, it is supposed at present that the optimal harvesting problem of cooperative systems and competitive systems are less studied. Therefore, this may also be a breakthrough of the optimal harvesting problem. Furthermore, we can investigate more complex models, such as random models with nonlinear functional responses (see [22]), Markov switching (see [30]), and Lévy jumps (see [31, 32]). In the following research, we hope to discuss these issues.