Long-time analysis of a stochastic chemostat model with instantaneous nutrient recycling

1.   Introduction

The chemostat is a laboratory apparatus used for the continuous culture of microorganisms. It plays an important role in biotechnology, microbial ecology, and population biology, and is the most simple idealization of a biological system for population studies ^[1]. In recent decades, chemostat modeling has drawn great attention from mathematics and ecology, see ^[2,3,4,5] and reference therein. One classic chemostat with single species and single substrate is expressed in the form of

where $S(t)$ and $x(t)$ represent the substrate concentration and the microbial concentration, respectively; $S^0$ is the original input concentration of the nutrient; $D$ is the dilution rate; $\delta$ is the growth yield constant reflect the conversion of nutrient to organism, obviously, $\frac{1}{\delta}\geq1$, see literature ^[6] for details; $p(S)$ is the per-capita growth rate of the species. The nutrient uptake function $p$ is generally assumed to satisfy: $p$ is continuously differentiable; $p(0) = 0$, $p(S) > 0$ for $S > 0$. Some kinds of the nutrient uptake function are often found in literature are as follows:

● Lotka-Volterra (see e.g., ^[7]): $p(S) = mS$.

● Michaelis-Menten (Monod) (see e.g., ^[8]): $p(S) = \frac{mS}{a+S}$, here $m$ is the maximal growth rate of the microbial species and $a$ is the half-saturation constant. The test result of the relationship between microbial proliferation rate and substrate concentration obtained by Monod is consistent with the form of Michaelis-Menten equation. So, function $p(S) = \frac{mS}{a+S}$ is also known as the Monod type function.

● Monod-Haldane (see e.g., ^[9]): $p(S) = \frac{mS}{a+S+bS^2}$, here the term $bS^2$ is an inhibition.

Butler et al. ^[10] established a multiple competing species chemostat model with a general class of functions describing nutrient uptake. If there is only one microbial species, we can obtain the dynamics of system (1.1) from ^[10].

The deterministic models can neglect the stochastic effects only at the macroscopic scale, thanks to the law of large numbers. However, ecosystem is inevitably affected by environmental stochasticity which cannot be neglected at microscopic scale. Also the accumulation of small perturbations in the context of multi-species could not be neglected (see Campillo et al. ^[11] for more details on this respect). Moreover, stochastic noises may change the behavior in a substantial manner by directly acting on the population densities or affecting the parameter values of ecosystems. In recent years, much research has considered the stochastic influence ^[12,13,14]. Persistence and extinction are the central question in ecosystems, including microorganism culture ^{[15,16,17,18]}. For example, Meng et al. ^[15] established sufficient criteria for extinction and weak persistence of a nonlinear impulsive stochastic chemostat system in a polluted environment. By using Lyapunov functions method, Lv et al. ^[17] obtained sufficient conditions for the existence of a unique ergodic stationary distribution of an impulsive stochastic chemostat model similar to that in ^[15]. While Nguyen ^[18] proved sufficient and almost necessary condition to determine the persistence and extinction of the model by defining a new threshold parameter. In addition to the research on the central question above, Yang et al. ^[19] constructed a stochastic chemostat model with degenerate diffusion using a discrete Markov chain. By solving the corresponding Fokker-Planck equation, they derived the explicit expression of the stationary joint probability density and investigated the effect of white noise on the variance and skewness of the concentration of microorganisms. Baratti et al. studied bioreactor models under the assumption of constant volume, and described the probability density function evolution by using the Fokker-Planck theory in ^[20,21]. Besides, Zhang and Yuan ^[22] investigated the existence of ergodic stationary distribution of a stochastic delayed chemostat model by using the stochastic Lyapunov analysis method. Subsequently, Zhang and Yuan ^[23] considered a stochastic chemostat model with a distributed delay and proved the sufficient condition for the extinction and the ergodicity of the solution. For other work related to the delay stochastic chemostat model, we can refer to ^[24,25,26].

In industrial applications and most ecological situations an undisturbed dilution rate cannot be expected. There are different possible approaches to introduce noise into stochastic differential equations, both from a biological and from a mathematical perspective. One traditional approach is analogous to that of Stephanopoulos et al. ^[27] who superimposed a one-dimensional white noise process on the dilution rate. In particular, the authors in ^[28] replaced the dilution rate $D$ by $D+\sigma_i\mathrm{d}B_i(t)$. Then, model (1.1) becomes the following system

where $B_i(t)$, $i = 1, 2$ are two mutually independent Brownian motions. Under the following assumptions:

(ⅰ) $p\in C^2([0, \infty), [0, \infty))$, $p(S)\leq cS$ for any $S\in(0, +\infty)$;

(ⅱ) $p''(S)S^3\geq m_0$, for any $S\in(0, +\infty)$,

the authors in ^[28] showed that system (1.2) admits a stationary distribution which is ergodic if

Obviously, Michaelis-Menten response function $p(S) = \frac{mS}{a+S}$ satisfies the above assumptions.

Moreover, nutrient recycling is an important factor among the many processes which influence ecosystem dynamics. Many researchers have studied the effect of nutrient recycling on ecosystem stability and persistence, see ^[29,30]. Usually, nutrient recycling is regarded as an instantaneous term ^[31]. So in this work, we consider a stochastic chemostat model with instantaneous nutrient recycling, and for convenience, we take $p(S) = \frac{m S}{a+S}$. Then, the following model is obtained

where, $\gamma$ is mortality rate of microorganism population, $b\in(0, 1)$ is the fraction of the nutrient recycled by bacterial decomposition of the dead microorganism. Obviously, if $\gamma = 0$, then system (1.3) becomes model (1.2) with the case $p(S) = \frac{mS}{a+S}$.

There exist some important questions concerning the system: (ⅰ) Under what conditions will microorganism population become extinct? (ⅱ) When do microorganism population is persistent? (ⅲ) How is dynamic behavior affected by the noise intensity? In the existing literature, one of the traditional method analogous to that in ^[28] to solve the above central question is to obtain the sufficient but not necessary conditions for the existence of ergodic stationary distribution of the system by using Khasmiskii theory. A different path taken by Imhofa and Walcherb ^[12] is to establish appropriate Lyapunov functions to analyze the long-term dynamic behavior of systems, including stochastic stability. In our paper, we develop a new treatment to deal with system (1.3) to prove the sufficient and almost necessary conditions of the extinction and persistence for the model and analyze the effect of noise intensity on system dynamics. Our results essentially improve the corresponding research in ^[28]. One of the main difficulties in dealing with this model stems from the fact that the stochastic comparison arguments do not work well for system (1.3) due to the nutrient cycling term. Therefore, some new techniques must require here. The rest of the paper is organised as follows. Section 2 presents some preliminary results and defines the threshold $\lambda$ which determines the persistence and extinction of the system. In Section 3, we focus on the condition for the extinction and persistence of system (1.3). Finally, some numerical examples and brief discussions are given in Sections 4 and 5.

2.   Preliminaries and the threshold $\lambda$

Let $(\mathrm{\Omega}, \mathcal{F}, \mathbb{P})$ be a complete probability space with a filtration $\left\{\mathcal{F}_{t}\right\}_{t \geq 0}$ satisfying the usual conditions (i.e., it is increasing and right continuous while $\mathcal{F}_{0}$ contains all $\mathbb{P}$-null sets). $\mathbb{R}_{+}^{2}: = \{(x, y): x \geq 0, y \geq 0\}$ and $\mathbb{R}_{+}^{2, \circ} = \{(x, y):x > 0, y > 0\}$. Denote $a \wedge b = \min \{a, b\}$, $a \vee b = \max \{a, b\}$.

Consider the deterministic system corresponding to (1.3)

We can derive from the monotonicity of the function $p(S) = \frac{mS}{a+S}$ that there exists a uniquely defined positive real number $\lambda_*$ such that $p(S) < D$ for $0 < S < \lambda^*$; $p(S) > D$ for $\lambda_* < S$; $p(S) = D$ for $\lambda_* = S$. Here, $\lambda_*$ represents the break-even concentration of the substrate for the species $x(t)$.

Define $\lambda_d = \frac{m S^0}{a+S^0}-D-\gamma$. Similar to ^[10], when $p(S^0) < D$, i.e., $\lambda_d < 0$, there exists a boundary equilibrium $E_0 = (S^0, 0)$ for system (2.1) which is asymptotically stable, while when $p(S^0) > D$, i.e., $\lambda_d > 0$, there exists an internal equilibrium $E_* = (s_*, x_*) = \left(\lambda_*, \frac{D\delta (S^0-\lambda_*)}{D+(1-b\delta)\gamma}\right)$ which is asymptotically stable.

The following theorem concerns the existence and uniqueness of positive solutions. The proof is standard, so we only present it without proof.

Theorem 2.1. For any given initial value $(S(0), x(0)) = (u, v)\in \mathbb{R}_{+}^{2}$, system (1.3) has a unique global solution $(S_{u, v}(t), x_{u, v}(t))$ on $t \geq 0$, and the solution will always remain in $\mathbb{R}_{+}^{2}$ with probability 1.

Lemma 2.1. For any initial value $(S(0), x(0))\in\mathbb{R}_+^2$, the solution of system (1.3) is stochastically ultimately bounded.

Proof. Let $V(t) = S(t)+\frac{1}{\delta} x(t)$, by making use of Itô formula, we have

where

Then, applying the argument of the proof of [23,Theorem 3.1], one can prove the theorem easily, we omit the details here and refer the reader to [23,Theorem 3.1].

Remark 2.1. Lemma 2.1 implies that for any $T > 0$, there exists positive $H_0$ such that for any $t > T$, $S(t)\leq H_0$, $x(t)\leq H_0$. $

Lemma 2.2. The following assertions hold:

(ⅰ) Let $p_1$ be a positive solution of

Then for any $0\leq q < p\leq1$ and $p < p_1$, there exists a constants $M$ such that,

(ⅱ) For any $\epsilon \in(0, 1)$, $T > 0$ and $H > 1$, there is $\bar{H} = \bar{H}(\varepsilon, T, H) > 1$ such that for any initial value $(u, v) \in[H^{-1}, H]\times[0, H]$,

Proof. Consider the Lyapunov function $V(S, x) = \left(S+\frac{1}{\delta} x\right)^{1+p}+S^{-q}$. By directly calculating the differential operator $\mathcal{L}V(S, x)$, we obtain

Let $q < \bar{p} < p$, Young's inequality says that

Therefore,

We choose $0 < H_{2} < (1+p) D-\frac{p(1+p)}{2}\left(\sigma_{1}^{2} \vee \sigma_{2}^{2}\right)$. Since $1+\bar{p} < 1+p$ and $\frac{q(1+\bar{p})}{\bar{p}} < 1+q$, we have

Hence,

Define the stopping time $\tau_{n} = \inf \{t \geq 0:\ V(S(t), x(t)) \geq n\}$. By Itô formula and (2.4) yield that

By letting $n \rightarrow \infty$, applying Fatou's lemma obtains the part (ⅰ) of Lemma 2.1,

The proof for the second part is standard, we omit the details here and refer the reader to [32,Lemma 2.1].

Consider the first equation of (1.3) on boundary when microorganism population is absent, i.e.,

We write $\widetilde{S}_{u}(t)$ for the solution of the Eq (2.5) with initial value $\widetilde{S}(0) = u$. By solving the Fokker-Planck equation, the Eq (2.4) has an ergodic stationary distribution $\nu_{0}$ with density function

where $\alpha = \frac{2 D}{\sigma_{1}^{2}}+1$, $\beta = \frac{2 D S^{0}}{\sigma_{1}^{2}}$. Then we have

Lemma 2.3. ^[33] The following assertions hold:

(ⅰ) $\lim\limits_{t \rightarrow \infty} \frac{\ln \widetilde{S}(t)}{t} = 0$, a.s.

(ⅱ) The unique stationary distribution $\nu_{0}$ of the Markov process $\widetilde{S}(t)$ satisfies

Define

The inequality $\frac{m s}{a+s} \leq \frac{m}{a} s$ implies that

As a result, $\lambda$ is well-defined.

Remark 2.2. The solution $\widetilde{S}(t)$ of the boundary equation (2.5) has a stationary distribution $\nu_{0}$. Let $S$ be a random variable with distribution $\nu_{0}$. Then, for the nonlinear convex response function $p(S) = \frac{m S}{a+S}$, there are parameters such that $\mathbb{E}_{\nu_0}[p(S)] > D+\frac{1}{2}\sigma_2^2 > D > p(S^0)$ for $\gamma = 0$, namely, $\lambda > 0 > \lambda_{d}$; further suppose that $\sigma_{2} = 0$, then there must be $\lambda > \lambda_{d}$ (see Figure 1).

Furthermore, we consider the perturbed equation of model (2.5)

According to the comparison theorem, $\widetilde{S}(t)\leq\widetilde{S}^{\theta}(t)$. Similar to (2.5), there exists a unique invariant probability measure $\nu_\theta$ to the solution $\widetilde{S}^{\theta}(t)$ of (2.7). The following key lemma reflects the relationship between the perturbation equation of (2.5) and itself.

Lemma 2.4. ^[34] We have

3.   Extinction and persistence

This section analyzes the long term properties of the system by using the parameter $\lambda$. Roughly speaking, if $\lambda < 0$, microorganism population will die out. If $\lambda > 0$, the system has an invariant probability measure and the transition probability of the solution process converges to the invariant measure and we refer to the cases as the system being persistent.

Lemma 3.1. Assume that $\lambda < 0$, for any $\epsilon > 0$, $H > 1$, $\rho > 0$, there exists $\delta_{0} = \delta_{0}(\epsilon, H)\in (0, H)$ such that for any initial value $(u^{\prime}, v^{\prime})\in[H^{-1}, H] \times(0, \delta_{0}]$,

Proof. By the ergidicity of $\widetilde{S}^{\theta}(t)$, we obtain

Consequently, for any $\epsilon > 0$, there exists a $T_{1} = T_{1}(\epsilon, H) > 0$ such that $\mathbb{P}_{H}\left(\Omega_{1}\right) \geq 1-\epsilon$, where

According to the uniqueness of solution, we have $\widetilde{S}_{u^{\prime}}^{\theta}(s) \leq \widetilde{S}_{H}^{\theta}(s)$, $\forall t \geq 0$, where the subscript of $\widetilde{S}_{u^{\prime}}^{\theta}(s)$ indicates the initial value $\widetilde{S}^{\theta}(0) = u^{\prime}$, which implies that $\mathbb{P}_{u^{\prime}}\left(\Omega_{1}\right) \geq 1-\varepsilon, \ \forall u^{\prime} \in\left[H^{-1}, H\right]$. From Lemma 2.4, we obtain that there is a $\theta_{0} \in[0, 1]$ such that $\mathbb{P}_{u^{\prime}}\left(\Omega_{2}\right) \geq 1-\epsilon$, where

The strong law of large numbers for Brownian motions

implies that there is a $T_{2} > 0$ such that $\mathbb{P}_{u^{\prime}, v^{\prime}}\left(\Omega_{3}\right) \geq 1-\epsilon$, where

Let $T: = T_{1} \vee T_{2}$. Using (2.3), there is $\bar{H} > 1$ such that $\forall\left(u^{\prime}, v^{\prime}\right) \in\left[H^{-1}, H\right] \times(0, H]$ and $H > 1$, $\mathbb{P}_{u^{\prime}, v^{\prime}}\left(\Omega_{4}\right) \geq 1-\varepsilon$, where

To proceed, we will establish the difference of $(\widetilde{S}^{\theta}(t)-S(t))$. It follows from (1.3) and (2.7) that

Denoting

Then, by constant variation formula we have

Choosing positive constants $\rho$ and $\theta$ satisfying

such that for $\forall\left(u^{\prime}, v^{\prime}\right) \in\left[H^{-1}, H\right] \times(0, \theta]$, $\mathbb{P}_{u^{\prime}, v^{\prime}}\left(\Omega_{5}\right) \geq 1-\epsilon$ where

Indeed, let $\varrho = \tau^{\theta} \wedge \inf \{t: S(t) \geq \bar{H}\}$, we can obtain by the inequality $\left(\sum\limits_{i = 1}^{n} a_{i}\right)^{2} \leq 2^{n} \sum\limits_{i = 1}^{n} a_{i}^{2}$, Young's inequality and Burkholder-Davis-Gundy inequality that for any $t\in[0, T]$,

where $M_{1} = 16\left(m^{2}+b^{2} \gamma^{2}\right) T^{2}, M_{2} = 16\left(D^{2} T+4 \sigma_{1}^{2}\right)$. Then further using Gronwall's inequality we have

Hence, for arbitrarily small $\theta$ as defined above,

By virtue of

we obtain

Define a stoping time

Applying exponential martingale inequality, there is $\delta_{1} \in(0, \theta)$ such that $\mathbb{P}_{u^{\prime}, v^{\prime}}\left(\Omega_{6}\right) \geq 1-\epsilon$, where

As a consequence, for $\omega \in \bigcap_{i = 1}^{6} \Omega_{i}$, we have $\tau^{\theta} \geq T$. Applying comparison argument, one has $S(t) \leq \widetilde{S}^{\theta}(t)$ for $\omega \in \bigcap_{i = 1}^{6} \Omega_{i}$ and $t \leq \tau^{\theta}$. Combining (3.2), (3.3) and (3.5), one has that for $\forall\left(u^{\prime}, v^{\prime}\right) \in\left[H^{-1}, H\right] \times\left(0, \delta_{1}\right]$, $\omega \in \bigcap_{i = 1}^{6} \Omega_{i}$ and $t \geq T$,

Therefore, for $\forall\left(u^{\prime}, v^{\prime}\right) \in\left[H^{-1}, H\right] \times\left(0, \delta_{1}\right], \omega \in \bigcap_{i = 1}^{6} \Omega_{i}, t \in[T, \tau]$,

In addition, by the law of iterated logarithm, there exits $q_{\epsilon}$ such that $\mathbb{P}_{u^{\prime}, v^{\prime}}\left(\Omega_{7}\right) \geq 1-\epsilon$ where

Denote $q_{\epsilon}(t) = q_{\epsilon} \sqrt{t(|\ln t|+1)}$. Hence, for all $\left(u^{\prime}, v^{\prime}\right) \in\left[H^{-1}, H\right] \times\left(0, \delta_{1}\right], \omega \in \bigcap_{i = 1}^{7} \Omega_{i}$ and $t \geq T$, one has

where

Let $\delta_{0} \in\left(0, \delta_{1}\right)$ satisfying

Then, for all $\left(u^{\prime}, v^{\prime}\right) \in\left[H^{-1}, H\right] \times\left(0, \delta_{0}\right], \ \omega \in \bigcap_{i = 1}^{7} \Omega_{i}$ and $t\geq T$, we obtain

It follows that $t \wedge \tau \leq \tau^{\rho}, \forall t \geq T$. Therefore,

Since $\tau = \tau^{\theta} \wedge \tau^{\rho}$, we have

Hence, for all $\left(u^{\prime}, v^{\prime}\right) \in\left[H^{-1}, H\right] \times\left(0, \delta_{0}\right], \omega \in \bigcap_{i = 1}^{7} \Omega_{i}, \ t \geq T$,

This means $t \wedge \tau^{\theta} < \tau^{\theta}, \ \forall t \geq T$. As a result

Hence, for any $\left(u^{\prime}, v^{\prime}\right) \in\left[H^{-1}, H\right] \times\left(0, \delta_{0}\right]$ and $\omega \in \bigcap_{i = 1}^{7} \Omega_{i}$,

where $\mathbb{P}_{u^{\prime}, v^{\prime}}\left\{\omega \in \bigcap_{i = 1}^{7} \Omega_{i}\right\} > 1-7 \varepsilon$.

Remark 3.1. 3.1 shows that when $\lambda < 0$, the component $x(t)$ of the solution $(S(t), x(t))$ starting from any initial value in the designated area $[H^{-1}, H]\times(0, \delta_0]$ finally enters the interval $(0, \delta_0]$ and tends to zero exponentially. Biologically, it implies that microorganism population is eventually extinct.

Theorem 3.1. If $\lambda < 0$, for any initial value $(u, v) \in \mathbb{R}_{+}^{2, \circ}$, we have $\mathbb{P}\left\{\lim\limits _{t \rightarrow \infty} \frac{\ln x_{u, v}(t)}{t} = \lambda\right\} = 1$, i.e., the number of microorganism population in the chemostat will tend to zero with exponential rate.

Proof. We shall use the manners established in ^[34] to prove this theorem. By virtue of Lemma 2.2 and Cheybshev's inequality, for any initial value $(u, v) \in \mathbb{R}_{+}^{2, \circ}$, one has

where

In view of Lemma 3.1, the process $(S(t), x(t))$ is not recurrent in $\mathbb{R}_{+}^{2, \circ}$. Since $B_{i}(t), \ i = 1, 2, 3$ are independent, $(S(t), x(t))$ must be transient ^[35]. Denote

which is a compact subset of $\mathbb{R}_{+}^{2, \circ}$. It follows from the definition of transience that

Combining (3.11) and (3.12), one has

Therefore, there exits $T^{\prime} > 0$ such that

The Markov property of the process, (3.1) and (3.13) deduce that

Since $\epsilon$ is arbitrary, one has

Theorem 3.2. If $\lambda > 0$, for any initial value $(u, v) \in \mathbb{R}_{+}^{2, \circ}$, the solution of the system (1.3) has a unique invariant probability measure $\pi$ with support $\mathbb{R}_{+}^{2, \circ}$. Moreover, for any $\pi$-integrable $f(x, y): \mathbb{R}_{+}^{2, \circ} \rightarrow \mathbb{R}$, we have

and

where $P(t, (u, v), \cdot)$ is the transition probability of $(S(t), x(t))$ and $\|\cdot\|$ is the total variation norm.

Proof. Let $(S(0), x(0)) = (u, v) \in \mathbb{R}_{+}^{2, 0}$. It follows from Lemma 2.1 that

and

From the first equation of system (1.3), we have

In view of (3.15), (3.17) and large numbers for martingales, we have

which follows from (2.6) of Lemma 2.3 that

Then,

From the second equation of system (1.3), we have

By (3.16), (3.19) and the strong large numbers for martingales, we have

Hence,

Therefore,

By combing (3.18) and (3.20) we obtain

For $0 < \eta < \Theta < \kappa < \infty$, Hölder's inequality yields that

By Fatou's lemma, we get

Applying (2.2), we obtain

Choosing sufficiently small $\eta$ and sufficiently large $\kappa$. It follows from (3.22) and (3.23) that

where $\mathcal{B} = \{(x, y): y \geq \eta, x+y \leq \kappa\}$. By the invariance of $\mathcal{M} = \{x \geq 0, y > 0\}$ under system (1.3), we consider the process $(S(t), x(t))$ on $\mathcal{M}$. Obviously, $(S(t), x(t))$ is a Feller Markov process. Thus, (3.24) and compactness of $\mathcal{B}$ in $\mathcal{M}$ implies that there is a invariant probability measure $\pi$ on $\mathcal{M}$. Then, by similar arguments of [33,Theothem 2.4] we can derive that $\pi$ is a unique invariant probability measure of $(S(t), x(t))$ on $\mathbb{R}_{+}^{2, \circ}$ and the strong of large numbers holds. The proof is complete.

Remark 3.2. Theorem 3.2 shows that when $\lambda > 0$, the positive solution of system (1.3) has a ergodic stationary distribution. Biologically, it implies that microorganism population is persistent. $

4.   Numerical simulations

In this section, we will introduce several numerical examples to illustrate our theoretical results. For the numerical simulation, we use Milstein's higher order method ^[36] to obtain the discretization equations of system (1.3)

where the time increment $\triangle t > 0$, $\xi_k$ are the independent Gaussian random variables which follow the distribution $N(0, 1)$ for $k = 1, 2, \ldots, n$.

Example 4.1. Consider the small noises. Without loss of generality, let $\sigma_{1} = \sigma_{2} = 0.08$. Set $S^{0} = 1.4$, $D = 0.2$, $m = 0.4$, $a = 0.45$, $\delta = 0.3$, $\gamma = 0.03$, $b\gamma = 0.02$. By calculation, we have $\lambda\approx0.0135 > 0$, which satisfies the condition of Theorem 3.2. Figure 2a, b show time series diagrams of the stochastic model (1.3) and the corresponding deterministic model. Figure 2c, d depict histograms of $S$ and $x$, respectively; Figure 2e shows the phase diagram of the stochastic model (1.3). Biologically, Figure 2 shows that the microorganism population will be persistent with lower noise intensity.

Although the theorem shows the existence of the stationary distribution $\pi$ of system (1.3) with small noise, its analytic formula determined by Kramers-Fokker-Planck (KFP) equation appears to be technically a very difficult problem. This means that it is not easy to obtain the threshold value of noise intensity containing (or causing) microbial persistence (or extinction). Under this circumstance, if $\sigma_1 = \sigma_2$, an asymptotic of the the stationary distribution in Gaussian form can be proposed to estimate the noise intensity ^[37]:

where $z = (S, x)$, $W$ is the stochastic sensitive function (SSF) of deterministic equilibrium $E_* = (0.6088, 0.2119)$ determined by the matrix equation:

where

Hence, the confidence ellipse of system (1.3) can be presented by SSF technique ^[37]:

where $P$ is a fiducial probability. By calculation, we obtain that

and the corresponding confidence ellipse equation is

Further if the fiducial probability $P = 0.95$, then the confidence ellipse for model (1.3) is shown in Figure 2e. One can see from Figure 2e that the random states of system are distributed around the corresponding deterministic equilibrium, and belong to the interior of the confidence ellipse with large probability $0.95$. $

Example 4.2. Consider the large noise $\sigma_{2}$. Let $\sigma_{2} = 0.55$ and keep other parameters the same as in Example 4.1. Then $\lambda\approx-0.1353 < 0$, the condition of Theorem 3.2 is satisfied. Figure 3 depicts the solutions $S(t)$ and $x(t)$ of system (1.3) and the corresponding deterministic model. We can see that microorganism population go to extinction in the stochastic system (1.3) while it is still persistent in the deterministic model. Thus, larger white noise $\sigma_2$ is harmful to microorganism population.

Moreover, it can be seen from the expression of $\lambda$ that $\sigma_{1}$ also has an important effect on microorganism population. From Figure 4 we can see that the peak of the density for $\nu_0$ becomes smaller with the increase of the noise $\sigma_1$, Resonating with the convexity of $p(S) = \frac{mS}{a+S}$, $\lambda$ as a function of $\sigma_1$ is not a monotonic function, while $\lambda_s$ is a decreasing function about $\sigma_1$. Remark 2.2 points out for the convex function $p(S) = \frac{mS}{a+S}$ and $\gamma = 0$, $\mathbb{E}_{\nu_0}[p(S)] > p(S^0)$. Hence, there are parameter sets such that $\lambda > 0 > \lambda_s$. If $\sigma_2 = 0$, $\sigma_1\neq0$, then $\lambda > \lambda_s$. $

Example 4.3. Consider the effect of the instantaneous nutrient recycling on microorganism. Let $S^{0} = 2$, $D = 0.1$, $m = 0.6$, $a = 0.1$, $\delta = 0.3$, $\gamma = 0.4$, $\sigma_1 = 0.08$ and $\sigma_2 = 0.07$. Take $b = 0, \ 0.75$. Although the nutrient recovery rate $b$ is not a key factor affecting the persistence and extinction of microorganism population, numerical simulations show that $b$ may affect the concentration of microorganism. From Figure 5 we can see that the concentration of microorganism population will become larger with the increase of $b$.

5.   Conclusions

In this paper, we proposed a stochastic chemostat model with instantaneous nutrient recycling. Then, we constructed a key threshold value $\lambda$ which enabled us to obtain the asymptotic behavior of the system. Concretely, if $\lambda < 0$, then $\mathbb{P}\left\{\lim \limits_{t \rightarrow \infty} \frac{\ln x(t)}{t} = \lambda\right\} = 1$, biologically, microorganism population will go extinct with exponential rate; if $\lambda > 0$, there exists an invariant probability measure such that the distribution of $(S(t), x(t))$ approximate to this invariant probability at infinite time, biologically, microorganism population is persistent. Although it failed to give an explicit analytic formula of the invariant probability, we constructed the stochastic sensitivity matrix by using the SSF technique to describe the scatter of trajectories in the neighborhood of stable equilibrium. It should be emphasized that if $\gamma = 0$, then system (1.3) becomes the (1.2) with the case $p(S) = \frac{mS}{a+S}$. We propose a new method to define the threshold value $\lambda$ instead of $\lambda_{s}$ in ^[28]. Moreover, different from ^[28], the condition in Theorem 3.1 (Theorem 3.2) is sufficient and almost necessary condition of the exponential extinction (persistence) for the model. However, the case of $\lambda = 0$ remains unsolved. We will continue this interesting work with the concluding of this paper.

Acknowledgments

This work was supported by Basic Research Project (Free Exploration) of Shanxi Province (No.20210302124256), Shanxi Province Department of Finace under Grant No. Z24179 and Talent Special Project of Research Project of Guangdong Polytechnic Normal University under Grant No.2021SDKYA053. We thank the three referees for the careful reading of our paper and all of the insightful suggestions and comments that greatly improved the presentation of the paper.

Conflict of interest

The authors declare that they have no conflict of interest.