Stochastic assessment of temperature–CO₂ causal relationship in climate from the Phanerozoic through modern times

1.   Intorduction

Although coral reefs only cover less than $0.2\%$ of ocean, they have close ties with fisheries, humans as well as the economy. They are the habitats and breeding ground for around $30\%$ of marine life, including nearly $25\%$ fish species ^[1]. In addition, they provide millions of employment opportunities for local fishers, coral reef fishing accounts for at least $6$ million metric tons of fish catches around the world each year ^[2]. Coral reefs could also create economic value through tourism, with net income estimated at close to $\$30$ billion annually ^[3]. However, during the last $50$ years, coral reefs around the world have suffered substantial deterioration due to the anthropogenic impacts such as overfishing and excess nutrients input from agriculture and domestic sewage ^[4], together with an increased severity and frequency of coral bleaching and mortality associated with climate change ^[1].

Considering today's rapidly varying environmental conditions and increasing human activity, many researchers investigated the mechanisms of coral reef ecosystems aiming to make efforts in coral reef conservation ^{[5,9,11,12,13,14,15,16,17,18]}. For example, Bellwood et al. ^[11] found out that the degradation of coral reefs usually manifested by phase shifts between two alternate stable states, specifically, from coral-dominated state to macroalgae-dominated state, with the decreasing of corals and the increasing of macroalgae abundance. The phase shift in coral reefs is mainly due to the fact that rapidly growing macroalgae dominate the competition for light and space by shading and reducing the available space for successful colonization of coral larvae ^[9,12,13,14]. Based on this mechanism, Mumby et al. ^[17] simplified the ecosystem into a three-state analytical model developed in terms of the cover (the fraction of available seabed in a given area) of macroalgae, corals and algal turfs respectively. Blackwood et al. ^[5] extended this model by explicitly including parrotfish grazing dynamics, demonstrated the effects of varying levels of fishing control measures on coral recovery time. Pal et al. ^[6] considered a mathematical model of interactions between coral, toxic seaweeds and herbivores, investigated how seeweed toxicity and overfishing negatively affect the ecological resilience of coral reefs through the trophic cascade. In this paper, we assumed a particular region of the seabed is covered entirely by macroalgae, coral, and algal turfs, the fraction of seabed they occupy is represented as $M$, $C$ and $T$, respectively. Obviously, $M+C+T = 1$ at any given time. We also suppose macroalgae could survive in the system irrespective of the abundance of corals. Algal turfs arise as a result of the natural mortalities of macroalgae and coral, as well as the grazing loss of macroalgae. Then the dynamics of macroalgae, corals and algal turf competing for available area in the particular region can be established as follows ^[5,6]:

where $0 < M(0), C(0), T(0) < 1$. Parrotfishes graze macroalgae and algal turfs without discrimination at a rate $g$ and $\frac{M}{M+T}$ is the proportion of grazing that affects macroalgae. More details of the parameters with biologically meaningful values are given in Table 1.

It is essential to take environmental stochasticity into account when investigating the dynamics of coral-reef system, since aquatic ecosystems are inevitably influenced by environmental fluctuations since many physical factors, such as nutrient availability, acidity, water temperature and so on, that embedded in aquatic ecosystems are usually unpredictable ^{[19,20,21,22,23,24,25,26,27,28,29,30]}. By the sensitivity analysis of parameters in system (1.1), the vegetative spread rate of macroalgae over algal turfs $r$ and the coral recruitment rate on turf algae $b$ are highly sensitive and could be easily affected by environmental stochasticity (more details can be seen in the Appendix). We then let $r\rightarrow r+\sigma_{1}\dot{B}_{1}(t)$, $b\rightarrow b+\sigma_{2}\dot{B}_{2}(t)$, the stochastic version of system (1.1) can be expressed as follows,

where $M(t)+C(t)+T(t) = 1$, for any $t > 0$. $\sigma_{i}^{2} > 0, \ i = 1, 2$ are the intensities of white noises, $B_{1}(t)$, $B_{2}(t)$ denote independent standard Brownian motions, which are defined in a complete probability space $(\Omega, \{\mathscr {F}_{t}\}_{t\geq0}, P)$, with a filtration $\{\mathscr {F}_{t}\}$ satisfying the usual normal conditions (right continuous and increasing while $\mathscr{F}_{0}$ contains all $P$-null set). The meaning of the other parameters is consistent with the deterministic system (1.1).

In this paper, we devote our main attention to investigate the dynamics of the coral reef benthic system with grazing intensity and immigrated macroalgae under the influence of environmental fluctuations. For the deterministic system (1.1), we investigate the existence and stability of equilibria, as well as the existence of bifurcations in Section 2. The results show that system (1.1) could possess more than one positive equilibrium and exhibits saddle-node bifurcation as well as bistability phenomenon, which implies that even a small environmental fluctuation may destroy the dynamics of the system. Consequently, we further consider the dynamics of the corresponding stochastic model (3.5) in Section 3. We investigate the existence of the unique positive solution, the existence of ergodic stationary distribution of system (3.5), as well as the extinction of corals. Moreover, the noise-induced transition from a macroalgae-coral coexistence state to one that without coral is also investigated by using the stochastic sensitivity functions (SSF) method. Finally, a brief conclusion summarizes our study in Section 4.

2.   Dynamics of the deterministic system

Denote

We first show $\tilde{\Omega}$ is a positively invariant of system (1.1).

Lemma 2.1. Solutions to system (1.1) with initial conditions in the set $\tilde{\Omega}$will remain there for all forward times.

Proof. To prove the positive invariance of $\tilde{\Omega}$, we need to examine the direction fields on the boundary of $\tilde{\Omega}$. When $M = 0$, $C > 0$, $T > 0$, $\frac{\mathrm{d}M}{\mathrm{d}t} > 0$. When $C = 0$, $M > 0$, $T > 0$, $\frac{\mathrm{d}C}{\mathrm{d}t} = 0$. When $T = 0$, $M > 0$, $C > 0$, $\frac{\mathrm{d}T}{\mathrm{d}t} = g+d_{1}M+d_{2}C > 0$. Therefore, all orbits starting from $\tilde{\Omega}$ cannot escape $\tilde{\Omega}$ from the boundaries. Besides, $M+C+T = 1$ holds for all $t\geq 0$, then $\tilde{\Omega}$ is a positively invariant of system (1.1). This concludes the proof. Without any loss of generality, system (1.1) can be simplified to the following form:

where $G_{1}(M, C) = aC-\frac{g}{1-C}+(r+\frac{\alpha}{M})(1-M-C)-d_{1}$, $G_{2}(M, C) = b(1-M-C)-d_{2}-aM$, and it is obvious that

is a positively invariant of system (2.2). We always assume that the initial value $(M(0), C(0))\in \Omega$ in the following analysis.

We then investigate the equilibria of system (2.2) and their stabilities. From system (2.2), it is obvious that $M\neq 0$. The equilibrium is the intersection of nullclines, where

$\bf{M-nullcline}$: $G_{1}(M, C) = 0$, where $\frac{\partial G_{1}(M, C)}{\partial M} < 0$ and $\frac{\partial^{2} G_{1}(M, C)}{\partial M^{2}} > 0$. Denote $l_{1}:\ G_{1}(M, C) = 0$;

$\bf{C-nullcline}$: $C = 0$ and $G_{2}(M, C) = 0$. Denote $l_{2}:\ G_{2}(M, C) = 0$.

The curve $l_{1}$ always has an intersection with $C = 0$, and its coordinate of the intersection $M_{1}$ satisfying $rM_{1}-\frac{\alpha}{M_{1}}-(r-\alpha-d_{1}-g) = 0$. $l_{2}$ is a line whose slop is $-\frac{a+b}{b}$ and the intersection values with M-axis and C-axis are $\frac{b-d_{2}}{a+b}$ and $1-\frac{d_{2}}{b}$, respectively, which are both proportional to the value $b$. With a small $b$, $l_{1}$ and $l_{2}$ do not intersect, see Figure 1(a). As the increasing of $b$, they will have at least one intersection once $b$ goes beyond a threshold $b^{*}$, when $b = b^{*}$, the curve $l_{1}$ is tangent to $l_{2}$ at $E_{3}(M_{3}^{*}, C_{3}^{*})$, and $b^{*}$ satisfying $f(b^{*}) = 0$, where

Moreover, as $b$ continues to increase, once $\frac{b-d_{2}}{a+b} > M_{1}$, $l_{1}$ and $l_{2}$ has only one intersection, see Figure 1(g). To sum up, system (2.2) possesses the following equilibria:

(i): a coral-free equilibrium $E_{1} = \left(M_{1}, 0\right)$, $M_{1} = \frac{r-\alpha-d_{1}-g+\sqrt{(r-\alpha-d_{1}-g)^{2}+4r\alpha}}{2r}$, $E_{1}$ always exists;

(ii): a unique interior equilibrium $E_{1}^{*} = (M_{1}^{*}, C_{1}^{*})$, when $b > \frac{M_{1}a+d_{2}}{1-M_{1}}: = \hat{b}$;

(iii): two different positive equilibria $E_{1}^{*} = (M_{1}^{*}, C_{1}^{*})$ and $E_{2}^{*} = (M_{2}^{*}, C_{2}^{*})$, provided that $b^{*} < b < \hat{b}$;

(iv): a double positive equilibrium $E_{3}^{*} = (M_{3}^{*}, C_{3}^{*})$ (i.e., $E_{1}^{*}$ and $E_{2}^{*}$ coincide), if $b = b^{*}$. Where $C_{i}^{*} = \frac{b-d_{2}-(a+b)M_{i}^{*}}{b}$, $M_{i}^{*}$ $(i = 1, 2, 3)$ satisfying the following polynomial equation $\sum_{j = 1}^{4}k_{i}M^{4-j} = 0$, $k_{1} = a(a+b)(a+b-r)$, $k_{2} = -a^{2}\alpha-a^{2}b+2a^{2}d_{2}-a\alpha b-ab^{2}+abd_{1}+2abd_{2}-2ad_{2}r+b^{2}d_{1}-bd_{2}r$, $k_{3} = -2a\alpha d_{2}-abd_{2}+ad_{2}^{2}-\alpha bd_{2}+b^{2}g+bd_{1}d_{2}-d_{2}^{2}r$, $k_{4} = -\alpha d_{2}^{2}$. More details can be seen in Figure 1.

We then investigate the stability of each equilibrium and verify the existence of saddle-node bifurcation. At $E_{1}$, the eigenvalues of the Jacobian matrix of system (2.2) are $b-d_{2}-M_{1}(a+b)$ and $-\sqrt{(r-\alpha-d_{1}-g)^{2}+4r\alpha}$, then

Lemma 2.2. System (2.2) is locally asymptotically stable at $E_{1}$, provided that $b < \hat{b}$. If $b > \hat{b}$, $E_{1}$ is unstable. System (2.2) undergoes a transcritical bifurcation at $E_{1}$, when $b$ crosses $\hat{b}$, provided that $\hat{b}\neq \frac{a[\alpha-M_{1}(a-g-r)]}{\sqrt{(r-\alpha-d_{1}-g)^{2}+4r\alpha}+M_{1}(a-g-r)-\alpha}$.

Proof. When $b = \hat{b}$, the Jacobian matrix $J|_{E_{1}}$ of system (2.2) at $E_{1}$ is as follows,

which exists a simple zero eigenvalue. Let $\textbf{U}_{1}$ and $\textbf{W}_{1}$ be the eigenvectors corresponding to the zero eigenvalue for $J|_{E_{1}}$ and $J|_{E_{1}}^{T}$, respectively. Then

Rewrite system (2.2) into a matrix form:

where $X = (M, C)^{T}\in R_{+}^{2}$ and

Then $\textbf{W}_{1}^{T}F_{b}\left(E_{1}, \hat{b}\right) = 0$, which implies that no saddle-node bifurcation occurs at $E_{1}$.

Also,

then

If $\hat{b}\neq \frac{a[\alpha-M_{1}(a-g-r)]}{\sqrt{(r-\alpha-d_{1}-g)^{2}+4r\alpha}+M_{1}(a-g-r)-\alpha}$, $\textbf{W}_{1}^{T}\left(D^{2}F\left(E_{1}, \hat{b}\right)(\textbf{U}_{1}, \textbf{U}_{1})\right)\neq 0$, then according to the Sotomayor theorem ^[31], system (2.2) undergoes a transcritical bifurcation at $E_{1}$ when $b$ crosses $\hat{b}$.

Lemma 2.3. When $b > \hat{b}$ holds, $E_{1}^{*}$ is locally asymptotically stable provided that $A_{2E_{1}^{*}} > 0$ and unstable if $A_{2E_{1}^{*}} < 0$. The expression of $A_{2E_{1}^{*}}$ can be seen in Eq (2.7).

Proof. The Jacobian $J_{1}^{*} = J_{1}|_{E_{1}^{*}}$ of the system (2.2) evaluated at the interior equilibrium $E_{1}^{*}$ is

The characteristic equation of the Jacobian $J_{1}^{*}$ is

where

Obviously, $A_{1E_{1}^{*}} > 0$, then Eq (2.5) could not admit a pair of pure imaginary roots, which implies that Hopf bifurcation does not exists. Besides, Eq (2.5) has at least one root with negative real part. If $A_{2E_{1}^{*}} > 0$ holds, both two roots have negative real part. Then $E^{*}$ is locally asymptotically stable. When $A_{2E_{1}^{*}} < 0$, Eq (2.5) must has a root with positive real part. The assertion is thus proved.

Similarly, we have the following conclusion.

Lemma 2.4. When $b^{*} < b < \hat{b}$, $E_{1}^{*}$ is locally asymptotically stable provided that $A_{2E_{1}^{*}} > 0$ and unstable if $A_{2E_{1}^{*}} < 0$. $E_{2}^{*}$ is a saddle point.

Lemma 2.5. When $b = b^{*}$, system (2.2) undergoes a saddle-node bifurcation at $E_{3}^*$ when $b$ crosses $b^{*}$, provided the following condition holds:

Proof. The characteristic equation of the Jacobian $J_{3}^{*}$ at $E_{3}^{*}$ of the system (2.2) is

where

Obviously, $A_{1E_{3}^{*}} > 0$, when $b = b^{*}$, $A_{2E_{3}^{*}} = C_{3}^{*}f(b^{*}) = 0$, the Jacobian $J_{3}^{*}$ of the system (2.2) has a simple zero eigenvalue. Let $\textbf{U}^{*}$ and $\textbf{W}^{*}$ are the eigenvectors corresponding to the zero eigenvalue for $J_{3}^{*}$ and $J_{3}^{*T}$, respectively. Then we have $\textbf{U}^{*} = \left(-\frac{b^{*}}{a+b^{*}}, \ 1\right)^{T}$, $\textbf{W}^{*} = \left(-\frac{(a+b^{*})C_{3}^{*}}{rM_{3}^{*}+\frac{\alpha (1-C_{3}^{*})}{M_{3}^{*}}}, \ 1\right)^{T}$. By simple calculation, we obtain $F_{b}(E_{3}^{*}, b^{*}) = \left(0, \ C_{3}^{*}(1-M_{3}^{*}-C_{3}^{*})\right)^{T}$ and

thus no transcritical bifurcation and pitchfock bifurcation occur at $E_{3}^{*}$ when $b$ crosses $b^{*}$. Also,

then when $\Theta\neq 0$,

It then follows from Sotomayor's theorem ^[31] that system (2.2) undergoes a saddle-node bifurcation at $E_{3}^{*}$ when $b$ crosses $b^{*}$.

To sum up, except the boundary equilibrium which is always exist, the system also can have one, two or a degenerated positive equilibrium points, see Figure 2. When system (2.2) has two positive equilibria, bistability phenomenon between a positive equilibrium and a boundary equilibrium occurs. At this situation, the eventual dynamics of the system is determined by the initial value. Moreover, once the two positive equilibria coincide each other and become a degenerated double positive equilibrium, it is a saddle-node at which the system undergoes a saddle-node bifurcation. At this situation, even a small fluctuation of the parameters near the bifurcation point may cause substantial changes on the dynamical behaviors of the system.

3.   Dynamics of the stochastic system

In Section (2), we studied the occurrence of the bistability and bifurcation phenomena of the deterministic system (1.1), which implies that the dynamics of this coral benthic system is highly sensitive to the parameters, and therefore could be significantly affected by the fluctuation of environment. Next, we are going to investigate how environmental fluctuations affect the population dynamics, including the global existence and uniqueness of positive solutions, the existence of the unique ergodic stationary distribution as well as the extinction of coral population. For starters, we establish the following lemma to illustrate that the stochastic system system (1.2) admits a unique global positive solution.

3.1.   The existence of the unique positive solution

Lemma 3.1. For any initial value $(M(0), C(0), T(0))\in R_{+}^{3}$, thereexists a unique solution $(M(t), C(t), T(t))$ to model (1.2) on $t\geq 0$, and the solution will remain in$R_{+}^{3}$ with probability one, namely, $(M(t), C(t), T(t))\in R_{+}^{3}$ for all $t\geq 0$ almost surlely (a.s.).

Proof. Since the coefficients of system (1.2) are locally Lipschitz in $R_{+}^{3}$, then there exists a unique local solution $(M(t), C(t), T(t))\in R_{+}^{3}$ of system (1.2) on the interval $(0, \tau_{e})$, where $\tau_{e}$ denotes the explosion time. To show the solution is global, we only need to prove that $\tau_{e} = \infty$. By the positivity and boundedness of initial value, there is a number $n_{0} > 0$ large enough such that $M(0), C(0), T(0)\geq\frac{1}{n_{0}}$. For each integer $n\geq n_{0}$, define the stopping times as

Obviously, $\tau_{n}$ is increasing as $n\rightarrow +\infty$. To show $(M(t), C(t), T(t))\in R_{+}^{3}$, we only need to show that $\tau_{\infty} = \infty$ a.s., then $\tau_{e} = \infty$, a.s. We prove it by contradiction. If $\lim_{n\rightarrow +\infty}\tau_{n} < \infty$, then there exists a pair of constants $T_{1} > 0$ and $\varepsilon_{1}\in (0, 1)$ such that $P\{\tau_{n}\leq T_{1}\} > \varepsilon_{1}$. Therefore, there exists an integer $n_{1}\geq n_{0}$ such that for all $n\geq n_{1}$,

Let

where $m_{1} = \frac{\min\{d_{1}, d_{2}\}}{g}$, then according to the It$\hat{\mathrm{o}}$'s formula ^[32], we have

Then,

Integrating the inequality from $0$ to $\tau_{n}\wedge T_{1}$ and taking the expectations of the above inequality leads to

Let $\Omega_{n} = \{\omega\in \Omega_{n}: \tau_{n} = \tau_{n}(\omega)\leq T_{1}\}$ for $n\geq n_{1}$ and in view of Eq (3.1), we know that $P(\Omega_{n}) > \varepsilon_{1}$. Note that for every $\omega\in \Omega_{n}$, at least one of $M(\tau_{n}, \omega)$, $C(\tau_{n}, \omega)$, $T(\tau_{n}, \omega)$ equals $\frac{1}{n}$. Hence, $V_{1}(M(\tau_{n}, \omega), C(\tau_{n}, \omega), T(\tau_{n}, \omega))$ is no less than

It follows from Eq (3.2) that

where $I_{\Omega_{n}}$ denotes the indicator function of $\Omega_{n}$. Letting $n\rightarrow \infty$ in Eq (3.2), we have

which leads to the contradiction and thus we have $\lim_{n\rightarrow +\infty}\tau_{n} = \infty$, a.s. This completes the proof.

Lemma 3.2. The population of macoralgae is weakly persistent, i.e, $\limsup_{t\rightarrow \infty }M(t) > 0$, a.s.

Proof. We prove it by contraction. Assume $\limsup_{t\rightarrow \infty }M(t)\leq0$, combining with the fact that $\liminf_{t\rightarrow \infty }M(t)\geq0$, we have $\lim_{t\rightarrow \infty }M(t) = 0$. Then there exists a sufficient small $\varepsilon_{2}$ and a positive $T_{2} > 0$, such that when $t > T_{2}$, $M(t)\leq \varepsilon_{2}$. By using the It$\hat{\mathrm{o}}$'s formula to system (1.2) we have,

Let $\varepsilon_{2}$ small enough, such that

then

Integrating Eq (3.4) from $T_{2}$ to $t$ and dividing by $t$ on both sides yields,

Moreover, from Lemma 3.1 we know that, $\limsup_{t\rightarrow \infty}\frac{\ln M(t)}{t}\leq 0$ and $\limsup_{t\rightarrow \infty}\frac{\ln T(t)}{t}\leq 0$. Let $t\rightarrow \infty$ and making use to the strong law of large numbers for martingales ^[32], we have,

This leads to the contradiction and thus we have $\limsup_{t\rightarrow \infty}M(t) > 0$, a.s.

According to Lemma 3.1, $\tilde{\Omega}$ as shown in Eq (2.1) is a positive invariant set of the stochastic system (1.2). Without any loss of generality, we only need to consider the following system in the rest of this paper,

It is obvious that $\Omega$ as shown in Eq (2.3) is a positive invariant set of system (3.5).

3.2.   Existence of stationary distribution

Before proving the main result in this subsection, we first present the following lemma from ^[33]. Suppose that $X(t)$ is a homogeneous Markov process in $n$-dimension Euclidean space $R^{n}$, satisfying the following stochastic differential equation:

where $\sigma_{r}(X) = (\sigma_{r}^{1}(X), \sigma_{r}^{2}(X)\ldots, \sigma_{r}^{n}(X))^{T}$, $A(X) = (a_{ij}(X))_{n\times n}$ is the diffusion matrix of $X(t)$ with $a_{ij}(X) = \sum _{r = 1}^{k}\sigma_{r}^{i}(X)\sigma_{r}^{j}(X)$.

Lemma 3.3. If there exists a bounded open domain $U\subset R^{n}$ with regularboundary, satisfying the following properties ^[33]:

$(\textbf{H1)}$ The diffusion matrix $A(x)$ is strictly positive definite for all $x\in U$;

$(\textbf{H2})$ There exists a non-negative $\mathcal{C}^{2}$-function $V(X)$ and a positive constant $\zeta$ such that $\mathscr{L}V(X)\leq-\zeta$ on $X\in R^{n}\setminus U$.

Then the Markov process $X(t)$ of the stochastic model (3.6) admits a unique stationary distribution $\nu(\cdot)$, and forany integrable function $f(\cdot)$ with regard to the measure$\nu$, the following equation holds,

In the following, we will apply Lemma 3.3 to prove the existence of a unique ergodic stationary distribution for system (3.5). We first make some notations and assumptions as follows:

Denote

Assumption 1. $\lambda_{1}, \lambda_{2} > 0$.

Assumption 2. $m_{2}\lambda_{1}+\lambda_{2} > (m_{3}+m_{4})(r+\alpha)+\alpha$,

where $m_{2} = \frac{\min\{d_{1}, d_{2}\}}{lg}$, the constant $l > 1$, $m_{3} = \frac{m_{2}(r+\alpha) +(r-g-d_{1}+\frac{\sigma_{1}^{2}}{2})}{r+\alpha+d_{1}+g}$, $m_{4} = \frac{a+b}{r+\alpha+d_{1}+g}$.

Theorem 3.1. If Assumptions 1, 2 hold, then for any initial value $(M(0), C(0))\in \Omega$, the stochastic system (3.5) admits a unique stationary distribution and it has ergodicproperty.

Proof. Define a $\mathcal{C}^{2}$-function $V_{2}: R_{+}^{2}\rightarrow R_{+}$ by

where $V_{21} = -\ln C+m_{4}M$, $V_{2}^{*}$ is the minimum value point of $V_{2}$.

By using the It$\hat{\mathrm{o}}$'s formula we have,

Then

Denote

where $k > 1$ is sufficient large. In the following, we will prove $\mathscr {L} V_{2} < -\zeta$ for any $(M, C)\in \Omega\backslash D_{k}$. To achieve that, we divide $D_{k}^{\mathcal{C}} = \Omega\backslash D_{k}$ into the following three parts,

Notice that $D_{k}^{\mathcal{C}} = \Omega\backslash D_{k} = D_{k}^{1}\cup D_{k}^{2}\cup D_{k}^{3}$, then we only need to prove $\mathscr {L} V_{2} < -\zeta$ respectively on the above three domains.

(1) When $(M, C)\in D_{k}^{1}$, i.e., $M+C\geq 1-\frac{1}{k}$, then $1-M-C\leq\frac{1}{k}$,

where $D_{2} = \min\{d_{1}, d_{2}\} +(m_3 + m_{4}+1)(r+\alpha)+m_{2}g+m_{2}d_{1} +2b+(m_3 + m_{4})a+\frac{(m_{2}+1)\sigma_{1}^{2}+\sigma_{2}^{2}}{2}$. It is not hard to find a sufficiently large $k$ such that $\mathscr {L}V_{2} < -\zeta$ for all $(M, C)\in D_{k}^{1}$.

(2) When $(M, C)\in D_{k}^{2}$, i.e., $M\leq \frac{1}{k^{2}}$ and $M+C\leq 1-\frac{1}{k}$, then $1-M-C\geq\frac{1}{k}$,

then there exists a sufficiently large $k$, such that $\mathscr {L}V_{2} < -\zeta$ for all $(M, C)\in D_{k}^{2}$.

(3) When $(M, C)\in D_{k}^{3}$, i.e., $C\leq \frac{1}{k}$,

Since $-m_{2}\lambda_{1}-\lambda_{2}+(m_3 + m_{4})(r+\alpha)+\alpha < 0$. Then there exists a sufficiently large $k$, such that $\mathscr {L}V_{2} < -\zeta$ holds for all $(M, C)\in D_{k}^{3}$. To sum up, we can conclude that

for all $(M, C)\in D_{k}^{\mathcal{C}}$ as long as $k$ is sufficient large. This confirms the condition (H2) in Lemma 3.3.

Moreover, system (3.5) can be rewritten as the following form:

with the diffusion matrix $A(X) = \mathrm{diag}\big(\sigma_{1}^{2}M^{2}(1-M-C)^{2}, \sigma_{2}^{2}C^{2}(1-M-C)^{2}\big)$. Besides, there exists a

such that

for all $(M, C)\in D_{k}$, $\xi = (\xi_{1}, \xi_{2})\in R^{2}$, which shows that condition (H1) of Lemma 3.3 is also satisfied. Then, based on Lemma 3.3, system (3.5) has a unique stationary distribution and it has ergodic property.

We then make some numerical simulations to verify our results. Let $b = 1$, $r = 0.8$, $l = 15$, $\sigma_{1} = 0.02$, $\sigma_{2} = 0.05$, by calculation, $\lambda_{1} > 0$, $\lambda_{2} > 0$ and $m_{2}\lambda_{1}+\lambda_{2}-(m_3 + m_{4})(r+\alpha)-\alpha = 0.0054 > 0$, then the conditions of Theorem 3.1 hold, system (3.5) has a unique stationary distribution and it has ergodic property, see Figure 3.

3.3.   Extinction of the coral population

According to Theorem 3.1, both macoralgae and coral populations will be persistent in mean, provided that the environmental fluctuations are small enough such that Assumptions 1 and 2 hold. Moreover, from Lemma (3.2), we know that the macroalgae population is always weakly persistent, no matter what the noise intensities are. In this section, we are going to investigate under what conditions the population of coral will go to extinction. Our main result is as follows.

Theorem 3.2. Let $(M(t), C(t))$ be the solution of system (3.5) with initial value $(M(0), C(0))\in \Omega$, then the coral population will go to extinction, i.e., $\lim_{t\rightarrow \infty}C(t) = 0, \ a.s.$, provided that one of the following conditions holds:

$(1)$ $\sigma_{2}^{2} > \frac{b^{2}}{2d_{2}}.$

$(2)$ $\sigma_{2}^{2} < b$ and $\lambda_{2} < 0$.

Proof. (1) We first prove when condition (1) is satisfied, coral population will go to extinction. By using the It$\hat{\mathrm{o}}$'s formula to the second equation of system (3.5) results in

When $\sigma_{2}^{2} > \frac{b^{2}}{2d_{2}}$, integrating Eq (3.7) from $0$ to $t$ and dividing by $t$ on both sides leads to

i.e., $\lim_{t\rightarrow \infty}C(t) = 0, \ a.s$. This completes the first part of Theorem 3.2.

(2) Now we are going to prove when condition (2) is satisfied, $\lim_{t\rightarrow \infty}C(t) = 0, \ a.s.$ also holds. Again, by using the It$\hat{\mathrm{o}}$'s formula we have,

Integrating both sides from $0$ to $t$ and dividing by $t$ leads to,

i.e., $\lim_{t\rightarrow \infty}C(t) = 0, \ a.s$. This completes the proof.

From the above analysis, $\lambda_{1}$ and $\lambda_{2}$ could be approximated as the indicators of net growth rate for macroalgae and coral population. When $m_{2}\lambda_{1}+\lambda_{2} > (m_3 + m_{4})(r+\alpha)+\alpha$, system (3.5) will admit a unique stationary distribution and it has ergodic property. When $\lambda_{2} < 0$, i.e., when the average recruitment rate of corals could not overcome the negative influence caused by its natural mortality and environmental noise, the coral population will go to extinction, see Figure 4.

Moreover, comparing Figure 4 to Figure 1, when $b = 0.365 > \hat{b} = 0.364$, the deterministic system will be persistent in mean, while the existence of environmental fluctuations will cause the extinction of corals. This implies that, the survival conditions for coral would be more rigorous with the existence of environmental fluctuations, and corals need to possess a larger recruitment rate to survive in such stochastic environment.

3.4.   Analysis of noise-induced transitions

According to the analysis in Section 2, the deterministic system (2.2) could have two different stable equilibria, namely, the coral extinction equilibrium $E_{1}$ and the macroalgae-coral coexistence equilibrium $E_{1}^{*}$. At this scenario, the eventual dynamics of system (2.2) is completely determined by its initial value. We then use the knowledge of SSF technique ^[34,35] to theoretically investigate the influence of environmental fluctuations in the bistable zone.

Choosing the parameters as in Figure 1(e), (f), and $\sigma_{1} = \sigma_{2} = \sigma$, denote the deterministic coexistence equilibrium $E_{1}^{*} = (M_{1}^{*}, C_{1}^{*})$, then applying the methods in ^[35], the confidence ellipse equation of $E_{1}^{*}$ can be represented by

When $\sigma = 0.005$, the random states of system (3.5) fluctuate around the coexistence equilibrium $E_{1}^{*}$ of the corresponding deterministic system, and they will locate in the interior of the confidence ellipse with probability $P = 0.95$, see Figure 5(a) and Figure 6(a), (b). As the noise intensity increases, the confidence ellipse starts to grow and once the noise intensity goes beyond a threshold value $\sigma^{*}\approx 0.0385$, which can be seen as the intensity corresponding to the confidence ellipse intersects the separatrix, it ultimately arrives the attraction basin of the coral extinction equilibrium $E_{1}$, see Figure 5(b) and Figure 6(c), (d).

4.   Conclusions

In this paper, we investigated a coral reef benthic system in which macroalgae and corals compete to occupy algae turfs with grazing intensity and immigrated macroalgae in both deterministic and stochastic environments. In deterministic environment, besides the coral-free equilibrium, system (2.2) can also have one, two or a degenerated macroalgae-coral coexistence equilibria, whose existence and stability have been clearly investigated. When system (2.2) exists two different coexistence equilibria, bistability between a macroalgae-coral coexistence equilibrium and the coral-free equilibrium occurs. Once the two macroalgae-coral coexistence coincide each other and become a degenerated double coexistence equilibrium, it becomes a saddle-node at which the model undergoes a saddle-node bifurcation. The complex dynamics of the deterministic system illustrate that the dynamics of the model could be highly sensitive to the parameters, and therefore can be significantly affected by the fluctuation of environment.

We then formulate the stochastic system (1.2), by adding environmental stochasticity into two highly sensitive parameters which are selected by the sensitivity analysis of system (1.1). By investigating the dynamics of the stochastic system (3.5), we obtained the conditions for the existence of the unique ergodic stationary distribution as well as coral extinction. Our results reveal that when $\lambda_{2} < 0$, i.e., the average recruitment rate of corals on turf algae could not overcome its own natural loss and the negative impact of environmental noises on its populations, it will go to extinction. When the average net growth indicators of coral and macroalgae populations, $\lambda_{1}$ and $\lambda_{2}$, both go beyond certain levels such that Assumptions 1, 2 hold, then system (3.5) would be persistent in mean and admits a unique ergodic stationary distribution. Moreover, the existence of environmental fluctuation makes the survival conditions for coral more harsher. Coral population that could survive in the deterministic system might be driven to extinction due to the existence of environmental fluctuation, and they need to possess larger recruitment rate to survive in such stochastic environment. For the bistable scenario between a coexistence equilibrium and the coral-free equilibrium in the absence of environmental fluctuations, the existence of environmental noises could cause state transition from macroalgae-coral coexistence to coral-extinction, and this transition is irreversible. Moreover, we can observe from Figure 1 and Figure 3 that, when coral and macroalgae could coexist, the persistent levels of coral are usually much higher than the corresponding levels of macroalgae. This could explain the phenomenon that coral reefs are widely exhibit at least two common stable states in the species-poor Caribbean ^[11], one such state corresponds to high levels of coral cover, and the other state corresponds to coral depletion ^[5].

As a summary, ignoring the impact of environmental stochasticity may lead to bias in the prediction and management of coral reef ecosystem status, and our analysis in this paper may provide theoretical help for reef managers to face the changing environment. There are other factors that we have not explicitly included that could clearly affect the coral reef dynamics. For example, coral reefs are occasionally suffered from hurricane ^[11], which could also significantly affect the coral reef ecosystem but can not be described by Gaussian white noise that introduced in this paper. Therefore, it is meaningful to further consider the influence of other type noises, for example, L$\acute{e}$vy noise ^[36,37] in coral reef ecosystem. Furthermore, macroalgae are highly seasonal in their occurrence, growth, and reproduction ^[38,39], it will be more reasonable to further include the seasonal influence on bioactivity of macroalgae in coral reefs. We leave these for future investigation.

Acknowledgments

This work was supported by the National Natural Science Foundation of China [grant numbers 12071293, 11671260].

Conflicts of interest

The authors declared that they have no conflict of interest.

Appendix: Sensitivity analysis of parameters in system (1.1)

In order to measure the relative changes of state variables when model parameter changes, we then perform the sensitivity analyses on system (1.1) at the interior equilibrium. The normalized forward sensitivity index of variable $\mu$ about parameter $\rho$ is as follows:

The higher the sensitivity index value, the more sensitive the variable is to the parameter.

Take $b = 0.365$, other parameters are same as in Figure 1. At this situation, system (1.1) exists a unique stable interior equilibrium $E^{*}\approx(0.0117, 0.3276, 0.6607)$. Applying the implicit function derivative rule, the values of the sensitivity indices for each component at $E^{*}$ are presented in Table A1 and Figure A1. Our result shows that the sensitivity index sign of the parameters for $M^{*}$ and $T^{*}$ remain consistent, while the sign for $C^{*}$ is opposite. Specifically, the overgrowth rate of macroglgae on coral $a$, the vegetative spread rate of macroalgae over algal turfs $r$, the colonization rate of newly immigrated macroalgae on algal turf $\alpha$ and natural mortality rate of coral reefs $d_{2}$ have positive impact on $M^{*}$ and $T^{*}$, and have negative impact on $C^{*}$. Other parameters have opposite effects on each component of $E^{*}$. Moreover, the most sensitive parameter for $M^{*}$, is coral recruitment rate on turf algae $b$ and coral reefs natural mortality rate $d_{2}$, other important parameters are the maximal macroalgae-grazing rate of parrotfish on macroalgae $g$ and the vegetative spread of macroglgae over algal turfs rate $r$. For $T^{*}$ and $C^{*}$, the most sensitive parameter is $b$, other important parameter is $d_{2}$.