Stochastic modeling of algal bloom dynamics with delayed nutrient recycling

1.   Introduction

It is known that many shallow lakes in China have been experiencing a problem of algal bloom resulting in shortage of drinking water supply and degradation of lake ecosystem ^[1]. For example, it is reported that in the spring of 2007, the accumulation of algal blooms in Tai Lake, the third largest freshwater lake in China, left approximately 2 million people of Wuxi City without drinking water for at least a week, and led to massive and rapid fish deaths ^[2,3,4]. The bloom of algal is a complex process ^[5] and is due to the increased input of nutrients, such as nitrogen and phosphorus, which directly affect the growth of algal ^[6] and nutrient-rich sediment (dead algae) ^[7], which sinks to the bottom of the lake forming detritus, and then is converted into nutrients due to the decomposition by bacteria ^[8,9].

It should be noted from ^[8,9]. that the conversion of detritus into nutrients is not instantaneous and requires time for the completion of bacterial decomposition resulting a delay. To describe the conversion of detritus into nutrients, authors of ^[8,9] proposed a nutrient-phytoplankton model

where $S(t)$ and $x(t)$ denote the concentrations of a limiting nutrient (nitrogen or phosphorus) and algae density at time $t$, respectively; $S^0$ is the input concentration of the nutrient, $D$ is the washout rate, $m$ is the maximum specific uptake rate of nutrient, $D_1$ is the death rate of the algae, $\gamma\in(0, 1)$, is the fraction of nutrient conversion, the delay kernel $f(s)$ is a non-negative bounded function defined on $[0, \infty)$ describing the contribution of the biomass dead in the past to the nutrient recycling at time $t$, $\mu\in(0, 1)$ is the fraction of the nutrient recycled by bacterial decomposition of the dead algae, $U(S)$ is the specific growth function. It is found in ^[9]that the delay should be relatively small to have global attractivity of positive equilibrium. Recently, Misra and Chandra ^[11] studied the effect of delayed nutrient recycling on a nutrient-algae-dissolved oxygen model. They found that if the delay is lager than some threshold then the concentration of algae fluctuates and may increase drastically; meanwhile, the concentration of dissolved oxygen may reduce drastically, resulting massive death of fish population in the lake. To avoid the massive death of fish population and other losses, the detritus from the lake should be regularly removed before the threshold value. To provide sound scientific advices, various studies using mathematical models have been conducted for algal bloom caused by nutrients in the lakes, where both nutrient and algae population dynamics are described by a system of coupled deterministic continuous processes ^{[12,13,14,15]}.

It is known that nutrient inputs from runoff vary not only in quantity (influenced by rainfall and other environmental factors), but also in composition (based on the form of fertilizer in use) ^[16]. These variations are random in nature and generally affect the nutrient balances, influence mass transfer efficiencies along food chains ^[17], and therefore play a key role in bloom initiation and maintenance ^[18]. To better understand algal bloom phenomenon, stochastic models are needed, describing the variation of nutrients and environmental forcing. Several studies have already addressed the stochastic modelling of the algal bloom. For example, according to the experimental data from Bohai Bay of China, Huang et al. ^[19] constructed a nonlinear stochastic model on densities of two kinds of typical HAB (harmful algae blooms) algae: diatom and dianoflagellate, and analyzed their characters of stability and Hopf bifurcation; Das et al. ^[20] proposed a stochastic model with the main purpose of considering the severity and duration of algal blooms in the ecological arena; Mandal et al. ^[21] established stochastic models of allelopathic interactions between two competing phytoplankton species as a continuous time Markov chain model and as an Itô stochastic differential equation model, in which approximate extinction probabilities for both species were obtained analytically for the continuous time Markov chain model.

Our goals of this paper are to extend the delayed model (1.1) into a stochastic one for algal bloom and to explore how the random fluctuations and the delay affect the dynamics of algae population. Recall that system (1.1) has a washout equilibrium $E_0 = (S^0, 0)$ when

and a positive equilibrium $E^*(S^*, x^*)$ when

with

and here $U^{-1}$ is the inverse function of $U$ ^[8,9]. Under random fluctuations, one naturally concerns whether or not there exists a steady state for the stochastic system, and how the dynamics of algae population can be affected.

To derive a reasonable stochastic analogue of the deterministic model (1.1), it generally needs the help of Markov chain ^[22]. But this approach is not justified for (1.1) because the delayed stochastic system is non-Markovian. So, in this paper, we assume that $f(s)$ takes the following special family of memory functions:

where $\alpha>0$ is a constant, $n$ is an nonnegative integer. In particular, we call (1.4) the weak kernel when $n = 1$ and strong kernel when $n = 2$. Notice also that in the limit case when $n\rightarrow\infty$, the delay kernel $f^n_\alpha(s)$ converges to a Dirac delta function $f(s) = \delta(s-\tau_f)$. Then (1.1) becomes a model with discrete time delay $\tau_f$. To properly propose our model, we first convert (1.1) into a system of ordinary differential equations by using chain-trick, and then derive a stochastic analogue by means of Markov chain. Based on the derived stochastic model we will investigate the effects of delay and environmental random fluctuations on the dynamics of the model.

The paper is organized as follows. In Section 2, a detailed derivation of the stochastic model is performed by using discrete Markov chain. In Section 3, we show the uniqueness and global existence of a positive solution of the stochastic model. In Section 4, we carry out the long term behavior analysis of the model: we first investigate the sufficient conditions for the stochastic stability of the washout equilibrium; then the spectral densities of the nutrient and the algae population are estimated by using Fourier transform method. Our study shows that the algae population can be extinct if experiencing a sufficiently large noise. Finally, numerical simulations to illustrate the results obtained and a brief discussion are presented in Section 5.

2.   Model formulation and some preliminaries

By taking into account of random effects, we consider a stochastic analogue of the deterministic model (1.1)

with initial value

where $\varphi_1(\theta), ~\varphi_2(\theta)\in \mathcal{C}((-\infty, 0], \mathbb{R}_+)$, the families of continuous functions from $(-\infty, 0]$ to $\mathbb{R}_+$. $(\xi_1(t), \xi_2(t))$ is a two dimensional Gaussian white noise process satisfying

where $\delta_{ij}$ is the Kronecker symbol and $\delta$ is the Dirac delta function. Since $\xi_i(t)$ is delta-correlated, so $\xi_i(t)dt$ can be written as $dB_i(t)$, where $B(t) = (B_1(t), B_2(t)$ is an independent Brownian motion or Wiener process. Thus (2.1) can be rewritten as the form

We now show in detail that model (2.2) or (2.1) is a reasonable stochastic analogue of the deterministic model (1.1). To this end, introduce

then (1.1) becomes the following system of coupled ordinary differential equations

Follow the idea and techniques in ^[22], next we derive a reasonable stochastic analogue of (2.3) using a discrete time Markov chain. For a fixed time increment $\Delta t>0$ and $t = 0, \Delta t, 2\Delta t, \cdots$ define a process

Here $S^{\Delta t}(t)$ denotes the nutrient concentration, $x^{\Delta t}(t)$ denotes the concentration of algae population, and $y_i^{\Delta t}(t)$, $i = 1, 2, \ldots, n$, are auxiliaries. Let the initial value be

and $\Big\{R_j^{\Delta t}(k)\Big\}_{k = 0}^\infty$, $j = 1, 2, \ldots, n+2$ denote the $n+2$ sequences of random variables. Suppose that these variables are jointly independent and that within each sequence the variables are identically distributed such that

where $\sigma_j\geq0$ are constants reflecting the size of the stochastic effects, and

Variable $R_2^{\Delta t}(k)$ is supposed to capture the effect of random influences (due to the external factors such as variation of nutrients, environmental forcing, or inherent factors, e.g. mutations) on the concentration of the algae during the period $[k\Delta t, (k+1)\Delta t)$. We assume that $x^{\Delta t}$ grows within the time period according to the deterministic system (2.3) and by the random amount $R_2^{\Delta t}(k)x^{\Delta t}(k\Delta t)$. Random effects on $S^{\Delta t}$ can be similarly modelled as $R_1^{\Delta t}(k)S^{\Delta t}(k\Delta t)$. Specifically, for $k = 0, 1, ...$, we set

and

We next show that $X^{\Delta t}$ converges to a diffusion process as $\Delta t\rightarrow0$. To this end, we first determine the drift coefficients of the diffusion. Let $\mathcal{P}^{\Delta t}(u, dv)$ denote the transition probabilities of the homogeneous Markov chain $\{X^{\Delta t}(k\Delta t)\}_{k = 0}^\infty$, that is

for all $u = (u_1, \ldots, u_{n+2})\in \mathbb{R}^{n+2}$ and all Borel sets $A \subset \mathbb{R}^{n+2}$. Let $u_1 = S, \ u_2 = x, \ u_3 = y_1, \ldots, u_{n+2} = y_n$. Then, by (2.4) and (2.5), we have

and

To determine the diffusion coefficients, we consider the moments

It follows from (2.4) and (2.5) that

Thus,

for all $K\in(0, \infty)$. Similarly, for all $K\in(0, \infty)$, we can obtain

and

Assuming that $E[R_i^{\Delta t}(k)]^4 = o(\Delta t)$ for $i = 1, 2$, one may verify that for all $K\in(0, \infty)$,

Extend the definition of $X^{\Delta t}(t)$ to all $t\geq 0$ by setting $X^{\Delta t}(t) = X^{\Delta t}(k\Delta t)$ for $t\in[k\Delta t, (k+1)\Delta t)$ and let $X(t)$ be the solution of system

with initial condition $X(0) = (S(0), x(0), y_1(0), \ldots, y_n(0))^T$. Then according to Theorem 7.1 of ^[23], we can from (2.6)-(2.15) obtain the following lemma.

Lemma 2.1. Given initial condition $X(0) = (S(0), x(0), y_1(0), \ldots, y_n(0))^T$, assume $X(t)$ is the unique solution of (2.15). Then $X^{\Delta t}(t)$ converges weakly to $X(t)$ as $\Delta t\rightarrow 0$.

Notice that if we take $n\rightarrow\infty$ for the kernel $f^n_\alpha(s)$, then it converges to a Dirac delta function. Thus (2.1) becomes a stochastic model with discrete time delay. Then, conditional probability density can be used to determine the "conditional average drift" and "conditional average diffusion" of the process (see Chapter 6 in ^[24]).

We should point out that (2.1) and (2.2) are two equivalent stochastic counterparts of deterministic model (1.1). In the sequel, we will select the specific form of the random model according to the specific needs. In the following sections, we will first show the uniqueness of positive solution of system (2.1), and then study the dynamics of the system on a complete probability space $(\Omega, \mathscr{F}, \{\mathscr{F}_t\}_{t\geq0}, P)$ with a filtration $\{\mathscr{F}_t\}_{t\geq0}$ satisfying the usual conditions, i.e., it is right continuous and $\mathscr{F}_0$ contains all P-null sets. Before that, we introduce a differential operator $L$ associated with a general $n$-dimensional stochastic functional differential equation ^[25]:

with initial condition $x_0 = \varphi\in\mathscr{H},$ where $\mathscr{H}$ is the space of $\mathscr{F}_0$-adapted random variables $\varphi$, with $\varphi(s)\in \mathbb{R}^n$ for $s\leq0$, and

$B(t)$ denotes $m$-dimensional standard Brownian motions. The differential operator $L$ is defined as

3.   Existence and uniqueness of the global positive solution of system (2.1)

The existence and uniqueness theorem of solutions for stochastic functional differential equations had been studied by Mao ^[25], see also the related results in ^[26,27]. However, these results are all in the case that the delay is finite. Recently, Wei and Wang ^[28] generalized the theorem to the case of infinite delay under the linear growth condition and local Lipschitz condition. For our model (2.1), obviously, the coefficients are locally Lipschitz continuous, but they do not satisfy the linear growth condition. So, the solution of system (2.1) may explode at a finite time. In this section, we shall show that the solution of system (2.1) with any positive initial value remains positive for all $t\geq0$.

Theorem 3.1. Given any initial value $(\varphi_1(\theta), \varphi_2(\theta))\in \mathcal{C}((-\infty, 0], \mathbb{R}_{+}^2)$, model (2.1) has a unique solution $(S(t), x(t))$ for all $t\geq0$; furthermore, the solution remains positive for all $t>0$ with probability $1$, namely $S(t)>0$ and $x(t)>0$ for all $t\geq0$ almost surely, if the specific growth function $U(S)$ satisfies

Proof. It is easy to verify that for any given initial value $(\varphi_1(\theta), \varphi_2(\theta))\in \mathcal{C}((-\infty, 0], \mathbb{R}_{+}^2)$, system (2.1) has a unique local solution $(S(t), x(t))$ on $t\in[0, \tau_{e})$, where $\tau_{e}$ is the explosion time. To prove this theorem, we only need to show this solution is also global, that is $\tau_{e} = \infty \ a.s.$

Let $k_{0}>0$ be sufficiently large so that $\frac{1}{k_0}\leq\varphi_1(0), \varphi_2(0)\leq k_0$, and for each integer $k\geq k_{0}$ define

which is known as the stopping time and increasing as $k\rightarrow\infty$. Set $\tau_{\infty} = \displaystyle\lim_{k\rightarrow\infty}\tau_{k}$, whence $\tau_{\infty}\leq\tau_{e}\ a.s.$. We next show $\tau_{\infty} = \infty\ a.s.$ by contradiction, which implies $\tau_{e} = \infty\ a.s.$.

If $\tau_{\infty}\neq\infty$, then there is a pair of constants $T>0$ and $\delta\in(0, 1)$ such that

Hence there is an integer $k_{1}\geq k_{0}$ such that for all $k\geq k_{1}$

Consider a Lyapunov function

where $C_1 = \frac{D+D_1-\gamma\mu D_1}{\gamma mU'(0)}$. By Itô's formula, one has

where

Notice for $S>0$

Then from (3.4), we deduce that

Therefore,

Integrating both sides of the above inequality from 0 to $\tau_{k}\wedge T$, and taking expectation, yields

One then has

Set $\Omega_{k} = \{\tau_{k}\leq T\}(k\geq k_{1})$ and by (3.1), $P(\Omega_{k})\geq\delta$. Note that for every $\omega\in\Omega_{k}$, there is at least one of $S(\tau_{k}, \omega)$, $x(\tau_{k}, \omega)$ which equals either $k$ or $\frac{1}{k}$. Then

It follows from (3.5) that

where $1_{\Omega_{k}}$ is the indicator function of $\Omega_{k}$. Let $k\rightarrow\infty$, then

which leads to a contradiction. So $\tau_{\infty} = \infty$. This completes the proof.

4.   Dynamical behavior of system (2.1)

In this section, we will study the long term dynamical behavior of model (2.1), particularly near the equilibria of the corresponding deterministic model (1.1). Obviously, when $\sigma_1 = 0$ and $\sigma_2\neq0$, $E_0$ is still an equilibrium of stochastic model (2.1), but $E^*$ is not; when $\sigma_1$ and $\sigma_2$ are positive, neither $E_0$ or $E^*$ is the equilibrium of model (2.1). Therefore, we will investigate

(a) the stochastic stability of $E_0$ when $\sigma_1$ is zero; and

(b) the spectral densities of the nutrient and the algae population when $\sigma_i>0, i = 1, 2.$

4.1.   Stability of system (2.1) at $E_0$ in case of $\sigma_1=0$

When $\sigma_1 = 0$, substituting $u_1 = S-S^0,$ $u_2 = x$ into model (2.2), one obtains

which has the linearized system

For the trivial solution of system (4.2), one has the following stability result. The method of constructing Liyapunov functions in its proof will help us to construct suitable Liyapunov functions in investigating the stability of the trivial solution of nonlinear system (4.1).

Proposition 4.1. Assume $C_2 = \mu D_1[D+D_1-\gamma mU(S^0)]$ and $\sigma_1 = 0$. Then the trivial solution of system (4.2) is asymptotically mean square stable if

where $T_f$ is the average delay with value $T_f = \frac{n}{\alpha}$.

Proof. Consider the function

It follows from (4.2) and Itô's formula that

where

Then we have

Define function

which is well defined since $\int_0^\infty s^2f(s)ds = \frac{n(n+1)}{\alpha^2}$. Then by Itô's formula, we have

It then follows from (4.7) and (4.9) that

We now consider a function

where $C_3 = \frac{2D}{2(D+D_1)-2\gamma mU(S^0)-\sigma_2^2}$. By Itô's formula, we have

For function

we have from (4.10) and (4.11) that

By (4.3), we choose $\varepsilon>0$ such that

Let $T = T(\varepsilon)>0$ such that $\int_t^\infty f(s)ds < \varepsilon$ for all $t\geq T$. Then for all $t\geq T$, one has

Integrating both sides of the above from $T$ to $t\geq T$ yields

Using the similar discussion as that in ^[9] and the Barbǎlat lemma, we conclude that $E(u_1^2(t)+u_2^2(t))\rightarrow 0$ as $t\rightarrow \infty.$ Applying the definition of mean square stability of the solution ^[25], we obtain the conclusion.

Next, we give the result about the stability of the trivial solution of non-linear system (4.1), that is, the stability of $E^0(S^0, 0)$ of system (2.1).

Theorem 4.1. Let $\sigma_1 = 0$ and assume conditions (4.3) and (4.4) hold. Then the trivial solution $(0, 0)$ of system (4.1) or the equilibrium $E^0(S^0, 0)$ of system (2.1) is stochastically stable.

Proof. Let $(u_1, u_2)$ be any solution of system (4.1). Define the stopping time

For the Lyapunov function $V_1(u_1, u_2)$ defined in (4.5), we obtain

where $h(t)$ is defined as in (4.6). For $V_2(u_1, u_2)$ in (4.8), one then has

Now define

where $C_4$ is a constant to be determined later. We have that

Therefore, for the function

it follows from (4.13) and (4.14) that

By (4.4), one can find a constant $\delta>0$ such that

Choose

Then, it follows from (4.15) and (4.16) that

Integrating both sides of the above formula from $0$ to $t\wedge T_{\varepsilon_1}$, yields

where $p = \frac{1}{2}\mu D_1(D+D_1-\gamma mU(S^0))\int_0^\infty s^2f(s)ds+C_4+\mu D_1T_f$. Now for $\varepsilon_1, ~\varepsilon_2\in (0, 1)$, let

Then, if $\|\varphi_1\|^2+\|\varphi_2\|^2 < \delta^2$, it follows that

On the other hand, we have

Hence, we have $P\{T_{\varepsilon_1}\leq t\}\leq \varepsilon_2$. Letting $t\rightarrow \infty$ gives

which is equivalent to

Then the definition of the stochastic stability of the solution implies the conclusion.

4.2.   Spectral density analysis of (2.1) in case of $\sigma_i>0$

In this subsection, we study the dynamics of system (2.1) near $E^*$, the positive equilibrium of (1.1). To this end, we assume condition (1.3) holds. Thus $E^*$ is a stable positive stable positive equilibrium for system (1.1) but not the equilibrium for system (2.1) when $\sigma_i\neq 0$. We aim to investigate the spectral densities, denoting the intensities of fluctuations, of the nutrient and the algae population by Fourier transform method.

To perform the spectral density analysis on model (2.1) with a general delay kernel function is very difficult, so we turn to consider the limit case when $n\rightarrow\infty$ in (1.4), i.e., $f(s) = \delta(s-\tau_f)$. Introducing $S(t) = \exp(u_1(t))$ and $x(t) = \exp(u_2(t))$, one converts (2.1) into

Substituting

into (4.18), yields

where $(u_1^*, u_2^*) = (\ln S^*, \ln x^*)$. The linearized system of (4.19) is

where

Given continuous function $v(t)$ over the interval $-\frac{T}{2} \leq t \leq \frac{T}{2}\ (T>0)$, define function

Since $\tilde{v}(\omega)$ is the Fourier transform of $v(t)$, we know

which implies that $\frac{1}{2\pi}\tilde{v}(\omega)$ is the amplitude density of the components of $v(t)$ in the angular frequency interval $\omega$ to $\omega+d\omega$. Thus $\frac{1}{2\pi}\tilde{v}(\omega)d\omega$ can be considered as an estimate of the amplitude of the component of $v(t)$ with angular frequency $\omega$. Taking Fourier transform of system (4.20), we obtain

where

Write (4.23) in the matrix form

where

$\tilde{\eta}(\omega) = (\tilde{\eta_1}(\omega), \tilde{\eta_2}(\omega))^T = (\sigma_1\tilde{\xi}_1(\omega), \sigma_2\tilde{\xi}_2(\omega))^T, \ \tilde{v}(\omega) = (\tilde{v}_1(\omega), \tilde{v}_2(\omega))^T$. Let $|A(\omega)| = \det(A(\omega))$ and we assume $|A(\omega)|\neq 0$. Then from (4.24) we have

where

Then

If the function $v(t)$ has zero mean value then the fluctuation intensity (variance) of the components in the frequency band $\omega$ and $\omega+d\omega$ is $S_v(\omega)d\omega$, where the spectral density $S_v(\omega)$ is formally defined by

Hence,

From (4.27) and (4.28), we have

since $\overline{\xi_i(t)} = 0$ and $\overline{\xi_i(t)\xi_j(t')} = \delta_{ij}\delta(t-t')$. Therefore, the fluctuation intensity (variance) in $v_i$ is

It follows from (4.26) that

where

When $\tau_f = 0$, we have

Following Gradshteyn and Ryzhik ^[29], the general integral encountered in calculations of fluctuation is of the type

where

When $n = 2$, the integral is given by

Then from the relation between $g_2(\omega)$ and $P_i(\omega) (i = 1, 2)$ and $h_2(\omega)h_2(-\omega)$ with $M(\omega)$ we obtain

Hence,

Now, the variances have been obtained in (4.31) for $\tau_f>0$ and in (4.36) for $\tau_f = 0$. But we should point out that the computed variance $\sigma_{v_i}^2$ is for the linear system (4.20), while for the nonlinear system (4.19), it is difficult to obtain. So, $\sigma_{v_i}^2$ in (4.31) is an estimate of the variance for system (4.19). Notice that $v_1 = \ln S-\ln S^*$ and $v_2 = \ln x-\ln x^*$, so the fluctuation intensities of the nutrient and the algae population can be estimated as $(S^*)^2\sigma_{v_1}^2$ and $(x^*)^2\sigma_{v_2}^2$, respectively. Thus, the results obtained so far in this subsection can be summarized as follow.

Theorem 4.2. Assume (1.3) holds. Then the fluctuation intensities of the nutrient and the algae population can be estimated as $(S^*)^2\sigma_{v_1}^2$ and $(x^*)^2\sigma_{v_2}^2$, respectively, where $\sigma_{v_i}^2$ is given in (4.31) for $\tau_f>0$ and in (4.36) for $\tau_f = 0$.

4.3.   Extinction

The following theorem shows that a sufficiently large noise can make the algae population extinct exponentially with probability one.

Theorem 4.3. For any given initial value $(S(0), x(0))\in \mathbb{R}^2_+$, the solution $(S(t), x(t))$ of system (2.1) has the following property:

In particular, if $\sigma_2^2\geq 2\gamma m-2(D+D_1)$, then

Proof. Define function $V(x) = \ln x$. Then by Itô's formula, we have

It follows that

integrating both sides of which from $0$ to $t$ yields

where $M(t) = \sigma_2B_2(t)$. Obviously, $M(t)$ is a locally continuous martingale and

which implies that

By Strong Law of Large Numbers, we obtain

Dividing $t$ on the both sides of (4.37) and letting $t\rightarrow \infty$, we have

This completes the proof of the Theorem.

Theorem 4.3 reveals that a large noise may induce the extinction of algae population even though its corresponding deterministic model (1.1) has a stable positive equilibrium $E^*$.

5.   Numerical simulations and discussions

In this paper, a delayed stochastic model (2.1) for algal bloom with nutrient recycling has been proposed and investigated. Such a stochastic model is useful in investigating and understanding various ecologically realistic features. We have focused our attention on two aspects:

(ⅰ) the random influences incorporated through perturbation on the nutrient and the algae population, and

(ⅱ) the conversion delay of detritus into nutrients.

For model (2.1), we first showed its reasonability by means of an approximate Markovian system of it under the assumption that the delay kernel $f(s)$ takes the family of generic delay kernel (1.4). Then we carried out the analysis of the uniqueness and the global existence of its positive solution with the help of the result in ^[28], since the incorporated delay in the system is infinite. Next, we analyzed its long time behaviors around the various equilibria of its corresponding deterministic model (1.1). Our findings in Theorem 4.1 reveal that $E_0$ is stochastically stable provided the intensity of the noise $\sigma_2$ and the average delay $T_f = \frac{n}{\alpha}$ are small. Though $E^*$ is not an equilibrium of model (2.1) when $\sigma_i\neq 0$, it is shown in Theorem 4.4 that the fluctuations intensities of the nutrient and the algae population can be estimated in a neighbourhood of $E^*.$ In addition, our result in Theorem 4.4 reveals that sufficiently large noise can make the algae population extinct exponentially with probability one, even if its corresponding deterministic model (1.1) has a stable positive equilibrium.

To illustrate the theoretic results obtained, numerical simulations are carried out by using Milstein scheme ^[30]. Here we assume that the specific growth function $U(S)$ is of Michaelis-Menten type

where $a_1$ is the half-saturation constant.

The first given example below concerns the effects of the random influence and the delay on the long time behavior around the washout equilibrium. Here we take $n = 1$ in (1.4), that is, the weak delay kernel $f(s) = \alpha e^{-\alpha s}$. Then the discretization of model (2.1) for $t = 0, \Delta t, 2\Delta t, \ldots, n \Delta t$ takes the form

where time increment $\Delta t>0$ and $\xi_i$ are $N(0, 1)-$distributed independent random variables which can be generated numerically by pseudorandom number generators.

Example 5.1. Consider model (2.1) with $D = 0.3, D_1 = 0.1, S^0 = 3, m = 0.54, a_1 = 0.4, \mu = 0.3, \gamma = 0.8$, $\sigma_1 = 0$, $\sigma_2 = 0.05$ and $\alpha = 10$.

It is easy to see that $T_f = 0.1$. Simple computations show that

and

Then by Theorem 4.1, the equilibrium $E_0$ of model (2.1) is stochastically stable. Our simulation supports this result as shown in Figure 1. To examine the effect of the delay in nutrient recycling, increasing the value of $T_f$ from $0.1$ to $10$ (i.e., decreasing the value of $\alpha$ from $10$ to $0.1$), we find that the equilibrium $E_0$ continues to be stochastically stable, but the levels of the nutrient and the algae population will be lower in the beginning of time (see Figure 1). It is because that a large delay in nutrient recycling can make the algae grow slowly.

To study effects of the random influence and delay on $E^*$, we give the following example with different values of noise intensities and delay. For the kernel function $f(s) = \delta(s-\tau_f)$, the discretization of model (2.1) for $t = 0, \Delta t, 2\Delta t, \ldots, n \Delta t$ takes the form

where time increment $\Delta t>0$ and $\xi_i$ are $N(0, 1)-$distributed independent random variables which can be generated numerically by pseudorandom number generators.

Example 5.2. Consider model (2.1) with $D = 0.3, ~D_1 = 0.1, ~S^0 = 5, ~m = 0.7, ~a_1 = 0.4, ~\mu = 0.3$ and $\gamma = 0.8$.

It is easy to see example 5.2 satisfies condition (1.3). By Theorem 4.2, the fluctuation intensities of the nutrient and the algae population can be estimated as $\sigma_S^2$ and $\sigma_x^2$, respectively, where $\sigma_S^2 = (S^*)^2\sigma_{v_1}^2$ and $\sigma_x^2 = (x^*)^2\sigma_{v_2}^2$. The variations of $\sigma_S^2$ and $\sigma_x^2$ with respect to the delay $\tau_f$ and noise intensities $\sigma_i$ are drawn in Figure 2 for different $\tau_f$ and $\sigma_i$, and the corresponding trajectories of the nutrient and the algae around $E^*$ are drawn in Figure 3.

Figures 2(a) and (b) show the effect of the delay on $\sigma_S^2$ and $\sigma_x^2$, which are drawn by taking $\sigma_1 = 0.01$ and $\sigma_2 = 0.01$. From Figure 2(a), $\sigma_S^2$ has a little decrease with the increase of the delay in the beginning, then it increases as the delay increases and finally it remains unchanged when the delay increases to $30$. From Figure 2(b), $\sigma_x^2$ first decreases as the delay increases, then it remains unchanged when the delay increases to $30$. These two sub-figures reveal that the delay has a different effect on $\sigma_S^2$ and $\sigma_x^2$ when it is small, while it does not affect the fluctuation intensities when it is large. The corresponding trajectories of the nutrient and the algae around $E^*$ are drawn by taking $\tau_f = 0.1$ and $\tau_f = 10$, respectively, please see Figure 3(a) and (b). These two sub-figures reveal that with a different value of $\tau_f$, trajectories of the nutrient and the algae oscillate ultimately around $E^*$ with a certain amplitude, but with a larger value of $\tau_f$, the levels of the nutrient and the algae will be lower in the beginning.

Figures 2(c) and (d) show the effect of noise intensities on $\sigma_S^2$ and $\sigma_x^2$, which are drawn by taking $\tau_f = \sigma_2 = 0$ in Figure 2(c) and $\tau_f = \sigma_1 = 0$ in Figure 2(d). From these two sub-figures, the fluctuation intensities of both the nutrient and the algae increase as $\sigma_1$ or $\sigma_2$ increases, and $\sigma_S^2$ increases faster in the absence of $\sigma_2$ and $\sigma_x^2$ grows faster in the absence of $\sigma_1$. The corresponding trajectories of the nutrient and the algae around $E^*$ are drawn by taking two different values of $\sigma_1$ (Figure 3(c)-(d)) and of $\sigma_2$ (Figure 3(e)-(f)).

From Figures 2 and 3, one can find that the fluctuation intensities $\sigma_S^2$ and $\sigma_x^2$ are more sensitive to the noise intensities $\sigma_i$ than the delay $\tau_f$, since the fluctuation intensities change within a narrow range with respect to $\tau_f$ (Figures 2 and 3(a)-(b)), and they change within a wide range with respect to $\sigma_S^2$ and $\sigma_x^2$ (Figures 2(c)-(d) and 3(c)-(f)). Furthermore, the fluctuation intensity of the nutrient is more sensitive to $\sigma_1$ than $\sigma_2$, since in the absence of $\sigma_2$, $\sigma_S^2$ has a faster increase with the increase of $\sigma_1$. Rather, the fluctuation intensity of the algae population is more sensitive to $\sigma_2$.

To further study the effects of the random influence on the extinction of system (2.1), we give the following example with a sufficiently large noise. Here the discretization of model (2.1) takes the form as in (5.1).

Example 5.3. Consider model (2.1) with $\alpha = 10$ (i.e., $T_f = 0.1$), $\sigma_1 = 0$ and $\sigma_2 = 0.5$ and all other parameters have the same values as in Example 5.2.

Simple computation shows that

By Theorem 4.3, the algae population will go to extinction (see Figure 4).

Examples 5.2 and 5.3 reveal that the algae population is more sensitive to $\sigma_2$. Under certain parametric conditions, there is fundamentally different behavior for the algae with different value of $\sigma_2$. If $\sigma_2 = 0.01$, the algae population will survive. Whereas, when $\sigma_2$ increase to $0.5$, the extinction of the algae will occur.

To sum up, this paper presents an investigation on the effect of the environmental noise and the delay occurred in the nutrient recycling on a nutrient-algae system. Our findings are useful for a better understanding of the dynamics of algal blooms. We should point out there are other interesting topics meriting further investigation, for example, the stationary distribution of the system. It is also very interesting to study the long time behavior of the multi-nutrient multi-algae system with noise. We leave these for future considerations.

Acknowledgments

Research is supported by the National Natural Science Foundation of China (11671260) and Shanghai Leading Academic Discipline Project (No. XTKX2012).

Conflict of interest

The authors declare no conflict of interest in this paper.