Stationary distribution and probability density function analysis of a stochastic Microcystins degradation model with distributed delay

1.   Introduction

Microcystins (MCs) represent a highly prevalent form of cyanotoxins derived from harmful cyanobacterial blooms. These toxins can introduce toxicity into vital organs including liver, kidney, heart etc., posing major threats to human health ^[1]. Biodegradation has emerged as an economically and environmentally favorable approach for MCs degradation, garnering considerable research attention ^[2,3]. Recently, Song et al. ^[4] proposed a degradation model as follows:

where $x_1(t), x_2(t)$ and $x_3(t)$ denote the concentrations of MCs, MC-degrading bacteria and degrading enzymes produced by MC-degrading bacteria at time $t$, respectively. $a_{10} > 0$ is the input concentration of MCs; $D$ is the washout rate; $a_{12}\geq 0$ is the consumption rate of MCs; $a_{21}\geq 0$ is the maximal growth rate of MC-degrading bacteria; $a_{20}\geq 0$ is the consumption rate of MC-degrading bacteria; $a_{30}\geq 0$ is the maximal growth rate of degrading enzymes; $a_{13}\geq 0$ and $a_{31}\geq 0$ are the degradation rate of MCs and the consumption rate of degrading enzymes, respectively; $d_i > 0\, (i = 1, 2, 3)$ represents the death rates of MCs, MC-degrading bacteria and degrading enzymes, as distinguished by, i = 1, 2, 3, respectively. Re-scaling model (1.1) by

the authors of ^[5] obtained the following dimensionless model:

Time delays cannot be ignored in models for MCs degradation, especially the delay between MC-degrading bacteria and the production of new degrading enzymes particles. There has been extensive discussion of the presence of such delays and their impact on microbial population dynamics (see ^[6,7,8]). Song et al. ^[9] considered the following model with time delay:

Here the constant $\tau \geq 0$ represents the time consumed by MC-degrading bacteria to convert MCs, once consumed, into viable biomass and $e^{-\delta \tau}$ represents the survival probability of degrading bacteria during the conversion process. However, a more realistic delayed MC-degrading model should be an equation with distributed delay, which describes the cumulative effect of the past history of a system. In this paper, we will mainly consider the following degradation model with distributed delay

According to MacDonald ^[10], the distributed delay function $f(t)$ can be described by the Gamma distribution

as a kernel, where $\alpha > 0$ and $n$ is a nonnegative integer. In this article, we mainly consider the weak kernel, that is $K(t) = \alpha e^{-\alpha t}\; \; (n = 0)$ with $\alpha > 0.$

On the other hand, in the real world, population systems are inevitably influenced by the environmental noise (see ^{[11,12,13,14,15]}). According to May ^[16], the birth rate, carrying capacity, competition coefficient, and other parameters should exhibit random fluctuations as a result of environmental noise. In order to investigate the impact of environmental noise on population dynamics, various scholars have added white noise into deterministic systems (see ^{[17,18,19,20,21]}).

Therefore, based on the model (1.2), we further consider the following stochastic MCs degradation model with distribution delay

where $B_i(t)$ denotes mutually independent standard Brownian motions and $\sigma_i^2 > 0$ denotes the intensity of white noises $i = 1, 2, 3.$ Set

Then, by using the linear chain technique, system (1.3) is transformed into the following equivalent system:

In this paper, we propose a Microcystins degradation model with distributed delay and stochastic perturbation. As is known to us, there is no work concerning the existence of stationary distribution in system (1.3). The main difficulty is to find appropriate Lyapunov functions. It should be emphasized that there are relatively few studies focusing on the explicit expression of probability density functions due to the difficulty of solving the corresponding Fokker-Planck equation.

The paper is organized as follows. In Section 2, we prove that the stochastic Microcystins degradation system (1.4) with distributed delay has a unique global positive solution. In Section 3, a unique global solution of system (1.4) is proved to be stationary by constructing several suitable Lyapunov functions. Numerical simulations are presented to illustrate our analytical results in Section 4.

2.   The positivity of the solution

To study the dynamical behavior of Microcystins degradation model, we first need to consider whether the solution is global and positive. In this section, we will prove the existence and uniqueness of a global solution to system (1.4) with any positive initial value.

Theorem 2.1: For any initial value $(x(0), y(0), z(0), w(0)) \in \mathbb{R }^4_+$, there is a unique solution of the system (1.4) on $t\ge 0$ and the solution will remain in $\mathbb{R }^4_+$ with probability one.

Proof. The coefficients of system (1.3) satisfy the local Lipschitz condition, so there is a unique local solution $(x(t), y(t), z(t), w(t))$ on $t\in [0, \tau_e)$, where $\tau_e$ is the explosion time. Now we shall show that this solution is global, i.e., prove $\tau_e = \infty$ a.s. Let $k_0\geq 1$ be sufficiently large so that $x(0), y(0), z(0)\; \text {and}\; w(0)$ all lie within the interval $[\frac{1}{k_0}, k_0]$. For each integer $k\geq k_0$, define the stopping time:

where throughout this paper we set $\inf \emptyset = \infty.$ Evidently, $\tau_k$ is increasing as $k\rightarrow \infty.$ Denote $\tau_{\infty} = \lim_{k\rightarrow \infty}\tau_k,$ where $\tau_{\infty}\leq \tau_{e}$ a.s. If we can prove that $\tau_{\infty} = \infty$ a.s., then $\tau_{e} = \infty$ and $(x(t), y(t), z(t), w(t))\in \mathbb{R }_+^4$ a.s. for all $t\geq 0.$ In other words, to complete the proof all we need to prove is that $\tau_{\infty} = \infty$ a.s. If this statement is false, then there is a pair of constants $T > 0$ and $\varepsilon \in (0, 1)$ such that

Hence there is an integer $k_1\geq k_0$ such that

Define the nonnegative $C^2-$ Lyapunov function $V_0$ as follows

where $a, b, c, d$ are positive constants which will be determined later.

Based on the basic inequality $u-1-\ln u \geq 0$, we have

Applying Itô's formula, we have

Choose $b = \frac{{a_2c_1}}{ \alpha D_3}, c = \frac{a_2}{D_3}, d = \frac{{b+D_1}}{c}$ and $a = \min\{\frac{a_1}{b_1}, \frac{{a_1D_3+c_1a_2}}{D_2D_3}\}$ such that $ab_1-a_1\leq 0\, \, \text{and } \, a_1-aD_2+\alpha b\leq 0,$ then, we have

where $K_1$ is a positive constant. Integrating equation (2.2) from $0 \, \text{to }\, \tau_k\wedge T$ and then taking the expectation on both sides, we have

Set $\Omega_k = \{\tau_{k}\leq T\}\; \; \text{for }k\geq k_1$ and in view of Eq (2.1), we have $\mathbb{P}(\Omega_k)\geq \epsilon.$ Note that for every $\omega \in \Omega_k,$ there is $x(\tau_k, \omega)\, \, \text{or } \, y(\tau_k, \omega)\, \, \text{or } \, z(\tau_k, \omega)\, \, \text{or } \, w(\tau_k, \omega)$ that is equal to either $k\, \, \text{or}\, \, \frac{1}{k}.$ Thus $V_0\big(x(\tau_k, \omega), y(\tau_k, \omega), z(\tau_k, \omega), w(\tau_k, \omega)\big)$ is no less than either

Therefore, we obtain

It follows from Eq (2.3) that

where $I_{\Omega_k(\omega)}$ is the indicator function of $\Omega_k$. Let $k\rightarrow \infty,$ we obtain

which is a contradiction and so we must have that $\tau_{\infty} = \infty$. This completes the proof.

3.   Stationary distribution

For the stochastic model, we mainly focus on the existence of a stationary distribution. Define a critical value

Consider the integral equation

Lemma 3.1. (^[22]) Suppose that the coefficient of (3.1) is independent of $t$ and satisfies the following condition for some constant $B$:

in $U_R\subset \mathbb{R }^d_+$ for every $R > 0$, also suppose that there exists a nonnegative $C^2$-function $V(x) \; \in \; \mathbb{R }^d_+$ such that

Then system (3.1) has a solution, which is a stationary distribution.

Remark 3.1. According to the proof of Lemma 2.1 in ^[23], we know that condition (3.2) in Lemma 3.1 can be replaced by the global existence of the solutions. Hence, Theorem 2.1 indicates that (3.2) in Lemma 3.1 is satisfied.

Theorem 3.1. For any initial value $(x(0), y(0), z(0), w(0)) \in \mathbb{R }^4_+$, if $R^s_0 > 1$, then system (1.4) has at least one stationary solution $(x(t), y(t), z(t), w(t))\in \mathbb{R }^4_+.$

Proof. By Lemma 3.1, we only need to construct a nonnegative $C^2$-function $\widetilde{V}(x, y, z, w)$ and a closed set $U\subset \mathbb{R }^4_+$ such that

Define

where $e_1$ is a positive constant to be chosen later. Then applying Itô's formula to $V_1$, we have

Taking $e_1 = \frac{b_1}{(D_2+\frac{1}{2}\sigma^2_2)^2}$, then in view of (3.3) one can see that

where

Define

Then applying Itô's formula to $V_2$, one gets

Select

where $0 < \theta < \min\{ \frac{D_1-\frac{\sigma_1^2}{2}}{D_1+\frac{\sigma_1^2}{2}}, \frac{D_2-\sigma_2^2}{D_2+\sigma_2^2}, \frac{D_3-\frac{\sigma_3^2}{2}}{D_3+\frac{\sigma_3^2}{2}}\}$. Applying the same method to $V_3$, we have

where

Consequently, we have

where $M_0 > 0$ satisfies that $-M_0R^s_0+c_0 \leq-2.$ Also,

It is obvious that

and $K_n = (\frac{1}{n}, n)\times(\frac{1}{n}, n)\times(\frac{1}{n}, n)\times(\frac{1}{n}, n).$ It is evident that $V(x, y, z, w)$ is a continuous functions, so there exists a minimum point $(x_0, y_0, z_0, w_0)\in \mathbb{R }^4_+.$ Thereafter, we formulate the nonnegative $C^2$ function $\widetilde{V}$ as follows:

Combining with (3.4), (3.5), (3.6) and (3.7), we obtain that

Next, we define the following bounded closed set

For convenience, we can divide $\mathbb{R }^4_+ \setminus U_\varepsilon$ into eight regions:

Case 1. If $(x, y, z, w)\in U_1$, in view of (3.8) we get

here, $\varepsilon$ is a positive constant small enough to satisfy the following conditions:

where

Case 2 : If $(x, y, z, w)\in U_2$ according to (3.9), it implies that

Case 3 : If $(x, y, z, w)\in U_3,$ by (3.10) we have

Case 4 : If $(x, y, z, w)\in U_4$, according to (3.11), we deduce that

Case 5 : If $(x, y, z, w) \in U_5$, by condition (3.12), we conclude that

Case 6 : If $(x, y, z, w) \in U_6$, and the condition (3.13) is satisfied, it follows that

Case 7 : If $(x, y, z, w) \in U_7$, in view of (3.14) we get

Case 8 : If $(x, y, z, w) \in U_8$, from (3.15), it is deduced that

Clearly, there exists a small enough constant $\varepsilon$ to make the conclusion

Remark 3.2. In the proof of the above theorem, the construction of $V_3$ means to obtain the constant $E$ and high-order terms of $x, \, y\, \text{and}\, z.$

4.   Probability density function analysis

Let $u_1(t) = \ln x, \, \, u_2(t) = \ln y, \, u_3(t) = \ln z\; \text{and}\; u_4(t) = \ln w.$ Using the Itô's formula, an equivalent equation of system (1.4) is given as follows

Assume that $R_0^s > 1,$ then, there is the quasi-endemic equilibrium $E^* = (x^*, y^*, z^*, w^*) = (e^{u_1^*}, e^{u_2^*}, e^{u_3^*}, e^{u_4^*})$ where

Before proving Theorem 4.1, we need to firstly introduce two important standard matrices in the following Lemmas 4.1 and 4.2.

Lemma 4.1 (^[24]). For the algebraic equation $H_0^2+B_0\Sigma_0+\Sigma_0B_0^T = 0$, where $H_0 = diag(1, 0, 0, 0)$ and $\Sigma_0$ is a real symmetric matrix, we have the standard matrix

If $\gamma_1 > 0, \; \gamma_3 > 0, \; \gamma_4 > 0$ and $\gamma_1\gamma_2\gamma_3-\gamma_3^2-\gamma_1^2\gamma_4 > 0$, then $\Sigma_0$ is a positive definite matrix, where

Here $B_0$ in this form is called the standard $R_1$ matrix.

Lemma 4.2 (^[24]). For the algebraic equation $H_0^2+E_0\Omega_0+\Omega_0E_0^T = 0$, where $H_0 = diag(1, 0, 0, 0)$ and $\Omega_0$ is a real symmetric matrix, we have the standard matrix

If $\tau_1 > 0, \; \tau_3 > 0\; \text{and}\; \tau_1\tau_2-\tau_3 > 0,$ then the matrix $\Omega_0$ is a semi-positive definite matrix, which means that

Denote $P = (p_1, p_2, p_3, p_4)^T = (u_1-u_1^*, u_2-u_2^*, u_3-u_3^*, u_4-u_4^*)^T \, \, (k = 1, 2, 3, 4).$ Therefore, the linearized equation system of Eq (4.1) is obtained as follows

where

Theorem 4.1. Assume that $R_0^s > 1$; for any initial value $(x(0), y(0), z(0), w(0))\in \mathbb{R }_+^4,$ then the solution $(x, y, z, w)$ of system (1.4) will have a log-normal probability density $\Phi (x, y, z, w)$ around $E^*$, which is given by

with $\Sigma = \Sigma_1+\Sigma_2+\Sigma_3$ as a positive definite matrix; $\Sigma_1, \Sigma_2\; \text{and}\; \Sigma_3$ are defined as follows

where $q_1 = \frac{ b_{21}\alpha\big[ b_{13} b_{33}-(b_{21}+ b_{33})\alpha\big]}{ b_{31}^2 b_{33}} b_{33}, \, q_2 = \frac{\alpha(- b_{11}+\alpha)}{ b_{12}}, q_3 = -\frac{ b_{13}\alpha}{ b_{12}}-\frac{ b_{12} b_{33}q_2}{ b_{12} b_{31}+ b_{33}\alpha}+\frac{ b_{12}q_2(- b_{12} b_{33}q_2+ b_{12} b_{31}\alpha+ b_{33}\alpha^2)}{(b_{12} b_{31}+ b_{33}\alpha)^2},$ and the matrices $T_1, \ldots T_{12}, \; \Omega_1, \overline{\Omega}_2, \widetilde{\Omega}_3$ are defined in the following proof.

Proof. For convenience and simplicity, let $B(t) = (B_1(t), B_2(t), B_3(t), 0)^\tau$,

$G = diag(\sigma_1, \sigma_2, \sigma_3, 0)$ and

the linearized system (4.2) can be equivalently rewritten as

Since G is a constant matrix, according to Gardiner ^[27], we have

In view of the independence of Brownian motions $B_1(t), B_2(t)$ and $B_3(t)$, by the principle of finite independent superposition, it is easy to know that Eq (4.4) can be equivalently developed into the sum of the solution to the following three algebraic sub-equations:

where $G^2_1 = diag(\sigma^2_1, 0, 0, 0), G^2_2 = diag(0, \sigma^2_2, 0, 0), G^2_3 = diag(0, 0, \sigma^2_3, 0)$ and $\Sigma = \Sigma_1+\Sigma_2+\Sigma_3.$

Next it will be proved that $M$ has all negative real-part eigenvalues. The characteristic polynomial of $M$ is defined as

where

By calculation, we can obtain

which yields that $r_2$ and $r_3 > 0$. Furthermore, we can verify that $r_1r_2-r_3 > 0\; \text{and}\; r_1r_2r_3-r^2_3-r^2_1r_4 > 0$

So, the matrix $M$ is a Hurwitz matrix. Next, we will prove the definiteness of $\Sigma$ by following three steps.

Step 1. Consider the algebraic equation

Let $M_1 = T_1 M T^{-1}_1$, where

Take $M_2 = T_2M_1T_2^{-1}$, where

where $q_1 = \frac{ b_{21}\alpha \left[ b_{31} b_{33}-(b_{21}+ b_{31})\alpha\right]}{a^2_{31} b_{33}}.$

Case 1.1: If $q_1 = 0$, the standard transformation matrix $T_3$ is given by

Then, we calculate

where

Therefore the equivalent equation of (4.4) can be written as $(T_3T_2T_1)G^2_1(T_3 T_2T_1)^T$+$Q_1 \left[ (T_3 T_2T_1)\Sigma_1 (T_3 T_2T_1)^T\right]$+$\left[ (T_3 T_2T_1) \Sigma_1(T_3 T_2T_1)^T\right] Q^T_1 = 0.$ From Lemma 4.2, we have

where

is a symmetric positive semi-definite matrix. Hence,

Case 1.2: If $q_1\neq0$, we can find standardized transformation $T_4$ such that $Q_2 = T_4 M_2T^{-1}_4$, where

where $j_1 = -\frac{ b_{31} b_{33}q_1(b_{33}+\alpha)}{ b_{21}}$, $j_2 = \frac{q_1 \left[b^2_{31} b_{33}(b_{33}+q_1)+ b_{31}(- b_{21}+ b_{31}) b_{33}\alpha+(b^2_{21}+ b_{21} b_{31}+b^2_{31}) \alpha^2\right]}{b^2_{31}}$, $j_3 = - b_{33}q_1(b_{33}+\alpha)-\frac{(b_{21}+ b_{31})\alpha \left[b^2_{31} b_{33}q_1+(b_{21}+ b_{31})^2 \alpha^2\right]}{b^3_{31}}.$

Then, we have

Likewise, (4.4) can be transformed into the following form: $(T_4T_2 T_1) G^2_1 (T_4T_2 T_1)^T+Q_2 \left[ (T_4T_2 T_1) \Sigma_1 (T_4T_2 T_1)^T\right]+ \left[ (T_4T_2 T_1)\Sigma_1 (T_4T_2 T_1)^T\right] Q^T_2 = 0.$ From Lemma 4.1, we have

where

is a positive definite symmetric matrix; hence,

is also a positive definite matrix.

Step 2. Consider

Let $M_3 = T_5 MT^{-1}_5$, where

Take $M_4 = T_6M_3T_6^{-1}$, where

with $q_2 = \frac{\alpha(- b_{11}+\alpha)}{ b_{12}}.$

Case 2.1: If $q_2 = 0$, we let $Q_3 = T_7 M_4 T^{-1}_7,$ where the standard transformation matrix $T_7$ is given by

where $j_4 = -\frac{ b_{13}\alpha(b_{11}+ b_{33}+\alpha)(b_{12} b_{31}+ b_{33}\alpha)}{b^2_{12}}, j_5 = -\frac{ b_{13}\alpha(b_{13}b_{31}+b^2_{33}+ b_{33}\alpha+\alpha^2)}{ b_{12}}, \, j_6 = -\alpha^3+\frac{ b_{13} b_{33}\alpha(b_{33}+2\alpha)}{ b_{12}}.$ Therefore, we have

where

Hence

is a positive define matrix.

Case 2.2: Consider that $q_2\neq0$. Let $M_5 = T_8 M_4T^{-1}_8$, where

where $q_3 = -\frac{ b_{13}\alpha}{ b_{12}}-\frac{ b_{12} b_{33}q_2}{ b_{12} b_{31}+ b_{33}\alpha}+\frac{ b_{12}q_2(- b_{12} b_{33}q_2+ b_{12} b_{31}\alpha+ b_{33}\alpha ^2)}{(b_{12} b_{31}+ b_{33}\alpha)^2}.$

Case 2.2.1. If $q_2\neq0\; \text{and}\; q_3 = 0,$ let $Q_4 = T_9M_5T^{-1}_9$, where

where $j_7 = b_{11} b_{31} + b_{33} (b_{31} + q_2) +\frac{ b_{33} \alpha(b_{11}+ b_{33})}{ b_{12}}$, $j_8 = b_{13} (b_{31} +\frac{ b_{33} \alpha}{ b_{12}}) + (b_{33} +\frac{ b_{12} b_{33} q_2}{ b_{12} b_{31} + b_{33} \alpha})^2$, $j_9 = - b_{33} (b_{33} + \alpha), \tau _{21} = b_{11}+ b_{33}(1+\frac{ b_{12}q_2}{ b_{12} b_{31}+ b_{33}\alpha})$, $\tau _{22} = -\big(- b_{12} b_{21} + b_{13} b_{31} - b_{11} b_{33} +\frac{ b_{13} b_{33} \alpha}{ b_{12}} - \frac{ b_{11} b_{12} b_{33} q_2}{ b_{12} b_{31} + b_{33}\alpha}\big)$, $\tau _{23} = \frac{b_{12} b_{21} b_{33}[b_{12} (b_{31} + q_2) + b_{33} \alpha]}{b_{12} b_{31} + b_{33 }\alpha}$, $\tau _{24} = -b_{33} \bigg\{\alpha (- b_{11}+ \alpha) + b_{12}\big[ b_{21} +\frac{ b_{33} q_2 (b_{11} b_{12} b_{31} + b_{12} b_{33} q_2 - 2 b_{12} b_{31}\alpha+ b_{11} b_{33} \alpha- 2 b_{33} \alpha^2)}{(b_{12} b_{31} + b_{33} \alpha)^2} \big]\bigg\}.$ Therefore, we conclude that

where

Beside, it follows from Lemmas 4.2 that the specific form of the positive semi-definite matrix $\widetilde{\Omega_2}$ is

Therefore, we conclude that

is semi-positive definite.

Case 2.2.2. If $q_2\neq0, q_3\neq0$, we let $Q_5 = T_{10} M_5T^{-1}_{10}$, where

with

Therefore, we have

where

Hence

is also a positive definite matrix.

Step 3. Consider

We let $M_6 = T_{11} MT^{-1}_{11}$, where

Let $Q_6 = T_{12} M_6 T^{-1}_{12}$, where the standard transformation matrix $T_{12}$ is given by

By direct calculation, we obtain that

Thus (4.6) can be transformed into the following from:

where

Thus,

is a positive definite matrix.

Summing up the above steps, we have $\Sigma = \Sigma_1+\Sigma_2+\Sigma_3$ is a positive definite matrix. The solution $(x(t), y(t), z(t), w(t))$ of system (1.4) has a log-normal probability density function

5.   Extinction

Theorem 5.1. Let $(x(t), y(t), z(t), w(t))$ be the solution of the system (1.4) with any initial value $(x(0), y(0), z(0), w(0))\in \mathbb{R }^4_+.$ If $R_0^E = \frac{b_1}{D_1(D_2+\frac{1}{2}\sigma_2^2)} < 1,$ then

which means the extinction of $y(t)$.

Proof. Consider the following SDE:

with the same initial value $\hat{x} (0) = x(0) > 0,$ making use of the stochastic comparison theorem, we have

Moreover,

where $h(s) = \bigg[\big(\frac{2}{\sigma_1^2}\big)^{1+\frac{2D_1}{\sigma_1^2}}\Gamma^{-1}\big(1+\frac{2D_1}{\sigma_1^2}\big)\bigg]s^{-2(1+\frac{2D_1}{\sigma_1^2})}e^{-\frac{2}{\sigma_1^2s}}, \, \, s > 0.$ Applying Itô's formula to $\ln y(t),$ we have

Integrating (5.1) from $0$ to $t$ on both sides, one can see that

Applying the strong law of large numbers yields

Next, we take the superior limit on both sides of (5.2)

which implies that

6.   Numerical simulation

In this section, employing the Milstein higher-order method, we will present numerical simulations to verify the theoretical results. The discretization equations of model (1.4) are given by

where $\eta_{i, k}(i = 1, 2, 3;\, k = 1, 2\ldots, n)$ denotes independent Gaussian random variables which follow the distribution $N(0, 1).$

Example 6.1. Let us choose $a_1 = 1, \, a_2 = 4, \, D_1 = 1.01, \, D_2 = 8, \, D_3 = 1.02, \, b_1 = 20, c_1 = 5, \, c_2 = 5, \, \alpha = 0.09, \, \sigma_1 = \sigma_2 = 0.01, \, \sigma_3 = 0.05$ and the initial value $(x(0), y(0), z(0), w(0)) = (0.9, 0.1, 0.1, 0.1).$ We conclude that $R_0 = \frac{b_1}{D_1D_2} = 2.4752 > 1$ and $R_0^s = \frac{b_1}{(D_1+\frac{1}{2}\sigma_1^2)(D_2+\frac{1}{2}\sigma_2^2)} = 2.4751 > 1.$ Theorem 4.1 shows that there exists an ergodic stationary distribution of stochastic model (1.4). Noting that $(x^*, y^*, z^*, w^*) = (0.4, 0.1955, 0.3236, 0.1955)$, we have following the covariance matrix

Thus, we have the corresponding probability density function

As a result, $\Phi(x, y, z, w)$ has the following four marginal density functions:

we can conclude that system (1.2) admits a global positive stationary solution on $\mathbb{R } _+^4$; see Figure 1 and Figure 2.

Example 6.2. Let us choose $b_1 = 20, \, D_1 = 2.3, \, D_2 = 8\; \text{and}\; \sigma_2 = 1.8.$ Then the condition $R_0^E = 0.693 < 1$ is satisfied. Theorem 5.1 shows that MC-degrading bacteria of system (1.4) will go to extinction with probability one, which is numerically confirmed by Figure 3.

7.   Conclusions

In the current paper, we consider a stochastic Microcystins degradation model with distributed delay. We have established sufficient conditions for the existence of an ergodic stationary distribution of the positive solutions to system (1.4) by constructing a suitable stochastic Lyapunov function. The result shows that a small amount of white noise can guarantee the existence of an ergodic stationary distribution of the positive solutions to system (1.4). In addition, we have obtained the exact probability density function around a quai-equilibrium point.

Some interesting topics deserve further consideration. On the one hand, the paper focuses on the dynamical evolution and stability of the system (1.4). It should be noted that model (1.4) can potentially be solved analytically by using the Lie algebra method ^[25,26]. On the other hand, because our model is autonomous and only disturbed by white noise, it would be interesting to introduce telegraph noise, such as a continuous-time Markov chain, into system (1.3). Moreover, it would also be interesting to study more complicated MCs degradation models, such as multi-group MCs degradation models. We will leave these problems as our future work.

Use of AI tools declaration

The authors declare that they have not used artificial intelligence tools in the creation of this article.

Acknowledgments

Thanks to the reviewers for their helpful comments and suggestions. This work was supported by the Hainan Provincial Natural Science Foundation of China (No. 122RC679,121RC554), Talent Program of Hainan Medical University (No. XRC2020030) and Northeast Petroleum University Special Research Team Project (No. 2022TSTD-05).

Conflict of interest

The authors declare that there is no conflict of interest.