Dynamics of a stochastic eutrophication-chemostat model with impulsive dredging and pulse inputting on environmental toxicant
Jianjun Jiao ORCID: orcid.org/0000-0002-4208-95981 &
Qiuhua Li2
In this paper, we present a stochastic eutrophication-chemostat model with impulsive dredging and pulse inputting on environmental toxicant. The sufficient condition for the extinction of microorganisms is obtained. The sufficient condition for the investigated system with unique ergodic stationary distribution is also obtained. The results show that the stochastic noise, impulsive dredging, and pulse input on the environmental toxicant play important roles in the extinction of microorganisms. The results also indicate the effective and reliable controlling strategy for water resource management. Finally, numerical simulations are employed to illustrate our results.
The chemostat is a device for continuous and impulsive cultures of microorganisms in laboratory [1–3]. Impulsive differential equations are found in almost every domain of applied science and have been studied in many investigations [4, 5]. With the development of society, the increasing amount of toxicants and contaminants have entered ecological systems. Environmental pollution has become one of the most important society-ecological problems. Therefore, it is very important to study the effects of toxicants on a population or community. Specially, the toxicant and abundant microorganisms in the water pollution environment are also a threat to the water resource management. Consequently, it is important to discuss chemostat models in a polluted environment [6, 7]. Zhou et al. [8] considered that reservoir dredging is the main and effective way to improve water quality by using a physical method. However, it is well known that many real-world systems may be disturbed by stochastic factors. Population systems are often subjected to various types of environmental noise. In ecology, it is critical to discover whether the presence of this noise has significant effects on population systems. Mao [9, 10] investigated stochastic differential equations and their applications. Lv et al. [11] presented an impulsive stochastic chemostat model with nonlinear perturbation.
Inspired by the above discussion, we consider a stochastic eutrophication-chemostat model with impulsive dredging and pulse inputting on environmental toxicant:
$$ \left \{ \textstyle\begin{array}{@{}l} \left . \textstyle\begin{array}{@{}l} dx(t)=[D(x_{0}-x(t)) \\ \hphantom{dx(t)={}}{}-\frac{ \beta x(t)y(t)}{k(A+x(t)+By(t))}]\,dt \\ \hphantom{dx(t)={}}{}+x(t)(\sigma _{11}+\sigma _{12}x(t))\,dB_{1}(t), \\ dy(t)=[\frac{\beta x(t)y(t)}{A+x(t)+By(t)} \\ \hphantom{dy(t)={}}{}-Dy(t)-rc_{o}(t)y(t)]\,dt \\ \hphantom{dy(t)={}}{}+y(t)(\sigma _{21}+\sigma _{22}y(t))\,dB_{2}(t), \\ dc_{o}(t)=(fc_{e}(t)-(g+m)c_{o}(t))\,dt, \\ dc_{e}(t)=(-hc_{e}(t))\,dt, \end{array}\displaystyle \right \}\quad t \neq (n+l)\tau , t\neq (n+1)\tau , \\ \left . \textstyle\begin{array}{@{}l} \triangle x(t)=0, \\ \triangle y(t)=-h_{1}y(t), \\ \triangle c_{o}(t)=0, \\ \triangle c_{e}(t)=-h_{2}c_{e}(t), \end{array}\displaystyle \right \}\quad t= (n+l)\tau ,n\in Z^{+}, \\ \left . \textstyle\begin{array}{@{}l} \triangle x(t)=0, \\ \triangle y(t)=0, \\ \triangle c_{o}(t)=0, \\ \triangle c_{e}(t)=\mu , \end{array}\displaystyle \right \}\quad t= (n+1)\tau ,n\in Z^{+}, \end{array}\displaystyle \right . $$
where \(x(t)\) is the concentration of the nutrient in a lake at time t. \(y(t)\) is the concentration of the microorganism in a lake at time t. \(c_{o}(t)\) is the concentration of the toxicant in the organism of the microorganism in a lake at time t. \(c_{e}(t)\) is the concentration of the toxicant in a lake at time t. D denotes the input rate from the lakes containing the nutrient and the wash-out rate of nutrients and microorganisms from the lake. \(\beta >0\) is the uptake constant of the nutrient. \(\frac{x(t)}{A+x(t)+By(t)}\) is a functional response of the Beddington–DeAngelis type. \(k>0\) is the yield of the microorganism y per unit mass of the nutrient. \(A>0\) and \(B>0\) are the saturating parameters of the Beddington–DeAngelis functional response. \(r>0\) is the depletion rate coefficient of the microorganism y due to the microorganism organismal toxicant. \(f>0\) is the coefficient of the population organism's net uptake of toxicant from the environment in a lake. \(-g<0\) and \(-m<0\), respectively, represent coefficients of the elimination and depuration rates of the toxicant in the organism in a lake. \(-h<0\) is the coefficient of the totality of toxicant losses from the system environment in a lake, including processes such as biological transformation, chemical hydrolysis, volatilization, microbial degradation, and photosynthetic degradation. τ is the period of impulsive dredging or the pulse input environmental toxin. \(0< h_{1}<1\) is the effect of impulsive dredging microorganism at time \(t=(n+l)\) (\(0< l<1\)). \(0< h_{2}<1\) is the effect of impulsive dredging environmental toxicant at time \(t=(n+l)\) (\(0< l<1\)). \(\mu \geq 0 \) is the amount of pulse input of environmental toxin concentration in a lake at \(t=(n+1)\tau \), \(n\in Z^{+}\), and \(Z^{+}=\{1,2,\ldots \}\).
The lemmas
In this paper, \((\varOmega , \mathscr{F},\mathscr{F}_{t\geq 0},P)\) stands for a complete probability space with filtration \(\mathscr{F}_{t\geq 0}\) satisfying the usual conditions. Define \(f^{l}=\inf_{t\in R_{+}}f(t)\), \(f(t)\) is a bounded function on \([0,+\infty )\), \(\langle f(t)\rangle =\frac{1}{t}\int _{0}^{t} f(s) \,ds\), where \(f(t)\) is an integrable function on \([0,+\infty )\).
Consider the subsystem of system (2.1) as follows:
$$ \left \{ \textstyle\begin{array}{@{}l} \left . \textstyle\begin{array}{@{}l} dc_{o}(t)=(fc_{e}(t)-(g+m)c_{o}(t))\,dt, \\ dc_{e}(t)=(-hc_{e}(t))\,dt, \end{array}\displaystyle \right \}\quad t \neq (n+l)\tau , t\neq (n+1)\tau , \\ \left . \textstyle\begin{array}{@{}l} \triangle c_{o}(t)=0, \\ \triangle c_{e}(t)=-h_{2}c_{e}(t), \end{array}\displaystyle \right \}\quad t= (n+l)\tau ,n\in Z^{+}, \\ \left . \textstyle\begin{array}{@{}l} \triangle c_{o}(t)=0, \\ \triangle c_{e}(t)=\mu , \end{array}\displaystyle \right \}\quad t= (n+1)\tau ,n\in Z^{+}. \end{array}\displaystyle \right . $$
With regard to system (3.1), we have the following equations with integrating and solving the first two equations of system (3.1) between pulses:
$$ \left \{ \textstyle\begin{array}{@{}l} c_{o}(t)= \textstyle\begin{cases} c_{o}(n\tau ^{+})e^{-(g+m)(t-n\tau )} + \frac{ fc_{e}(n\tau ^{+})(1-e^{-(h-g-m)(t-n\tau )})}{(h-g-m)}, \\ \quad t\in (n\tau , (n+l)\tau ], \\ c_{o}((n+l)\tau ^{+})e^{-(g+m)(t-(n+l)\tau )} \\ \quad {}+ \frac{ fc_{e}((n+l)\tau ^{+})(1-e^{-(h-g-m)(t-(n+l)\tau )})}{(h-g-m)}, \\ \quad t\in ((n+l)\tau ,(n+1)\tau ], \end{cases}\displaystyle \\ c_{e}(t)= \textstyle\begin{cases} c_{e}(n\tau ^{+})e^{-h(t-n\tau )},& t\in (n\tau ,(n+l) \tau ], \\ c_{e}((n+l)\tau ^{+})e^{-h(t-(n+l)\tau )},&t\in ((n+l) \tau ,(n+1)\tau ]. \end{cases}\displaystyle \end{array}\displaystyle \right . $$
The stroboscopic map of system (3.1) is obtained by the last two equations of system (3.1):
$$ \left \{ \textstyle\begin{array}{@{}l} c_{o}((n+1)\tau ^{+})=c_{o}(n\tau ^{+})e^{-(g+m)\tau } \\ \hphantom{c_{o}((n+1)\tau ^{+})={}} {}+ \frac{ c_{e}(n\tau ^{+})f(e^{-(g+m)(1-l)\tau }-e^{-(h-g-m)l\tau -(1-l)\tau })}{(h-1)} \\ \hphantom{c_{o}((n+1)\tau ^{+})={}} {}+ \frac{ (1-h_{2})c_{e}(n\tau ^{+})f(e^{-hl\tau }-e^{-(h-g-m)\tau )})}{(h-g-m)}, \\ c_{e}((n+1)\tau ^{+})=(1-h_{2})e^{-h\tau }c_{e}(n\tau ^{+})+ \mu . \end{array}\displaystyle \right . $$
We can easily have a unique fixed point \((c_{o}^{\ast },c_{e}^{\ast })\) of system (3.3) as follows:
$$ \left \{ \textstyle\begin{array}{@{}l} c_{o}^{\ast }=\frac{\mu f}{1-e^{-(g+m)\tau }}\times [ \frac{ (e^{-(g+m)(1-l)\tau }-e^{-(h-g-m)l\tau -(1-l)\tau })}{(h-g-m)(1-(1-h_{2})e^{-h\tau })} \\ \hphantom{c_{o}^{\ast }={}}{}+ \frac{ (1-h_{2})(e^{-h_{l}\tau }-e^{-(h-g-m)\tau )})}{(h-g-m)(1-(1-h_{2})e^{-h\tau })}], \\ c_{e}^{\ast }=\frac{\mu }{1-(1-h_{2})e^{-h\tau }}. \end{array}\displaystyle \right . $$
The unique fixed point \((c_{o}^{\ast },c_{e}^{\ast })\) of system (3.3) is globally asymptotically stable for the eigenvalues of the coefficient matrix of system (3.3)
$$ \begin{pmatrix} e^{-\tau } & \star \\ 0& (1-h_{2})e^{-h\tau } \end{pmatrix} $$
are less 1, there is no need for calculating ⋆.
Similar to Lemma 3.3 in reference [6], we can obtain the following lemma.
Lemma 3.1
System (3.1) has a unique positive τ-periodic solution \((\widetilde{c_{o}(t)},\widetilde{c_{e}(t)})\), which is also globally asymptotically stable, \(\widetilde{c_{o}(t)}\)and \(\widetilde{c_{e}(t)}\)are defined as follows:
$$ \left \{ \textstyle\begin{array}{@{}l} \widetilde{c_{o}(t)}= \textstyle\begin{cases} c_{o}^{\ast }e^{-(g+m)(t-n\tau )} + \frac{ fc_{e}^{\ast }(1-e^{-(h-g-m)(t-n\tau )})}{(h-g-m)}, \\ \quad t\in (n\tau , (n+l)\tau ], \\ c_{o}^{\ast \ast }e^{-(g+m)(t-(n+l)\tau )}+ \frac{ fc_{e}^{\ast \ast }(1-e^{-(h-g-m)(t-(n+l)\tau )})}{(h-g-m)}, \\ \quad t\in ((n+l)\tau ,(n+1)\tau ], \end{cases}\displaystyle \\ \widetilde{c_{e}(t)}= \textstyle\begin{cases} c_{e}^{\ast }e^{-h(t-n\tau )},& t\in (n\tau ,(n+l)\tau ], \\ c_{e}^{\ast \ast }e^{-h(t-(n+l)\tau )},&t\in ((n+l)\tau ,(n+1) \tau ], \end{cases}\displaystyle \end{array}\displaystyle \right . $$
where \(c_{o}^{\ast }\), \(c_{e}^{\ast }\)are defined as (3.4), and
$$ \left \{ \textstyle\begin{array}{@{}l} c_{o}^{\ast \ast }=c_{o}^{\ast }e^{-(g+m)l\tau }+ \frac{ fc_{e}^{\ast }(1-e^{-(h-g-m)l\tau })}{(h-g-m)}, \\ c_{e}^{\ast \ast }=(1-h_{2})e^{-hl\tau }c^{\ast }_{e}. \end{array}\displaystyle \right . $$
Remark 3.2
For any positive solution \((c_{o}(t),c_{e}(t))\) of system (3.1) with the initial value \((c_{o}(0),c_{e}(0))\in R^{+}_{2}\), we can obtain
$$ \begin{aligned}[b] &\lim_{t\rightarrow +\infty }\bigl\langle c_{o}(t)\bigr\rangle \\ &\quad =\frac{c^{\ast }_{o}(1-e^{-(g+m)l\tau })}{(g+m)\tau } \\ &\qquad {} +\frac{flc^{\ast }_{e}}{h-g-m}- \frac{fc^{\ast }_{e}(1-e^{-(h-g-m)l\tau })}{(h-g-m)^{2}\tau } \\ &\qquad {}+ \frac{c^{\ast \ast }_{o}(e^{-(g+m)l\tau }-e^{-(g+m)\tau })}{\tau } \\ &\qquad {} +\frac{f(1-l)c^{\ast \ast }_{e}}{h-g-m}- \frac{fc^{\ast \ast }_{e}(e^{-(h-g-m)l\tau }-e^{-(h-g-m)\tau })}{(h-g-m)^{2}\tau } \stackrel{\Delta }{=} \widetilde{c_{0}}, \end{aligned} $$
where \(c^{\ast }_{o}\) and \(c^{\ast }_{e}\) are defined as (3.4), and \(c^{\ast \ast }_{o}\) and \(c^{\ast \ast }_{e}\) are defined as (3.6).
For convenience, we consider the following notation:
$$ \tau _{n}=n\tau ,\qquad \tau _{n+l}=(n+l)\tau ,\qquad h_{n+l}=h_{1}. $$
Define \((x(t),y(t))\) and \((w(t),z(t))\) are the solutions of the subsystem of system (2.1), respectively:
$$ \left \{ \textstyle\begin{array}{@{}l} \left . \textstyle\begin{array}{@{}l} dx(t)=[D(x_{0}-x(t)) \\ \hphantom{dx(t)={}}{}-\frac{\beta x(t)y(t)}{k(A+x(t)+By(t))}]\,dt \\ \hphantom{dx(t)={}}{} +x(t)(\sigma _{11}+\sigma _{12}x(t))\,dB_{1}(t), \\ dy(t)=[\frac{\beta x(t)y(t)}{A+x(t)+By(t)} \\ \hphantom{dy(t)={}}{} -Dy(t)-rc_{o}(t)y(t)]\,dt \\ \hphantom{dy(t)={}}{} +y(t)(\sigma _{21}+\sigma _{22}y(t))\,dB_{2}(t), \end{array}\displaystyle \right \} \quad t\neq (n+l)\tau , \\ \left . \textstyle\begin{array}{@{}l} \triangle x(t)=0, \\ \triangle y(t)=-h_{1}y(t), \end{array}\displaystyle \right \}\quad t= (n+l)\tau ,n\in Z^{+}, \end{array}\displaystyle \right . $$
and the following SDE without impulsive perturbations:
$$ \left \{ \textstyle\begin{array}{@{}l} dw(t)=[D(x_{0}-w(t)) \\ \hphantom{dw(t)={}}{}- \frac{\beta w(t)\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z(t)}{k(A+w(t)+B\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z(t))}]\,dt \\ \hphantom{dw(t)={}}{} +w(t)(\sigma _{11}+\sigma _{12}w(t))\,dB_{1}(t), \\ dz(t)=[ \frac{ \beta w(t)z(t)}{A+w(t)+B\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z(t)} \\ \hphantom{dz(t)={}}{} -Dz(t)-rc_{o}(t)z(t)]\,dt \\ \hphantom{dz(t)={}}{} +z(t)[\sigma _{21}+\sigma _{22}\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z(t)]\,dB_{2}(t), \end{array}\displaystyle \right . $$
with the initial value \(w(0)=x(0)\) and \(z(0)=y(0)\).
The solutions \((x(t),y(t))\)of the subsystem of system (2.1) can also be expressed as follows:
$$ \left \{ \textstyle\begin{array}{@{}l} x(t)=w(t), \\ y(t)=\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z(t), \end{array}\displaystyle \right . $$
where \((w(t),z(t))\)is the solution of (3.10).
One can find that \((x(t),y(t))\) is continuous on the interval \((\tau _{n},\tau _{n+l})\), and for \(t\neq \tau _{n+l}\),
$$ \left \{ \textstyle\begin{array}{@{}l} dx(t)=dw(t) \\ \hphantom{dx(t)}=[D(x_{0}-w(t)) - \frac{ \beta w(t)\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z(t)}{k(A+w(t)+B\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z(t))}]\,dt \\ \hphantom{dx(t)={}}{} +w(t)(\sigma _{11}+\sigma _{12}w(t))\,dB_{1}(t) \\ \hphantom{dx(t)} =[D(x_{0}-w(t)) -\frac{\beta w(t)y(t)}{k(A+w(t)+By(t))}]\,dt \\ \hphantom{dx(t)={}}{} +w(t)(\sigma _{11}+\sigma _{12}w(t))\,dB_{1}(t), \\ dy(t)=\prod_{0< \tau _{n+l}< t}(1-h_{n+l})dz(t) \\ \hphantom{dy(t)} =\prod_{0< \tau _{n+l}< t}(1-h_{n+l})\{ [ \frac{\beta w(t)z(t)}{A+w(t)+B\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z(t)} \\ \hphantom{dy(t)={}}{} -Dz(t)-rc_{o}(t)z(t)]\,dt \\ \hphantom{dy(t)={}}{} +z(t)[\sigma _{21}+\sigma _{22}\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z(t)]\,dB_{2}(t) \} \\ \hphantom{dy(t)} =[ \frac{\beta w(t)\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z(t)}{A+w(t)+B\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z(t)} \\ \hphantom{dy(t)={}}{} -D\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z(t)-rc_{o}(t)\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z(t)]\,dt \\ \hphantom{dy(t)={}}{} +\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z(t)[\sigma _{21}+\sigma _{22} \prod_{0< \tau _{n+l}< t}(1-h_{n+l})z(t)]\,dB_{2}(t) \\ \hphantom{dy(t)} =[\frac{\beta w(t)y(t)}{A+w(t)+By(t)} -Dy(t)-rc_{o}(t)y(t)]\,dt \\ \hphantom{dy(t)={}}{} +y(t)[\sigma _{21}+\sigma _{22}y(t)]\,dB_{2}(t). \end{array}\displaystyle \right . $$
For every \(n\in N\), and \(\tau _{n+l}\in [0,+\infty )\),
$$ \begin{aligned}[b] y\bigl(\tau ^{+}_{n+l} \bigr)&=\lim_{t\rightarrow \tau ^{+}_{n+l}}\prod_{0< \tau _{j}< t}(1-h_{j)}z(t) \\ &=\prod_{0< \tau _{j}\leq \tau _{n+l}}(1-h_{\tau _{j})}z\bigl(\tau ^{+}_{n+l}\bigr) \\ &=(1-h_{\tau _{n+l})}\prod_{0< \tau _{j}< \tau _{n+l}}(1-h_{\tau _{j})}z( \tau _{n+l}) \\ &=(1-h_{\tau _{n+l})}y(\tau _{n+l}), \end{aligned} $$
$$ \begin{aligned}[b] y\bigl(\tau ^{-}_{n+l} \bigr)&=\lim_{t\rightarrow \tau ^{-}_{n+l}}\prod_{0< \tau _{j}< t}(1-h_{j)}z(t) \\ &=\prod_{0< \tau _{j}< \tau _{n+l}}(1-h_{\tau _{j})}z\bigl(\tau ^{-}_{n+l}\bigr) \\ &=\prod_{0< \tau _{j}< \tau _{n+l}}(1-h_{\tau _{j})}z(\tau _{n+l})=y( \tau _{n+l}). \end{aligned} $$
Assumption 3.4
There exists a bounded domain \(U\subset E_{d}\) with regular boundary, then
(\(A_{1}\)):
In the open domain U and some neighborhood thereof, the smallest eigenvalue of the diffusion matrix \(A(x)\) is bounded away from zero;
If \(x\in E_{d}\setminus U\), the mean time τ at which a path issuing from x reaches the set U is finite, and \(\sup_{x\in K}E_{x}\tau <\infty \) for every compact subset \(K\subset U_{d}\).
Assumption 3.4 is a general assumption which is the condition for Lemma 3.6.
If Assumption 3.4holds, the Markov process \(X(t)\)has a stationary distribution \(\mu (\cdot )\), and
$$ \mathbb{P} \biggl\{ \lim_{T\rightarrow \infty }\frac{1}{T} \int ^{T}_{0}f \bigl(x(t) \bigr)\,dt= \int _{E_{d}}f(x)\mu (dx) \biggr\} =1, $$
where f is an integrable function with respect to the measure μ.
The dynamics
In the following theorem, we devote ourselves to investigating system (3.10).
Theorem 4.1
If \(\frac{\beta }{A}\int ^{\infty }_{0}w\phi (w)\,dw< D+r\widetilde{c_{o}}+ \frac{\sigma ^{2}_{21}}{2}\)holds, then
$$ \lim_{t\rightarrow +\infty }z(t)=0\quad \textit{a.s.}, $$
where for \(x\in (0,+\infty )\)
$$ \begin{aligned} \phi (x)={}&Cx^{-2- \frac{2(2Dx_{0}\sigma _{12}+D\sigma _{11})}{\sigma ^{3}_{1}}}\times ( \sigma _{11}+\sigma _{12}x)^{-2+ \frac{2(2Dx_{0}\sigma _{12}+D\sigma _{11})}{\sigma ^{3}_{11}}} \\ &{} \times e^{-\frac{2}{\sigma _{11}(\sigma _{11}+\sigma _{12}x)}{( \frac{Dx_{0}}{x}+ \frac{2Dx_{0}\sigma _{12}+D\sigma _{11}}{\sigma _{11}}})}, \end{aligned} $$
and constant C satisfies that \(\int ^{\infty }_{0}\phi (x)\,dx=1\).
Constructing the following auxiliary differential equation:
$$ dW(t)=\bigl[D\bigl(x_{0}-W(t)\bigr)\bigr]\,dt+W(t) \bigl(\sigma _{11}+\sigma _{12}W(t)\bigr) \,dB_{1}(t), $$
with the initial value \(W(0)=x(0)>0\), we assume that \(W(t)\) is the solution of (4.3). Obviously, the following inequality can be obtained by the comparison theorem for stochastic differential equations:
$$ w(t)\leq W(t)\quad \mbox{a.s.} $$
We set
$$ a(w)=D\bigl(x_{0}-w(t)\bigr),\qquad \sigma (w)=w( \sigma _{11}+\sigma _{12}w),\quad w \in (0,+\infty ), $$
and compute the following indefinite integral:
$$ \begin{aligned}[b] \int \frac{a(t)}{\sigma ^{2}(t)}\,dt={}& \int \frac{D(x_{0}-t)}{t^{2}(\sigma _{11}+\sigma _{12}t)^{2}}\,dt \\ ={}&\frac{2Dx_{0}\sigma _{12}+D\sigma _{11}}{\sigma ^{3}_{11}}\ln \frac{\sigma _{11}+\sigma _{12}t}{t} \\ &{}-\frac{Dx_{0}}{\sigma _{11}t(\sigma _{11}+\sigma _{12}t)}- \frac{2Dx_{0}\sigma _{12}+D\sigma _{11}}{\sigma ^{2}_{11}t(\sigma _{11}+\sigma _{12}t)}+C. \end{aligned} $$
$$ e^{\int \frac{a(t)}{\sigma ^{2}(t)}\,dt}=e^{C}\biggl( \frac{\sigma _{11}+\sigma _{12}t}{t}\biggr)^{ \frac{2Dx_{0}\sigma _{12}+D\sigma _{11}}{\sigma ^{3}_{11}}} e^{- \frac{1}{\sigma _{11}(\sigma _{11}+\sigma _{12}t)}(\frac{Dx_{0}}{t}+ \frac{2Dx_{0}\sigma _{12}+D\sigma _{11}}{\sigma _{11}})}. $$
$$ \begin{aligned}[b] \int ^{\infty }_{0}\frac{1}{\sigma ^{2}(w)}e^{\int ^{w}_{0} \frac{2a(s)}{\sigma ^{2}(s)}\,ds} \,dw={}& \int ^{\infty }_{0} w^{-2}(\sigma _{11}+ \sigma _{12}w)^{2}\biggl( \frac{\sigma _{11}+\sigma _{12}w}{w}\biggr)^{ \frac{2(2Dx_{0}\sigma _{12}+D\sigma _{11})}{\sigma ^{3}_{11}}} \\ &{}\times e^{-\frac{2}{\sigma _{11}(\sigma _{11}+\sigma _{12}w)}( \frac{2(2Dx_{0}\sigma _{12}+D\sigma _{11})}{\sigma ^{3}_{11}})}\,dw< \infty . \end{aligned} $$
This indicates that SDE (4.3) has the ergodic property. By the ergodic theorem, we have
$$ \lim_{t\rightarrow +\infty }\frac{1}{t} \int ^{t}_{0}w(s)\,ds= \int ^{\infty }_{0}w\phi (w)\,dw \quad \mbox{a.s.} $$
Applying Itô's formula, we have
$$ \begin{aligned}[b] d\ln z(t)={}&\biggl[ \frac{\beta w(t)}{A+w(t)+B\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z(t)}-D-rc_{0}(t) \\ &{}-\frac{1}{2}\biggl(\sigma _{21}+\sigma _{22}\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z(t) \biggr)^{2}\biggr]\,dt \\ &{}+\biggl(\sigma _{21}+\sigma _{22}\prod _{0< \tau _{n+l}< t}(1-h_{n+l})z(t)\biggr)\,dB_{2}(t) \\ \leq{}& \biggl[\frac{ \beta w(t)}{A}-D-rc_{o}(t)- \frac{\sigma ^{2}_{21}}{2} \\ &{}-\sigma _{21}\sigma _{22}\prod _{0< \tau _{n+l}< t}(1-h_{n+l})z(t)- \frac{(\prod_{0< \tau _{n+l}< t}(1-h_{n+l})\sigma _{22})^{2}}{2}z^{2}(t) \biggr]\,dt \\ &{} +\biggl(\sigma _{21}+\sigma _{22}\prod _{0< \tau _{n+l}< t}(1-h_{n+l})z(t)\biggr)\,dB_{2}(t). \end{aligned} $$
Integrating with respect to t from 0 to t on both sides of (4.10), we have
$$ \begin{aligned}[b] \ln z(t)\leq{}& \frac{\beta }{A} \int ^{t}_{0}w(s)\,ds-\biggl(D+ \frac{\sigma ^{2}_{21}}{2}\biggr)t- \int ^{t}_{0}rc_{o}(s)\,ds \\ &{} - \frac{(\sigma _{22}\prod_{0< \tau _{n+l}< t}(1-h_{n+l}))^{2}}{2} \int ^{t}_{0}z^{2}(s)\,ds+ \sigma _{21}B_{2}(t)+M(t)+\ln y(0), \end{aligned} $$
where \(M(t)=\sigma _{22}\prod_{0<\tau _{n+l}<t}(1-h_{n+l})\int ^{t}_{0}z(s)\,dB_{2}(t)\) and its quadratic variation is given by
$$ \langle M,M\rangle (t)=\biggl(\sigma _{22}\prod _{0< \tau _{n+l}< t}(1-h_{n+l})\biggr)^{2} \int ^{t}_{0}z^{2}(s)\,ds. $$
According to the exponential martingales inequality, for any positive τ, α, β,
$$ P\biggl\{ \sup_{0\leq t\leq T}\biggl[M(t)- \frac{a}{2\langle M,M\rangle }(t)\biggr]>b\biggr\} \leq e^{-ab}. $$
Let \(T=k_{1}\), \(a=1\), \(b=\ln k_{1}\), then
$$ P\biggl\{ \sup_{0\leq t\leq T}\biggl[M(t)- \frac{1}{2\langle M,M\rangle }(t)\biggr]>\ln k_{1}\biggr\} \leq \frac{1}{k_{1}}. $$
There exists random \(k^{0}_{1}\in k_{1}(\omega )\) such that \(k_{1}>k^{0}_{1}\) for almost all \(\omega \in \varOmega \). We can obtain the following by the Borel–Cantelli lemma:
$$ \sup_{0\leq t\leq T}\biggl[M(t)- \frac{1}{2\langle M,M\rangle }(t)\biggr]\leq \ln k_{1}. $$
$$ \begin{aligned}[b] M(t)&\leq \ln k_{1}+ \frac{1}{2}\langle M,M\rangle (t) \\ & =\ln k_{1}+ \frac{(\sigma _{22}\prod_{0< \tau _{n+l}< t}(1-h_{n+l}))^{2}}{2} \int ^{t}_{0}z^{2}(s)\,ds, \\ &\quad \mbox{for all } 0< t< k_{1}, k_{1}>k^{0}_{1} \mbox{ a.s.} \end{aligned} $$
Considering (4.11) and (4.16), we have
$$ \ln z(t)\leq \frac{\beta }{A} \int ^{t}_{0}w(s)\,ds-\biggl(D+ \frac{\sigma ^{2}_{21}}{2}\biggr)t- \int ^{t}_{0}rc_{o}(s)\,ds+\sigma _{21}B_{2}(t)+ \ln k_{1}+\ln y(0). $$
Then, for \(0\leq k_{1}-1\leq t\leq k_{1}\), we have
$$ \begin{aligned}[b] \frac{\ln z(t)}{t}\leq{}& \frac{\beta }{At} \int ^{t}_{0}w(s)\,ds-\biggl(D+ \frac{\sigma ^{2}_{21}}{2}\biggr) \\ &{}-\frac{1}{t} \int ^{t}_{0}rc_{o}(s)\,ds+ \frac{\sigma _{21}B_{2}(t)}{t}+ \frac{\ln k_{1}}{k_{1}-1}+\frac{\ln y(0)}{t}. \end{aligned} $$
Taking the superior limit on both sides of (4.18), note that
$$ \lim_{t\rightarrow +\infty }B_{2}(t)=0 \quad \mbox{a.s.} $$
and \(t\rightarrow +\infty \Rightarrow k_{1}\rightarrow +\infty \), we have
$$ \lim_{k_{1}\rightarrow +\infty }\frac{\ln k_{1}}{k_{1}}=\lim _{k_{1} \rightarrow +\infty }\frac{1}{k_{1}}=0. $$
Then we obtain
$$ \lim \sup_{t\rightarrow +\infty }\frac{\ln z(t)}{t}\leq \frac{\beta }{A} \int ^{+\infty }_{0}w\phi (w)\,dw-\biggl(D+r \widetilde{c_{o}}+ \frac{\sigma ^{2}_{21}}{2}\biggr). $$
This implies that if \(\frac{\beta }{A}\int ^{\infty }_{0}w\phi (w)\,dw< D+r\widetilde{c_{o}}+ \frac{\sigma ^{2}_{21}}{2}\) holds, then
$$ \lim_{t\rightarrow +\infty }z(t)=0\quad \mbox{a.s.} $$
This completes the proof. □
According to Lemma 3.3 and Theorem 4.1, one can easily obtain
$$ \lim_{t\rightarrow +\infty }y(t)=0\quad \mbox{a.s.} $$
If \(\frac{\beta }{A}\int ^{\infty }_{0}w\phi (w)\,dw< D+r\widetilde{c_{o}}+ \frac{\sigma ^{2}_{21}}{2}\)holds, the distribution of \(x(t)\)converges weakly to the measure which has the density \(\pi (x)\).
For any small \(\varepsilon >0\), there exist \(t_{0}\) and a set \(\varOmega _{\varepsilon }\subset \varOmega \) such that \(\mathbb{P}(\varOmega _{\varepsilon })>1-\varepsilon \) and \(\frac{ xy}{k(A+x+By)}\leq \varepsilon x\) for \(t\geq t_{0}\) and \(\omega \in \varOmega _{\varepsilon }\). Then
$$\begin{aligned}& \bigl[ \bigl(x_{0}-x(t) \bigr)-\varepsilon x(t) \bigr]\,dt+x(t) \bigl( \sigma _{11}+\sigma _{12}x(t) \bigr) \,dB_{1}(t) \\& \quad \leq dx(t)\leq \bigl[ \bigl(x_{0}-x(t) \bigr) \bigr]\,dt+x(t) \bigl( \sigma _{11}+\sigma _{12}x(t) \bigr) \,dB_{1}(t). \end{aligned}$$
This shows that the distribution of the process \(x(t)\) converges weakly to the measure with density \(\pi (x)\). □
If \(Dx_{0}\beta >(D+\sigma _{11}^{2}+ \frac{2\sigma _{12}Dx_{0}}{\sigma _{11}})(D+rc^{u}_{o}+\frac{1}{2} \sigma _{21}^{2})\)holds, system (3.10) admits a unique stationary distribution and it has ergodic property for initial \((w(0),z(0))\in \mathbb{R}^{2}_{+}\).
$$\begin{aligned}& V^{\ast }(w,z) \\& \quad =M \biggl[-c_{1}\ln w-c_{2}l\ln z+ \frac{2c_{1}}{\theta (1-\theta )\sigma _{11}^{\theta }}(\sigma _{11}+ \sigma _{12}w)^{\theta } \biggr]+\frac{1}{p}w^{p}+\frac{1}{p}z^{p} \\& \quad \stackrel{\Delta }{=}MV_{1}^{\ast }+V_{2}^{\ast }, \end{aligned}$$
where \(\theta ,p\in (0,1)\), and M is a sufficiently large constant satisfying the following condition:
$$ M \biggl[-2Dw_{0} \biggl(\sqrt{ \frac{Dx_{0}\beta }{(D+\sigma _{11}^{2}+\frac{2\sigma _{12}Dx_{0}}{\sigma _{11}})(D+rc^{u}_{o}+\frac{1}{2}\sigma _{21}^{2})}}-1 \biggr) \biggr]+f^{u}+g^{u} \leq -2, $$
$$\begin{aligned}& c_{1}= \frac{Dw_{0}}{1+\sigma ^{2}_{11}+\frac{2\sigma _{12}Dw_{0}}{(D-\theta )\sigma _{11}}}, \\& c_{2}=\frac{Dw_{0}}{D+rc^{u}_{o}+\frac{1}{2}\sigma ^{2}_{21}}, \\& f^{u}=\sup_{w\in R_{+}} \biggl\{ Dw_{0}w^{p-1}-Dw^{p}- \frac{1-p}{2}\sigma ^{2} _{12}w^{p+2} \biggr\} , \\& g^{u}=\sup_{w\in R_{+}} \biggl\{ -Dz^{p}-rc^{l}_{o}z^{p}- \frac{1-p}{2} \biggl[ \sigma _{22}\prod _{0< \tau _{n+l}< t}(1-h_{n+l}) \biggr]^{2}y^{p+2} \biggr\} . \end{aligned}$$
It shows that
$$ \liminf_{\varepsilon \rightarrow 0, (w,z)\in R_{+}\setminus D}V^{ \ast }(w,z)=+\infty , $$
where \(D=(\varepsilon ,\frac{1}{\varepsilon })\times (\varepsilon , \frac{1}{\varepsilon })\) and \(V^{\ast }(w,z)\) is a continuous function. Therefore, \(V^{ \ast }(w,z)\) has a minimum point \((w_{0},z_{0})\) in the interior of \(\mathbb{R}^{2}_{+}\). Thus,we can define a nonnegative \(C^{2}\)-function \(V:\mathbb{R}^{2}_{+}\rightarrow \mathbb{R}_{+} \)
$$ V(w,z)=V^{\ast }(w,z)-V^{\ast }(w_{0},z_{0}). $$
We can obtain the following equation by Itô's formula:
$$ LV=LV^{\ast }=MLV^{\ast }_{1}+LV^{\ast }_{2}. $$
$$\begin{aligned} LV^{\ast }_{1} =&-\frac{c_{1}Dw_{0}}{w}+c_{1}D+ \frac{c_{1}\beta (\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z)}{K[A+w+B(\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z)]} \\ &{}+\frac{2c_{1}\sigma _{12}Dw_{0}}{(1-\theta )\sigma ^{\theta }_{11}}( \sigma _{11}+\sigma _{12}w)^{\theta -1}- \frac{2c_{1} \sigma _{12}w}{(1-\theta )\sigma ^{\theta }_{11}}(\sigma _{11}+\sigma _{12}w)^{ \theta -1} \\ &{}- \frac{2c_{1}\sigma _{12} \beta w (\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z)}{(1-\theta )\sigma ^{\theta }_{11}k[A+w+B(\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z)]}( \sigma _{11}+\sigma _{12}w)^{\theta -1} \\ &{}+\frac{c_{1}}{2}(\sigma _{11}+\sigma _{12}w)^{2}- \frac{c_{1}\sigma _{12}^{2}}{\sigma _{11}^{\theta }}(\sigma _{11}+ \sigma _{12}w)^{\theta }w^{2} \\ &{}-\frac{c_{2}\beta w}{A+w+B(\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z)}+c_{2}D+c_{2}rc_{o}(t) \\ &{}+\frac{c_{2}}{2} \biggl[\sigma _{21}+\sigma _{22} \prod_{0< \tau _{n+l}< t}(1-h_{n+l})z \biggr]^{2} \\ \leq& -\frac{c_{1}Dw_{0}}{w}-c_{1} \beta w+c_{2}\beta w +c_{1}D+ \frac{c_{1} \beta w}{kA}+ \frac{2c_{1}\sigma _{12}Dw_{0}}{(1-\theta )\sigma 11} \\ &{}+c_{1}\sigma _{11}^{2}- \frac{c_{1}}{2}(\sigma _{11}-\sigma _{12}w)^{2}+c_{2}D+c_{2}rc^{u}_{o}+ \frac{c_{2}}{2} \biggl(\sigma _{21}+\sigma _{22} \prod_{0< \tau _{n+l}< t}(1-h_{n+l})z \biggr)^{2} \\ \leq& -2\sqrt{c_{1}c_{2}Dw_{0}\beta }+c_{1} \biggl[D+\sigma ^{2}_{11}+ \frac{2\sigma _{12}Dw_{0}}{(1-\theta )\sigma _{11}} \biggr]+c_{2} \biggl[D+rc^{u}_{o}+ \frac{\sigma ^{2}_{21}}{2} \biggr] \\ &{}+ \biggl[\frac{c_{1}\beta }{kA}+c_{2}\sigma _{21} \biggl( \sigma _{22}\prod_{0< \tau _{n+l}< t}(1-h_{n+l}) \biggr) \biggr]w \\ &{}+ \frac{c_{2}(\sigma _{22}\prod_{0< \tau _{n+l}< t}(1-h_{n+l}))^{2}z^{2}}{2}+c_{2} \beta w \\ =&-2w_{0} \biggl(\sqrt{ \frac{w_{0}}{(1+\sigma _{11}^{2}+\frac{2\sigma _{12}w_{0}}{\sigma _{11}})(1+rc^{u}_{0}+\frac{1}{2}\sigma _{21}^{2})}}-1 \biggr) \\ &{}+ \biggl[\frac{c_{1}\beta }{kA}+c_{2}\sigma _{21} \biggl( \sigma _{22}\prod_{0< \tau _{n+l}< t}(1-h_{n+l}) \biggr) \biggr]z \\ &{}+ \frac{c_{2}(\sigma _{22}\prod_{0< \tau _{n+l}< t}(1-h_{n+l}))^{2}}{2}z^{2}+c_{2} \beta w, \end{aligned}$$
where \(c_{1}\) and \(c_{2}\) are such that
$$ c_{1} \biggl[D+\sigma _{11}^{2}+ \frac{2\sigma _{12}Dw_{0}}{(1-\theta )\sigma _{11}} \biggr]=c_{2} \biggl[D+rc^{u}_{o}+ \frac{1}{2}\sigma ^{2}_{21} \biggr]=Dw_{0}. $$
The function \(\frac{Dw_{0}\beta }{(D+\sigma _{11}^{2}+\frac{2\sigma _{12}Dw_{0}}{\sigma _{11}})(D+rc^{u}_{o}+\frac{1}{2}\sigma _{21}^{2})}>1\) is continuous. Choose \(\varepsilon >0\) sufficiently small such that
$$\begin{aligned} LV^{\ast }_{1} \leq& -2Dw_{0} \biggl(\sqrt{ \frac{Dw_{0}\beta }{(D+\sigma _{11}^{2}+\frac{2\sigma _{12}Dw_{0}}{\sigma _{11}})(D+rc^{u}_{o}+\frac{1}{2}\sigma _{21}^{2})}}-1 \biggr) \\ &{}+ \biggl[\frac{c_{1}\beta }{kA}+c_{2}\sigma _{21} \biggl( \sigma _{22}\prod_{0< \tau _{n+l}< t}(1-h_{n+l}) \biggr) \biggr]z \\ &{}+ \frac{c_{2}(\sigma _{22}\prod_{0< \tau _{n+l}< t}(1-h_{n+l}))^{2}z^{2}}{2}+c_{2} \beta w. \end{aligned}$$
We can also have
$$\begin{aligned} LV^{\ast }_{2} \leq& D w_{0}w^{p-1}-Dw^{p}- \frac{\beta \prod_{0< \tau _{n+l}< t}(1-h_{n+l})zw^{p}}{k(A+w+B(\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z)])} \\ &{}-\frac{1-p}{2}w^{p-2}w^{2}(\sigma _{11}+\sigma _{12}w)^{2}+ \frac{\beta w[\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z]^{p}}{A+w+B(\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z} \\ &{}-Dw^{p}-rc_{o}(t)w^{p}- \frac{1-p}{2} \biggl[\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z \biggr]^{p}\biggl( \sigma _{21}+\sigma _{22} \biggl[ \prod_{0< \tau _{n+l}< t}(1-h_{n+l})z \biggr]^{2}\biggr) \\ \leq& D w_{0}w^{p-1}-Dw^{p}- \frac{1-p}{2}\sigma _{12}^{2}w^{p+2}-Dw^{p}-rc^{l}_{o}z^{p} \\ &{}-\frac{1-p}{2}\sigma ^{2}_{22} \biggl[\prod _{0< \tau _{n+l}< t}(1-h_{n+l})z \biggr]^{p+2}+ \frac{\beta w (\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z)^{p}}{A}. \end{aligned}$$
$$\begin{aligned} LV \leq& M \biggl\{ -2Dw_{0} \biggl(\sqrt{ \frac{Dw_{0}\beta }{(D+\sigma _{11}^{2}+\frac{2\sigma _{12}Dw_{0}}{\sigma _{11}})(D+rc^{u}_{o}+\frac{1}{2}\sigma _{21}^{2})}}-1 \biggr) \\ &{}+ \biggl[\frac{c_{1}\beta }{kA}+c_{2}\sigma _{21} \sigma _{22}\prod_{0< \tau _{n+l}< t}(1-h_{n+l}) \biggr]z+ \frac{c_{2}(\sigma _{22}\prod_{0< \tau _{n+l}< t}(1-h_{n+l}))^{2}}{2}z^{2}+c_{2} \beta w \biggr\} \\ &{}+Dw_{0}x^{p-1}-Dx^{p}-rc^{l}_{o}z^{p} \\ &{}-\frac{1-p}{2}\sigma ^{2}_{22} \biggl[\prod _{0< \tau _{n+l}< t}(1-h_{n+l})z \biggr]^{p+2}+ \frac{\beta x (\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z)^{p}}{A}. \end{aligned}$$
$$\begin{aligned} H(w,z) =& M \biggl\{ -2Dw_{0} \biggl(\sqrt{ \frac{Dw_{0}\beta }{(D+\sigma _{11}^{2}+\frac{2\sigma _{12}Dw_{0}}{\sigma _{11}})(D+rc^{u}_{o}+\frac{1}{2}\sigma _{21}^{2})}}-1 \biggr) \\ &{}+ \biggl[\frac{c_{1}\beta }{kA}+c_{2}\sigma _{21} \sigma _{22}\prod_{0< \tau _{n+l}< t}(1-h_{n+l}) \biggr]z \\ &{}+ \frac{c_{2}(\sigma _{22}\prod_{0< \tau _{n+l}< t}(1-h_{n+l}))^{2}}{2}z^{2}+c_{2} \beta w \biggr\} \\ &{}+f(x)+g(z)+ \frac{ \beta x (\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z)^{p}}{A}, \end{aligned}$$
where \(f(x)=Dw_{0}w^{p-1}-Dx^{p}-\frac{1-p}{2}\sigma _{12}^{2}x^{p+2}\) and \(g(z)=-Dx^{p}-rc^{l}_{o}z^{p} -\frac{1-p}{2}\sigma ^{2}_{22}[\prod_{0< \tau _{n+l}<t}(1-h_{n+l})z]^{p+2}\). Then we can get
$$ H(w,z)\leq \textstyle\begin{cases} H(+\infty ,z)\rightarrow -\infty ,\quad \mbox{as } w\rightarrow + \infty , \\ H(w,+\infty )\rightarrow -\infty ,\quad \mbox{as } z\rightarrow + \infty , \\ M[-2Dw_{0}(\sqrt{ \frac{Dw_{0}\beta }{(D+\sigma _{11}^{2}+\frac{2\sigma _{12}Dw_{0}}{\sigma _{11}})(D+rc^{u}_{o}+\frac{1}{2}\sigma _{21}^{2})}}-1)] \\ \quad {}+f^{u}+g^{u}\leq -2,\quad \mbox{as } w\rightarrow 0^{+},z\rightarrow 0^{+}. \end{cases} $$
Therefore, there exists sufficiently small \(\varepsilon >0\) such that
$$ LV\leq -1, \quad \mbox{for any } (w,z)\in \mathbb{R}^{2}_{+} \setminus \mathbb{D}, $$
where \(\mathbb{D}=(\varepsilon ,\frac{1}{\varepsilon })\times (\varepsilon , \frac{1}{\varepsilon })\).
On the other hand, the diffusion matrix of system (3.10) is given by
$$\begin{aligned} \sum_{i,j=1}^{2}a_{ij}(w,z) \xi _{i}\xi _{j} =& \biggl( \bigl(\sigma _{11}w+ \sigma _{12}w^{2} \bigr) \xi _{1}, \biggl( \sigma _{21}z+\sigma _{22} \prod_{0< \tau _{n+l}< t}(1-h_{n+l})z^{2} \biggr) \xi _{2} \biggr) \\ &{}\times \begin{pmatrix} (\sigma _{11}w+\sigma _{12}w^{2})\xi _{1} \\ (\sigma _{21}z+\sigma _{22}\prod_{0< \tau _{n+l}< t}(1-h_{n+l})z^{2}) \xi _{2} \end{pmatrix} \\ =& \bigl(\sigma _{11}w+\sigma _{12}w^{2} \bigr)\xi ^{2}_{1}+ \biggl(\sigma _{21}z+ \sigma _{22} \prod_{0< \tau _{n+l}< t}(1-h_{n+l})z^{2} \biggr)\xi ^{2}_{2} \\ \geq& G \Vert \xi \Vert ^{2} \quad \mbox{for any } (w,z)\in D_{e}\subset \mathbb{D}^{2}_{+}, \end{aligned}$$
where \(\xi =(\xi _{1}, \xi _{2})\in \mathbb{D}^{2}_{+}\), \(G=\min_{(w,z)\in D_{e}}\{(\sigma _{11}w+\sigma _{12}w^{2})\xi ^{2}_{1}+( \sigma _{21}z+\sigma _{22}\prod_{0<\tau _{n+l}<t}(1- h_{n+l})z^{2}) \xi ^{2}_{2}\}\), and \(D_{e}=[\frac{1}{e},e]\times [\frac{1}{e},e]\).
From Theorem 4.7 in reference [9], it can be known that system (3.11) is ergodic and has a unique stationary distribution, and the distribution of the process converges weakly to the measure with density. □
From Lemma 3.3 and Theorem 4.4, we can easily know that if \(Dx_{0}\beta >(D+\sigma _{11}^{2}+ \frac{2\sigma _{12}Dx_{0}}{\sigma _{11}})(D+rc^{u}_{o}+\frac{1}{2} \sigma _{21}^{2})\) holds, system (2.1) admits a unique stationary distribution and it has ergodic property for initial \((x(0),y(0))\in \mathbb{R}^{2}_{+}\).
If it is assumed that \(x(0)=0.8\), \(y(0)=0.3\), \(c_{o}(0)=0.9\), \(c_{e}(0)=0.9\), \(D=1\), \(\beta =0.1\), \(k=0.1 A=0.05\), \(B=0.002\), \(r=0.1\), \(f=0.1\), \(g=0.5\), \(m=0.5\), \(h=0.1\), \(h_{1}=0.02\), \(h_{2}=0.2\), \(\sigma _{11}=0.01\), \(\sigma _{12}=0.01\), \(\sigma _{21}=0.01\), \(\sigma _{22}=0.01\), \(l=0.25\), \(\tau =4\), when \(\mu =0.1\), the microorganism \(y(t)\) will survive (it can be seen in (a) of Fig. 1), when \(\mu =0.9\), the microorganism \(y(t)\) will be extinct(it can be seen in (b) of Fig. 1). From Theorem 4.1, Remark 4.5, and the computer simulations in Fig. 1, we conjecture that there must exist a threshold \(\mu ^{\ast }\), if \(\mu >\mu ^{\ast }\), the microorganism \(y(t)\) will be extinct, if \(\mu <\mu ^{\ast }\), the microorganism \(y(t)\) will survive. If it is assumed that \(x(0)=0.8\), \(y(0)=0.3\), \(c_{o}(0)=0.9\), \(c_{e}(0)=0.9\), \(D=1\), \(\beta =0.1\), \(k=0.1\), \(A=0.05\), \(B=0.002\), \(r=0.1\), \(f=0.1\), \(g=0.5\), \(m=0.5\), \(h=0.1\), \(\mu =0.1\), \(h_{2}=0.2\), \(\sigma _{11}=0.01\), \(\sigma _{12}=0.01\), \(\sigma _{21}=0.01\), \(\sigma _{22}=0.01\), \(l=0.25\), \(\tau =4\), when \(h_{1}=0.01\), the microorganism \(y(t)\) will survive (it can be seen in (c) of Fig. 2), when \(h_{1}=0.1\) the microorganism \(y(t)\) will be extinct(it can be seen in (d) of Fig. 2). From Theorem 4.1, Remark 4.5, and the computer simulations in Fig. 2, we conjecture that there must exist a threshold \(h_{2}^{\ast }\), if \(h_{2}< h_{2}^{\ast }\), the microorganism \(y(t)\) will be extinct, if \(h_{2}>h_{2}^{\ast }\), the microorganism \(y(t)\) will survive. If it is assumed that \(x(0)=0.8\), \(y(0)=0.3\), \(c_{o}(0)=0.9\), \(c_{e}(0)=0.9\), \(D=1\), \(\beta =0.1\), \(k=0.1\), \(A=0.05\), \(B=0.002\), \(r=0.1\), \(f=0.1\), \(G=0.5\), \(M=0.5\), \(h=0.1\), \(\mu =0.1\), \(h_{1}=0.07\), \(\sigma _{11}=0.01\), \(\sigma _{12}=0.01\), \(\sigma _{21}=0.01\), \(\sigma _{22}=0.01\), \(l=0.25\), \(\tau =4\), when \(h_{2}=0.9\), the microorganism \(y(t)\) will survive (it can be seen in (e) of Fig. 3), when \(h_{2}=0.1\), the microorganism \(y(t)\) will be extinct (it can be seen in (f) of Fig. 3). From Theorem 4.1, Remark 4.5, and the computer simulations in Fig. 3, we conjecture that there must exist a threshold \(h_{1}^{\ast }\), if \(h_{1}< h_{1}^{\ast }\), the microorganism \(y(t)\) will be extinct, if \(h_{1}>h_{1}^{\ast }\), the microorganism \(y(t)\) will survive. If it is assumed that \(x(0)=0.8\), \(y(0)=0.3\), \(c_{o}(0)=0.9\), \(c_{e}(0)=0.9\), \(D=1\), \(\beta =0.1\), \(k=0.1\), \(A=0.05\), \(B=0.002\), \(r=0.1\), \(f=0.1\), \(g=0.5\), \(m=0.5\), \(h=0.1\), \(\mu =0.1\), \(h_{1}=0.02\), \(h_{2}=0.1\), \(\sigma _{12}=0.1\), \(\sigma _{21}=0.05\), \(\sigma _{22}=0.05\), \(l=0.25\), \(\tau =4\), and when \(\sigma _{11}=0.3\), the microorganism \(y(t)\) will survive (it can be seen in (g) of Fig. 4), when \(\sigma _{11}=0.6\), the microorganism \(y(t)\) will be extinct(it can be seen in (h) of Fig. 4). From Theorem 4.1, Remark 4.5, and the computer simulations in Fig. 4, we conjecture that there must exist a threshold \(\sigma _{11}^{\ast }\), if \(\sigma _{11}<\sigma _{11}^{\ast }\), the microorganism \(y(t)\) will be extinct, if \(\sigma _{11}>\sigma _{11}^{\ast }\), the microorganism \(y(t)\) will survive.
Threshold analysis of parameter μ in system (2.1) with \(x(0)=0.8\), \(y(0)=0.3\), \(c_{o}(0)=0.9\), \(c_{e}(0)=0.9\), \(D=1\), \(\beta =0.1\), \(k=0.1\), \(A=0.05\), \(B=0.002\), \(r=0.1\), \(f=0.1\), \(g=0.4\), \(m=0.5\), \(h=0.1\), \(h_{1}=0.02\), \(h_{2}=0.2\), \(\sigma _{11}=0.01\), \(\sigma _{12}=0.01\), \(\sigma _{21}=0.01\), \(\sigma _{22}=0.01\), \(l=0.25\), \(\tau =4\), (a): \(y(t)\) survival with parameter \(\mu =0.1\); (b): \(y(t)\) extinction with parameter \(\mu =0.9\)
Threshold analysis of parameter \(h_{1}\) in system (2.1) with \(x(0)=0.8\), \(y(0)=0.3\), \(c_{o}(0)=0.9\), \(c_{e}(0)=0.9\), \(D=1\), \(\beta =0.1\), \(k=0.1\), \(A=0.05\), \(B=0.002\), \(r=0.1\), \(f=0.1\), \(g=0.5\), \(m=0.5\), \(h=0.1\), \(\mu =0.1\), \(h_{2}=0.2\), \(\sigma _{11}=0.01\), \(\sigma _{12}=0.01\), \(\sigma _{21}=0.01\), \(\sigma _{22}=0.01\), \(l=0.25\), \(\tau =4\), (c): \(y(t)\) survival with parameter \(h_{1}=0.01\); (d): \(y(t)\) extinction with parameter \(h_{1}=0.1\)
Threshold analysis of parameter \(h_{2}\) in system (2.1) with \(x(0)=0.8\), \(y(0)=0.3\), \(c_{o}(0)=0.9\), \(c_{e}(0)=0.9\), \(D=1\), \(\beta =0.1\), \(k=0.1\), \(A=0.05\), \(B=0.002\), \(r=0.1\), \(f=0.1\), \(g=0.5\), \(m=0.5\), \(h=0.1\), \(\mu =0.1\), \(h_{1}=0.07\), \(\sigma _{11}=0.01\), \(\sigma _{12}=0.01\), \(\sigma _{21}=0.01\), \(\sigma _{22}=0.01\), \(l=0.25\), \(\tau =4\), (e): \(y(t)\) survival with parameter \(h_{2}=0.9\); (f): \(y(t)\) extinction with parameter \(h_{2}=0.1\)
Threshold analysis of parameter \(\sigma _{11}\) in system (2.1) with \(x(0)=0.3\), \(y(0)=0.3\), \(c_{o}(0)=0.3\), \(c_{e}(0)=0.3\), \(D=1\), \(\beta =0.1\), \(k=0.1\), \(A=1.5\), \(B=1\), \(r=0.2\), \(f=0.1\), \(g=0.5\), \(m=0.5\), \(h=0.1\), \(\mu =0.2\), \(h_{1}=0.3\), \(h_{2}=0.2\), \(\sigma _{12}=0.01\), \(\sigma _{21}=0.01\), \(\sigma _{22}=0.01\), \(l=0.25\), \(\tau =4\). (g): \(y(t)\) survival with parameter \(\sigma _{11}=0.03\); (h): \(y(t)\) extinction with parameter \(\sigma _{11}=0.4\)
In this work, we consider a stochastic eutrophication-chemostat model with impulsive dredging and pulse inputting on environmental toxicant. The sufficient condition for the extinction of microorganisms is obtained. The sufficient condition for the investigated system with unique ergodic stationary distribution is also obtained by the Lyapunov functions method. The results of mathematical analysis and numerical analysis show that the stochastic noise, impulsive diffusion, and pulse input on environmental toxicant play important roles in the extinction and survival of the microorganisms. These results indicate the effective and reliable controlling strategy for water resource management with eutrophication.
Smith, H., Waltman, P.: The Theory of the Chemostat: Dynamics of Microbial Competition. Cambridge University Press, Cambridge (1995)
Li, Z., Chen, L., Liu, Z.: A chemostat model with variable yield and impulsive state feedback control is considered. Appl. Math. Model. 36, 1255–1266 (2012)
Zhang, X., Yuan, R.: Sufficient and necessary conditions for stochastic near-optimal controls: a stochastic chemostat model with non-zero cost inhibiting. Appl. Math. Model. 78, 601–626 (2020)
Lakshmikantham, V.: Theory of Impulsive Differential Equations. World Scientific, Singapore (1989)
Jiao, J., et al.: Dynamics of a stage-structured predator–prey model with prey impulsively diffusing between two patches. Nonlinear Anal., Real World Appl. 11, 2748–2756 (2010)
Jiao, J., et al.: Dynamical analysis of a five-dimensioned chemostat model with impulsive diffusion and pulse input environmental toxicant. Chaos Solitons Fractals 44, 17–27 (2011)
Lv, X., Meng, X., Wang, X.: Extinction and stationary distribution of an impulsive stochastic chemostat model with nonlinear perturbation. Chaos Solitons Fractals 110, 273–279 (2018)
Zhou, P.: Reservoir desilting is an effective way to alleviate water crisis. Energy Conserv. Environ. Prot. 4, 2014S3919 (2014)
Mao, X.: Stochastic Differential Equations and Applications. Elsevier, Amsterdam (2007)
Wang, K.: Stochastic Biomathematical Models. Science Press, Beijing (2010) (in Chinese)
Liu, Q., Jiang, D.: Stationary distribution and extinction of a stochastic SIR model with nonlinear perturbation. Appl. Math. Lett. 73, 8–15 (2017)
Data sharing not applicable to this article as no datasets were generated or analysed during the current study.
This paper is supported by the National Natural Science Foundation of China (11761019, 11361014), the Science Technology Foundation of Guizhou Education Department (20175736-001, 2008038), the Science Technology Foundation of Guizhou (2010J2130),the Project of High Level Creative Talents in Guizhou Province (No. 20164035), and the Joint Project of Department of Commerce and Guizhou University of Finance and Economics (No. 2016SWBZD18).
School of Mathematics and Statistics, Guizhou University of Finance and Economics, Guiyang, 550025, P.R. China
Jianjun Jiao
Key Laboratory for Information System of Mountainous Area and Protection of Ecological Environment of Guizhou Province, Guizhou Normal University, Guiyang, 550001, P.R. China
Qiuhua Li
All authors contributed equally to the writing of this paper. All authors read and approved the final manuscript.
Correspondence to Jianjun Jiao.
The authors declare that there is no conflict of interests.
Jiao, J., Li, Q. Dynamics of a stochastic eutrophication-chemostat model with impulsive dredging and pulse inputting on environmental toxicant. Adv Differ Equ 2020, 447 (2020). https://doi.org/10.1186/s13662-020-02905-5
Stochastic eutrophication-chemostat model
Impulsive dredging
Pulse input