Dynamical analysis of a stochastic toxin-producing phytoplankton–fish system with harvesting

In this paper, we analytically and numerically study the dynamics of a stochastic toxin-producing phytoplankton–fish system with harvesting. Mathematically, we give the existence and stability of the positive equilibrium in the deterministic system (i.e., the system without environmental noise fluctuations). In the case of the stochastic system (i.e., the system with environmental noise fluctuations), in addition to the existence and uniqueness of the positive solution, we provide the properties of the stochastic dynamics including the stochastic extinction and persistence in the mean, almost sure permanence and uniform boundedness, and the existence of ergodic stationary distribution for the phytoplankton and fish. Ecologically, via numerical analysis, we find that (1) the small random environmental fluctuations can ensure the persistence of phytoplankton and fish, but the larger one can result in the extinction of these populations; (2) an appropriate increase in harvest rate can reduce the irregular random variation of phytoplankton and fish; (3) the increase of toxin liberate rate is capable to decrease the height of probability density function of phytoplankton. These results may help us to better understand the phytoplankton–fish dynamics.


Introduction
In recent years, phytoplankton blooms, especially harmful algal blooms, occur frequently in lakes, reservoir, wetland and marine, such as Erie Lake [1], Taihu Lake and Asia-Pacific region [2]. Phytoplankton blooms can threaten drinking water, human health, and even the ecological balance [3]. For example, when harmful algal blooms emerge, a large amount of toxic substances will be released into water body, and then some fishes or shellfish may die. At present, the unpredictability of phytoplankton blooms still remains, so the control and management of phytoplankton blooms are mainly algae removal. However, there is a time lag between the occurrence of phytoplankton blooms and algae removal, so phytoplankton blooms are still the threat to the public health and ecosystem. Hence, it is very important to understand the dynamic mechanisms of phytoplankton growth, by which we can predict the evolution of phytoplankton in time and space.
Mathematical models, as powerful tools to give an important insight into understanding of the dynamic mechanisms of phytoplankton blooms, have drawn increasing attention in recent years [4][5][6][7][8]. By introducing functional responses, we can use different modifications of the classical prey-predator models to model nutrient-phytoplankton-herbivore dynamics for various aquatic environments. Due to the eutrophication, nutrient is abundant enough to support uptake of phytoplankton in water body, so nutrient may not be the main factor limiting phytoplankton growth under this case. Then the phytoplanktonherbivore model becomes the focus to study dynamics of phytoplankton growth. In nature, many phytoplankton population (e.g. Cyanobacteria) have the capacity to produce toxic substances, and then the toxin released by phytoplankton not only affects themselves but may also restrains the growth of their herbivore. The work from Chattopadhayay et al. [9] shows that toxin-producing phytoplankton may act as a biological control for the termination of phytoplankton blooms based on field observe and mathematical modelling. The results by Banerjee et al. [10] indicate that toxin released by phytoplankton may affect the growth of zooplankton, but it does not drive the zooplankton to go extinct. Currently, many researchers have been taking interest in the ecotoxicological effects of toxicant released by phytoplankton population, one may refer to [10][11][12][13][14][15][16][17][18][19][20].
In the present paper, we consider a phytoplankton-fish model where the growth of fish population is affected by some toxicant released by phytoplankton. Additionally, the nonselective harvesting of both populations is taken into consideration. Recently, many research results related to how harvesting efforts affect plankton-fish dynamics or predatorprey dynamics have been reported [21][22][23][24][25][26]. Panja et al. [22] indicate that the harvesting effort plays an important role in the stability of the phytoplankton-zooplankton-fish system. Upadhyay et al. [23] summarized that a variation of harvesting term can generate rich dynamics such as limit cycles and chaos in plankton system. Then the following dynamical system proposed as a simple model of phytoplankton-fish interaction: where P(t) and F(t) are the densities of toxin-producing phytoplankton and herbivorous fish at time t, respectively; r is the intrinsic growth rate of phytoplankton; the terms d 1 P 2 and d 2 F 2 represent intraspecific competition of phytoplankton and fish, respectively, where d 1 and d 2 are the corresponding intraspecific competition coefficient; the term αPF model phytoplankton-fish interaction, which implies that the density (or biomass) of phytoplankton consumed by fish per unit time is given by αP; β denotes the rate of biomass conversion (0 < β < α ≤ 1); μ 1 is the death rate of fish; the term ρPF γ +P [9] represents the distribution of toxic substances released by phytoplankton that ultimately contributes to the death of fish population, where γ is a half-saturation constant and ρ is the toxin liberation rate. Also here both the phytoplankton and the fish are subjected to a combined harvesting effort E, and q 1 is the catchability coefficient of the phytoplankton, and q 2 is the catchability coefficient of the fish and E is the harvesting effort. Let h 1 = Eq 1 and h 2 = Eq 2 . All the parameters are assumed to be positive.
Additionally, as most aquatic ecosystems are exposed within the open environment, the randomly fluctuating environmental factors affecting phytoplankton growth may exist, such as nutrition loading, light availability, water temperature variation [20]. Especially, Gard [27] has indicated that these factors can be described by the process of white noise. Actually, May has also pointed out the fact [28,29] that many parameters involved in the system are not constant, but exhibit random fluctuation to a greater or lesser extent because of the effects of white noise, such as the birth rates, carrying capacity, competition coefficient, and so on. Recently, some authors have taken white noise perturbations into account in approaches on the dynamics of plankton systems [30][31][32][33][34]. Mandal et al. [31] indicated that the environmental white noise perturbations play an important role in maintaining the coexistence of phytoplankton populations. Valenti et al. [32] suggested that the increase in the strength of random fluctuations is capable to cause the rapid extinction of picophytoplankton populations. Hence, under the randomly fluctuating environment, it is reasonable that the deterministic model to describe a phytoplankton-fish system is extended to a stochastic version by introducing white noise. The application of stochastic mathematical models in other fields, such as infectious disease dynamics, can be found in the literature [35][36][37][38][39].
Motivated by all this, in this paper we introduce the white noise into system (1) using the method in [40], that is, it is assumed that the environmental noise is proportional to these variables. For small t, it is appreciate to model X(t) = (P(t), F(t)) T as a Markov process with the following specifications [41]: and Then system (1) can be rewritten as the following form: where B j (t) are mutually independent standard Brownian motions defined on the probability space (Ω, F, {F } t≥0 , P) with a filtration {F } t≥0 satisfying the usual conditions (i.e., it is increasing and right continuous with F 0 contains all P-null sets) and δ 2 j denotes the intensity of the white noise, j = 1, 2.
The rest of this paper is organized as follows: In Sect. 2, we analyze the existence and stability of positive equilibrium in system (1). And then we state the existence of the positive solutions with respect to system (2) in Sect. 3. Section 4 is devoted to studying the existence of a unique ergodic stationary distribution of system (2). In Sect. 5, based on geometric structure of invariant set, we investigate the almost sure permanence and uniform boundedness of the system. In Sect. 6, the sufficient conditions, which guarantee stochastic extinction and persistence in the mean of each population, are explored. In Sect. 7, some numerical simulations are carried out to analyze dynamics of system (2) in depth. Finally, the paper ends with conclusions in Sect. 8.

Existence and stability of positive equilibrium in system (1)
In this section, it is stated first that the first quadrant is positive invariant in system (1), and then the following lemma holds.
Proof For system (1), we have P(t) > 0 and F(t) > 0 for all t ∈ [0, τ ], where τ is any positive real number. Suppose that is not true. Then there exists 0 < t τ < τ such that, for all t ∈ [0, t τ ], P(t) > 0, F(t) > 0 and either P(t τ ) = 0 or F(t τ ) = 0. Define Then, from system (1), we have P(t) = P(0) exp( This completes the proof. From system (1), it is not difficult to find that there is always a trivial equilibrium E 0 (0, 0) and the boundary equilibrium E 1 ( r-h 1 d 1 , 0) exists under r > h 1 . Obviously, the existence of E 1 depends on the harvesting level. When the condition, r ≤ h 1 , holds, the phytoplankton population will become extinct from the first equation of system (1), and then the fish population also will become extinct. Hence, it is always assumed that the condition, r > h 1 , holds in the rest of this paper.
Moreover, a unique positive equilibrium, E * (P * , F * ), can be obtained when r > h 1 and d 1 P * < rh 1 , where F * = (rh 1d 1 P * )/α and P * is the positive root of the following equation: From Eq. (3), it is obvious that there must be a positive root which is unique when the condition, r > h 1 , holds. Hence, F * > 0 when the condition, d 1 P * < rh 1 , holds, that is, the unique positive equilibrium E * exists. Then we can obtain the following theorem.
Theorem 1 When the condition, r > h 1 , holds, there exists a unique positive equilibrium E * (P * , F * ) in system (1) if the condition, d 1 P * < rh 1 , holds, which is locally asymptotically stable.
Proof From the above analysis, the existence and uniqueness of positive equilibrium E * is clear. Now, we will prove its locally asymptotical stability. By system (1), we can get the Jacobian matrix at E * , as follows: And then we can obtain the corresponding characteristic equation, as follows: where . From the second equation of system (1), it is not difficult to find that the condition, Hence, the positive equilibrium E * is locally asymptotically stable by Eq. (4).
This completes the proof. (1), then it is globally asymptotically stable.

Theorem 2 If the positive equilibrium E * exists in system
Proof In order to prove the globally asymptotical stability of E * , we consider the following Dulac function: and let (L 1 , Obviously, the function B ∈ C 1 (R 20 + ) and B > 0, where R 20 + is the interior of R 2 + . Furthermore, we have Due to the positivity of P(t) and F(t), it can be got that ∂(BL 1 ) + , and hence, by virtue of Dulac-Bendixon criterion, system (1) has no limit cycle in the positive plane. Because there exists a unique positive equilibrium of system (1) in this positive plane, and all the positive solutions of phytoplankton and fish tend to E * . Considering the local asymptotic stability of again (see Theorem 1), it can be concluded that the positive equilibrium is globally asymptotically stable if it exists. This completes the proof.

Existence and uniqueness of the positive solution in system (2)
In this section, we will discuss the existence and uniqueness of the positive solution in system (2), which can be presented using the following theorem.
Theorem 3 For any given initial value (P(0), F(0)) ∈ R 2 + , there exists a unique positive solution (P(t), F(t)) on t ≥ 0 in system (2), and the solution will remain in R 2 + with probability one.
Proof Since the coefficients of system (2) satisfy the local Lipschitz condition, then, for any initial values (P(0), F(0)) ∈ R 2 + , there exists a unique local solution (P(t), F(t)) on t ∈ [0, τ e ), where τ e denotes the explode time [42]. To show that this solution is global in R 2 + , we need to prove τ e = ∞ a.s. Hence, we choose a sufficiently large non-negative number 0 such that P(0) and F(0) lie within the interval [ 1 0 , 0 ]. For each integer , we can define the stopping time: In the following, we need to show that τ e = ∞ a.s. If this statement is violated, there exist two constants T > 0 and σ ∈ (0, 1) such that P{τ ∞ ≤ T} > σ . Hence, we can find an integer Considering that z + 1lnz ≥ 0 for all z > 0, the function V (·) is positive definite for all (P, F) ∈ R 2 + . Calculating the differential of V along the solution trajectories of system (2) using Itô's formula, we can get where LV (P, F) where K is a positive constant. The remainder of this proof is similar to reasoning in another study (e.g., see [43]) and hence is omitted here. Therefore, the solution with the initial value (P(0), F(0)) ∈ R 2 + is positive and will remain in R 2 + with probability one. This completes the proof.

Existence of ergodic stationary distribution for system (2)
In this section, we mainly study the existence of a unique ergodic stationary distribution for system (2).
Let X(t) be a regular temporally homogeneous Markov process in E l ⊂ R l described by the stochastic differential equation and the diffusion matrix is defined as follows:

Lemma 2 ([44])
Assume that there exists a bounded domain U ⊂ E l with regular boundary, which has the following properties: (a) In the domain U and some neighborhood thereof, the smallest eigenvalue of the diffusion matrix A(x) is bounded away from zero. (b) If x ∈ E l \ U, the mean time τ at which a path issuing from x reaches the set U is finite, and sup x∈Q E x τ < ∞ for every compact subset Q ⊂ E l . If assumptions (a) and (b) hold, the Markov process X(t) has a unique stationary distribution π(·) and it has ergodic property.
Then we have the following result.
there exists a unique stationary distribution π(·) for system (2) with any initial value (P 0 , F 0 ) ∈ R 2 + , and it is ergodic.
Proof In the deterministic system (1), there exists a positive equilibrium E * (P * , F * ) and therefore Set V : where θ is positive constant will be determined in later. Taking the derivative of (8) and making use of (7), we have Then we have We define c 1 = 0.5δ 2 1 P * + 0.5δ 2 2 θ F * , and if c 1 ≤ min{d 1 (P * ) 2 , d 2 θ (F * ) 2 }, then the ellipsoid lies entirely in R 2 + . One can take U as any neighborhood of the ellipsoid such thatÛ ⊂ R 2 + , whereÛ is a closure of U. Thus, we have LV < 0 for any (P, F) ∈ R 2 + \ U, which implies that the condition (b) in Lemma (2) holds. On the other hand, there is a constant M 0 > 0 such that for all (P, F) ∈Û, ξ ∈ R 2 ,which indicates that the condition (a) in Lemma (2) is also satisfied. Hence, system (2) has a unique stationary distribution π(·) and it is ergodic. This completes the proof.
Remark 1 Making use of Itô's formula on (8), we can get Integrating Eq. (10) from 0 to t on both sides, and then taking the expectation gives Thus, let = min{d 1 , αd 2 β } and then we can derive lim sup That is, the positive solution of the stochastic system (2) oscillates around the positive equilibrium E * in the deterministic system (1). Moreover, the small noise intensity can derive the solution of the stochastic system (2) to approach the positive equilibrium E * in the deterministic system (1).

Almost sure permanence and uniform boundedness of system (2)
Theorem 4 signifies that there exists a unique stationary distribution of the solutions in system (2). In this section, we shall show that the stationary distribution of the solutions lie in the interior of the first quadrant by using the geometric structure of invariant set, which implies that the solutions of the system maintain almost sure permanence and uniform boundedness.
Based on the proof of Theorem 4, we have By introducing the following notations: Θ := 0.5δ 2 1 P * + 0.5δ 2 2 F * θ , we get The derivation goes through two steps.
Step 1: Solution (P(t), F(t)) to the stochastic system (2) starts from the exterior of ξ 1 (P, F) will almost surely enter the interior of ξ 1 (P, F) in finite time.

Lemma 3 For any initial value
Proof For any ε ∈ (0, Θ), we will show a similar conclusion for ξ ε 1 (x, y) at first, i.e.
P inf Then this conclusion holds because of the arbitrariness of ε. Now, we are position to show Eq. (11). Define the first arrive time at ξ ε 1 as τ a 1 , i.e.
Step 2: Solution (P(t), F(t)) to the stochastic system (2) starts from the interior of B n 2 (P, F) will almost surely lie in the interior of B n 2 (P, F) in any finite time.
Then Eqs. (18), (19), and (20) imply LV P(s), F(s) ds P(χ) Let t → ∞, this leads to a contradiction, since V (P(η e ), F(η e )) = V (P(η a ), F(η a )) = n 2 for any χ ∈Ω. This completes the proof From Lemma 3 and Lemma 4 and based on the fact that B n 2 (P, F) lies in the interior of the first quadrant we conclude to the following theorem.
Theorem 5 Under conditions of Theorem 4, for any initial value (P 0 , F 0 ) ∈ R 2 + , system (2) is almost sure permanent, that is, there two positive constants m and M such that for any 0 ≤ t ≤ ∞. This also indicates that system (2) is uniformly ultimately bounded almost surely.

Stochastic extinction and persistence in the mean for system (2)
This section is devoted to investigating the stochastic extinction and persistence in the mean of each population for system (2).
Before presenting the main results, we first give the following notations, lemma and definition: If there exist three positive constants T, τ and τ 0 such that, for all t > T, Now we are in the position to give our main result in this section.
Remark 2 From (1) of Theorem 6, one can find that ϕ 1 is a threshold that determines whether phytoplankton will die out or survive in the future. From (2) of Theorem 6, it is worth noting that if the phytoplankton population is extinct, the fish population that live on the phytoplankton will eventually go to extinction with probability one.

Numerical simulations
In this part, we numerically study the resulting dynamics of system (1) and system (2) obtained in previous sections, in order to further explore the phytoplankton-fish dynamics.
As an example, we choose such set of parameter values as follows: α = 0.8, β = 0.6, γ = 4, ρ = 0.5, d 2 = 0.8, d 1 = 0.2, μ 1 = 0.8, r = 1.5, h 1 = 0.3 and h 2 = 0.8, where some of the parameters are derived from the literature [47]. These parameters are used for the following numerical simulations. We first analyze the dynamics of system (1) in the absence of noise perturbations. For the parameters selected above, based on Theorem 1, dy direct calculations, it is not difficult to derive E * = (3.8, 0.55). Taking ρ and h 1 as control parameters, the positive equilibrium point E * in system (1) with respect to parameters ρ and h 1 are displayed in Fig. 1. From Fig. 1(a), it is evident that the parameter ρ can promote the growth of phytoplankton population. The growth of fish population, by contrast, is restrained by the parameter ρ (see Fig. 1(b)). This result indicates that toxin released by phytoplankton is advantage for the growth of phytoplankton but is disadvantage for the growth of fish. For parameter h 1 , Figs. 1(c) and 1(d) indicate that both the growth of phytoplankton and fish are restrained by parameter h 1 . Obviously, the growth of phytoplankton and fish in system (1) can be significantly influenced by the toxin liberation rate and harvesting level. From Theorem 2, we see that all the trajectories of P(t) and F(t) in system (1) converge to the positive equilibrium E * = (3.8, 0.55) as time goes on, which implies that it is globally asymptotically stable, as is shown in Fig. 2. Next, we focus on the effects of white noise on the dynamics of system (2) using the method in [48]. Taking δ 1 and δ 2 as control parameters, the param- eter plane δ 1δ 2 can be divided into three parts according to Theorem 6 (see Fig. 3(a)). In region I of Fig. 3(a), the phytoplankton population is persistent, but extinction of fish population occurs. However, phytoplankton population and fish population are both persistent in the mean in region II. When noise intensity δ 1 is beyond the critical value (i.e. entering region III), both phytoplankton population and fish population will be extinct Figure 4 The existence of stationary distribution for system (2). (a) Bifurcation diagram with respect to δ 1 and δ 2 , where I indicates that there is no stationary distribution for system (2) and II denotes that there is a stationary distribution for system (2). (b) The histogram of probability density function of phytoplankton and fish for system (2), where the gray smoothed curve signifies the probability density function of phytoplankton or fish with probability one. More specifically, we take (δ 1 , δ 2 ) = (0.1, 1.4), (0.1, 0.2), (1.3, 1.4), then the corresponding numerical solutions of system (2) are given in Fig. 3(b).
When we analyze the existence of a stationary distribution, it is shown from Fig. 4(a) that the parameter plane δ 1δ 2 can be divided into two regions, where the stationary distribution does not exist in region I but exists in region II. For further investigation, in the persistent case of both phytoplankton and fish populations, we further adopt (δ 1 , δ 2 ) = (0.1, 0.2), (0.3, 0.2), (0.5, 0.2) and (0.7, 0.2). By Theorem 4, the existence conditions for the stationary distribution are satisfied. Consequently, the stationary distribution of phytoplankton P(t) and fish F(t) are obtained from 10,000 numerical simulations at time 200. The corresponding distributions are presented in Fig. 4(b). From (i) to (iv), it can be seen that the significant change in the stationary distribution with the increasing magnitude of noise intensities, which indicates that the mean value and skewness of the distribution for P(t) and F(t) are varying as the noise intensities increase. More precisely, for the phytoplankton population, when (δ 1 , δ 2 ) = (0.1, 0.2), the distribution appears closer to the normal distribution (see (i) in Fig. 4(b)), but as noise intensities increase to (δ 1 , δ 2 ) = (0.7, 0.2), the distribution is positively skewed (see (ii) in Fig. 4(b)).
For the case (δ 1 , δ 2 ) = (0.1, 0.2) in Fig. 4, we only vary ρ = 0.2, 0.4 and 0.6, it is found that the height of the probability density function of phytoplankton in system (2) is decreasing as the toxin liberation rate ρ increases (see Fig. 5). Furthermore, when the parameter h 1 varies, it is observed that the biomass of phytoplankton and fish in system (2) decreases with the increase of h 1 , which are clearly shown in Fig. 6(a) and Fig. 6(b), respectively; if the value of h 1 continues to increase, both phytoplankton and fish populations will eventually go to extinction, as are demonstrated in Fig. 6(c) and Fig. 6(d), respectively. Comparing Fig. 6(a) and Fig. 6(c) (or Fig. 6(b) and Fig. 6(d)), it can be seen that the irregularity of random variation and the range of the fluctuations of P(t) (or F(t)) decrease as harvest rate h 1 increases properly (see red and green curves). This indicates that the appropriate harvesting strategy may be beneficial to the persistence of phytoplankton and fish populations. Comparing Fig. 6(c) and Fig. 6(d), one can find that with the increase of h 1 , the phytoplankton P(t) is persistent but fish F(t) goes extinct (see the green curve in Fig. 6(c) and Fig. 6(d)), which suggests that if the phytoplankton on which herbivorous fish lives is harvested in large quantities, the fish population will not have enough food to survive.
In order to explore the relationships between the positive solution of system (2) and the positive equilibrium E * of system (1), we choose two sets of noise intensity values (δ 1 , δ 2 ) = (0.01, 0.02) and (δ 1 , δ 2 ) = (0.2, 0.4), respectively, and then the corresponding results are displayed by Fig. 7, which show that the positive solutions of the stochastic sys- tem (2) are oscillating around the steady state E * = (3.8, 0.55) of system (1) (see Fig. 7). Furthermore, it can be seen that the irregularity of stochastic variation of the solutions in Fig. 7(b) is stronger than that in Fig. 7(a), and the range of fluctuation of the solutions is larger than that in Fig. 7(a). By comparison, it can be asserted that the smaller the random environmental fluctuation is, the closer the solutions are to the steady state E * , which agrees well with Remark 1.

Conclusions
In this paper we proposed a phytoplankton-fish model under random environmental fluctuations (i.e. system (2)), where the effects of harvesting and toxin-producing phytoplankton are both considered. In order to discuss the effect of random environmental fluctuations on phytoplankton-fish dynamics, the phytoplankton-fish model without stochastic effect was presented firstly (i.e. system (1)), and it was found by analysis that the positive equilibrium is globally asymptotically stable if it exists in system (1), which implies that phytoplankton population and fish population can coexist at the positive equilibrium eventually. In system (1), when the positive equilibrium does not exist, fish population is always extinct, but the permanence of phytoplankton population depends upon the harvesting level.
However, when we take random environmental fluctuations into consideration, population extinction can occur from stochastic dynamic viewpoint even if the positive equilibrium exists in system (1). By constructing stochastic Lyapunov function, it has been derived that the system (2) has a unique stationary distribution, and the unqiue stationary distribution lies in the interior of the first quadrant using the geometric structure of invariant set, which means that the solutions of the system remain the almost sure permanence and uniform boundedness. Furthermore, we obtained the sufficient conditions guaranteeing the stochastic extinction and persistence in the mean of each population.
Additionally, numerical analysis indicated that the small random environmental fluctuations can ensure the existence of a unique stationary distribution denoting the persistence of phytoplankton and fish in system (2), while the large random environmental fluctuations can lead to the extinction of the two populations, even if the unique positive equilibrium was globally asymptotically stable in system (1). Furthermore, numerical simulations shows that the increase of harvest rate can reduce the irregularity of random variation and the range of fluctuation of the phytoplankton and fish, and the larger harvest rate could lead these populations to go extinct eventually (see Fig. 6). Consequently, an appropriate harvesting strategy may be beneficial to the balance of ecosystems.
Moreover, it is worth noting that system (2) almost preserves the property of the global stability when the white noise intensities are much smaller. In this case, we can ignore the effects of white noise, and phytoplankton-fish dynamic can be presented approximately using of system (1). However, the larger intensities of white noise can force the solution of system (2) to oscillate strongly around the positive steady state E * , even the extinction of population can emerge. Therefore, in these cases, we cannot ignore the effects of random environmental fluctuations.
The numerical simulations are agreement with theoretical results very well, which indicates that the mathematical models proposed in this paper are feasible to study phytoplankton-fish dynamics. Nevertheless, some other factors that can affect the growth of phytoplankton, such as cell size [49], deserve further investigation. These will be left for our future work.