A stochastic mussel-algae model under regime switching

1.   Introduction

The invasion of alien species has caused significant economic and ecological problems in a diverse array of ecosystems. One of the most instructive examples was the zebra mussel, a highly invasive bivalve in the Great Lakes in 1986, which has high fecundity and strong ability, not only causing ecological modifications and economic damages, but also bringing great inconvenience to human life ^[1,2], such as blocking pipes, polluting water sources and crowding out native species. In recent years, governments in many countries have invested lots of capital in mussel control and damage repair ^[3], and an increasing number of scholars have begun to study this hot issue. Many significant existing literatures indicate that the growth and survival of mussel are mainly dependent on algae, which means that the food supply of algae will limit the intake of mussel (see ^[4,5,6,7]). To date, quite a few ecologists and applied mathematical researchers have paid close attention to the development of mussel with various types of mathematical mussel-algae models, and numerous results have been presented about mussel-algae models on account of their extensive applications ^{[8,9,10,11,12]}. For example, Turing pattern ^[8], Turing-Hopf bifurcation ^[9,10], wavetrain solution ^[11] and so on.

Reviewing many of the above investigations, we can see that deterministic model can simplify complex system, which is conducive to theoretical research and analysis. However, when the external disturbance is large, the population number may be significantly reduced, the law of large numbers is not on trial, and the deterministic model is not objective. If the noises are ignored, the model will be inaccurate. So in order to describe the species in the real world more objectively, the environmental fluctuations should be taken into account^{[13,14,15,16]}. It has been found that environmental stress, salinity and so on have a remarkable effect on the microbiota associated with zebra mussel, therefore influence the survival analysis of mussel ^[17]. Usually, one of the most common and indispensable environmental fluctuations in an ecosystem, white noises, which can be used to describe these continuous and small fluctuations, have recently become quite attractive (see ^{[18,19,20,21,22]}).

However, there are several moderate fluctuations which can not be characterised by white noises ^[23], let alone some sudden environmental abrupt fluctuations. Take temperature as an example, organisms have a certain average range of optimal growth temperature, and the change of temperature during individuals of the same species or different physiological periods of the same species, may have an effect on population separation and behaviour differences ^[24]. With the exception of birds and mammals, which are thermostats, all other animals are poikilotherms, namely, their body temperature always varies in response to the changes in external temperature ^[25]. When the effect of temperature change exceeds the threshold of environment or ecosystem, results can range from mild individual death to species turnover or ecosystem change ^[25]. There is one vivid example in Finland during the 1950s, the pikeperch population crumbled just because of a trifling cooling together with the grown fishing mortality ^[26], manifesting that even a slight change of mean temperature, like $0.5-1.0^ {\circ}C$, may appreciably affect populations living nearby ^[27]. Although the temperature in annual specific season is relatively constant in general, the stochastic effects on temperature still may lead to tremendous variations in larval fish survival ^[28], and hence none of the white noises or the period describing form may be applicable for this fluctuation. Coincidentally, it has been observed that temporal patterns in the life cycle of the zebra mussel vary on increasing temperature in experiment ^[29]. Thus, when the system jumps back and forth caused by moderate fluctuations or sudden environmental distributions, a further step to investigate telephone noises, another common noises type, can be used to elucidate a switching between different (two or more) environmental regimes, is noteworthy ^[30]. Because this switching is memoryless and the waiting time of two different switches obeys an exponential distribution ^[31], using state switching of a continuous-time Markov chain to drive the changes in the main parameters of system can satisfy our needs. Numerous literatures focus on these issues, for example ^[32,33,34] and the reference therein.

The necessary environmental fluctuations exist everywhere in various forms. In recent years, as two most common and representative noises, white noises and telephone noises are widely used in stochastic population systems, which has been studied with the efforts of many researchers, see ^{[35,36,37,38,39]}. Especially, a stochastic phytoplankton allelopathy model formulated by Zhao et al. ^[35], which derives the stationary distribution and its statistical characteristics (i.e., the mean and variance). Considering a Holling-II stochastic predator-prey model, Jiang et al. ^[36] established a set of sufficient criteria for the strong persistence in the mean, extinction as well as stationary distribution. Very recently, Cai et al. ^[37] dealt with the long-time properties of a stochastic delay foraging arena predator-prey model, including stochastic ultimate boundedness, pathwise estimation and extinction.

As is well known, there are few papers related to stochastic mussel models with white noises and telephone noises. There is a doubt whether or not we can derive that the extinction and the existence of unique stationary distribution of such regime-switching differential system. Fortunately, the answer is affirmative and the related detailed proofs will be given later.

Focusing on two essentials, the extinction of mussel species and ergodic stationary distribution, the paper is bent on exploring the long-time properties of model (2.4). The present paper has the following structure: in Section 2 we introduce a stochastic mussel-algae model simultaneously affected by white noises and telephone noises. Section 3 ensures the existence and positivity of global solution. Some dynamic behaviors about extinction are carried out in Section 4 and ergodic stationary distribution are proved in Section 5, respectively. Several numerical examples are set to demonstrate our mathematical results in Section 6. Section 7 ends with some conclusions and future directions. A necessary theory listed in Appendix. Last but not least, several literatures inspired the present work are listed.

2.   Model and preliminaries

For the final export of the model we will discuss, let us first introduce the following nondiffusive mussel-algae model which corresponds to a nondiffusive version proposed by Koppel et al.^[12].

The biological meanings of all variables and parameters concerned model (2.1) are explained as below.

With the idea of the influence of intraspecific competition, the above model (2.1) is recently extended by Zhou et al. ^[40] to the following deterministic version

where $g$ represents the intraspecific competition strength of mussel.

The white noises are taken into account model (2.2) by Zhou et al. ^[40] and a following stochastic version is derived

where $B_1(t)$ and $B_2(t)$ are the standard independent Brownian motions, $\sigma_1^2$ and $\sigma_2^2$ represent their intensities. Thanks to the analysis of the dynamical behaviors of model (2.3), not only the survival of mussel as well as ESD (ergodic stationary distribution) are obtained, but an interesting result is discovered, namely the lower amplitude noises could cause mussel species outbreak, while the higher amplitude noises could make mussel species become extinct.

As an extension of model (2.3), motivated by the above studies we further introduce the telephone noises into model (2.3), and suppose that Markov chain $r(t)$ take values in state space $\mathbb{S} = \{1, 2, ..., N\}$ controlling the switching between the environmental regimes ^[41]. Then the resulting stochastic mussel-algae model simultaneously affected by white noises and telephone noises can be expressed as

$r(0)\in\mathbb{S}$, $X(0) > 0$, $Y(0) > 0$. Here, suppose that Markov chain $r(t)$ is independent of the Brownian Motions $B_1(t)$ and $B_2(t)$. The coefficients $f(k)$ are all nonnegative constants for any $k\in\mathbb{S}$, where $f = g$, $b$, $c$, $\beta$, $l$, $\sigma_1$, $\sigma_2$, $Y_{up}$.

For the convenience of subsequent sections, we list the following preliminaries. Define probability space by $(\Omega, \mathcal{F}, \{\mathcal{F}_t\}_{t\geq0}, \mathbb{P})$ with a filtration $\{\mathcal{F}_t\}_{t\geq 0}$ satisfying the usual conditions. Assign $\mathbb{R}_+^n = \{z\in\mathbb{R}^n|z_i > 0, 1\leq i\leq n\}$. Denote right-continuous Markov chain $r(t)\ (t\geq0)$ on the complete probability space whose domain is a finite state space $\mathbb{S}$. Let $\hat{h} = \min_{k\in\mathbb{S}}\{h(k)\}$ and $\check h = \max_{k\in\mathbb{S}}\{h(k)\}$, where $h = (h(1), ..., h(N))$. The generator $\tilde{\Gamma} = (\mu_{ij})_{N\times N}$ is given, for $\triangle > 0, \mu_{ij}\geq 0$, by

where $\mu_{ij}$ implies the transition rate from state $i$ to state $j$, and keeps nonnegative if $i\neq j$ satisfying $\sum_{i = 1}^N\mu_{ij} = 0$. When $r(t)$ is irreducible, then it admits a unique stationary distribution $\pi = \{\pi_1, \pi_2, ..., \pi_N\}$ such that $\pi\tilde\Gamma = 0$ satisfying $\sum_{k = 1}^N\pi_k = 1 \; (\pi_k > 0)$ for any $k\in\mathbb{S}$.

By the help of $r(\cdot)$, define the diffusion process $(q(t), r(t))$ as below:

where $B(\cdot)$ signifies Brownian motion. Assign $\varsigma(\cdot, \cdot):\mathbb{R}^n\times\mathbb{S}\rightarrow \mathbb {R}^n, \mathcal{G}(\cdot, \cdot):\mathbb{R}^n\times\mathbb S\rightarrow \mathbb{R}^{n\times n}$ such that $\mathcal{G}(q, k)\mathcal{G}^ T(q, k) = (d_{ij}(q, k))$. For each $k\in\mathbb{S}$, let $\mathcal{W}(q, k)$ be twice continuously differentiable with respect to $q$, and define the operator $L$ as:

3.   Global positive solution

We here dedicate to proving an essential property of model (2.4).

Lemma 3.1. Given arbitrary initial data $(X(0), Y(0), r(0))\in\mathbb{R}_{+}^2\times\mathbb{S}$, model $(2.4)$ owns a unique global positive solution with probability one.

Proof. Notice the following model

$(\widetilde{X}(0), Y(0)) = (\ln X(0), Y(0))$. It's clear that the coefficients of model (3.1) obey locally Lipschitz continuous conditions. Thus, a unique solution $(\widetilde{X}(t), Y(t))$ on $[0, \tau_e)$ to model (3.1) which belongs to $\tau_e\leq+\infty$ exists. In other words, the presence of a unique local positive solution of model (2.4) is guaranteed on $[0, \tau_e)$, denoted by $(X(t), Y(t)) = (e^{\widetilde{X}(t)}, Y(t))$. The only thing we need to do is to verify $\tau_e = +\infty \; a.s.$ Denote $n_0 > 0$ be enough large such that $1/n_0\leq X(0), Y(0)\leq n_0$. For each $n\geq n_0$, assign

Set $\tau_\infty = \lim_{n\rightarrow +\infty}{\tau_n}$, then $\tau_\infty\leq\tau_e$. Testify $\tau_\infty = +\infty$ can finish this proof. If this assertion is false, then assume that $\tau_e < \infty$. There are two constants $T > 0$ and $\epsilon\in(0, 1)$ such that $P\{\tau_\infty\leq T\} > \epsilon$. Consequently there exists an integer $n_1\geq n_0$ satisfying

Define

from which we conclude, together with Itô's formula, that

where

Integrating the above inequation on $[0, \tau_n\wedge T]$ and taking expectation lead to

Set $\Omega_n = \{\tau_n\leq T\}$ for $n\geq n_1$, which together with Eq (3.2) leads to, $P(\Omega_n)\geq\epsilon$. For any $\theta\in\Omega_n$, at least one of $X(\tau_n, \theta), Y(\tau_n, \theta)$ equals to $n$ or $1/n$. Thereby,

Recalling Eq (3.3) and one has

where $1_{\Omega_n}$ representatives the indicator function of $\Omega_n$. Letting $n\rightarrow +\infty$ gives that

which causes a contradiction. So we get $\tau_\infty = +\infty$ which completes this proof.

4.   Extinction

We now concern with the extinction of mussel species.

Theorem 4.1. Given initial value $(X(0), Y(0), r(0))\in\mathbb{R}_{+}^2\times\mathbb{S}$, the solution $(X(t), Y(t), r(t))$ to model $(2.4)$ satisfy

Particularly, if $\sum_{k = 1}^{N}\pi_k\big(\frac{a\check b}{\hat l}l(k)Y_{up}(k)-\frac{1}{2}\sigma_1^2(k)\big) < 0$ holds, then

Proof. Define

Utilizing the Itô's formula to Eq (4.1), one obtains

Integrating both sides of $dV(X, Y, k)$ on the interval $[0, t]$ and dividing it by $t$ yields

Reordering the above equation, we have

After integrating the first equation of model (2.4) from $0$ to $t$, we divide it by $t$, then

Substituting Eq (4.2) into the above inequality, we further have that

where

Note that $\lim_{t\rightarrow +\infty}t^{-1}\varphi(t) = 0$ in the light of the strong law of large numbers. Then, letting $t\rightarrow +\infty$ gives

which implies the required assertion.

5.   Ergodic stationary distribution

To investigate the ESD of model (2.4) with Markov switching, we recall an important lemma (Hasminskii's Theory ^[42]) in Lemma A.1.

Define

Theorem 5.1. If $R_0 > 1$ and $\hat l > \check\sigma_2^2$, thenmodel $(2.4)$ admits a unique ESD.

Proof. Above all, the condition $\mu_{ij} > 0$ for $i\neq j$ (see Appendix) guarantees the condition (i) of Lemma A.1.

Next, the diffusion matrix of model (2.4) has the following form

Choose $\Lambda = \min_{(X, Y)\in{\mathcal{D}^c\times\mathbb{S}}}\{\sigma_1^2(k)X^2, \sigma_2^2(k)Y^2\} > 0$, where $\mathcal{D} = [\varepsilon, \frac{1}{\varepsilon}]\times[\varepsilon, \frac{1}{\varepsilon}]$, we have

which ensures that the condition (ii) of Lemma A.1 is valid.

At last, we prove that the condition (iii) holds. Consider the function

where

Using the Itô's formula to Eq (5.1) and element inequality $a+b\geq2\sqrt{ab}$, one computes

where

Now, we define $R_0 = \big(R_0(1), R_0(2), \ldots, R_0(N)\big)^{T}$. Due to the generator matrix is $\Gamma$, there is a solution of the Poisson system can be expressed as

such that

Then we have

which combines with the definitions of $m_1$ and $m_2$ leads to

where $D_1 = 2(\sqrt{R_0}-1)$.

We further define

The application of Itô's formula to $V_2$ yields

where $D_2 = \check l+\check\sigma_2^2/2$ and

Since $\hat{l} > \check{\sigma}_2^2$, we have $D_3 < +\infty$.

Now assign

where $\lambda_1 > 0$ is sufficiently large. According to Eqs (5.2)–(5.3), we know that

where $\mathcal{U}_{\ell_1} = (\frac{1}{\ell_1}, \ell_1)\times(\frac{1}{\ell_1}, \ell_1)$, $\ell_1 > 1$ is a sufficiently large number. Continuous function $W(X, Y, k)$ guarantees that it has a minimum point $(X_0, Y_0, k)$ in $\mathbb{R}_{+}^2\times\mathbb{S}$, we further define the following $C^2$-function

Therefore, we can naturally get

Let $0 < \varepsilon < 1$ is sufficient small and define a bounded close set

We can split $\mathbb{R}_{+}^2\setminus\mathcal{U}$ into the following four ranges:

Case 1 If $(X, Y, k)\in \mathcal{U}_1$, then we derive that

Case 2 If $(X, Y, k)\in \mathcal{U}_2$, then one can see that

where

Case 3 If $(X, Y, k)\in \mathcal{U}_3$, then

where

Case 4 If $(X, Y, k)\in \mathcal{U}_4$, then we obtain

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

According to Eqs (5.4)–(5.8), we have

Summarizing the above discussions, the conditions (i)-(iii) of Lemma A.1 have been checked, which means that model (2.4) has a ESD. The proof is completed.

6.   Numerical simulations

Here, it's time to verify our analytical results by some numerical simulations. Give the initial values of model (2.4) by $(X(0), Y(0)) = (0.3, 0.2)$. Suppose that $r(t)$ ranges in state space $\mathbb{S} = \{1, 2\}$ with the generator

and the corresponding stationary distribution $\pi = (\pi_1, \pi_2) = (0.4, 0.6)$.

Example 6.1 We fix $a = 0.2$, $p = 0.1$. When $r(t) = 1$, we choose a set of parameters with $c(1) = 0.02$, $g(1) = 0.01$, $b(1) = 0.1$, $Y_{up}(1) = 1$, $l(1) = 0.4$, $\beta(1) = 150$, $\sigma_1(1) = 0.3$, $\sigma_2(1) = 0.3$. When $r(t) = 2$, we select the other group parameters with $c(2) = 0.3$, $g(2) = 0.02$, $b(2) = 0.2$, $Y_{up}(2) = 1.5$, $l(2) = 0.6$, $\beta(2) = 200$, $\sigma_1(2) = 0.4$, $\sigma_2(2) = 0.4$. By calculating, we get $\sum_{k = 1}^{N}\pi_k\big(\frac{a\check b}{\hat l}l(k)Y_{up}(k)-\frac{1}{2}\sigma_1^2(k)\big) = -0.005 < 0$, which together with Theorem 4.1 shows that mussel $X(t)$ of model (2.4) tends to be extinct, depicted in Figure 1.

Example 6.2 Fix $a = 0.7$, $p = 0.1$. When $r(t) = 1$, select a set of parameters with $c(1) = 0.03$, $g(1) = 0.02$, $b(1) = 0.1$, $Y_{up}(1) = 1.2$, $l(1) = 0.6$, $\beta(1) = 150$, $\sigma_1(1) = 0.05$, $\sigma_2(1) = 0.05$. And choose the other set of parameters with $c(2) = 0.04$, $g(2) = 0.01$, $b(2) = 0.2$, $Y_{up}(2) = 2.2$, $l(2) = 0.7$, $\beta(2) = 200$, $\sigma_1(2) = 0.1$, $\sigma_2(2) = 0.1$ when $r(t) = 2$. We then compute $R_0\approx6.18 > 1$ and $\check\sigma_2^2-\hat l = -0.59 < 0$. As a result, we know from Theorem 5.1 that model (2.4) has a unique ESD. Simulations shown in Figure 2 can confirm this clearly.

7.   Discussion and conclusions

Based on the previous work ^[40], this article focuses on the long-time behaviors of a stochastic mussel-algae model incorporating both white noises and telephone noises. To the best of our knowledge, it is the first attempt to analyze the above stochastic mussel-algae model with two types of environmental noises. We establish two sets of sufficient criteria to verify the extinction and the existence of a unique ESD. And some crucial influences of environmental noises on mussel evolution are revealed. Based on the above theoretical analysis, this paper draws the followings conclusions.

● We can conclude, by Theorem 4.1, that the mussel will be extinct if

After further observation, if the conversion rate of ingested algae to mussel production $a$ and the rate of consumption of the algae by a mussel $b(k)$ are suitably small while the amplitude noises $\sigma_1(k)$ is relatively large, then the mussel species die out (see Figure 1), which means environmental fluctuations can be used as a strategy to control the growth of mussel.

● We have demonstrated, according to Theorem 5.1, that model (2.4) has a unique ESD if $R_0 > 1$ and $\hat l > \check\sigma_2^2$ which can be used to estimate the probability of mussel outbreaks. This fact implies that large conversion rate of ingested algae to mussel production $a$, large rate of consumption of the algae by a mussel $b(k)$, while small maximal per capita mussel death rate $c(k)$ and small noises amplitude $\sigma_{i}(k)(i = 1, 2)$ are advantage for the persistence of species mussel (see Figure 2).

Looking forward, some challenging topics deserve further research studies. We may investigate more complex model by introducing different noises perturbations into mussel-algae model (2.4), such as Ornstein-Uhlebeck process, seim-Markov, Lévy jump^[43]. Taking into account the long-term memory of population systems, further consideration of fractional-order modeling methods has important theoretical and practical implications^{[44,45,46,47,48]}. We will leave these tasks in future.

Acknowledgements

We would like to thank the referees and editor for their careful reading of the original manuscript and their valuable comments and suggestions that greatly improved the presentation of this work. The work is supported by the NNSFs of China (Nos. 11871201, 12101211).

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

The necessary theory listed as below, which is fundamental for our analysis result.

Lemma A.1 System $(1.5)$ is ergodic and positive recurrent if the following assumptions hold (^[42])

(i) $\mu_{ij} > 0$ for any $i\neq j$;

(ii) given $k\in\mathbb{S}$, $P(z, k) = (p_{ij}(z, k))$ is symmetric satisfying

with some constant $\varrho\in(0, 1]$;

(iii) there exists a non-negative function $\mathcal{W}(\cdot, k):\mathcal{D}^c\rightarrow \mathbb{R}$, for each $k\in\mathbb{S}$, such that $\mathcal{W}(\cdot, k)$ is twice continuously differentable and for some $\alpha > 0$,

where $\mathcal{D}^c$ is the complement of a bounded open subset $\mathcal{D}\in\mathbb{R}^n$ with a smooth boundary. Moreover, for any Borel measurable function $\varsigma(\cdot, \cdot):\mathbb{R}^n\times\mathbb{S}\rightarrow \mathbb{R}$, there is a unique stationary density $\pi(\cdot, \cdot)$ such that

which implies