Long-time behaviors of two stochastic mussel-algae models

1.   Introduction

A nail-cap-sized zebra mussel was first discovered in the US waters in 1988 and has a powerful reproductive capacity ^[1]. The invasion of zebra mussel has caused great inconvenience to people, such as blocking pipes, polluting water sources, crowding out local species and causing serious economic losses. According to an estimation from the Center for Invasive Species Research at UC Riverside ^[1], the US spends as much as 500 million dollars every year to manage mussel in the Great Lakes. As a result, many biologists, ecologists, and mathematicians have studied the invasion of mussel from different perspectives.

The growth and survival of mussel depend heavily on the availability of food sources for algae. A lot of literatures have revealed that the food supply of algae can limit mussel intake ^[2,3,4]. In order to uncover the relationships between mussel and algae, Koppel et al. ^[5] proposed a diffusive mussel-algae model, considering the corresponding nondiffusive form:

where $M(t)$ and $A(t)$ respectively denote the size of mussel and algae, $\beta$ represents the conversion rate of ingested algae to mussel production, $c$ is the consumption constant, $\mu$ is the maximal per capita mussel death rate, $k$ stands for the value of $M(t)$ at which mortality is half maximal, $A_{up}$ denotes the concentration of algae in the upper water layer, $f$ describes the rate of exchange between the lower and upper water layers, $h$ is the height of the lower water layer. All the parameters are positive.

The research on model (1.1) has attracted much attention. For example, based on model (1.1), Koppel et al. ^[5] analyzed the scale-dependent feedback and regular spatial patterns of young mussel beds, and uncovered that the self-organization patterns would affect the emergent properties of ecosystems in large-scale space. Cangelosi et al. ^[6] established a mussel-algae model with Turing patterns and carried out a series of stability analyses. Song et al. ^[7] dissected Turing-Hopf bifurcation of model (1.1) with reaction-diffusion. Similar diffusive models to study the spatial dynamics of mussel-algae can be found in ^[8,9,10,11]. In addition, quite a few researchers pay attention on the control of mussel-algae. A model describing mussel bed appearance was proposed in ^[12] to explore the habitat suitability analysis for littoral mussel beds in the Dutch Wadden Sea. The effects of 19 macroalgal species on the settlement and metamorphosis of the mussel were investigated ^[13].

Considering that the growth of mussel is affected by intraspecific competition, we transform model (1.1) into the following model:

where $a$ is the intraspecific competition strength of mussel and positive. Other parameters are defined in the same as in model (1.1).

Note that the above studies are all deterministic models. However, environmental uncertainties are ubiquitous in aquatic ecosystems, the populations are always inevitably influenced by environmental noises, which is a momentous element in ecosystems ^[14]. Environmental stochasticity may involve water temperature, noise, salinity, depth and predators, which might affect the growth and evolution of the populations. Accordingly, stochastic models are usually more realistic, and it is essential to bring environmental stochasticity into model (1.2). Quite a few existing literatures focus on this and obtain excellent results, e.g., survival analysis ^[15], asymptotic stability ^[16], stationary distribution ^[17], optimal harvesting ^[18,19] and so on. However, as we know, a very little bit of work has been done with stochastic mussel-algae models, especially the corresponding stochastic version of model (1.2).

For $M(t)$ and $A(t)$ in model (1.2), given $\Delta t > 0$ is a fixed step size. Define $\Gamma^{\Delta t}(p\Delta t) = (M^{\Delta t}(p\Delta t)$, $A^{\Delta t}(p\Delta t))$, $p = 0, 1, 2, \ldots$. Let a normal distribution random variable sequence $\{\Theta_i^{\Delta t}(p)\}_{p = 0}^\infty$ satisfy $\mathbb{E}[\Theta_i^{\Delta t}(p)] = 0$, $\mathbb{E}[\Theta_i^{\Delta t}(p)]^2 = \sigma_i^2\Delta t$, $i = 1, 2$, where the constants $\sigma_1^2$ and $\sigma_2^2$ reflect the size of the random perturbations. In each time period $[p\Delta t, (p+1)\Delta t]$, we hypothesize that $\Gamma^{\Delta t}$ grows in the light of the discrete modification of model (1.2) as well as a stochastic amount $(M^{\Delta t}(p\Delta t)\Theta_1^{\Delta t}(p)$, $A^{\Delta t}(p\Delta t)\Theta_2^{\Delta t}(p))$, then we get

On the basis of ^[20] (Theorem 7.1 and Lemma 8.2), as $\Delta t\rightarrow0$, $\Gamma^{\Delta t}$ converges weakly to the solution of the following stochastic differential equation:

where $B_1(t)$ and $B_2(t)$ are independent standard Brownian motions defined on a complete probability space $(\Omega, F, \{F\}_{t\geq0}, P)$.

The effects of a periodically varying environment are important as populations evolve influenced by external effects, for example, seasonal changes, food supply, living habits and other factors, which changes significantly through the whole life of populations. This idea has found much attention and is incorporated into dynamical models ^{[21,22,23,24,25]}. Till date, to investigate whether these models will exist period solutions or not is still worth noting. Keeping given this fact, model (1.3) may need to be extended into the following periodic version:

where the coefficients $c(t)$, $a(t)$, $\mu(t)$, $k(t)$, $A_{up}(t)$, $f(t)$ are positive continuous $\mathrm{T}$-periodic functions.

It is well known that stability is one of the key topics in mathematical biology. For autonomous stochastic population models, scholars are concerned with the stable "stochastic positive equilibrium"---stationary distribution. For periodic stochastic population models, positive periodic solution is an attractive concept. To the best of our knowledge, however, both the stationary distribution of model (1.3) and the existence of periodic solution of model (1.4) have not been considered. The objectives of this paper are to test these two issues. The rest arrange of this paper is as follows. In the next section, the existence and uniqueness of the global positive solution are testified. In Section 3, the extinction, nonpersistent in the mean and weak persistence of model (1.3) are probed. Section 4 provides the conditions under which model (1.3) possesses a unique ergodic stationary distribution. In Section 5, we explore the existence of $\mathrm{T}$-periodic solution of model (1.4). To illustrate the theoretical findings, some numerical simulations are given in Section 6. A few biological meanings of conditions and results are discussed to end Section 7.

2.   Existence and uniqueness of the global positive solution

Theorem 2.1. For arbitrary initial data $(M(0), A(0))\in \mathbb{R}_{+}^2$, model $(1.3)$ has a unique global positive solution with probability one.

Proof. Recalling model $(1.3)$, assign $M(t) = e^{\widetilde{M}(t)}$, we obtain

with $(\widetilde{M}(0), A(0)) = (\ln M(0), A(0))$. One can see that the coefficients of model (2.1) obey the locally Lipschitz continuous conditions, as a result, it possesses a unique solution $(\widetilde{M}(t), A(t))$ on $[0, \tau_e)$, where $\tau_e\leq+\infty$. Accordingly, model (1.3) possesses a unique positive solution $(M(t), A(t)) = (e^{\widetilde{M}(t)}, A(t))$ on $[0, \tau_e)$. To finish the proof, we only need to testify that $\tau_e = +\infty$ a.s. Choose an integer $n_0 > 0$ which obeys that $1/n_0\leq M(0), A(0)\leq n_0$. For every $n\geq n_0$, define

Set $\tau_\infty = \lim_{n\rightarrow +\infty}{\tau_n}$. As a result, $\tau_\infty\leq \tau_e$. Now we only need to testify that $\tau_\infty = +\infty$. If it is not true, then one can find two constants $T > 0$ and $\epsilon\in(0, 1)$ such that $P\{\tau_\infty \leq T\} > \epsilon$. As a result, one can set an integer $n_1\geq n_0$ which satisfies

Define

Taking advantage of Itô's formula, one has

where

Integrating both sides of Eq (2.3) from $0$ to $\tau_n\wedge T$ yields

Taking expectation on both sides results in

Set $\Omega_n = \{\tau_n\leq T\}$ for $n\geq n_1$. According to Eq (2.2), $P(\Omega_n)\geq\epsilon$. For any $\theta\in\Omega_n$, at least one of $M(\tau_n, \theta), A(\tau_n, \theta)$ equals to $n$ or $1/n$. Thus, we derive

Therefore, Eq (2.4) implies that

where $1_{\Omega_n}$ denotes the indicator function of $\Omega_n$. Letting $n\rightarrow +\infty$ causes the contradiction:

This finishes the proof.

Remark 2.1. Similar to the proof of Theorem 2.1, one can testify that model $(1.4)$ has a unique global positive solution with probability one, and the details are left out.

3.   Extinction and persistence of model (1.3)

Lemma 3.1. Given $(M(0), A(0)) \in\mathbb{R}_{+}^2$, model $(1.3)$ admits $\limsup_{t\rightarrow +\infty}{[M(t)+\beta hA(t)]} < +\infty$ and

Proof. Denote $Z(t) = M(t)+\beta hA(t)$. From model (1.3), we have

where $\delta = \min\{2a, f\} > 0$. Consider

The solution of model (3.1) is

where

is a local martingale satisfying $N(0) = 0$ a.s. Thus

where

with $Q(0) = P(0) = 0$. Clearly, $Q(t)$ and $P(t)$ are continuous increasing functions. By ^[26], we have $\lim_{t\rightarrow +\infty}{Y(t)} < +\infty$ a.s., then by stochastic comparison theorem, one has $\limsup_{t\rightarrow +\infty}{Z(t)} < +\infty$ a.s.

Let $N_1 = \int_0^t\sigma_1M(s)\mathrm{d}B_1(s)$ and $N_2 = \int_0^t\sigma_2A(s)\mathrm{d}B_2(s)$. Through calculation, we obtain

then

In light of ^[27], $\lim_{t\rightarrow +\infty}t^{-1}N_1(t) = 0$ a.s. Similarly, we have $\lim_{t\rightarrow +\infty}t^{-1}N_2(t) = 0$.

Theorem 3.1. If $\lambda_0 = \beta cA_{up}-\mu-\sigma_1^2/2 < 0$ and $a > \mu/k$, then $M(t)$ is extinct a.s.

Proof. We deduce from model (1.3) that

which implies that

where

satisfying $\lim_{t\rightarrow +\infty}{\varphi_1(t)/t} = 0$. In light of Eq (3.2), we have

By the first equation of model (1.3) and using Itô's formula, we obtain

then together with Eq (3.3), one has

Since the strong law of numbers implies that

Thus, it follows from Eqs (3.4) and (3.5) that

which implies the required assertion.

Theorem 3.2. If $\lambda_0 = 0$ and $a > \mu/k$, then $M(t)$ is nonpersistent in the mean a.s., namely, $\lim_{t\rightarrow +\infty}t^{-1}\int_0^tM(s)\mathrm{d}s = 0$ a.s.

Proof. Let $\rho > 0$ be a constant which satisfies $\rho < a-\mu/k$. From Eq (3.4), we have

Note that for any $\varepsilon > 0$, there is $T > 0$ such that for $t\geq T$,

Substituting Eq (3.7) into Eq (3.6), we have

Set $\varrho(t) = \int_0^tM(s)\mathrm{d}s$, then we get

Hence for $t > T$, we have

Integrating this inequality from $T$ to $t$, one can derive that

That is,

Taking the logarithm of both sides of Eq (3.8) results in

Note that $\varrho(t) = \int_0^tM(s)\mathrm{d}s$, then one can obtain that

Applying L'Hospital's rule leads to

It then follows from the arbitrariness of $\varepsilon$ that $\limsup_{t\rightarrow +\infty}t^{-1}\int_0^t M(s)\mathrm{d}s\leq0$. This proof is complete.

Theorem 3.3. If $\lambda_0 > 0$, then $M(t)$ is weakly persistent a.s., namely, $\limsup_{t\rightarrow +\infty} M(t) > 0$ a.s.

Proof. We first testify that

By Itô's formula,

Integrating the both sides from $0$ to $t$, we have

where $W(t) = \int_0^te^s\sigma_1\mathrm{d}B_1(s)$ is a local martingale with the quadratic form

By the exponential martingale inequality (see ^[26] on page 44), for arbitrary positive constants $T_0$, $\iota$ and $\nu$, one has

Choose $T_0 = \vartheta r$, $\iota = e^{-\vartheta r}$ and $\nu = \varpi e^{\vartheta r}\ln r$, then we obtain

where $\varpi > 1$, $\vartheta > 0$. By the Borel-Cantalli lemma (see ^[26] on page 7), for almost all $\zeta\in\Omega$, there exists a $r_0(\zeta)$ such that for $r\geq r_0(\zeta)$,

Combining Eq (3.10) with Eq (3.11), we obtain

Since $\ln M(t)+\beta cA(t)-aM(t)-\frac{1}{2}\sigma_1^2+0.5e^{t-\vartheta r}\sigma^2_1$ is bounded, for any $0\leq s\leq\vartheta r$, there is a constant $C$ independent of $r$ such that

Substituting Eq (3.13) into Eq (3.12), we obtain

Dividing the both sides of Eq (3.14) by $e^t$ leads to

Consequently, if $\vartheta (r-1)\leq t\leq\vartheta r$ and $r\geq r_0(\zeta)$, then one can observe that

which is the needed assertion Eq (3.9) by letting $r\rightarrow +\infty$.

Now let us testify $\limsup_{t\rightarrow +\infty}M(t) > 0$ a.s. If not, then we denote $S = \{\limsup_{t\rightarrow +\infty}M(t) = 0\}$, $P(S) > 0$. In light of Eq (3.4), one has

For all $\zeta\in S$, we have $\lim_{t\rightarrow +\infty}{M(t, \zeta)} = 0$, and the law of large numbers for local martingales indicates that $\lim_{t\rightarrow +\infty}{\frac{1}{t}\int_0^t\sigma_1\mathrm{d}B_1(s)} = 0$. Thus we have $\limsup_{t\rightarrow +\infty}t^{-1}\ln M(t, \zeta) = \lambda_0 > 0$. By Eq (3.9), a contradiction arises.

4.   Ergodic stationary distribution (ESD) of model (1.3)

Now we dissect the stationary distribution for model (1.3) by taking advantage of Has'minskii's results ^[28]. Denote by $X(t)$ a time-homogeneous Markov process in $\mathbb{R}^n$ which obeys

Let $I(x) = (a_{ij}(x))$ be the diffusion matrix of $X(t)$, where

For any $C^2$- function $V_1(x)$, define

Lemma 4.1. If there is a bounded domain $U\subset\mathbb{R}^d$ with regular boundary such that (^[28])

● there is a positive number $\Lambda$ which obeys

● there is a nonnegative $C^2$- function $V_2$ such that $LV_2(x) < -1$ for any $x\in\mathbb{R}^d\setminus U$,

then $X(t)$ admits a unique ESD.

Define

Theorem 4.1. If $\lambda_0 > 0$ and $\sigma_2^2$ is sufficiently small such that

then model $(1.3)$ admits a unique ESD.

Proof. Considering the function $V_3(M, A) = -m_1\ln M-m_2\ln A$, and $m_1$, $m_2$ are positive constants to be chosen later, we obtain

where

Define

we have

Notice that $h^2\beta^2f-h\beta\sigma_2^2 > 0$, hence

In addition,

where $D_3 = f+\sigma_2^2/2$.

Therefore,

Now define $V_{5}(M, A) = \lambda V_3(M, A)+V_4(M, A)$, where $\lambda > 0$ is sufficiently large. Hence,

where $U_{q_1} = (\frac{1}{q_1}, q_1)\times(\frac{1}{q_1}, q_1)$, $q_1$ is a sufficiently large number. Notice that $V_{5}(M, A)$ is continuous. Thus $V_{5}(M, A)$ has a minimum point $(M_0, A_0)$ in $\mathbb{R}_{+}^2$. Define

Thus, we can get

Define a bounded close set:

where $0 < \varepsilon < 1$ is sufficient small. We can split $\mathbb{R}_{+}^2\setminus U$ into the following four ranges,

Case 1. If $(M, A)\in U_1$, then we have

Case 2. If $(M, A)\in U_2$, then one can see that

where

Case 3. If $(M, A)\in U_3$, then one has

where

Case 4. If $(M, A)\in U_4$, then we obtain

In $\mathbb{R}_{+}^2\setminus U$, let $\varepsilon$ be sufficiently small which satisfies

It follows from Eqs (4.1)–(4.5) that

The diffusion matrix of model (1.3) has the form

Choosing $\Lambda = \min_{(M, A)\in U_q}\{\sigma_1^2M^2, \sigma_2^2A^2\} > 0$, we have

According to Eqs (4.6), (4.7) and Lemma 4.1 that we complete the proof.

5.   Existence of $\mathrm{T}$-periodic solution of model (1.4)

Consider the stochastic periodic equation

where $v(t)$ and $g(t)$ are $\mathrm{T}$-periodic functions in $t$.

Lemma 5.1. If there exists a function $V_{7}(t, x)\in C^2$ which is $\mathrm{T}$-periodic and satisfies the conditions (^[28])

● $\inf_{|x| > \Theta}V_{7}(t, x)\rightarrow \infty$ as $\Theta\rightarrow \infty$,

● $LV_{7}(t, x)\leq-1$ on the outside of some compact set,

then there exists a periodic solution to Eq $(5.1)$.

Define

Define $\langle g\rangle_\mathrm{T} = \frac{1}{\mathrm{T}}\int _ 0^\mathrm{T} g(s)\mathrm{d}s$, where $g(t)\in[0, \infty)$ is an integrable function.

Define $g^u = \max_{t\in [0, +\infty)}{g(t)}$, $g^l = \min_{t\in [0, +\infty)}{g(t)}$, where $g(t)\in[0, +\infty)$ is a bounded function.

Theorem 5.1. If $R_1 > 1$ and $(\sigma^2_2)^u < h\beta f^l$, then model $(1.4)$ admits a positive $\mathrm{T}$-periodic solution.

Proof. Define

and $b_1$, $b_2$ are positive constants to be chosen later. By Itô's formula, we have

where

Let $\bar{\omega}(t)$ be the solution of the following equation

Then $\bar{\omega}(t)$ is a $\mathrm{T}$-periodic function. On the basis of Eqs (5.2) and (5.3), we can obtain

where $\alpha_1 = 2(R_1-1) > 0$.

Define

Applying Itô's formula, one has

where

Similarly, one deduces

According to Eqs (5.5) and (5.6) one can get

Define

where $H$ is positive constant. Clearly,

where $U_{q_2} = (\frac{1}{q_2}, q_2)\times(\frac{1}{q_2}, q_2)$, $q_2$ is a sufficiently large number. Combining with Eqs (5.4) and (5.7), we have

Define a bounded close set:

where $0 < \varepsilon < 1$ is a sufficient small number. We divide $\mathbb{R}_{+}^2\setminus \mathcal {K}$ into the following four ranges

Case 1'. If $(M, A)\in \mathcal {K}_1$, then we get

Case 2'. If $(M, A)\in \mathcal {K}_2$, then we have

where

Case 3'. If $(M, A)\in\mathcal{K}_3$, then we derive

where

Case 4'. If $(M, A)\in\mathcal{K}_4$, then one has

In the set $\mathbb{R}_{+}^2\setminus\mathcal{K}$, we choose $\varepsilon$ sufficiently small such that

It then follows from Eqs (5.8)–(5.15) that $LV_{10}(t, M, A) < -1$ for all $(M, A)\in\mathbb{R}_{+}^2\setminus\mathcal{K}$.

6.   Discussion and simulations

In this section, we take advantage of some real data (see Table 1) and the Euler-Maruyama method ^[34] to illustrate the above results. For model (1.3), we pay attention to the discretization equation:

where $\zeta_{1n}$, $\zeta_{2n}$ mean independent Gaussian random variable.

From Theorems 3.1–3.3, one can find out that under the assumption $a > \mu/k$, $\lambda_0$ is sufficient conditions determining the persistence and extinction of the mussel species. More precisely, under the assumption $a > \mu/k$, if $\lambda_0 < 0$, then the mussel species dies out; if $\lambda_0 = 0$, then the mussel species is nonpersistent in the mean. If $\lambda_0 > 0$, then the mussel species is weakly persistent. Note that $\lambda_0 = \beta cA_{up}-\mu-\frac{1}{2}\sigma_1^2$, which suggests that white noise can greatly influence the survival of mussel: when the intensity of the noise is large enough, it could make mussel become extinct.

Theorem 4.1 suggests that if $\lambda_0 > 0$ and $\sigma_2^2$ is sufficiently small such that

then model (1.3) possesses a unique ESD on $\mathbb{R}_{+}^2$. This ESD can be used to estimate the outbreak possibility of mussel.

Figure 1(a)–(c) characterize the persistence and extinction of the mussel species in model (1.3) with different $\sigma_1$. We choose $\sigma_2 = 0.3$, initial value $M(0) = 0.1$ and $A(0) = 0.12$. Figure 1(a) is with $\sigma_1 = 0.3$, which reflects that the mussel species dies out with probability one; Figure 1(b) is with $\sigma_1 = 0.1$, which suggests that the mussel species is nonpersistent in the mean; Figure 1(c) is with $\sigma_1 = 0.01$, which shows that the mussel species is weakly persistent. Comparing Figure 1(c) with Figure 1(a), one can observe that with the increasing value of $\sigma_1$, the mussel species tends to go to extinction. In other words, large noise could lead to the extinction of mussel.

Figure 2 plots the probability density function (PDF) of the stationary distribution of model (1.3) with $\sigma_1 = 0.005$ and $\sigma_2 = 0.0048$.

For model (1.4), we focus on the following discretization equation:

where $\zeta_{1n}$, $\zeta_{2n}$ mean independent Gaussian random variable.

Theorem 5.1 provides sufficient conditions (i.e., $R_1 > 1$ and $(\sigma^2_2)^u < h\beta f^l$) under which model (1.4) admits a positive periodic solution. This offers some insightful understanding on how environmental fluctuations affect the survival of mussel and algae.

The initial values $M(0) = 0.1$ and $A(0) = 0.12$ are kept the same as in Fig. 1, and we choose $\beta = 0.2$, $h = 1$, $A_{up}(t) = 10+0.1\sin(\pi t)$, $a(t) = 0.2+0.1\sin(\pi t)$, $f(t) = 1+0.1\sin(\pi t)$, $c(t) = 0.5+0.1\sin(\pi t)$, $\mu(t) = 0.1+0.1\sin(\pi t)$, $k(t) = 10+0.1\sin(\pi t)$, $\sigma_1(t) = 0.03+0.01\sin(\pi t)$, $\sigma_2(t) = 0.04+0.01\sin(\pi t)$. It follows from Theorem 5.1 that model (1.4) has a T-periodic solution, see Figsure 3(a), (c). Moreover, Figsure 3(a), (b) show that $M(t)$ and $A(t)$ fluctuate periodically, that is, mussel and algae will not die out. In particular, the effect of environmental noises can be easily found by comparing Figsure 3(a), (c) with Figsure 3(b), (d).

7.   Conclusions

Understanding the effect of random perturbations on the evolution of the mussel is useful for managing this species. This paper proposed two stochastic mussel-algae models (one is autonomous, and the other is non-autonomous) to test the effect of environmental fluctuations on the evolution of mussel. For the autonomous model, the critical value between extinction and weak persistence was obtained. In addition, sufficient conditions for the existence of an ESD were established. For the non-autonomous model, the existence of a positive periodic solution was examined. Some vital impacts of environmental fluctuations on the evolution of mussel were uncovered.

In comparison with the existing papers, this research has the following contributions:

●   Our models consider the environmental fluctuations which are more reasonable. Actually, to the best of our knowledge, this research is the first attempt to dissect the stochastic mussel-algae models.

●   We obtain the critical value between extinction and weak persistence for the mussel, and uncover that environmental fluctuations can significantly affect the extinction/persistence of the mussel.

●   We give some conditions under which model (1.3) has an ESD. This ESD is useful to estimate the outbreak probability of the mussel.

●   We provide sufficient conditions for existence of a positive periodic solution of model (1.4). This positive periodic solution is helpful for the understanding how environmental fluctuations affects the survival of mussel and algae.

Some studies on mussel-algal models are worth further investigations. Actually, Theorem 3.1 and Theorem 3.2 have an assumption $a > \mu/k$, what happens when $a < \mu/k$ is still unclear. In addition, Theorem 4.1 testifies that if $\lambda_0 > 0$ and

then model (1.3) possesses a unique ESD on $\mathbb{R}_{+}^2$. It is interesting to relax the restriction on $\sigma_2^2$. Finally, one may put forward some more realistic and meaningful models, such as considering the effects of Lévy jump ^[35,36], impulsive perturbations ^[37,38], time delay ^[39,40] or fractional order ^[41,42]. We will leave these for future works.

Acknowledgments

The authors thank the editor and referees for their careful reading and valuable comments. This work is supported by the National Natural Science Foundations of China (Nos.11771174, 11871201).

Conflict of interest

The authors declare there is no conflict of interest.