Stationary Probability Density Analysis for the Randomly Forced Phytoplankton–Zooplankton Model with Correlated Colored Noises

Abstract

In this paper, we propose a stochastic phytoplankton–zooplankton model driven by correlated colored noises, which contains both anthropogenic and natural toxins. Using Khasminskii transformation and the stochastic averaging method, we first transform the original system into an Itô diffusion system. Afterwards, we derive the stationary probability density of the averaging amplitude equation by utilizing the corresponding Fokker–Planck–Kolmogorov equation. Then, the stability of the averaging amplitude is studied and the joint probability density of the original two-dimensional system is given. Finally, the theoretical results are verified by numerical simulations, and the effects of noise characteristics and toxins on system dynamics are further illustrated.

1. Introduction

Phytoplankton are tiny floating plants, usually single-celled algae, which sit at the bottom of the food chain and thus support secondary and tertiary productivity in the ocean [1,2,3]. As the major primary producers in the marine ecosystem, phytoplankton are not only a basic food source for zooplankton, but also provide large amounts of oxygen for humans and other animals after absorbing nearly half of the universal carbon dioxide through photosynthesis [4,5,6]. However, unfortunately, the rapid and sustained growth of phytoplankton biomass can alter energy flows and disrupt the normal functioning of marine ecosystem [7]. Especially in the stage of algal bloom demise, the decomposition process after the death of high concentration algae can absorb a large amount of dissolved oxygen in the water, resulting in the death of marine organisms asphyxiation, which in turn affects the water quality, tourism, and fishery resources [8,9].

There are many factors that contribute to algal blooms, such as nutrient level [10], temperature [11] and light availability [12], but the key cause for the formation of algal blooms is still not fully understood. Perhaps the main reason behind population succession and algal blooms is the toxins produced by toxin-producing phytoplankton (TPP). For this reason, Chattopadhayay et al. [13] proposed a general form of phytoplankton–zooplankton interaction model with TPP, and concluded that TPP may act as a biological control for planktonic blooms through the field-collected samples and theoretical analysis. Since then, a large number of plankton models with TPP have been developed, with similar results—including, but not limited to, Refs. [14,15,16,17,18,19]. Notably, most studies only consider the impact of toxic substances released by TPP on the grazing pressure of zooplankton, while ignoring the impact of anthropogenic toxins in the environment. In fact, chemical toxins widely exist in earth’s ecosystem as a result of multitudes of human activities, seriously affecting the survival of plankton and human health [20]. Assuming that prey and predator have different rates of exposure to anthropogenic toxins, Das et al. [21] investigated the bioeconomic harvesting of a prey–predator model in the presence of toxicity. Later, Chakraborty and Das [22] extended this model to three-dimensional case, which includes two-zooplankton and one-phytoplankton. To study the direct and indirect effects of anthropogenic toxins on competitive species, Shan and Huang [23] established a toxin-dependent competition model and found that the level of toxins and the distinct vulnerabilities of two species to toxins can influence competition outcomes in many counterintuitive ways.

Note that the aforementioned models are based on a deterministic approach. Species concentrations, however, are susceptible to fluctuations in the real environment that can not be neglected when seeking a better understanding of the dynamics of complex living systems [24,25,26,27,28,29,30,31,32]. This is especially true for marine organisms due to the unpredictability of water temperature, nutrient availability, photosynthetically active radiation and other physical factors embedded in aquatic ecosystems [33,34,35,36]. Consequently, nonlinear stochastic ecosystem with noise term has recently attracted the attention of scholars [37,38,39,40,41]. For example, authors of [42,43] established different Itô type stochastic plankton models and studied their asymptotic behavior and stability by Lyapunov function method. Significantly, Gauss white noise in the sense of Itô is regarded as ideal white noise with infinite power, which does not exist in reality. As a result, researchers [44] used Stratonovich type stochastic differential equation instead of Itô type to characterize the interaction between plankton and marine environment. In particular, Spagnolo and his co-workers [45,46,47] developed different types of stochastic population models with multiplicative white Gaussian noise to simulate the complex dynamic behavior of ecosystems and observed that theoretical results can effectively reproduce experimental data. Huang et al. [48] explored bifurcation dynamics of two competing algal species by constructing a stochastic nonlinear model with multiplicative and additive noises. Except for biomathematics, multiplicative and additive noises are also widely used in other many fields, such as medicine [49], physics [50], neurology [51], and signal processing [52]. It is important to note that the above stochastic models were established based on the assumption that noises have different sources, i.e., they are independent of each other. Meanwhile, the researchers [53] uncovered that in some cases noises may have a common source and therefore can be correlated. Although the correlation time of actual noises may be small physically, it can never be strictly zero [54]. Especially, it is important to consider the non-zero correlation time when the fluctuation is large and the driving process relaxation time scale is long [55]. Noise is generally considered harmful and can cause disturbance of the dynamical system. However, in recent years, the constructive and counterintuitive role of noise in nonlinear system, such as noise-enhanced stability [56], noise-induced resonance [57], noise-induced transport [58], etc., has been widely studied theoretically and experimentally. For more information, please refer to Refs. [59,60,61,62].

Inspired by the above discussion, it is more reasonable to incorporate natural and anthropogenic toxins and environmental noises into the plankton model. To this end, a randomly forced phytoplankton–zooplankton model driven by correlated colored noises, which contains both anthropogenic and natural toxins, is constructed in Section 2. In Section 3, to explore the stochastic dynamic properties of the model, we first transform the original model into an averaging Itô diffusion system by using the Khasminskii transformation [63] and stochastic averaging method [64,65]. Then, the stationary probability density of the averaging amplitude equation by utilizing the corresponding Fokker–Planck–Kolmogorov equation is derived. Afterwards, the stability of the averaging amplitude is studied and the joint probability density of the proposed two-dimensional system is given. In Section 4, some numerical simulations are carried out to illustrate the obtained theoretical results and the effects of noise characteristics and toxins on system dynamics are further illustrated. Finally, we conclude the paper by a brief discussion in Section 5.

2. Model Formulation

Let $P\equiv P\left(t\right)$ and $Z\equiv Z\left(t\right)$, respectively, denote the TPP and zooplankton populations at time t. Now some assumptions for establishing the model are presented.

(i)

In the absence of zooplankton, the growth of TPP population follows the logistic law with intrinsic growth rate r and environmental carrying capacity K.

(ii)

In the absence of limiting factors, the chance of an individual zooplankton encountering prey is proportional to its abundance, so the predation rate is assumed to obey the simple law of mass action. On the other hand, no matter how large the TPP population is, each zooplankton individual has a maximum consumption rate. Therefore, in the presence of toxic algae, the more common choice is the Holling type II functional response to describe this grazing phenomena [13].

(iii)

TPP population is directly affected by anthropogenic toxins, while zooplankton population feeding on contaminated TPP is indirectly affected by toxins [21].

(iv)

Some external factors such as temperature and climate affect the biological individual, which result in a multiplicative noise $\xi \left(t\right)$[52], while internal factors such as competition between individuals for food and living environment may alter the population directly, leading to an additive noise $\eta \left(t\right)$ [53].

Based on the above assumptions, the stochastic phytoplankton–zooplankton model with TPP is given as follows:

$$\left\{\begin{array}{c}\frac{dP}{dt}=rP\left(1-\frac{P}{K}\right)-{\alpha }_{1}PZ-{\gamma }_1{P}^3+P\xi \left(t\right)+\eta \left(t\right),\hfill \\ \frac{dZ}{dt}={\alpha }_{2}PZ-\mu Z-\frac{bPZ}{a+P}-{\gamma }_2{Z}^2+Z\xi \left(t\right)+\eta \left(t\right).\hfill \end{array}\right.$$

(1)

Worth noting that phytoplankton and zooplankton have different exposure rates to anthropogenic toxins released by external environment due to their density and size. Hence, the effect of anthropogenic toxins on zooplankton population is less than that on TPP population by assumption (iii), namely, the coefficients of toxicity satisfy $0<{\gamma }_2<{\gamma }_1$ [22]. In addition, both phytoplankton and zooplankton are assumed to be affected by the same multiplicative noise source $\xi \left(t\right)$ and the same additive noise source $\eta \left(t\right)$, where $\xi \left(t\right)$ and $\eta \left(t\right)$ are correlated Gaussian colored noises with zero mean and satisfy

$$\begin{array}{cc}\hfill & \langle \xi \left(t\right)\xi \left(s\right)\rangle =\frac{D_1}{{\tau }_1}exp(-\frac{|t-s|}{{\tau }_1}),\hfill \\ \hfill & \langle \eta \left(t\right)\eta \left(s\right)\rangle =\frac{D_2}{{\tau }_2}exp(-\frac{|t-s|}{{\tau }_2}),\hfill \\ \hfill & \langle \xi \left(t\right)\eta \left(s\right)\rangle =\frac{D_3}{{\tau }_3}exp(-\frac{|t-s|}{{\tau }_3}).\hfill \end{array}$$

(2)

Here, ${\tau }_i$ ($i=1,2,3$) denote the correlation time, $D_1$ and $D_2$ are the intensity of the two colored noises, and $D_3=\lambda \sqrt{D_1{D}_2}$, where $\lambda$ represents the correlation degree between the two noises with $\left|\lambda \right|<1$. Biological explanations of other parameters are shown in

Table 1. Biological explanation of parameters in model (1).

Parameter	Description
r	Intrinsic growth rate of TPP population
K	Environmental carrying capacity of TPP population
${\alpha }_1$	Rate of predation of zooplankton on TPP population
${\alpha }_2$	Ratio of biomass consumed by zooplankton for its growth
$\mu$	Mortality rate of zooplankton
a	Half saturation constant
b	Rate of toxin liberation by TPP population
${\gamma }_1$	Coefficient of toxicity to phytoplankton
${\gamma }_2$	Coefficient of toxicity to zooplankton

. In the following, the stationary probability density of model (1) will be presented.

3. Stochastic Analysis of The Model

Without considering the impact of noise, model (1) degenerates into a deterministic one. Then the equilibria for the deterministic model is $E_0(0,0)$, $E_1(P_1,0)$, $E_{*}(P_{*},Z_{*})$, where

$$P_1=\frac{1}{2{\gamma }_1}(-\frac{r}{K}+\sqrt{{\left(\frac{r}{K}\right)}^2+4r{\gamma }_1}),Z_{*}=\frac{1}{{\alpha }_1}(r-\frac{r}{K}{P}_{*}-{\gamma }_1{P}_{*}^2),$$

and $P_{*}$ satisfies

$$d_0{P}_{*}^3+d_1{P}_{*}^2+d_2{P}_{*}+d_3=0,$$

(3)

here

$$\begin{array}{cc}\hfill & d_0={\gamma }_1{\gamma }_2,\hfill \\ \hfill & d_1={\alpha }_1{\alpha }_2+{\gamma }_2(a{\gamma }_1+r{K}^{-1}),\hfill \\ \hfill & d_2=ar{\gamma }_2{K}^{-1}-r{\gamma }_2+a{\alpha }_1{\alpha }_2-\mu {\alpha }_1-b{\alpha }_1,\hfill \\ \hfill & d_3=-a\mu {\alpha }_1-ar{\gamma }_2.\hfill \end{array}$$

It is obvious from (3) that the necessary and sufficient condition for the existence of positive equilibrium $E_{*}$ is

$$r-r{P}_{*}{K}^{-1}-{\gamma }_1{P}_{*}^2>0.$$

(4)

Next, the stationary probability density and stability of model (1) are discussed in detail. We first introduce a transformation of variables, $X=P-P_{*}$, $Y=Z-Z_{*}$, and substitute this into model (1), which will lead to the following form by ignoring the higher powers:

$$\left\{\begin{array}{c}\frac{dX}{dt}=a_{11}X+a_{12}Y+(X+P_{*})\xi \left(t\right)+\eta \left(t\right),\hfill \\ \frac{dY}{dt}=a_{21}X+a_{22}Y+(Y+Z_{*})\xi \left(t\right)+\eta \left(t\right),\hfill \end{array}\right.$$

(5)

where

$$\begin{array}{cc}\hfill a_{11}=& -r{K}^{-1}{P}_{*}-2{\gamma }_1{P}_{*}^2,a_{12}=-{\alpha }_1{P}_{*},\hfill \\ \hfill a_{21}=& {\alpha }_2{Z}_{*}-\frac{ab{Z}_{*}}{{(a+P_{*})}^2},a_{22}=-{\gamma }_2{Z}_{*}.\hfill \end{array}$$

(6)

It follows from the Khasminskii transformation $X=\rho cos\theta$ and $Y=\rho sin\theta$ [63] that

$$\left\{\begin{array}{c}\frac{d\rho }{dt}=f_1(\rho ,\theta )+g_{11}(\rho ,\theta )\xi \left(t\right)+g_{12}(\rho ,\theta )\eta \left(t\right),\hfill \\ \frac{d\theta }{dt}=f_2(\rho ,\theta )+g_{21}(\rho ,\theta )\xi \left(t\right)+g_{22}(\rho ,\theta )\eta \left(t\right),\hfill \end{array}\right.$$

(7)

where

$$\begin{array}{cc}\hfill & f_1(\rho ,\theta )=\rho (a_{11}{cos}^2\theta +(a_{12}+a_{21})sin\theta cos\theta +a_{22}{sin}^2\theta ),\hfill \\ \hfill & f_2(\rho ,\theta )=a_{21}{cos}^2\theta -a_{12}{sin}^2\theta +(a_{22}-a_{11})sin\theta cos\theta ,\hfill \\ \hfill & g_{11}(\rho ,\theta )=\rho +P_{*}cos\theta +Z_{*}sin\theta ,g_{12}(\rho ,\theta )=cos\theta +sin\theta ,\hfill \\ \hfill & g_{21}(\rho ,\theta )={\rho }^{-1}(Z_{*}cos\theta -P_{*}sin\theta ),g_{22}(\rho ,\theta )={\rho }^{-1}(cos\theta -sin\theta ).\hfill \end{array}$$

(8)

After applying the stochastic averaging method [64,65] for system (7), we obtain

$$\left\{\begin{array}{c}d\rho =m_{\rho }\left(\rho \right)dt+{\sigma }_{11}\left(\rho \right)d{w}_{\rho }+{\sigma }_{12}\left(\rho \right)d{w}_{\theta },\hfill \\ d\theta =m_{\theta }\left(\rho \right)dt+{\sigma }_{21}\left(\rho \right)d{w}_{\rho }+{\sigma }_{22}\left(\rho \right)d{w}_{\theta }.\hfill \end{array}\right.$$

(9)

Here, $w_{\rho }$ and $w_{\theta }$ stand for independent standard Wiener processes, and the drift and diffusion coefficients are

$$\begin{array}{cc}\hfill m_{\rho }\left(\rho \right)=& \frac{1}{2}\rho (a_{11}+a_{22})+\frac{D_1(2{\rho }^2{\tau }_1^2+2{\rho }^2+P_{*}^2+Z_{*}^2)}{2\rho ({\tau }_1^2+1)}\hfill \\ \hfill & +\frac{D_2}{\rho ({\tau }_2^2+1)}+\frac{D_3(P_{*}+Z_{*})}{\rho ({\tau }_3^2+1)},\hfill \\ \hfill m_{\theta }\left(\rho \right)=& \frac{1}{2}(a_{21}-a_{12}),{\sigma }_{12}^2\left(\rho \right)=0,{\sigma }_{21}^2\left(\rho \right)=0,\hfill \\ \hfill {\sigma }_{11}^2\left(\rho \right)=& \frac{D_1(2{\rho }^2{\tau }_1^2+2{\rho }^2+P_{*}^2+Z_{*}^2)}{{\tau }_1^2+1}+\frac{2D_2}{{\tau }_2^2+1}+\frac{2D_3(P_{*}+Z_{*})}{{\tau }_3^2+1},\hfill \\ \hfill {\sigma }_{22}\left(\rho \right)=& \frac{D_1(P_{*}^2+Z_{*}^2)}{{\rho }^2({\tau }_1^2+1)}+\frac{2D_2}{{\rho }^2({\tau }_2^2+1)}+\frac{2D_3(P_{*}+Z_{*})}{{\rho }^2({\tau }_3^2+1)}.\hfill \end{array}$$

For simplicity, we introduce some notations:

$$\begin{array}{cc}\hfill & {\phi }_1=\frac{1}{2}(a_{11}+a_{22})+D_1,{\phi }_2=2D_1,{\phi }_3=\frac{1}{2}(a_{21}-a_{12}),\hfill \\ \hfill & {\phi }_4=\frac{D_1(P_{*}^2+Z_{*}^2)}{2({\tau }_1^2+1)}+\frac{D_2}{{\tau }_2^2+1}+\frac{D_3(P_{*}+Z_{*})}{{\tau }_3^2+1}.\hfill \end{array}$$

(10)

Then the probability density function $p\left(\rho \right)$ for the diffusion process $(\rho ,\theta )$ satisfies the Fokker–Planck–Kolmogorov equation

$$\frac{\partial p}{\partial t}+\frac{\partial }{\partial \rho }\left(m_{\rho }p\right)+\frac{\partial }{\partial \theta }\left(m_{\theta }p\right)=\frac{1}{2}\left[\frac{{\partial }^2}{\partial {\rho }^2}\left({\sigma }_{11}^{2}p\right)+\frac{{\partial }^2}{\partial {\theta }^2}\left({\sigma }_{22}^{2}p\right)\right].$$

(11)

Note that the amplitude $\rho$ in (9) does not depend on the phase $\theta$, so Equation (11) degenerates into

$$\frac{\partial p}{\partial t}=-\frac{\partial }{\partial \rho }\left(m_{\rho }p\right)+\frac{1}{2}\frac{{\partial }^2}{\partial {\rho }^2}\left({\sigma }_{11}^{2}p\right),$$

(12)

i.e.,

$$\frac{\partial p}{\partial t}=-\frac{\partial }{\partial \rho }\left[({\phi }_1\rho +{\phi }_4{\rho }^{-1})p\right]+\frac{1}{2}\frac{{\partial }^2}{\partial {\rho }^2}\left[(2{\phi }_4+{\phi }_2{\rho }^2)p\right].$$

(13)

It then follows from $\frac{\partial p}{\partial t}=0$ that the stationary probability density function $p\left(\rho \right)$ satisfies

$$-\frac{\partial }{\partial \rho }\left[({\phi }_1\rho +{\phi }_4{\rho }^{-1})p\right]+\frac{1}{2}\frac{{\partial }^2}{\partial {\rho }^2}\left[(2{\phi }_4+{\phi }_2{\rho }^2)p\right]=0.$$

(14)

From Equations (6) and (10), we know that ${\phi }_2>0$ and ${\phi }_4>0$. Therefore, when ${\phi }_2>2{\phi }_1$, the following expression of stationary probability density function can obtained from (14)

$$p\left(\rho \right)=\frac{{\phi }_2}{{\phi }_4}\rho {\left(1+\frac{{\phi }_2}{2{\phi }_4}{\rho }^2\right)}^{\frac{{\phi }_1}{{\phi }_2}-\frac{3}{2}}\Gamma \left(\frac{3}{2}-\frac{{\phi }_1}{{\phi }_2}\right){\Gamma }^{-1}\left(\frac{1}{2}-\frac{{\phi }_1}{{\phi }_2}\right),$$

(15)

where $\Gamma \left(x\right)={\int }_0^{\infty }{t}^{x-1}{e}^{-t}dt$ is the Gamma function with $t>0$.

One can see from $\frac{dp\left(\rho \right)}{d\rho }=0$ that the extreme point of $p\left(\rho \right)$ is $\sqrt{\frac{{\phi }_4}{{\phi }_2-{\phi }_1}}$ (${\phi }_2>2{\phi }_1$). Further calculations show that $\frac{d^{2}p\left(\rho \right)}{d{\rho }^2}{|}_{\rho =\sqrt{\frac{{\phi }_4}{{\phi }_2-{\phi }_1}}}<0$, which implies that $\overline{\rho }=\sqrt{\frac{{\phi }_4}{{\phi }_2-{\phi }_1}}$ (${\phi }_2>2{\phi }_1$) is a maximum point of $p\left(\rho \right)$. According to Namachchivaya’s theory [66], the extreme value of an invariant measure contains the most important essence of nonlinear stochastic system, namely, the sample trajectory will stay near $\overline{\rho }$ with bigger probability, which means that $\overline{\rho }$ is stable in the sense of probability.

According to $\rho ={(X^2+Y^2)}^{\frac{1}{2}}$, we can further obtain the following joint probability density of variables X and Y:

$$\Psi (X,Y)=\frac{{\phi }_2}{{\phi }_4}{(X^2+Y^2)}^{\frac{1}{2}}{\left[1+\frac{{\phi }_2}{2{\phi }_4}(X^2+Y^2)\right]}^{\frac{{\phi }_1}{{\phi }_2}-\frac{3}{2}}\Gamma \left(\frac{3}{2}-\frac{{\phi }_1}{{\phi }_2}\right){\Gamma }^{-1}\left(\frac{1}{2}-\frac{{\phi }_1}{{\phi }_2}\right).$$

(16)

Hence, $\Psi (X,Y)$ exists a maximum value at the points of the stable limit cycle

$$X^2+Y^2=\frac{{\phi }_4}{{\phi }_2-{\phi }_1}.$$

For system (1), we can obtain similar results, namely, the probability density of variables P and Z is

$$\Psi (P,Z)=\frac{{\phi }_2}{{\phi }_4}\sqrt{\mathsf{\Theta }(P,Z)}{\left[1+\frac{{\phi }_2}{2{\phi }_4}\mathsf{\Theta }(P,Z)\right]}^{\frac{{\phi }_1}{{\phi }_2}-\frac{3}{2}}\Gamma \left(\frac{3}{2}-\frac{{\phi }_1}{{\phi }_2}\right){\Gamma }^{-1}\left(\frac{1}{2}-\frac{{\phi }_1}{{\phi }_2}\right),$$

(17)

where $\mathsf{\Theta }(P,Z)={(P-P_{*})}^2+{(Z-Z_{*})}^2$.

4. Numerical Simulation

In this section, we will perform some numerical simulations to verify the theoretical results and further illustrate the effects of noise characteristics and anthropogenic and natural toxins on system dynamics. Unless otherwise specified, all numerical simulations in this paper are carried out by fixing parameters $r=0.2$, $K=108$, ${\alpha }_1=1$, ${\alpha }_2=0.15$, $\mu =0.075$, $a=5.7$ and changing $D_i$, ${\tau }_i$, ${\gamma }_j$ and b, $i=1,2,3$, $j=1,2$.

(i) Assuming that one noise intensity is fixed while the other two change continuously, and selecting ${\tau }_1=0.5$, ${\tau }_2=0.8$, ${\tau }_3=0.1$, ${\gamma }_1=0.057$, ${\gamma }_2=0.02$ and $b=0.5$, we plot Figure 1, which shows the effect of noise intensity on the maximum amplitude $\overline{\rho }$. We find from Figure 1 that the maximum amplitude $\overline{\rho }$ decreases with the increase in multiplicative noise intensity $D_1$, but the effect of $D_2$ and $D_3$ is just the opposite.

(ii) Let $D_1=0.1$, $D_2=0.05$, $D_3=0.15$, ${\gamma }_1=0.057$, ${\gamma }_2=0.02$, $b=0.5$. Then, the effect of correlation time on the maximum amplitude $\overline{\rho }$ is shown in Figure 2. Clearly, $\overline{\rho }$ tends to decrease due to the increase in ${\tau }_i$ ($i=1,2,3$), which means that the correlation time is not conducive to the stability of the system.

(iii) Figure 3 is plotted by choosing $D_1=0.1$, $D_2=0.05$, $D_3=0.15$, ${\tau }_1=0.5$, ${\tau }_2=0.8$, ${\tau }_3=0.1$ and $b=0.5$, which illustrates the effect of anthropogenic toxins. We see from Figure 3 that $\overline{\rho }$ decreases with increasing ${\gamma }_1$, while the effect of ${\gamma }_2$ is more complex. To be specific, $\overline{\rho }$ increases only when ${\gamma }_2$ is smaller and decreases when ${\gamma }_2$ is larger, so ${\gamma }_2$ exists a maximum point. This implies that the coefficient of anthropogenic toxicity to zooplankton ${\gamma }_2$ may be used as a biological control measure for algal blooms.

(iv) Let $D_1=0.1$, $D_2=0.05$, $D_3=0.15$, ${\tau }_1=0.5$, ${\tau }_2=0.8$, ${\tau }_3=0.1$, ${\gamma }_1=0.057$, ${\gamma }_2=0.02$. Then, the effect of TPP population on the maximum amplitude $\overline{\rho }$, i.e., the effect of natural toxins, is shown in Figure 4. Interestingly, Figure 4 exhibits two extreme values, a maximum and a minimum, and $\overline{\rho }\to \infty$ when $b\to \infty$. That is to say, toxin release rate by TPP population b may affect the intensity of the algal blooms, so it plays an important role in the formation of algal blooms.

(v) Let $D_1=D_2=D_3=0.01$, ${\tau }_1={\tau }_2={\tau }_3=0.02$, ${\gamma }_1=0.057$, ${\gamma }_2=0.02$, $b=0.5$. Then, we plot the joint probability density of system (5), as shown in Figure 5. Therefore, $\Psi (X,Y)$ exists a maximum value at the points of the stable limit cycle: $X^2+Y^2=\frac{{\phi }_4}{{\phi }_2-{\phi }_1}=0.63$.

5. Discussion

Considering internal and external random factors, in this paper we proposed a stochastic phytoplankton–zooplankton model with correlated colored noises. In order to explore the stochastic dynamics of the model, we first employed the Khasminskii transformation and stochastic averaging method to transform the original model into an averaging Itô diffusion system. Then, the stationary probability density of the diffusion process was obtained by utilizing the corresponding Fokker–Planck–Kolmogorov equation. Finally, we discussed the stability of the averaging amplitude with the help of maximum value and presented the joint probability density of the original two-dimensional system. Summarizing the theoretical and numerical results, we can draw the following interesting conclusions.

•

The noise intensity $D_1$, the correlation time ${\tau }_i$ ($i=1,2,3$) and the coefficient of anthropogenic toxicity ${\gamma }_1$ may reduce the level of $\overline{\rho }$, namely, the distribution range of population density will be more concentrated with the increase in $D_1$, ${\tau }_i$ ($i=1,2,3$), and ${\gamma }_1$. Ecologically, these parameters are favorable for maintaining a balanced plankton population, which may seem counterintuitive.

•

The noise intensities $D_2$ and $D_3$ can enhance the level of $\overline{\rho }$, which implies that the distribution range of population density will be enlarged with the increase in $D_2$ and $D_3$. As a result, they weaken the stability of the system.

•

The influence of anthropogenic toxicity coefficient ${\gamma }_1$ and the toxin release rate by TPP population b is more complicated, depending on the content of the two toxins. In other words, these two parameters can be used as a means of controlling algal blooms.

This paper is concerned with stationary probability density analysis for the randomly forced phytoplankton–zooplankton model with correlated colored noises. The obtained results may enrich the research of aquatic ecosystem dynamics. Note here that the model proposed in this paper is a two-dimensional model, which is based on top-down mechanism. How to construct a high-dimensional stochastic food chain model based on both bottom-up and top-down mechanisms to describe the interaction between nutrients and plankton, and to carry out its stationary probability density analysis is a problem worth discussing. We leave this for future investigations.