The relationship between macro-socioeconomics determinants and COVID-19 vaccine distribution

1.   Introduction

Plankton, the organisms that the freely floating and weakly swimming in aquatic environments, occupy the first tropic level and the second trophic level of any aquatic food chains. Phytoplankton are the photosynthetic microorganisms and commonly unicellular and microscopic in size, and zooplankton are the heterotrophic plankton that live on phytoplankton. In addition to recognizing the importance of plankton for the wealth of the aquatic ecosystems and ultimately for the planet itself ^[1], the variation of plankton biomass is an important factor influencing the real aquatic environments, and understanding of plankton dynamics can be helpful to estimate the productivity of aquatic ecosystems ^[2,3] and regulate the balance of plankton ecosystem. However, planktonic blooms can occur under some conducive environments, which may cause seriously environmental issues and threat to human health. But the processes underlying the formation of planktonic blooms are not yet well understood. In this respect, thus, the great effort has been made towards the understanding of the complex dynamics of plankton, and then mathematical models can be acted as a useful tool to investigate the dynamics of plankton ecosystem, which can provide a deeper understanding of the dynamic mechanisms of changes in plankton populations.

Actually, many mathematical models have been constructed to study the dynamical behaviors of plankton since the pioneering work of Riley et al. ^[4], and many physical and biological processes underlying the mechanisms of plankton dynamics in the aquatic environments have been investigated ^{[1,5,6,7,8,9,10]}. For example, in order to study how the nutrient affects the dynamics of phytoplankton blooms, Huppert et al. ^[1] presented a simple nutrient-phytoplankton model and identified an important threshold effect that a bloom will only be triggered when nutrients exceed a certain defined level using mathematical model analysis. Caperon ^[6] concluded that the time-lag effect exists in the growth process of phytoplankton, and further suggested that models play an important role in understanding the growth dynamics of phytoplankton characterized by time delays. Lin et al. ^[8] used a nutrient-phytoplankton-zooplankton model to examine the patterns and consequences of adaptive changes in plankton body size and suggested that evolutionary interactions between phytoplankton and zooplankton may have contributed to observed changes in phytoplankton sizes and associated biogeochemical cycle over geological time scales.

In recent years, the dynamical behaviors of phytoplankton-zooplankton systems with various biological factors, such as stability, bifurcation and spatiotemporal pattern, have been explored extensively ^[6,10,11,12]. Nevertheless, some phytoplankton species are harmful phytoplankton that can produce potent toxic or allelopathic substances during phytoplankton blooms ^[13], which can affect species interaction by suppressing the growth and establishment of other phytoplankton species ^[14]. Moreover, some laboratory experiments ^[15,16], as well as field observation ^[17] have suggested that the toxicity may be as a strong mediator in the zooplankton feeding rate. As a result, some researchers have taken this important factor of toxic production released by TPP into account when studying the phytoplankton-zooplankton systems ^[18,19,20]. For example, Scotti et al. ^[18] indicated that a toxic phytoplankton may destabilize the spatially homogeneous coexistence and trigger the formation of spatial pattern, and further concluded that local blooms more likely occur when the strength of the toxicity is of a certain level. Additionally, some results from field observations and model analysis concluded that the toxic substances can affect the interaction between phytoplankton and zooplankton and reduce the growth of zooplankton, indicating TPP may act as a biological control way for the termination of planktonic blooms ^{[21,22,23,24]}. Sarkar et al. ^[25] proposed the following mathematical model consisting of two harmful phytoplankton and zooplankton species:

where $P_{1}(t), \; P_{2}(t)$ and $Z(t)$ are the population densities of two harmful phytoplankton and zooplankton species at time $t$, respectively; $r_{1}$ and $r_{2}$ denote the intrinsic growth rates of two TPP, respectively; $K_{1}$ and $K_{2}$ are their corresponding environmental carrying capacities; $a_{1}$ and $a_{2}$ represent be the inhibitory effects of two harmful phytoplankton species; $\alpha$ and $\beta$ are the maximum zooplankton ingestion rates for both two TPP species, respectively; $m$ and $n$ are the maximum zooplankton conversion rates, respectively; $d$ is the natural death rate of zooplankton; $\gamma$ and $\delta$ are the rates of toxin liberation by two TPP species, respectively; $a$ and $b$ denote the half-saturation constants for two TPP species. In ^[25], the authors studied the asymptotic stability of the system (1) and claimed that the presence of two harmful phytoplankton has a positive impact for the termination of planktonic blooms.

In the real world, however, the unpredictability and ubiquity of environmental fluctuations in the natural aquatic ecosystems, for example, the necessary nutrient availability, water temperature, light and turbulence, can greatly cause the growths of plankton populations to experience random fluctuations. Systems with such kinds of environmental fluctuations can be described by stochastic differential equations, which play a significant role in the population dynamics as they can provide some additional degree of realism compared to their corresponding deterministic counterparts ^[26]. Thus, stochastic population systems, as an important application in ecological and biological systems, have attracted increasing attention ^{[27,28,29,30,31,32,33]}. Especially, stochastic plankton systems with white noise have been the common area of interest among researchers ^{[9,23,30,34,35,36,37,38,39,40,41,42]} in recent years and many interesting results have been shown. For example, Sarkar and Chattopadhayay ^[34] proposed a toxic phytoplankton-non-toxic phytoplankton-zooplankton with stochastic perturbation around the positive equilibrium, and they concluded that TPP and stochastic fluctuations can significantly affect the coexistence of species. Yu et al. ^[41] investigated a nutrient-phytoplankton system with TPP under environmental fluctuations, and they obtained some conditions for extinction, persistence and the existence of ergodic stationary distribution. All these works greatly stimulate researchers to explore the way how environmental stochasticity and toxin production affect the coexistence and survival prospect of plankton populations in the presence of harmful phytoplankton. Obviously, it is meaningful to further incorporate the environmental fluctuations into the underlying model (1). Moreover, there are few literatures to study the dynamics of the stochastic phytoplankton-zooplankton system with two harmful phytoplankton, and the dynamics of the stochastic phytoplankton-zooplankton system with two harmful phytoplankton is still not very clear currently. Hence, we mainly present the influence of the effects of environmental white noise and toxic liberation rates produced by two TPP on the dynamics of phytoplankton-zooplankton system in this paper. Motivated by the works above, we assume that the intrinsic growth rates of two harmful phytoplankton and the death rate of zooplankton are influenced by the environmental fluctuations effect, and thus introduce the white noise into underlying system (1), resulting in the following form:

where $B_{i}(t)$ are mutually independent standard Brownian motions with $B_{i}(0) = 0$ ^[43], and $dB_{i}(t)$ are standard white noise and $\sigma_{i}^{2}(t)$ are their intensities, $i = 1, 2, 3$.

By now, we have successful introduced a stochastic phytoplankton-zooplankton system with two toxic phytoplankton focusing on the effects of environmental stochasticity and TPP, and our research questions include: (i) How does environmental stochasticity affects the dynsmics of plankton populations? (ii) What influences the peak of the outbreaks of planktonic blooms in a fluctuating environment? The rest of this paper is organized as follows: Section 2 presents the basic assumptions firstly, and then we investigate the existence and uniqueness of global positive solutions, and apply the It$\hat{o}$'s formula to obtain the sufficient conditions for the extinction and persistence in the mean of system (2), and the existence of a unique ergodic stationary distribution by establishing a appropriate stochastic Lyapunov function. A series of numerical simulations are carried out to further study the dynamics of system (2) in section 3. In section 4, we summarize the results and present our conclusions.

2.   Main results

In this section, we investigate mainly the existence and uniqueness of global positive solutions, the extinction and persistence in the mean of system (2), and discuss the positive recurrence and ergodic property of system (2) as well.

2.1.   Preliminaries

Denote $\mathbb{R}_{+} = [0, +\infty)$ and $\mathbb{R}_{+}^{n} = \{(x_{1}, \cdots, x_{n})\in \mathbb{R}^{n}:x_{i} > 0, i = 1, 2, \cdots, n\}$, and $|x| = \sqrt{\sum_{i = 1}^{n}x_{i}^{2}}$. Throughout this paper, unless otherwise indicated, we always assume that $(\Omega, \mathscr{F}_{t}, \{\mathscr{F}_{t}\}_{t\geq0}, \mathbb{P})$ is a completed probability space with a filtration $\{\mathscr{F}_{t}\}$ satisfying the usual normal conditions (i.e., it is right-continuous and increasing while $\{\mathscr{F}_{0}\}$ contains all $\mathbb{P}$-null sets). For convenience, if $\varphi(t)$ is a integrable function on $\mathbb{R}_{+}$, we define $\langle\varphi\rangle = \frac{1}{T}\int_{0}^{T}\varphi(s)ds, T > 0$.

Generally, we consider the $n$-dimentional stochastic differential equation:

with initial value $x(t_{0}) = x_{0}\in \mathbb{R}^{n}$, where $B(t)$ denotes a $n$-dimentional standard Brownian motion defined on the completed probability space $(\Omega, \mathscr{F}_{t}, \{\mathscr{F}_{t}\}_{t\geq0}, \mathbb{P})$. Denote by $C^{2, 1}(\mathbb{R}^{n}\times \mathbb{R_{+}}; \mathbb{R_{}})$ the family of all non-negative functions $V(x, t)$ defined on $\mathbb{R}^{n}\times \mathbb{R_{+}}$ such that they are continuously twice differentiable in $x$ and once in $t$. Define a differential operator $\mathscr{L}$ associated with Eq (3) by ^[44] as follows:

Let $V(x, t)\in C^{2, 1}(\mathbb{R}^{n}\times \mathbb{R}, \mathbb{R})$, then

where

By It$\hat{o}$'s formula, if $x(t)\in \mathbb{R}^{n}$, then

Next, we introduce the criterion of positive recurrent and the ergodic properties. Before it, we consider the stochastic equation:

where $x(t)$ is a homogeneous Markov process in $n$-dimentional Euclidean space $\mathbb{R}^{n}$. The diffusion matrix is $D(x) = (d_{ij}(x))$ and $(d_{ij}) = \sum_{r = 1}^{k}\sigma_{r}^{i}(x)\sigma_{r}^{j}(x)$. Thus, a lemma which describes the criterion of stationary distribution can be given.

Lemma 2.1 (see ^[45]) Suppose that there exists a bounded open set $E\subset \mathbb{R}^{n}$ with a smooth regular boundary $\Theta$ satisfying the following conditions:

(i) the diffusion matrix $D(x)$ is strictly positive definite for all $x\in E$;

(ii) there exists a non-negative $C^{2}$-function $V(x)$ and a positive constant $M$ such that $\mathscr{L}V\leq-M$ for $\forall x\in \mathbb{R}^{n}/E$.

Then there exists a solution $x(t)$ of the system (4) which is a stationary Markov process with a stationary distribution $\mu(\cdot)$ and for any integrable function $g(\cdot)$ with respect the measure $\mu$, we have

2.2.   Existence and uniqueness of global positive solutions

Before investigating the stochastic dynamics of system (2), we should first guarantee whether the solution is global and positive. Therefore, based on the biological interpretation, in this subsection, we just take the non-negative solutions into account for system (2) and discuss the existence of global positive solutions in system (2). The following result can be presented.

Theorem 2.2 For any given initial value $(P_{1}(0), P_{2}(0), Z(0))\in \mathbb{R}_{+}^{3}$, system (2) exists a unique solution $(P_{1}(t), P_{2}(t), Z(t))$ on $\mathbb{R}_{+}$ and the positive solution will remain in $\mathbb{R}_{+}^{3}$ with probability one, that is, $(P_{1}(t), P_{2}(t), Z(t))\in \mathbb{R}_{+}^{3}$ for all $t > 0$ almost surely.

Proof. From the method of the Lemma 2.1 in ^[46], it is obvious to find that, for any given initial value $(P_{1}(0), P_{2}(0), Z(0))\in \mathbb{R}_{+}^{3}$, all coefficients of system (2) are locally Lipschitz continuous and the system admits a unique local solution $(P_{1}(t), P_{2}(t), Z(t))$ on $[0, \tau_{e})$, where $\tau_{e}$ represents the explosion time. In the following, thus, we need to illustrate the solution is global, that is, we only need to prove $\tau_{e} = \infty$ a.s. Let $n_{0}\geq 1$ enough large such that $(P_{1}(0), P_{2}(0), Z(0))\in [\frac{1}{n_{0}}, n_{0}]$, and for each integer $n\geq n_{0}$, we define the stopping time by the following form:

and the set $\inf \emptyset = \infty$ ($\emptyset$ denotes the empty set). Obviously, one can obtain that $\tau_{n}$ is increasing as $n\rightarrow \infty$. Thus, let $\tau_{\infty} = \lim_{t\rightarrow +\infty}\tau_{n}$, then $\tau_{\infty}\leq \tau_{e}$ a.s. If we can show that $\tau_{\infty} = \infty$ a.s., then $\tau_{e} = \infty$ and $(P_{1}(0), P_{2}(0), Z(0))\in \mathbb{R}_{+}^{3}$ a.s. for all $t\geq0$. In other words, to complete the proof, we only need to proof $\tau_{\infty} = \infty$ a.s. If the statement is false, then there exist two constants $T > 0$ and $\varepsilon\in (0, 1)$ such that $\mathbb{P}\{\tau_{\infty}\leq T\} > \varepsilon.$ Hence, for all $n\geq n_{1}$, there exists an integer $n_{1}\geq n_{0}$ such that $\mathbb{P}\{\tau_{\infty}\leq T\}\geq\varepsilon.$

Define a $C^{2}$-function $\bar{V}:\mathbb{R}_{3}^{+}\rightarrow \mathbb{R}_{+}$ by

Obviously, the function $\bar{V}(P_{1}, P_{2}, Z, )$ is non-negative since the inequality $x-1-\log x > 0$ holds for all $x > 0$. Applying It$\hat{o}$'s formula to $\bar{V}(P_{1}, P_{2}, Z)$ yields

where $\mathscr{L}\bar{V}:\mathbb{R}_{+}^{3}\rightarrow \mathbb{R}$ is defined by

where

Notice that

for all $Z > 0$, then one can obtain that

where

The remainder of the proof follows that in the Theorem 3.3 ^[37], here, we omit it. This completes the proof.

2.3.   Persistence in the mean and extinction

Based on Theorem 2.2, we will discuss the extinction and persistence in the mean of system (2) and derive sufficient conditions for them in this subsection.

Theorem 2.3 Suppose that $(P_{1}(t), P_{2}(t), Z(t))$ is the solution of system (2) with initial value $P_{1}(0), P_{2}(0), Z(0))\in \mathbb{R}_{+}^{3}$, then we have the following results:

I Let $A = r_{1}-\frac{\sigma_{1}^{2}}{2}, \theta_{1} = A-\frac{a_{1}K_{2}}{r_{2}}B-\alpha D$, the species $P_{1}(t)$ is

(i) extinct if $A < 0$;

(ii) persistence in the mean if $\theta_{1} > 0$.

II Let $B = r_{2}-\frac{\sigma_{2}^{2}}{2}, \theta_{2} = B-\frac{a_{2}K_{1}}{r_{1}}A-\beta D$, the species $P_{2}(t)$ is

(i) extinct if $B < 0$;

(ii) persistence in the mean if $\theta_{2} > 0$.

III Let $C = d+\frac{\sigma_{3}^{2}}{2}, D = \frac{mK_{1}}{r_{1}}A+\frac{nK_{2}}{r_{2}}B-C$, the species $Z(t)$ is

(i) extinct if $D < 0$;

(ii) persistence in the mean if

Proof. I. Applying the It$\hat{o}$'s formula to system (2) yields

and

Integrating the above from 0 to $t$ and dividing $t$ on both sides, we have

and

where

Moreover, the quadratic variation of $M_{i}(t)\; (i = 1, 2, 3)$ satisfy

By the strong law of large numbers for martingales ^[45] yields

Thus, according to (5), we have

which implies that $\lim_{t\rightarrow \infty}P_{1}(t) = 0\; \; a.s.$ This completes the proof of (i).

Next, we give the prove of (ii). By making some estimations of (5), we have

where $B$ and $D$ will be determined later. In addition, since the fact that

Thus, from the properties of the limit, for arbitrary $\epsilon_{1}$, there exists a constant $T_{1} > 0$ such that

Substituting above inequalities into (9) and for all $t\geq T_{1}$, we have

where $\theta_{1} = A-\frac{a_{1}K_{2}}{r_{2}}-\alpha D.$ From the Lemma 4 in ^[47], we have

This completes the proof of (ii).

II. The case II are similar to the case I, here we omit it.

III. From (I), (II), one can obtain that

According to Eq (7), we have

Combining with (8) and (10), and taking upper limit on both sides of (11) yields

which implies $\lim_{t\rightarrow \infty}Z(0) = 0 \; \; a.s.$ This completes the proof of (i).

Now, we show the proof of (ii). Computing $(5)\times \frac{mK_{1}}{r_{1}}+(6)\times \frac{nK_{2}}{r_{2}}+(7)$, we have

By the strong law of large numbers for martingales, we can derive that

From the Lemma 2.3 in ^[48], we can obtain that

Thus, taking the limit superior in (12) and from the lemma 4 in ^[49], one can see that

which implies

This completes the proof of (ii).

2.4.   Positive recurrence and ergodic property

In this subsection, by constructing a suitable Lyapunov function, we discuss the positive recurrence of system (2) and some sufficient conditions for the existence and uniqueness of stationary distribution are obtained.

Theorem 2.4 If the following condition

holds. Then for any given initial value $(P_{1}(0), P_{2}(0), Z(0))\in \mathbb{R}_{+}^{3}$, the solution $(P_{1}(t), P_{2}(t), Z(t))$ of system (2) admits a unique ergodic stationary distribution.

Proof. In order to prove this theorem, we only need to prove both two conditions of Lemma 2.1 one by one. Obviously, we can obtain that the diffusion matrix $d_{ij}(x, K) = diag\{\sigma_{1}^{2}, \sigma_{2}^{2}, \sigma_{3}^{2}\}$ of system (2) is positive definite, which implies that the condition (i) in Lemma 2.1 holds.

Next, we prove the condition (ii) in Lemma 2.1. We construct a $C^{2}$-function $\tilde{V}:\mathbb{R}_{+}^{3}\rightarrow \mathbb{R}$ as follows:

where $p$ is a positive constant satisfying $0 < p < \frac{d}{3^{p}2[\sigma_{1}^{2}\vee\sigma_{2}^{2}\vee\sigma_{3}^{2}]}$ and $M > 0$ will be determined later. Since the function $\tilde{V}(P_{1}, P_{2}, Z)$ is continuous, then there exists a unique point $(\tilde{P_{1}}, \tilde{P_{2}}, \tilde{Z})$ in $\mathbb{R}_{+}^{3}$ which is the minimum value of $\tilde{V}(P_{1}, P_{2}, Z)$. Thus, construct a non-negative $C^{2}$-function $V:\mathbb{R}_{+}^{3}\rightarrow \mathbb{R}$ by

Using the generalized It$\hat{o}$'s formula, one can see that

Noting that the inequality $|\sum_{i = 1}^{n}c_{i}|^{n}\leq k^{n-1}\sum_{i = 1}^{n}|c_{i}|^{n}$, we can obtain that

where

and

Thus,

Considering the following compact subset $E$:

where $\varepsilon$ is a sufficiently small constant satisfying the following conditions:

where

Then

where

Now, we need to prove the negativity of $\mathscr{L}V(P_{1}, P_{2}, Z)$ for any $(P_{1}, P_{2}, Z)\in \mathbb{R}_{+}^{3}/E$.

Case I: If $(P_{1}, P_{2}, Z)\in E_{1}$, then $P_{1}Z\leq \varepsilon Z \leq \varepsilon \left(1+Z^{2}\right)$, and one can obtain that

By (13), we have $\mathscr{L}V(P_{1}, P_{2}, Z)\leq-1$ for any $(P_{1}, P_{2}, Z)\in E_{1}$.

Case II: If $(P_{1}, P_{2}, Z)\in E_{2}$, we have

By (14), we have $\mathscr{L}V(P_{1}, P_{2}, Z)\leq-1$ for any $(P_{1}, P_{2}, Z)\in E_{2}$.

Case III: If $(P_{1}, P_{2}, Z)\in E_{3}$, the similar analysis to case I, we can obtain that

By (15), we have $\mathscr{L}V(P_{1}, P_{2}, Z)\leq-1$ for any $(P_{1}, P_{2}, Z)\in E_{3}$.

Case IV: If $(P_{1}, P_{2}, Z)\in E_{4}$, the similar analysis to case II, we have

which follows from (16), then we have $\mathscr{L}V(P_{1}, P_{2}, Z)\leq-1$ for any $(P_{1}, P_{2}, Z)\in E_{4}$.

Case V: If $(P_{1}, P_{2}, Z)\in E_{5}$, then $P_{1}Z\leq \varepsilon P_{1}\leq \varepsilon\left(1+P_{1}^{2}\right), P_{2}Z\leq \varepsilon P_{2}\leq \varepsilon\left(1+P_{2}^{2}\right)$, we have

According to (17), we have $\mathscr{L}V(P_{1}, P_{2}, Z)\leq-1$ for any $(P_{1}, P_{2}, Z)\in E_{5}$.

Case VI: If $(P_{1}, P_{2}, Z)\in E_{6}$, we have

From (18), we have $\mathscr{L}V(P_{1}, P_{2}, Z)\leq-1$ for any $(P_{1}, P_{2}, Z)\in E_{6}$.

Hence, the condition (ii) of Lemma 2.1 is verified. Thus, system (2) admits a unique ergodic stationary distribution. This completes the proof.

3.   Numerical simulations

In this section, based on the Milstein's Higher Order Method mentioned in ^[50], some numerical simulations are carried out to further study the effects of the environmental noise and toxin rate released by TPP on the dynamics of system (2). Due to in a phytoplankton-zooplankton system that takes toxic phytoplankton into consideration, the intrinsic growth rates of two toxic phytoplankton and the death rate of zooplankton are parameters that are most susceptible to environmental influences and are relatively important. Therefore, we only consider the intrinsic growth rates and death rate affected by white noise. In the following numerical simulations, unless otherwise specified, the following parameter values are always used: $r_{1} = 0.55$, $r_{2} = 0.5$, $a_{1} = 0.005$, $a_{2} = 0.004$, $\alpha = 0.15$, $\beta = 0.15$, $m = 0.09$, $n = 0.075$, $d = 0.09$, $a = 0.1$, $b = 0.12$, $K_{1} = 2$, $K_{2} = 5$, the initial value $(P_{1}(0), P_{2}(0), Z(0)) = (0.5, 0.5, 0.5)$, and other parameters are chosen as control parameters. In order to study how the white noise and two TPP affect the dynamics of system (2), we firstly consider system (2) does not experience the white noise, that is, system (2) becomes its corresponding deterministic system. According to ^[25], we can obtain that the system (1) possesses a unique positive interior equilibrium $E^{*}(0.8178, 1.8066, 2.1071)$ which is locally asymptotically stable, depicting the coexistence of all three species.

Next, the impact of white noise on the stochastic dynamics of system (2) will be shown. We first fix $\sigma_{1} = 0.1$ and $\sigma_{3} = 0.1$, and let $\sigma_{2}\; (0\leq \sigma_{2} \leq 1.2)$ vary to see how the white noise influences the survival of plankton populations. According to the Theorem 2.3, we can obtain that all three species of system (2) will undergo extinction when the white noise reaches some a critical value. Obviously, we can find from Figure 1(a) that the species $P_{2}(t)$ of system (2) is always persistence in the mean in the area of $I$ and persistence in the mean or extinct alternating in the area of $II$, but species $P_{2}(t)$ becomes die out rapidly in the space $III$ when $\sigma_{2}$ is beyond $\sigma_{2} = 1$. Figure 1(b) depicts that the stochastic dynamical behaviors of species $P_{2}(t)$ with respect to Figure 1 (a) for $\sigma_{2} = 0.1, \sigma_{2} = 0.8$, $\sigma_{2} = 1.1$ and $t$ on $[1000,2000]$, respectively. Moreover, we can observe from Figure 1(b) that, with the increase in the magnitude of the environmental fluctuations, the random variation of plankton density becomes more significant, which implies that white noise can accelerate the stochastic oscillation of plankton density. For example, let $\sigma_{2} = 0.1$, it is not difficult to find that the two TPP and zooplankton of system (2) can coexist at a relatively stable state and their densities exhibit oscillation around the deterministic steady state values $P_{1}^{*} = 0.8178$, $P_{2}^{*} = 1.8066$ and $Z^{*} = 2.1071$, respectively (see Figure 2(a)-(c)). Actually, following the Theorem 2.3 and Theorem 2.4, system (2) is persistence in the mean and has a unique ergodic stationary distribution, which are consistent with our numerical analysis. From the stationary distribution of all three species, it can be seen clearly that they are distributed normally around the values 0.8, 1.8 and 2.1, respectively, which illustrates that the standard deviation $\sigma_{1}, \sigma_{2}$ and $\sigma_{3}$ can keep processes $P_{1}(t)$, $P_{2}(t)$ and $Z(t)$ moving around the solution of system (1). In other words, system (2) can preserve some stability in the random sense when the intensities of white noise are relatively weak. Now let $\sigma_{2}$ vary within some a level, it can be concluded that the weaker the environmental fluctuations are, the closer the solutions of system (2) are to steady state $E^{*}$ (The Figures here are not given due to the similarity to Figure 2). However, when we increase the density of environmental forcing to $\sigma_{2} > 1$, for example $\sigma_{2} = 1.1$, we can easily get from Theorem 2.3 that $B = -0.105 < 0$, which implies that species $P_{2}(t)$ of system (2) tends to go rapid extinct, even if its corresponding deterministic system (1) still presents obvious stability, indicating a different phenomenon from its deterministic system (see Figure 3). This also shows that white noise intensity can help to control the density of toxic phytoplankton. Comparing Figures 1, 2 and 3, it is obvious to find that the intensity of white noise can not only aggravate the stochastic oscillation of plankton density, but also significantly change the dynamics of the plankton system. That is, a high white noise intensity can accelerate the extinction of the plankton populations, which implies that the white noise can help control the biomass of plankton populations and provide a guide for us to the termination of planktonic blooms. This is consistent with the results obtained by the work of Sarkar and Chattopadhayay ^[34], who demonstrated the controlling of plankonic blooms by artifical eutrophication or the intensity of white noise from their experimental and field observations. Thus, it is worth pointing out that the results from the Figures 1, 2 and 3 can support that the plankton systems incorporating white noise can better simulate planktonic blooms than its corresponding deterministic counterparts. Similarly, if we impose the intensities of white noise on species $P_{1}(t)$ and $Z(t)$, respectively, we can easily obtain the similar results, thus here we omit it.

On the one hand, in order to study how the effect of one toxin liberation rate on population density dynamic evolution trend under the environmental fluctuations, we choose $\gamma$ as a control parameter and all other parameters are the same as Figure 2. Clearly, we can observe from Figure 4 that the species $P_{1}(t)$ and $P_{2}(t)$ are persistence in the mean and their biomass will increase as the increasing value of $\gamma$, while species $Z(t)$ undergoes extinction when $\gamma$ beyond a certain value, here the colorbars denote the biomass of species $P_{1}(t)$, $P_{2}(t)$ and $Z(t)$, respectively. Actually, it is easy to obtain that $\Pi > 0$ under the condition of $0\leq\gamma < 0.3846$ and $\Pi = 0$ if and only if $\gamma \approx 0.3846$, which indicates that system (2) is persistence in the mean and exists the stationary distribution under the condition of $\gamma < 0.3846$. In contrast, species $Z(t)$ of system (2) will die out. Therefore, it can be asserted that TPP can significantly affect the coexistence of all the three species. For precisely, we take three different values of $\gamma\; (\gamma = 0.1, 0.2, 0.35)$, then system (2) has a unique stationary distribution. Figure 5 depicts the relative frequency density of $P_{1}(t), P_{2}(t)$ and $Z(t)$ with these different values, respectively, where the smoothed curves are the probability density functions of system (2). More importantly, we can obtain the result from the Figure 5 that with the increasing value of $\gamma$, the distributions of two TPP appear closer to the normal distribution, but the distribution of zooplankton becomes more skewness, which implies that the increase of the liberation rate $\gamma$ can increase the survival chance of two harmful phytoplankton but decrease the biomass of zooplankton. Additionally, it can be seen from Figure 5 that the peak value of the probability density functions of system (2) will be higher as $\gamma$ increases. All these results indicate that for the value of the toxin liberation rate $\gamma$ satisfying the condition of Theorem 2.4, its enhancement will contribute to the persistence in the mean of system (2) though the termination of planktonic blooms. In addition, from the Theorems 2.3 and 2.4, we want to know what happens if $\Pi < 0$? Selecting $\gamma = 0.39$, we can easily check $\Pi\approx -0.0054 < 0$ and $\Gamma\approx -0.0194 < 0$, which imply that the conditions of the Theorems 2.3 and 2.4 are not satisfied, but system (2) has a stationary distribution, indicating all the three species are persistence in the mean (see Figure 6). However, when we choose $\gamma = 0.42$, similarly, we can obtain that $\Pi\approx -0.0354 < 0$ and $\Gamma\approx -0.0494 < 0$. Obviously, the two toxic phytoplankton can coexist while the zooplankton tends to go extinct (see Figure 7). For the case of $\gamma = 0.06$ and changing the value of $\delta$, we can easily get the similar results, which are omitted here.

Finally, the combined effects of two toxin liberation rates on the dynamics of system (2) are studied as well. Figure 8 depicts how the combined role of $\gamma$ and $\delta$ affect the dynamics of system (2), where the red smoothed curves are probability density functions of system (2). By a straightforward computation, the condition of Theorem 2.4 can be verified, which means system (2) has a unique stationary distribution (see Figure 8). Comparing Figure 5 and Figure 8, one can see that the mean values of two harmful phytoplankton are larger than the case of $\gamma = 0.2, \; \delta = 0.07$, while that of zooplankton is smaller than that case, indicating the high abundance of toxic phytoplankton because of the low level or the absence of other toxic phytoplankton, while the peaks of zooplankton population come down. Thus, we can obtain that, by controlling two harmful phytoplankton, the mean values of both the harmful phytoplankton are larger than the value observed when considering in the case of any one of toxic phytoplankton, while that of zooplankton is smaller than the case of any one of toxic phytoplankton is present, which imply that a considerable decrease in the biomass accumulations of plankton as compared with the case of in the presence of a single toxic phytoplankton. Therefore, the introduction of two harmful phytoplankton is contribute to the persistence of system (2) and plays an important role in the termination of planktonic blooms.

4.   Conclusions

It is now well recognized that stochastic population dynamics play a significant role in population dynamics, since environmental fluctuations can affect the growth process of species, such as the growth rate and death rate, which can be described by white noise ^[51]. And the Gaussian white noise can been theoretically preferred to model environmental fluctuations because of its irregularity and thus a good approximation to the phenomena of rapid fluctuations ^[52]. The study of stochastic population dynamics goes back to the pioneering work by Haminskii ^[53], who introduced two white noise to stabilize a linear system. After that, lots of attention has been paid on stochastic population dynamics studies ^[54,55,56]. Mao et al. ^[54] pointed out that stochastic noise can suppress potential explosion in population dynamics. Wang et al. ^[55] showed that time-periodic forcing can lead to the transitions from a spatially homogeneous stationary state to a periodic oscillation in time. Additionally, lots of stochastic plankton growth systems have been derived by numerous researchers ^[40,42,57], and stochastic plankton systems involve toxin-producing phytoplankton have become a hot topic in ecological studies due to harmful phytoplankton can significantly affect the dynamics of plankton systems ^{[18,19,20,21,22,23,24,25,34,41]}.

In this paper, therefore, we first propose a stochastic phytoplankton-zooplankton system with two harmful phytoplankton populations, where the intrinsic growth rates of two harmful phytoplankton and the natural death rate of zooplankton are influenced by the environmental noise, and we then study the effects of TPP and white noise on the dynamics of system (2) theoretically and numerically. In order to ensure that the system is biologically meaningful, the existence and uniqueness of global positive solutions of system is discussed, and the results demonstrate that for any initial value $(P_{1}(0), P_{2}(0), Z(0))\in \mathbb{R}_{+}^{3}$, the solution will remain in $\mathbb{R}_{+}^{3}$ with probability one. Based on this situation, we derive some sufficient conditions for the extinction and persistence in the mean of the system. Obviously, those conditions are great significance to study the extinction and persistence in the mean for the phytoplankton-zooplankton system ^[30,34]. Significantly, when the system is persistence in the mean, we also investigate the existence and uniqueness of positive recurrent of solution for the system, which implies the system has a unique stationary distribution under some conditions. Numerical analysis illustrates our theoretical results and further indicates that both two TPP and environmental fluctuations have a significant effect on the controlling of planktonic blooms.

On the one hand, from Theorems 2.3 and 2.4, which follows that, when the low level intensity of white noise satisfies the conditions $\theta_{1} > 0, \theta_{2} > 0,$ and $\Pi > 0$, system (2) is persistence in the mean and exists a ergodic stationary distribution, indicating the coexistence of all those three species in the random sense for a long time (see Figure 2). However, when we increase the intensity of environmental forcing such that the condition of $\theta_{2} < 0$ holds, the species $P_{2}(t)$ undergoes rapid extinction while other two species are persistence in the mean, as it is shown in Figure 3. Comparing Figures 2 and 3, it can be asserted that white noise can aggravate the stochastic oscillation of plankton density and significantly change the dynamics of phytoplankton-zooplankton system. Especially, a high intensity of white noise can accelerate the extinction of plankton populations. Consequently, these results may be more realistic than that of in ^[25], which implies that the controlling of rapid environmental fluctuations may be a good way in the termination of planktonic blooms. Therefore, it is great ecological significance to consider environmental noise when studying phytoplankton-zooplankton interaction in the presence of harmful phytoplankton.

On the other hand, it is investigated how the dynamics of system (2) strongly depends on TPP. By controlling one toxin liberation rate, the dynamic behaviours of system (2) can be changed. That is, when the toxin liberation rate is beyond some a critical value, two harmful phytoplankton can coexist, while zooplankton tends to extinction (see Figure 4). Moreover, when controlling any one of the two TPP, it is obvious to survey from Figure 5 that the increasing value of one toxin liberation rate can reduce the biomass of zooplankton, while increase the survival chance of two phytoplankton. In addition, in the presence of both two TPP, it can be seen from Figures 5 and 8 that the combined effects of two liberation rates on the changes in plankton populations are stronger than that of controlling any one of the two TPP. Thus, the introduction of two harmful phytoplankton is conducive to the persistence of the system (2) through the termination of planktonic blooms. Therefore, TPP has a profound impact on the dynamics of phytoplankton-zooplankton systems and may be used as a biological way to control planktonic blooms.

There are some interesting topics waiting for us to further explore. For example, the zooplankton mortality will occur after some time lapse due to the bloom of toxic phytoplankton ^[58], it seems to more reasonable to study a stochastic toxic-producing phytoplankton-zooplankton system with time delay. Another problem of interest is to consider impulsive perturbations or regime switching into the system. We leave those for our future research goals.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant No.61871293), the National Key Research and Development Program of China (Grant No. 2018YFE0103700) and the National Natural Science Foundation of China (Grant No. 61901303).

Conflict of interest

The author declares that there are no conflicts of interest.