2 The stochastic differential equation system

Consider the dynamics model describing water quality in terms of BOD and DO interaction with rates of reaction kinetics of nonlinear BOD degradation and DO depletion of Michaelis–Menten type [Reference Pimpunchat12], through a first-order system of differential equations, where the BOD and DO concentrations are denoted by $x_1(t)$ and $x_2(t)$ , respectively, and t is the time. Then, the BOD–DO interaction system is given by

(2.1) $$ \begin{align} \begin{aligned} \frac{d x_1(t)}{dt}= & -k_1\frac{x_2(t)}{k+x_2(t)}x_1(t)+q, \\ \frac{d x_2(t)}{dt}= &-k_2\frac{x_2(t)}{k+x_2(t)}x_1(t)+\alpha(x_{2_c}-x_2(t)), \end{aligned} \end{align} $$

where $k_1$ is the BOD decay rate, $k_2$ is the DO deoxygenation rate, k is the half-saturated oxygen demand concentration, $\alpha $ is the reaeration rate, $x_{2_c}$ is the oxygen saturation constant and q is the source term for BOD.

The overall dynamics of BOD and DO (that is, the model in equation (2.1)) under the influence of random fluctuations can be described by the stochastic differential equation system:

(2.2) $$ \begin{align} \begin{aligned} d x_1(t)&= \bigg(-k_1\frac{x_2(t)}{k+x_2(t)}x_1(t)+q\bigg)\,dt+\sigma_1 x_1 \, dW_1(t), \\ d x_2(t)&= \bigg(-k_2\frac{x_2(t)}{k+x_2(t)}x_1(t)+\alpha(x_{2_c}-x_2(t))\bigg)\,dt+\sigma_2 x_2 \, dW_2(t), \end{aligned} \end{align} $$

where $W_1(t)$ and $W_2(t)$ are independent Gaussian white noises defined on complete probability space $(\Omega , \mathcal {F},(\mathcal {F}_t)_{t\geq 0}, \mathbb {P})$ , and nonnegative constants $\sigma _1$ and $\sigma _2$ represent the noise intensities. For convenience, it is assumed that the variables $x_1$ and $x_2$ are dimensionless, that is, the equations describe changes in relative concentration size.

3 Analysis for stochastic persistence

To study the stochastic persistence and the impacts of random noise, we first discuss the dynamics of the system in equation (2.2) in the absence of noise by determining the stability of the solution. Therefore, we examine the equilibria and perform a phase plane analysis.

For the equilibrium points, the steady-state solution is considered. The steady-state equilibrium point $E_s=(x_{1_s}, x_{2_s})$ is given by

$$ \begin{align*} E_s=\bigg(\bigg(1+\frac{k}{x_{2_s}}\bigg)\frac{q}{k_1}, x_{2_c}-\frac{k_2 q}{\alpha k_1}\bigg),\, k_1,\, \alpha,\, x_{2_s}\neq 0. \end{align*} $$

This interior equilibrium point exists whenever $x_{2_c}>k_2q/\alpha k_1$ and $q\neq 0$ . To analyse the stability of this equilibrium point, we determine the characteristic equation of the associated Jacobian matrix of equation (2.1) given by

$$ \begin{align*} J(x_1, x_2)= \begin{pmatrix} -\dfrac{k_1 x_2}{k+x_2} & -\dfrac{k_1 k x_1 x_2}{(k+x_2)^2} \\[.3cm] -\dfrac{k_2 x_2}{k+x_2} & -\dfrac{k_2 k x_1 x_2}{(k+x_2)^2}-\alpha \end{pmatrix}. \end{align*} $$

Using the Routh–Hurwitz criterion [Reference Chang and Chen2], which guarantees that all the roots of the characteristic equation $\lambda ^2+a\lambda +b =0$ have negative real part if and only if $a>0$ and $b>0$ , and considering that the model parameters are positive, the coefficients of the characteristic equation of the Jacobian matrix $J(x_{1_s}, x_{2_s})$ are

$$ \begin{align*} \begin{split} a = & \frac{k_1 x_{2_s}}{k+x_{2_s}} +\frac{k_2 k x_{1_s} x_{2-s}}{(k+x_{2_s})^2}+\alpha>0 , \\ b = & \frac{\alpha k_1 x_{2_s}}{k+x_{2_s}}>0. \end{split} \end{align*} $$

Hence, $E_s$ is a stable equilibrium point. To verify that the model has a stable interior equilibrium $E_s$ , we take parameter values as in [Reference Guaca and Poletti6]: $k_1=0.15, k_2=0.5, q=0.2, x_{2_c}=9.5, \alpha =0.35, k=3$ . The resulting phase portrait is shown in Figure 1.

Figure 1 This figure shows the phase portrait of equation (2.1) for given parameter values.

Then, to study the stochastic persistence and the impact of random noise, the model system in equation (2.2) will be used.

3.1 Analytical results

In this subsection, we present some analytical results concerning the stochastic persistence for the model system in equation (2.2). For this, we follow the method of Gray et al. [Reference Gray, Greenhalgh, Hu, Mao and Pan5].

Theorem 1. For any given initial values $x_i(0)=x_{i_0}, i=1,2$ , the solution of the model in equation (2.2) obeys

(3.1) $$ \begin{align} \limsup_{t\rightarrow \infty} x_i^{-1}(t)\geq \xi_i \quad a.s., \end{align} $$

and

(3.2) $$ \begin{align} \liminf_{t\rightarrow \infty} x_i^{-1}(t)\leq \xi_i \quad a.s., \end{align} $$

where $\xi _1$ and $\xi _2$ are solutions of the following equations:

$$ \begin{align*} \begin{aligned} f_1(\xi_1)&= \frac{q}{\xi_1}-\frac{k_1 x_{2_s}}{k}-\frac{\sigma_1^2}{2} =0, \\ f_2(\xi_2)&=\frac{\alpha x_{2_c}}{\xi_2}- \frac{\alpha k+k_2 x_{1_s}}{k}- \frac{\sigma_2^2}{2}=0, \end{aligned} \end{align*} $$

respectively. (Here a.s. means almost surely.) That is, $x_i^{-1}(t)$ will rise to or above the level $\xi _i$ infinitely often with probability one.

Proof. We consider approximations for $dx_1(t)$ and $dx_2(t)$ in equation (2.2) and prove the assertion in equation (3.1). If it is not true, then there is a sufficiently small $\epsilon \in (0,1)$ such that

$$ \begin{align*} \mathbb{P}(\Omega_1)>\epsilon, \end{align*} $$

where $\Omega _1=\{\limsup _{t\rightarrow \infty } x_i^{-1}(t)\leq \xi _i-2\epsilon \}$ . Hence, for every $\omega \in \Omega _1$ , there is a ${T=T(\omega )>0}$ such that

(3.3) $$ \begin{align} x_i^{-1}(t,\omega) \leq \xi_i-\epsilon \quad \text{whenever } \, t\geq T(\omega). \end{align} $$

It follows from equation (3.3) that

(3.4) $$ \begin{align} f(x_i(t,\omega))\geq f(\xi_i-\epsilon) \quad \text{whenever } t\geq T(\omega). \end{align} $$

Moreover, by the large number theorem for martingales, there is a $\Omega _2\subset \Omega $ with $\mathbb {P}(\Omega _2)=1$ such that for every $\omega \in \Omega _2$ ,

(3.5) $$ \begin{align} \lim_{t\rightarrow \infty}\frac{1}{t}\int_0^t \sigma_i \, dW_i(s,\omega)=0. \end{align} $$

Now, fix any $\omega \in \Omega _1\cap \Omega _2$ . It then follows from the Ito formula [Reference Øksendal10] and equation (3.4) that, for $t\geq T(\omega )$ ,

$$ \begin{align*} \log (x_i(t,\omega))\geq \log(x_{i_0})+\int_0^{T(\omega)} f_i(x_i(s,\omega))\, ds +f_i (\xi_i-\epsilon)(t-T(\omega))+\int_0^t \sigma_i \,dW_i(s,\omega). \end{align*} $$

This yields

$$ \begin{align*} \liminf_{t\rightarrow \infty} \frac{1}{t} \log (x_i(t,\omega))\geq f_i(\xi_i-\epsilon)>0, \end{align*} $$

whence

$$ \begin{align*} \lim _{t\rightarrow \infty} x_i^{-1}(t,\omega)=\infty. \end{align*} $$

However, this contradicts equation (3.3). We therefore must have the desired assertion in equation (3.1).

Next we prove the assertion in equation (3.2). If it is not true, then there is a sufficiently small $\delta \in (0,1)$ such that

$$ \begin{align*} \mathbb{P}(\Omega_3)>\delta, \end{align*} $$

where $\Omega _3=\{\liminf _{t\rightarrow \infty } x_i^{-1}(t) \geq \xi _i+2\delta \}$ . Hence, for every $\omega \in \Omega _3$ , there is a ${\tau =\tau (\omega )>0}$ such that

(3.6) $$ \begin{align} x_i^{-1}(t,\omega)\geq \xi_i+\delta \quad \text{whenever } \, t\geq \tau(\omega). \end{align} $$

Now, fix any $\omega \in \Omega _3\cap \Omega _2$ . It then follows from the Ito formula that, for $t\geq \tau (\omega )$ ,

$$ \begin{align*} \log (x_i(t,\omega))\leq \log(x_{i_0})+\int_0^{\tau(\omega)} f_i(x_i(s,\omega))\, ds +f_i (\xi_i+\delta)(t-\tau(\omega))+\int_0^t \sigma_i \, dW_i(s,\omega). \end{align*} $$

This, together with equation (3.5), yields

$$ \begin{align*} \limsup_{t\rightarrow \infty} \frac{1}{t} \log (x_i(t,\omega))\leq f_i(\xi_i+\delta)<0, \end{align*} $$

whence

$$ \begin{align*} \lim _{t\rightarrow \infty} x_i^{-1}(t,\omega)=0. \end{align*} $$

However, this contradicts equation (3.6). We therefore obtain the desired assertion in equation (3.2). Hence, the proof is complete.

Moreover, we have that $x_i$ converges to a unique stationary distribution with the probability density given by

$$ \begin{align*} P_i(x_i)=\frac{b_i^{a_i}}{\Gamma(a_i)}x^{-(a_i+1)} e^{-b_i/x_i},\quad x_i>0,\quad i=1,2, \end{align*} $$

where

$$ \begin{align*} a_1&=\frac{2 \beta_1+\sigma_1^2}{\sigma_1^2},\quad b_1=\frac{2 q}{\sigma_1^2}, \quad \beta_1=\frac{k_1x_{2_s}}{k}, \\ a_2&=\frac{2 \beta_2+\sigma_2^2}{\sigma_2^2},\quad b_2=\frac{2 \alpha x_{2_c}}{\sigma_2^2},\quad \beta_2=\frac{k\alpha+k_2 x_{1_s}}{k}, \end{align*} $$

and $\Gamma (.)$ is the Gamma function, $\Gamma (x)=\int _0^\infty t^{x-1} e^{-t}\,dt$ , and we have

$$ \begin{align*} \lim_{t\rightarrow \infty} \frac{1}{t} \int_0^t x_i(s) \, ds\approx \zeta_i,\quad i=1,2, \quad \text{a.s.}, \end{align*} $$

respectively, with $\zeta _1=kq/k_1 x_{2_s}$ and $\zeta _2=k \alpha x_{2_c}/(k\alpha +k_2 x_{1_s})$ .

3.2 Computer simulations

Here, we present some computer simulations by applying the Milstein method [Reference Higham7, Reference Kloeden and Platen9] to illustrate the behaviour of the model system in equation (2.2) and support the analytical results. The discrete equations are

(3.7) $$ \begin{align} \begin{aligned} x_{1_{j+1}} & =x_{1_j}+\bigg(-k_1\frac{x_{2_j}}{k+x_{2_j}}x_{1_j}+q \bigg)\Delta t \\ &\quad +\sigma_1 x_{1_j}\sqrt{\Delta t}W_{1_j}+\frac{\sigma_1^2}{2}x_{1_j}(W_{1_j}^2-1)\Delta t, \\ x_{2_{j+1}} & =x_{2_j}+\bigg(-k_2\frac{x_{2_j}}{k+x_{2_j}}x_{1_j}+\alpha (x_{2_c}-x_{2_j}) \bigg)\Delta t \\ &\quad +\sigma_2 x_{2_j} \sqrt{\Delta t}W_{2_j}+\frac{\sigma_2^2}{2} x_{2_j}(W_{2_j}^2-1)\Delta t, \end{aligned} \end{align} $$

where $\Delta t$ is time increment, and $W_{1_j}$ and $W_{2_j}\,(j=1,2,\ldots , n)$ are the independent Gaussian random variables.

We solve equation (3.7) with the time step $10^{-2}$ for computation of the stochastic solution. We consider the influence of noise for the same parameter values ${k_1=0.15, k_2=0.5, q=0.2, x_{2_c}=9.5, \alpha =0.35, k=3}$ , where the system in equation (2.1) possesses a stable equilibrium $E_s=(x_{1_s}, x_{2_s})=(1.89,7.52)$ . Resulting solutions for initial conditions $(x_1(0),x_2(0))=(12,6)$ are presented in Figures 2 and 3. In Figure 3, we show a plot of the sample paths of the solution $(x_1(t), x_2(t))$ for noise intensities $\sigma _1=\sigma _2= 0.1$ . Figure 3 also shows sample paths of $x_1(t)$ and $x_2(t)$ for large noise intensities $\sigma _1=\sigma _2= 0.4$ . We note that we obtain the same behaviour for different initial conditions. These figures clearly illustrate persistence and the effect of changing the noise intensities. As we can observe, for small noise intensities, the system shows small-amplitude noisy oscillations while for larger noise intensities, the system exhibits large-amplitude oscillations of complex form.

Figure 2 Sample paths of the system in equation (2.2) with small values of the noise intensities represented by red lines and its deterministic case represented by blue lines. Here, $\sigma _1=\sigma _2= 0.1$ . (Colour available online.)

Figure 3 Sample paths of the system in equation (2.2) with large values of the noise intensities and its deterministic case. Here, $\sigma _1=\sigma _2= 0.4$ . (Colour available online.)

The transition from small- to large-amplitude stochastic oscillations with increasing noise are accompanied by the significant changes of the form of the probability density. In Figures 4, 5, 6 and 7, we show the approximate probability density functions of the stationary distribution obtained by computer simulation of sample paths of $x_1(t)$ and $x_2(t)$ of equation (3.7) in the case of keeping the parameters the same but with different noise intensities $\sigma _1=\sigma _2= 0.1, 0.2, 0.3$ and $0.4$ , respectively. The simulations were run for $2$ million iterations with step size $10^{-2}$ . From these figures, one can see the changes of probability densities under increasing values of the noise intensities $\sigma _1$ and $\sigma _2$ . The probability densities reveal a maximum indicating the most probable concentration size and influence by the noise intensity. This maximum of the stationary distribution depends on the noise intensity. Clearly, the maximum is shifted for large noise intensities. These computer simulations of the solution also show that for large noise intensities, the distribution of the solution is skewed in the sense that the distribution of the solution becomes more symmetric for small noise intensities.

Figure 4 Plot of the probability density functions of the stationary probability distributions of the variables $x_1(t)$ and $x_2(t)$ for small noise intensities, corresponding to Figure 2, $\sigma _1=\sigma _2=0.1$ .

Figure 5 Plot of the probability density functions of the stationary probability distributions of the variables $x_1(t)$ and $x_2(t)$ for increasing the noise intensities $\sigma _1=\sigma _2=0.2$ .

Figure 6 Plot of the probability density functions of the stationary probability distributions of the variables $x_1(t)$ and $x_2(t)$ for increasing the noise intensities $\sigma _1=\sigma _2=0.3$ .

Figure 7 Plot of the probability density functions of the stationary probability distributions of the variables $x_1(t)$ and $x_2(t)$ for large noise intensities, corresponding to Figure 3, $\sigma _1=\sigma _2=0.4$ .

We note also that values of means of these probability distributions of $x_1(t)$ and $x_2(t)$ are $0.5266$ and $5.0379$ , respectively, which are useful in assessing the behaviour of the solution process.