Structure analysis of the attracting sets for plankton models driven by bounded noises

1.   Introduction

Plankton are located in the first trophic layer of the aquatic food chain and are the base of the aquatic ecosystem ^[1]. They not only generate organic compounds by absorbing carbon dioxide dissolved in the surrounding environment but also perform photosynthesis, which has an important impact on large-scale global processes such as the global carbon cycle, climate change and ocean-atmosphere dynamics ^[2].

Toxin producing phytoplankton (TPP) are a kind of harmful plankton with the ability to release toxic chemicals into the environment. The toxic chemicals may inhibit predation pressure from phytoplankton and other predator populations in planktonic systems ^[3] and then contribute to the formation of harmful algal blooms (HABs) ^[4]. For example, some experimental observations in ^[5] have indicated that the toxic dinoflagellate Alexandrium fundyense can negatively affect the growth rate of the copepod Acartia hudsonica, and the toxic effects may have profound implications on the ability of grazers to control the HABs. Over the past few decades, research on the complex dynamics of planktonic systems has attracted great interests of researchers; see previous studies ^{[6,7,8,9,10,11,12]} and the references therein.

For example, for a nutrient-phytoplankton system with TPP, Chakraborty et al. ^[1] established the following nonlinear mathematical model:

where $N(t)$ and $P(t)$ are the concentrations of nutrient and TPP at time $t$, respectively. Parameter $a$ is the external nutrient inflow rate, $b$ is the nutrient uptake rate of phytoplankton, $c\; (c < b)$ is the conversion rate of nutrients into phytoplankton, $h$ is the loss rate of nutrients, $d$ is the mortality rate of phytoplankton, $k\; (k < d)$ is the nutrient recycle rate due to the death of phytoplankton, $\mu$ is the half saturation constant, and $\theta$ represents the rate of release of toxic chemicals by the TPP population. All the parameters are assumed to be positive. The authors showed that, for a certain range of $\theta$, model (1.1) exhibits periodic solutions. They also observed that toxin produced by the TPP may act as a biological control in the termination of the planktonic bloom, which is in good agreement with some earlier findings.

Because the real environment is full of stochasticity, and every ecosystem is inevitably affected by environmental noise, it seems more appropriate to develop some stochastic ecological models by considering the influence of environmental noise ^{[13,14,15,16,17]}. For example, Ji et al. ^[13] established a stochastic Lotka-Volterra predator-prey model with white noise, and obtained some criteria for persistence and extinction of the species. Zhang et al. ^[14] proposed and studied a stochastic non-autonomous prey-predator model with impulsive effects. They showed that the stochastic noise and impulsive perturbations have crucial effects on the persistence and extinction of each species.

In fact, plankton systems are more susceptible to environmental fluctuations such as light, water temperature and water pH ^[18,19,20]. Therefore, based on deterministic model (1.1), Yu et al. ^[18] constructed the following plankton model with white noise:

where $W_{i}(t)$ are standard Wiener processes with intensities $\alpha_{i}$, $i = 1, 2$, and $W_{i}(t)$ are defined on a complete probability space $\left(\Omega, \mathcal{F}, \mathbb{P}\right)$ with filtration $\left\{\mathcal{F}_{t}\right\}_{t\ge 0}$. For stochastic model (1.2), the authors gave sufficient criteria for the existence of ergodic stationary distribution and investigated the extinction and persistence of the phytoplankton species. They also showed that the TPP and environmental fluctuations may have great influence on planktonic blooms.

Different from the standard Wiener process, the Ornstein-Uhlenbeck (O-U) process ^[21] can be used to model the bounded environmental fluctuation in a real ecosystem. The ecological model driven by O-U process is closer to reality, as stated by Caraballo et al. ^[22]: "The most common stochastic process that is considered when modeling disturbances in real life is the well-known standard Wiener process. Nevertheless, this stochastic process has the property of having continuous but not bounded variation paths, which does not suit to the idea of modeling real situations since, in most of cases, the real life is subjected to fluctuations which are known to be bounded." So, ecological models driven by O-U process have been proposed and analyzed by some scholars; see ^{[22,23,24,25,26,27,28]} and the references therein. For example, Caraballo et al. ^[22] used the O-U process to model the bounded noise perturbations in a logistic system and competitive Lotka-Volterra system. They present an example testing the theoretical result with real data and verified that this method is a realistic one. A general eco-epidemiological system, in which the birth rate of prey population is driven by O-U process, was considered in ^[24]. The authors proved the existence of a global random attractor and the persistence of susceptible prey, and provided some conditions for the simultaneous extinction of infective preys and predators. In ^[25], López-de-la-Cruz derived a random chemostat model driven by O-U process and investigated the existence and internal structure of the attracting set (or attractor) for the random model. In reality, the internal structure of the attracting set can reflect the survival of species in ecosystem ^[25,29]. For other systems with O-U process we refer also to ^{[30,31,32,33]}.

In view of the latest research and the advantages of O-U process, for the nutrient-phytoplankton system with TPP, we consider the environmental fluctuations to be bounded and model the bounded noise by using suitable O-U process in this paper. Based on deterministic model (1.1) and stochastic model (1.2), respectively, we first construct two random plankton models and then investigate the existence and internal structures of the attracting sets for these models.

This paper is organized as follows: In Section 2, we analyze a random plankton model corresponding to deterministic model (1.1) in which the external nutrient inflow rate $a$ is driven by an O-U process. In Section 3, we use the O-U process to transform stochastic plankton model (1.2) into a random one and investigate the attracting set for the random model. A simple discussion is given in Section 4. For completeness, some mathematical backgrounds of the O-U process and random dynamical system is given in the Appendix.

2.   Random plankton model related to deterministic system (1.1)

In this section, we will consider a suitable O-U process to perturb the external nutrient inflow rate $a$ in a deterministic plankton system in the same way as in ^[25,34]. Particularly, we are interested in replacing $a$ by the random term $a+\sigma z^{\ast}(\theta_{t}\omega)$ in deterministic model (1.1), where $z^{\ast}(\theta_{t}\omega)$ denotes the O-U process which will be introduced in the Appendix, and $\sigma > 0$ represents the intensity of perturbation. In such a way, the resulting random model is given by the following system of random differential equations:

We would also like to note that, thanks to the property $\lim_{\lambda\to\infty}z^{\ast}(\theta_{t}\omega) = 0$ shown in Proposition A.1, for every fixed $\omega\in\Omega$, it is possible to take $\lambda$ large enough such that $a+\sigma z^{\ast}(\theta_{t}\omega)\in(a_1, a_2)$ for every $t\in\mathbb{R}_+$, where $a_1$ and $a_2$ are positive constants.

We will introduce the main results for random model (2.1) in the following three subsections, including existence and uniqueness of the global positive solution and existence and internal structure of the attracting set.

2.1.   Existence and uniqueness of global positive solution

Theorem 2.1. For any initial value $S(0): = (N(0), P(0))\in\mathbb{R}_+^2$ and any $\omega\in\Omega$, model (2.1) possesses a unique global positive solution

with $S(0; 0, \omega, S(0)) = S(0)$.

Proof. Let $S\left(t; 0, \omega, S(0)\right) = \left(N(t; 0, \omega, N(0)), P(t; 0, \omega, P(0))\right)$, then model (2.1) can be rewritten as

where

and $F:\mathbb{R}_{+}^{2}\times \mathbb{R}_{+}\to \mathbb{R}^{2}$ is given by

where $\eta = (\eta_1, \eta_2)\in\mathbb{R}_{+}^{2}$.

One can find that $F(\eta, \theta_{t}\omega)$ is locally Lipschitz with respect to $(\eta_1, \eta_2)$, then model (2.1) possesses a unique local solution. To prove the local solution is a global one, we define the new variable $U = N+P$. It is easy to see that $U$ satisfies the following equation:

Notice that $a+\sigma z^{\ast}(\theta_{t}\omega)\in(a_1, a_2)$, and from Eq (2.2), we know

where $m_1 = \min\left\{h, d-k\right\}$. It is straightforward to check that $U$ does not blow up at any finite time, and the same happens to $N$ and $P$. Therefore, the unique local solution can be extended to a global one.

Moreover, from Eq (2.1), we know that

for all $P\ge0$, and

for all $N\ge0$. Thus, the unique global solution $S\left(t; 0, \omega, S(0)\right)$ of random model (2.1) remains in the positive quadrant $\mathbb{R}_+^2$ for every initial value $S(0)\in\mathbb{R}_+^2$.

Remark 2.1. Define a mapping $\varphi_{S}:\mathbb{R}_{+}\times\Omega\times\mathbb{R}_{+}^{2}\to \mathbb{R}_{+}^{2}$ given by

Since the function $F$ is continuous in $(S, t)$ and is measurable in $\omega$, we obtain the $(\mathcal{B}(\mathbb{R}_{+})\times\mathcal{F}\times\mathcal{B}(\mathbb{R}_+^2), \mathcal{B}(\mathbb{R}_+^2))$-measurability of this mapping, which defines a random dynamical system generated by the solution mapping of model (2.1).

2.2.   Existence of attracting set

Theorem 2.2. There exists a deterministic compact attracting set

for the solution of model (2.1), where $U_1 = \frac{a_1}{M_1}$, $U_2 = \frac{a_2}{m_1}$, $\underline{N}_{1} = \frac{a_1}{bU_2+h}$, and

Proof. According to inequality (2.3), we can obtain

On the other hand, it follows from equation (2.2) that

By setting $M_1 = \max\left\{(b-c)U_2+h, \; d-k+\frac{\theta}{2\mu}\right\}$, we get

and then

According to inequality (2.4), we know that, for every initial value $S(0)\in\mathbb{R}_+^2$ and any given $\varepsilon > 0$, there exists some time $T(\omega, S(0), \varepsilon) > 0$ such that $U(t)\le U_2+\varepsilon$ for all $t\ge T(\omega, S(0), \varepsilon)$. Therefore, we know $N(t)+P(t)\le U_2$ holds for every time $t$ large enough. It then follows from $a+\sigma z^{\ast}(\theta_{t}\omega) > a_1$ and $P\le U_2$ that

and

Thus, from inequalities (2.4), (2.5) and (2.6), we can obtain that

is a deterministic attracting set for the solution of model (2.1).

Remark 2.2. The existence of attracting set $B_0$ indicates that the inequalities

hold for every time $t$ large enough, where $S(t) = (N(t), P(t))$ is the solution of model (2.1). In what follows, we always believe that these inequalities are true, because the purpose of this paper is to explore the long-time behavior of the plankton species.

Taking parameters $a_1 = 1.3$, $a_2 = 2.7$, $b = 0.8$, $c = 0.7$, $k = 0.1$, $h = 0.4$, $d = 0.5$, $\theta = 0.1$ and $\mu = 2$ in model (2.1), we can calculate $U_1 = 1.21$, $U_2 = 6.75$ and $\underline{N}_{1} = 0.22$, and the simulation of the attracting set is shown in Figure 1.

2.3.   Internal structure of the attracting set

Theorem 2.3. For model (2.1), if the condition

holds, then the attracting set $B_0$ is reduced to a line segment on the coordinate axis. More precisely, it is

where $\underline{N}_{2} = \frac{a_1}{h}$ and $\overline{N} = \frac{a_2}{h}$.

Proof. According to Remark (2.2), we know $N < U_2$. It then follows from

that

If $cU_2-d < 0$, we know that

and then, for every time $t$ large enough, the first equation of model (2.1) can be written as

It follows from $a+\sigma z^{\ast}(\theta_{t}\omega)\in(a_1, a_2)$ that

and then

Therefore, the attracting set of model (2.1) will become

which is a line segment on the coordinate axis.

Taking parameters $a = 2$, $a_1 = 1.3$, $a_2 = 2.7$, $b = 0.8$, $c = 0.1$, $k = 0.1$, $h = 0.6$, $d = 0.6$, $\theta = 0.1$, $\lambda = 20$, $\mu = 2$ and $\sigma = 0.5$ in model (2.1), we can calculate $\underline{N}_{2} = 2.16$ and $\overline{N} = 4.5$. The simulation of the attracting set and three trajectories with different initial values is shown in Figure 2. One can see that the trajectory of model (2.1) eventually enters the line segment $B_0$ on the coordinate axis, which indicates that the phytoplankton species will go extinct, and only the nutrient can be persistent.

Theorem 2.4. For model (2.1), if the condition

holds, then the attracting set $B_0$ is reduced to a plane region in the first quadrant. More precisely, it is

where $\underline{N}_{3} = \frac{a_1+kP_1}{bU_2+h}$, $P_1 = \frac{cU_1-(d+\frac{\theta}{2\mu})}{c}$.

Proof. From Remark 2.2, we know $U_1-P\le N$. It follows from

that

If $cU_1-(d+\frac{\theta}{2\mu}) > 0$, we can conclude

For every time $t$ large enough, it follows from $a+\sigma z^{\ast}(\theta_{t}\omega) > a_1$ and $P_1\leq P\le U_2$ that

and then

Therefore, the attracting set of model (2.1) will become

At this time, the attracting set lies completely in the first quadrant plane.

Taking parameters $a = 2$, $a_1 = 1.3$, $a_2 = 2.7$, $b = 0.8$, $c = 0.7$, $k = 0.1$, $h = 0.4$, $d = 0.5$, $\theta = 0.1$, $\lambda = 20$, $\mu = 2$ and $\sigma = 0.5$ in model (2.1), we can calculate $\underline{N}_{3} = 0.23$, $P_1 = 0.45$, $U_1 = 1.2$ and $U_2 = 6.75$. The simulation of the attracting set and three trajectories with different initial values is shown in Figure 3. One can see that the trajectory of model (2.1) eventually enters the plane region $B_0$ in the first quadrant, which indicates that the phytoplankton species and nutrient can be simultaneously persistent.

3.   Random plankton model related to stochastic system (1.2)

In this section, we assume that the nutrient and phytoplankton species in plankton system are affected by the same white noise, and then model (1.2) is reduced to the following stochastic model in Itô's sense:

Due to the properties of Stratonovich integrals following the classical rules in calculus, with the help of the well-known conversion between Itô's and Stratonovich's senses, we further rewrite model (3.1) as the following stochastic model in Stratonovich's sense:

In what follows, we use the O-U process to transform stochastic model (3.2) into a random one. To this end, we first define two new variables $x(t)$ and $y(t)$ as follows:

For the sake of simplicity, we will write $z^{\ast}$ instead of $z^{\ast}(\theta_{t}\omega)$, $x$ instead of $x(t)$, and $y$ instead of $y(t)$. From Eq (3.2) and the Langevin equation shown in the Appendix, we know that variables $x$ and $y$ satisfy the following equations:

According to the property $\lim_{\lambda\to 0}\lambda z^{\ast}(\theta_{t}\omega) = 0$ shown in Proposition A.1, for every fixed $\omega\in\Omega$, it is possible to choose a suitable $\lambda$ such that $\frac{\alpha^{2}}{2}-\alpha\lambda z^{\ast}\in(l_1, l_2)$ for every $t\in\mathbb{R}_+$, where $l_1 < l_2 < \infty$, and so that both $h+l_1$ and $d-k+l_1$ are positive.

We will introduce the main results for random model (3.3) in the following three subsections, including existence and uniqueness of the global positive solution and existence and internal structure of the attracting set.

3.1.   Existence and uniqueness of global positive solution

Theorem 3.1. For any initial value $X(0): = (x(0), y(0))\in\mathbb{R}_+^2$ and any $\omega\in\Omega$, model (3.3) possesses a unique global positive solution

with $X(0; 0, \omega, X(0)) = X(0)$.

Proof. Let $X(t; 0, \omega, X(0)): = \left(x(t; 0, \omega, x(0)), y(t; 0, \omega, y(0))\right)$. Then, model (3.3) can be rewritten as

where

and $F:\mathbb{R}_{+}^{2}\times\mathbb{R}_{+}\to\mathbb{R}^{2}$ is given by

where $\eta = (\eta_1, \eta_2)\in\mathbb{R}_{+}^{2}$.

We can find that $F(\eta, \theta_{t}\omega)$ is locally Lipschitz with respect to $\eta = (\eta_1, \eta_2)$, and then model (3.3) possesses a unique local solution. To prove the local solution is a global one, we define the new state variable $V = x+y$. It is easy to see that $V$ satisfies the following equation:

Notice that $\frac{\alpha^{2}}{2}-\alpha\lambda z^{\ast}\in(l_1, l_2)$, and from Eq (3.4), we know

where $m_2 = \min\left\{h+l_1, d-k+l_1\right\}$. It is straightforward to check that $V$ does not blow up at any finite time, and the same happens to $x$ and $y$. Therefore, the unique local solution can be extended to a global one.

Moreover, from Eq (3.3), we know that

for all $y\ge 0$, and

for all $x\ge 0$. Thus, the unique global solution $X(t; 0, \omega, X(0))$ of random model (3.3) remains in the positive cone $\mathbb{R}_+^2$ for every initial value $X(0)\in\mathbb{R}_+^2$.

Remark 3.1. Define a mapping $\varphi_{X}: \mathbb{R}_+\times\Omega\times\mathbb{R}_+^2\to\mathbb{R}_+^2$ given by

Since the function $F$ is continuous in $(X, t)$ and is measurable in $\omega$, we obtain the $(\mathcal{B}(\mathbb{R}_{+})\times\mathcal{F}\times\mathcal{B}(\mathbb{R}_{+}^{2}), \mathcal{B}(\mathbb{R}_{+}^{2}))$-measurability of the mapping, which defines a random dynamical system generated by the solution mapping of model (3.3).

3.2.   Existence of deterministic attracting set

Theorem 3.2. There exists a deterministic compact attracting set

for the solution of model (3.3), where $V_1 = \frac{a}{M_2}e^\frac{l_1-\frac{\alpha^2}{2}}{\lambda}$, $V_2 = \frac{a}{m_2}e^\frac{l_2-\frac{\alpha^2}{2}}{\lambda}$, $\underline{x}_{1} = \frac{ae^{\frac{l_1-\frac{\alpha^2}{2}}{\lambda}}} {bV_2e^{\frac{\frac{\alpha^2}{2}-l_1}{\lambda}}+h+l_2}$ and

Proof. According to inequality (3.5), we can obtain

Also, from equation (3.4) and $\frac{\alpha^2}{2}-\alpha\lambda z^{\ast}\in(l_1, l_2)$, we can obtain

By setting $M_2 = \max\Big\{(b-c)V_2e^\frac{\frac{\alpha^2}{2}-l_1}{\lambda}+h\notag+l_2, \; d-k+\frac{\theta}{2\mu}+l_2\Big\}$, we can get

and find

According to inequalities (3.6), we know that, for every initial value $X(0)\in\mathbb{R}_+^2$ and any given $\varepsilon > 0$, there exists some time $T(\omega, X(0), \varepsilon) > 0$ such that $V(t)\le V_2+\varepsilon$ for all $t\ge T(\omega, X(0), \varepsilon)$. Therefore, we know $x+y\le V_2$ holds for every time $t$ large enough. It then follows from $\frac{\alpha^{2}}{2}-\alpha\lambda z^{\ast}\in(l_1, l_2)$ and $y\le V_2$ that

and

Therefore,

is a deterministic attracting set for the solution of model (3.3).

Remark 3.2. The existence of attracting set $B_0$ indicates that the inequalities

hold for every time $t$ large enough, where $X(t) = (x(t), y(t))$ is the solution of model (3.3). In what follows, we always believe that these inequalities are true, due to the purpose of this paper is to explore the long-time behavior of the plankton species.

Taking parameters $a = 2$, $b = 1.5$, $c = 1.4$, $k = 0.05$, $h = 0.7$, $d = 0.6$, $\theta = 0.1$, $\mu = 0.8$, $\alpha = 0.1$, $\lambda = 0.5$, $l_1 = -0.13$ and $l_2 = 0.13$ in model (3.3), we can calculate $V_1 = 0.93$, $V_2 = 6.11$ and $\underline{x}_{1} = 0.118$, and the simulation of the attracting set is shown in Figure 4.

3.3.   Internal structure of the attracting set

Theorem 3.3. For model (3.3), if the condition

holds, then the attracting set $B_0$ is reduced to a line segment on the coordinate axis. More precisely, it is

where $\underline{x}_{2} = \frac{a}{h+l_2}e^{\frac{l_1-\frac{\alpha^2}{2}}{\lambda}}$ and $\overline{x} = \frac{a}{h+l_1}e^{\frac{l_2-\frac{\alpha^2}{2}}{\lambda}}$.

Proof. According to Remark 3.2, we know that $x < V_2$. Then, it follows from

and $l_1\le\frac{\alpha^2}{2}-\alpha\lambda z^{\ast}$ that

If $cV_2e^\frac{\frac{\alpha^2}{2}-l_1}{\lambda}-(d+l_1) < 0$, we know that

and then, for every time $t$ large enough, the first equation of model (3.3) can be rewritten as

It follows from $\frac{\alpha^{2}}{2}-\alpha\lambda z^{\ast}\in(l_1, l_2)$ that

and then

Therefore, the attracting set of model (3.3) will become

which is a line segment on the coordinate axis.

Taking parameters $a = 0.9$, $b = 1.5$, $c = 0.1$, $k = 0.05$, $h = 0.7$, $d = 0.6$, $\theta = 0.1$, $\lambda = 0.5$, $\alpha = 0.1$, $\mu = 0.8$, $l_1 = -0.13$ and $l_2 = 0.13$ in model (3.3), we can calculate $\underline{x}_{2} = 0.82$ and $\overline{x} = 2.02$. The simulation of the attracting set and three trajectories with different initial values is shown in Figure 5. One can see that the trajectory of model (3.3) eventually enters the line segment $B_0$ on the coordinate axis, which indicates that the phytoplankton species will go to extinct, and only the nutrient can be persistent.

Theorem 3.4. For model (3.3), if the condition

holds, then the attracting set $B_0$ is reduced to a plane region in the first quadrant. More precisely, it is

where $\underline{x}_{3} = \frac{ae^{\frac{l_1-\frac{\alpha^2}{2}}{\lambda}}+ky_1} {bV_2e^{\frac{\frac{\alpha^2}{2}-l_1}{\lambda}}+h+l_2}$ and $y_1 = \frac{cV_1e^{\frac{\frac{\alpha^2}{2}-l_2}{\lambda}} -(d+l_2+\frac{\theta}{2\mu})}{ce^{\frac{\frac{\alpha^2}{2}-l_2}{\lambda}}}$.

Proof. From Remark 3.2, we know that $V_1-y\le x$. Then, it follows from

and $\frac{\alpha^{2}}{2}-\alpha\lambda z^{\ast}\le l_2$ that

If $cV_1e^{\frac{\frac{\alpha ^2}{2}-l_2}{\lambda}}-(d+l_2+\frac{\theta}{2\mu}) > 0$, then

For every time $t$ large enough, it follows from $\frac{\alpha^{2}}{2}-\alpha\lambda z^{\ast}\in(l_1, l_2)$ and $y_1\leq y\le V_2$ that

and then

Therefore, the attracting set of model (3.3) will become

In that case, the attracting set lies completely in the first quadrant plane.

Taking parameters $a = 2$, $b = 1.5$, $c = 1.4$, $k = 0.05$, $h = 0.7$, $d = 0.6$, $\theta = 0.1$, $\mu = 0.8$, $\alpha = 0.1$, $l_1 = -0.13$, $l_2 = 0.13$ and $\lambda = 0.5$ in model (3.3), we can calculate $\underline{x}_{3} = 0.12$, $y_1 = 0.2$, $V_1 = 0.93$ and $V_2 = 6.11$. The simulation of the attracting set and three trajectories with different initial values is shown in Figure 6. One can see that the trajectory of model (3.3) eventually enters the plane region $B_0$ in the first quadrant, which indicates that the phytoplankton species and nutrient can be simultaneously persistent.

Remark 3.3. From the expressions $V_1$ and $V_2$ shown in Theorem 3.2 and the expressions $\underline{x}_{3}$ and $y_1$ shown in Theorem 3.4, we can find that the values of $V_1$, $V_2$, $\underline{x}_{3}$ and $y_1$ will decrease with the increase of $\alpha$. That is to say, when the perturbation intensity $\alpha$ increases, the attracting set $B_0$ will move towards the origin of coordinates. Biologically speaking, the perturbation is adverse to the survival of the plankton system.

4.   Summary and discussion

We have considered two random plankton models for the plankton systems driven by bounded noise. To this end, we make use of the O-U process to ensure the random perturbations are bounded in some interval. The first random model (i.e., model (2.1)) is related to deterministic system (1.1) in which the external nutrient inflow rate $a$ is perturbed by the O-U process. The second one (i.e., model (3.3)) is related to stochastic system (1.2), which can be achieved by appropriate variable substitution associated with the O-U process.

We first proved, respectively, in Theorem 2.1 and Theorem 3.1 that the random models possess unique global solutions for any positive initial conditions. Then, we proved, respectively, in Theorem 2.2 and Theorem 3.2 the existence of attracting sets for the solutions of random model (2.1) and random model (3.3). In order to have more detailed information about the long-time behavior of the plankton species, we further investigated the internal structures of the attracting sets. Specifically, Theorem 2.3 and Theorem 3.3 state some conditions under which the attracting set is reduced to a line segment on the coordinate axis (biologically speaking, the phytoplankton species will go to extinct). Theorem 2.4 and Theorem 3.4 state some conditions under which the attracting set is reduced to a plane region in the first quadrant (biologically speaking, the phytoplankton species can be persistent).

It is important to point out that the attracting sets for the solutions of model (2.1) (show in Theorems 2.2, 2.3 and 2.4) do not depend on the intensity of the perturbation, but the attracting sets for the solutions of model (3.3) (show in Theorems 3.2, 3.3 and 3.4) will move towards the origin of coordinates when the perturbation intensity $\alpha$ increases. In Figure 7, by taking initial value (4, 2); parameters $a = 2$, $b = 1.5$, $c = 1.4$, $k = 0.05$, $h = 0.7$, $d = 0.6$, $\theta = 0.1$, $\mu = 0.8$, $l_1 = -0.13$, $l_2 = 0.13$, $\lambda = 0.5$; and different noise intensities $\alpha = 0.1$, $\alpha = 0.4$ and $\alpha = 0.7$, we show trajectories of model (3.3). One can see from Figure 7 that the region that the trajectory finally enters will move towards the origin of coordinates when the perturbation intensity $\alpha$ increases.

The results in the present paper seem to be able to help us better understand the dynamics of the plankton system in a stochastic sense. One can further use the O-U process to model the real bounded fluctuations existing in other ecological systems.

Acknowledgments

This work was jointly supported by grants from the National Natural Science Foundation of China (62173161, 12161005, 11801224) and Natural Science Foundation of Jiangsu Province (BK20180856).

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

In this section, we will recall briefly some useful definitions and results about the O-U process and random dynamical systems to make our presentation as complete as possible.

A.1.   O-U process

Let $W$ be a two sided Wiener process. Kolmogorov's theorem ensures that $W$ has a continuous version, which we will denote by $\omega$, whose canonical interpretation is as follows: Let $\Omega$ be defined by

$\mathcal{F}$ be the Borel $\sigma$-algebra on $\Omega$ generated by the compact open topology ^[35] and $\mathbb{P}$ be the corresponding Wiener measure on $\mathcal{F}$. We consider the Wiener shift flow given by

Then, $(\Omega, \mathcal{F}, \mathbb{P}, \left\{\theta_{t}\right\}_{t\in\mathbb{R}})$ is a metric dynamical system ^[35].

Now, let us introduce the following O-U process, defined on $(\Omega, \mathcal{F}, \mathbb{P}, \left\{\theta_{t}\right\}_{t\in\mathbb{R}})$ as the random variable given by

which solves the Langevin equation ^[35,36]

where $\lambda > 0$ is a mean reversion constant that represents the strength with which the process is attracted by the mean or, in other words, how strongly our system reacts under some perturbation. There are some important properties ^{[28,35,36,37]} of the O-U process:

Proposition A.1. If there exists a $\theta_{t}$-invariant set $\tilde{\Omega}\in\mathcal{F}$ of $\Omega$ of full $\mathbb{P}$-measure, then

● for a.e. $\omega\in\tilde{\Omega}$ and every $\lambda > 0$,

● for a.e. $\omega\in\tilde{\Omega}$ and all $t\in\mathbb{R}$,

A.2.   Random dynamical system (RDS)

Let $(\Omega, \mathcal{F}, \mathbb{P})$ be a probability space and $(\mathcal{X}, \left\|\cdot\right\|_{\mathcal{X}})$ be a separable Banach space. The following definitions about the RDS can be found in ^[35,38].

Definition A.1. An RDS on $\mathcal{X}$ consists of two ingredients: (a) a metric dynamical system $(\Omega, \mathcal{F}, \mathbb{P}, \{\theta_{t}\}_{t\in\mathbb{R}})$ with a family of mappings $\theta:\Omega\to\Omega$ such that

● $\theta_{0} = \mathrm{d}_\Omega$,

● $\theta_{t+s} = \theta_{t}\circ\theta_s$ for all $t, s\in\mathbb{R}$,

● the mapping $(t, \omega)\mapsto\theta_{t}\omega$ is measurable, and

● the probability measure $\mathbb{P}$ is preserved by $\theta_{t}$, i.e., $\theta_{t}\mathbb{P} = \mathbb{P}$;

and (b) a mapping $\psi:[0, +\infty)\times\Omega\times\mathcal{X}\to\mathcal{X}$ which is $(\mathcal{B}([0, +\infty))\times\mathcal{F}\times\mathcal{B}(\mathcal{X}), \mathcal{B}(\mathcal{X}))$ -measurable, such that for a.e. $\omega\in\Omega$,

● the mapping $\varphi(t, \omega):\mathcal{X}\to\mathcal{X}, \; x\mapsto\varphi(t, \omega)x$ is continuous for every $t\geq0$,

● $\varphi(0, \omega)$ is the identity operator on $\mathcal{X}$, and

● $\varphi(t+s, \omega) = \varphi(t, \theta_s\omega)\varphi(s, \omega)$ for all $t, s\geq0$.

Definition A.2. A random set $K$ is a measurable subset of $\mathcal{X}\times\Omega$ with respect to the product $\sigma$-algebra $\mathcal{B}(\mathcal{X})\times\mathcal{F}$. Moreover, $K$ will be called a closed or a compact random set if $K(\omega) = \left\{x:(x, \omega)\in K\right\}$, $\omega\in\Omega$, is closed or compact for $\mathbb{P}$-almost all $\omega\in\Omega$, respectively.

Definition A.3. A bounded random set $K(\omega)\subset\mathcal{X}$ is said to be tempered with respect to $\{\theta_{t}\}_{t\in\mathbb{R}}$ if for a.e. $\omega\in\Omega$ and all $\lambda > 0$,

Definition A.4. A random set $B(\omega)\subset\mathcal{X}$ is called a random absorbing set in $\mathcal{E}(\mathcal{X})$, if for any $E\in\mathcal{E}(\mathcal{X})$ and a.e. $\omega\in\Omega$, there exists $T_E(\omega) > 0$ such that for all $t\geq T_E(\omega)$,

Definition A.5. Let $\{\varphi(t, \omega)\}_{t\geq0, \omega\in\Omega}$ be an RDS over $(\Omega, \mathcal{F}, \mathbb{P}, \{\theta_{t}\}_{t\in\mathbb{R}})$ with state space $\mathcal{X}$, and let $A(\omega)$ be a random set. Then, $\mathcal{A} = \{A(\omega)\}_{\omega\in\Omega}$ is called a global random $\mathcal{E}$-attractor (or pullback $\mathcal{E}$-attractor) for $\{\varphi(t, \omega)\}_{t\geq0, \omega\in\Omega}$ if

● $A(\omega)$ is a compact set of $\mathcal{X}$ for a.e. $\omega\in\Omega$;

● $\varphi(t, \omega)A(\omega) = A(\theta_t\omega)$ holds for a.e. $\omega\in\Omega$ and all $t\geq0$;

● for a.e. $\omega\in\Omega$ and any $E\in\mathcal{E}(\mathcal{X})$,

where $\mathrm{dist}(G, H)_{\mathcal{X}} = \sup_{g\in G}\inf_{h\in H}\left\|g-h\right\|_{\mathcal{X}}$ is the Hausdorff semi-metric for $G, H\subseteq\mathcal{X}$.

Proposition A.2. ^[39,40] Let $B\in\mathcal{E}(\mathcal{X})$ be an absorbing set for the continuous RDS $\{\varphi(t, \omega)\}_{t\geq0, \omega\in\Omega}$ which is closed and satisfies the asymptotic compactness condition for a.e. $\omega\in\Omega$, i.e., each sequence $x_{n}\in\varphi(t_n, \theta_{-t_n}\omega)B(\theta_{-t_n}\omega)$ has a convergent subsequence in $\mathcal{X}$ when $t_{n}\to\infty$. Then, $\varphi$ has a unique global random attractor $\mathcal{A} = \left\{\mathcal{A}(\omega)\right\}_{\omega\in\Omega}$ with component subsets

Proposition A.3. ^[21] Let $\varphi_u$ be an RDS on $\mathcal{X}$. Suppose that the mapping $\mathcal{T}:\Omega\times\mathcal{X}\to\mathcal{X}$ possesses the following properties:

● for fixed $\omega\in\Omega$, the mapping $\mathcal{T}(\omega, \cdot)$ is a homeomorphism on $\mathcal{X}$;

● for fixed $x\in\mathcal{X}$, the mappings $\mathcal{T}(\cdot, x)$ and $\mathcal{T}^{-1}(\cdot, x)$ are measurable.

Then, the mapping

is a conjugated RDS.