Dynamical analysis of a toxin-producing phytoplankton-zooplankton model with refuge.

To study the impacts of toxin produced by phytoplankton and refuges provided for phytoplankton on phytoplankton-zooplankton interactions in lakes, we establish a simple phytoplankton-zooplankton system with Holling type II response function. The existence and stability of positive equilibria are discussed. Bifurcation analyses are given by using normal form theory which reveals reasonably the mechanisms and nonlinear dynamics of the effects of toxin and refuges, including Hopf bifurcation, Bogdanov-Takens bifurcation of co-dimension 2 and 3. Numerical simulations are carried out to intuitively support our analytical results and help to explain the observed biological behaviors. Our findings finally show that both phytoplankton refuge and toxin have a significant impact on the occurring and terminating of algal blooms in freshwater lakes.


1.
Introduction. The outbreaks of algal bloom in freshwater lakes have been becoming more and more frequent and popular on the globe. For example, Tai Lake (or Taihu), China's third-largest freshwater lake, has repeatedly suffered from the disastrous algae blooms, which dates back to at least 1987. The annual duration of algal blooms became longer and longer from 1987 to 2007 and the date of initial blooming was in approximately 11.42 days advancement per year since 1998 [1]. Recently, algae blooms frequently broke out in Lake Erie, the fourth largest lake of the five Great Lakes in north America, from 2008 to 2013 [2]. The harmful algal blooms bring a great economic loss to the nation and people. for example, in 2011, a sixth of the surface of Lake Erie was covered by thick of toxic blue-green algae, which was estimated to have brought a 2.4 million dollar loss to Ohios recreational fishing industry, and a 1.3 million dollar loss to the Maumee Bay State Park [3]. In a report to the Congress of USA [4], "harmful algal blooms were considered as one of the most scientifically complex and economically damaging issues challenging our ability to safeguard the health of our Nations aquatic and marine ecosystems." For algal bloom, there have been extensive modelling studies to understand bloom dynamics by considering the factors, including nutrient, temperature, light, viral disease, harvesting, and toxin released by some phytoplankton and so on ( [5,6,7,8,9,10,11,12] and references therein). However, the potential effect of toxin released by phytoplankton on algae bloom does not receive much attention as it should. Hallam et al. [13] observed that some plankton are able to produce and release toxic substances into lake environment which has an adverse effect on other species. Noctiluca scintillans, for example, is harmful to other planktonic organisms including fish by releasing more toxin when algal blooms occur. Turner et al. [14] pointed out that the mechanisms of toxin between toxic phytoplankton and their grazers are very complex. So far the related researches have been focused on considering the impact of toxin released by toxic-phytoplankton in phytoplanktonzooplankton models to study the dynamical behaviors of the toxic-phytoplankton. Based on both observations and mathematical models, Chattopadhyay et al. [15] concluded that toxin released by phytoplankton has a capacity to terminate the planktonic blooms by decreasing the grazing pressure of zooplankton and can serve as a biological controller. At the same time, by proposing and analyzing a simple phytoplankton-zooplankton model under three forms of the liberation of toxic substances: (a) a Holling type II, (b) a gamma distribution and (c) a discrete type, respectively, Chattopadhyay et al. [16] once again confirmed that toxin produced by phytoplankton can play a significant role in terminating the blooms. In a nutrientphytoplankton-zooplankton model, the similar finding was also observed by Pal et al. [17], they showed that the concentration of toxin that surpasses a threshold level dampens phytoplankton-zooplankton population oscillations and has a stabilizing effect on the lake ecosystem. In order to explore the dynamics of seasonally recurring bloom phenomena, Chakraborty et al. [18] proposed a simple nutrientphytoplankton model by considering the toxin liberation rate with a periodic function to reflect the seasonal changes. Their findings showed that with a changing liberation rate, toxin can contribute to the explanation for algal bloom and a wide range of complex dynamical behaviors, such as simple cyclical blooms, irregular chaotic blooms and skipping phenomenon, were also observed by varying the toxin liberation rate. By analyzing a phytoplankton-toxic phytoplankton-zooplankton model with a Monod-Haldane response function, Banerjee et al. [19] concluded that the zooplankton population can survive with the existence of toxic phytoplankton. Therefore, it is essential to take toxin into consideration when studying the interactions between phytoplankton and zooplankton in a freshwater lake.
In nature, there is another fascinating factor, the refuge provided to prey, can have an important impact upon the dynamics of the food web and the balance of prey-predator interaction. For a lake ecosystem, the concept of a refuge is worthy of attention since it might stabilize the biomass of plankton by preventing extinction of phytoplankton from the predation of zooplankton temporarily [20].
At present, the effect of refuge on some predator-prey interactions has been widely studied. Collings et al. [21] demonstrated that a refuge incorporated in mite predator-prey interactions can offer the prey some degree of protection from predation and help to prolong a predator-prey interaction by reducing the possibility of extinction due to predation, and they also showed that increasing the amount of refuge can raise prey densities and even trigger population outbreaks due to the presence of multiple stable states. So some prey-predator types of models with prey refuge have been proposed and analyzed. Based on the principle of biomass conversion, González-Olivares et al. [22] studied the dynamic consequences of a continuous models with the role of prey's refuge, they analytically proved that the local stability of equilibria and existence of limit cycle. In terms of the work of González-Olivares, Chen et al. [23] further gave a complete dynamical analysis for the same model in [22] by proving the global stability of the positive equilibrium and the uniqueness of limit cycle when exists. Kumar [24] dealt with a prey-predator model with Holling type II response function by taking the refuge of prey into consideration. By the qualitative analysis of ordinary differential equation, they demonstrated that a refuge of prey is an important biological control factor for a pest and the density of prey can increase and even reach a peak as the amount of refuge increased. With a Holling type III response function instead of a Holling type II, Huang et al. [25] illustrated that refuge has an important impact on the stability of prey-predator interactions by using the similar mathematical analysis as in [24]. In a lake ecosystem, the refuge for phytoplankton is also a common phenomenon. Schindler et al. [26] claimed that benthic sediments can provide the phytoplankton with a refuge to produce eggs or be dormant, contributing to temporal escape from the predation of zooplankton. Wiles et al. [27] showed that stratification of water column can also be used as a temporary refuge for phytoplankton to recover. However, the effect of such refuges of phytoplankton has not been well considered in modelling the interactions between phytoplankton and zooplankton.
In freshwater lakes, zooplankton as secondary consumers feed on phytoplankton, which obviously shows a prey-predation mechanism between phytoplankton and zooplankton. Meanwhile, some phytoplankton have the ability to utilize refuge and release toxin to protect themselves from grazing. Therefore, it is reasonable and interesting to establish a new model by considering the impact of the above factors on the interaction between phytoplankton and zooplankton simultaneously to explore the mechanism of algae bloom in lakes. Hence, our aim of this study is to explore how both refuge and toxin affect the forming and ending of algal bloom in a lake ecosystem and to better understand the mechanism of the blooms for the purpose of their prevention and control.
In section 2, we will first propose a two-dimensional phytoplankton-zooplankton model considering the impacts of toxin released by phytoplankton and refuge of phytoplankton. And then we will carry out a complete dynamical analysis of the model to study bifurcations and related dynamics, including Hopf and Bogdanov-Takens bifurcations, which reveals the complex dynamics and the reason behind the complexity. Finally, some numerical simulations will help us to have an intuitionistic view for the effects of refuge and toxin on the occurrence and termination of algal bloom.
2. The model with toxin-producing and refuge of phytoplankton. Based on the prey-predation relationship between phytoplankton-zooplankton, the following general model has been introduced and studied in [28,29,30], where P and Z denote the population size of phytoplankton and zooplankton respectively; φ 1 (P ) represents the predatory response function; r is the intrinsic growth rate of phytoplankton; K represents the environmental carrying capacity of phytoplankton; β 1 is the specific predation rate of zooplankton; β 2 denotes the conversion rate of consumed phytoplankton into zooplankton; d is regarded as the mortality rate of zooplankton due to natural death as well as higher predation. We also assume that total toxic phytoplankton population follows a logistic growth. When the influence of toxin is taken into consideration, the grazing pressure of zooplankton due to the existence of toxin will be decreasing according to the observation in [15]. In this case, an additional mortality of zooplankton can be observed due to insufficient food acquisition. Therefore, the following model framework adopted from (1) will be used to study the toxin liberation processes for planktonic dynamics, where θ denotes the rate of toxin production per phytoplankton species; ψ(P ) represents the distribution of toxic substances. We will consider a simple case by assuming that the total phytoplankton are responsible for the release of toxic substances. The amount of toxin increases with the growing amount of phytoplankton population. However, toxic substances is limited due to the saturation of phytoplankton, Therefore, it is reasonable to use the Holling type II function, ψ(P ) = P a2+P , to model the distribution of toxic substances, where a 2 is the half saturation constant.
The refuge is favorable for phytoplankton growth by preventing the grazing of zooplankton. For the effect of refuges, one naturally subdivides the total phytoplankton into two groups, one portion is in the refuge and the other portion is outside the refuge. Usually, one would consider that the population size of phytoplankton in refuge as one variables depends on different growth and predation pressure. However, the effectiveness of refuges for phytoplankton depends on various factors, for example, temperature, radiation and winds mainly affect the process of stratification which offers a refuge for phytoplankton in some ways [31]. Due to the complexity of the biological process, we restrict ourselves to study a special case that the refuge capacity denoted by m remains a constant and is taken as a parameter in the model analysis. Besides, we assume that phytoplankton can be able to take full advantage of refuges. Therefore, we consider that the predation of zooplankton on phytoplankton depends on the population size of phytoplankton outside refuges expressed by P − m and the predation of zooplankton on unprotected phytoplankton increases and can approach its maximum biological limit with the saturation of phytoplankton in water body. Motivated by the study in [22], we employ Holling type II functional response, i.e., φ 1 (P ) = P −m a1+(P −m) , to describe the predation of zooplankton on unprotected phytoplankton, where a 1 is the the half saturation constant.
From the above discussion, we will extend the previous model (2) for toxinproducing phytoplankton (P) with refuge and zooplankton (Z) as follows:  (3) that if β 2 ≤ d, the population size of zooplankton will continue to drop or remain unchanged over time, which means that the dynamical behaviors of zooplankton can not be described objectively in the actual situation when the conversion rate of phytoplankton is no more than the mortality rate of zooplankton. Therefore, it makes sense to consider only the case β 2 > d. Besides, when phytoplankton and zooplankton co-exist in lakes, the population size of (protected) phytoplankton and the refuge size can not reach the maximal carrying capacity. Thus we take 0 < m < K in the following study.
From the first equation of the system (3), it is easy to derive that lim inf Proof. Let (P (t), Z(t)) be a solution of system (3) with initial conditions P (0) > m, Z(0) > 0. It is easy to obtain from first equation of system (3) that lim sup hence, for sufficiently small ε > 0, there is a sufficiently large T > 0 such that P (t) < K + ε for all t > T . Define the function w(P, Z) = cP (t) + Z(t), where c = β2 β1 , then the time derivative of w(P, Z) along the solutions of system (3) is given by ≤ −min{r, d}w + 2cr(K + ε).
And then we have which follows from the Gronwall inequality [32]. Therefore, for sufficiently large t, w(t) ≤ 2cr(K+ε) min{r,d} holds. The proof is finished.
This lemma implies that the set Ω = (P, min{r,d} is an invariant set. For practical biological significance, we will focus our discussion about system (3) in Ω.
Obviously, E 0 (K − m, 0) is always a boundary equilibrium of system (8). Suppose (P * , Z * ) is a positive equilibrium of the system (8), by putting the right-hand side of system (8) to zero, except for E 0 , we have where P * satisfies the quadratic equation Of special note here is that P * also satisfies 0 < P * < K − m. The interior equilibria (P * , Z * ) of system (8) in Ω 2 exist if and only if the equation (10) has positive root(s) in the interval (0, K − m).
In order to explore the effect of toxin and refuge on the interactions between phytoplankton and zooplankton, we will only choose θ and m as two key parameters. We will discuss the number and distribution of the equilibria of system (8) for the parameters in the region Λ = {(θ, m)|0 < θ, 0 < m < K} of plane (θ, m).
For the convenience of expression, we first define where d 1 < d 2 < d 3 < d 0 and the meaning of this definition will become clear in the following statement. Since f (0) = c < 0, f (P ) will have positive root(s) in three cases. In this case, we will discuss the existence conditions of the positive root(s) of the equation (10) in details. It is easy to have that the equation (10) has positive roots only and only if any one of three cases is satisfied in Fig. 1(a)-(c).
One can verify that −b/(2a) = 0 defines a curve Γ 1 , By c < 0, if the positive solution of the equation (10) (13), we can get that if d < d 0 with a 2 > a 1 , the positive root of (10) might exist. Or else, the equation (10) has no any positive root. Therefore, the existence conditions of the positive root of the equation (10) will be further investigated in the following only when d < d 0 .
With d < d 0 , along the curve Γ 1 , the derivative of m 1 (θ) with respect to θ is In addition, it is easy to have that m 1 (θ) with respect to θ is monotone decreasing. Besides, θ = β 2 − d is a vertical asymptote of m 1 (θ) and m 1 (∞) → −a 1 . By calculating, the curve Γ 1 travels through two points , 0) and , K) in the region Λ. ∆ = 0 can define two curves. By writing the discriminant ∆ = b 2 − 4ac as a function of θ and m, a straightforward calculation leads to then solving the equation ∆(m, θ) = 0 for m, one can get two curves In order to ensure the existence of the root of the equation (10), we first need ∆ ≥ 0. From the above discussion, it is not difficult to obtain that ∆ > 0 if and only if m > m 2 (θ) or m < m 3 (θ) and ∆ = 0 if and only if m = m 2 (θ) = m 3 (θ). Combined with the above discussion on the curve Γ 1 , by comparing, we have m 3 (θ) < m 1 (θ) < m 2 (θ). Thus, (10) may have two positive roots if 0 < m < m 3 (θ) and an unique positive root if 0 < m = m 3 (θ), otherwise, there is no root. With these, one can see that the curve Γ 3 pays a key role in the existence of the positive root(s) of (10).
when m 4 = 0. Furthermore, the derivative of m 4 (θ) with respect to θ is given by the following form Finally, f (K − m) = 0 defines a curve Γ 5 , which implies that Since θ > β 2 − d, we get from the characters of the curve Γ 5 that equation (10) When d < d 2 , we have that only three curves Γ 3 , Γ 4 and Γ 5 are critical to determine the existence of the positive root(s) of (10).
By solving two equations (18) and (20) simultaneously, we can have that two curves Γ 4 and Γ 5 uniquely intersect in the region Λ when d < d 1 . Through further validation, the intersection point also satisfies the equation of Γ 3 . That is, when d < d 1 , three curves have an unique common intersection point (θ * , m * ). Otherwise, three curves can not intersect at the same point in the region Λ. Let (θ * , m * ) be the intersection point, by computing, then we get According the above analysis, the relative locations of Γ 3 , Γ 4 and Γ 5 can be approximately depicted in the region Λ (see Fig. 2(a) and (b)).
(1c) when d ≥ d 2 , then system (8) does not have any positive equilibrium.
Based on the definition of the curve Γ 5 in the region Λ, As a consequence, we reach a conclusion that the positive equilibrium is existent and unique for system (8) Case 3. a > 0, i.e., θ < β 2 − d By using c < 0, it is easy to see that the equation (10) has an unique positive root if and only if f (K − m) > 0. See Fig. 1(e).
According to the analysis of the curve Γ 5 , f (K − m) > 0 is also equivalent to 0 < m < m 5 (θ) in the case. By the monotonicity of m 5 (θ), then m 5 (0) > 0, i.e., d < d 3 , must be satisfied, which means there is an unique positive equilibrium for system (8)  From the discussion above, three curves Γ 3 , Γ 4 and Γ 5 are important to determine the existence and the numbers of the positive equilibrium points of system (8). Here we define C 1 , C 2 and C 1 as which divide the whole region Λ into the following subregions where the positive equilibrium points exist: therefore, we summarize the above discussion about the existence and number of the positive equilibria in the following theorem.
According to the definition of the curve Γ 5 , by simply rearranging the expression of f (K − m), we have , the ratio of growth to mortality of zooplankton at the critical situation. Based on the preceding assumption 0 < m < K, then f (K − m) > 0 if and only if R 0 > 1 and f (K − m) < 0 if and only if R 0 < 1.
For d < d 1 and m ∈ (0, m * ), the curve Γ 3 is also a threshold line which determines the change of the number of positive equilibria from 1 to 2 or from 2 to 1, which is clearly seen from Fig. 2(c). The curve Γ 3 is defined by 4a(a1+K−m)(da2+dK+θK) and 0 <R 0 < 1 due to θ > β 2 − d. We summarize the above discussion into the following proposition. (d) if R 0 <R 0 < 1, system (8) has none positive equilibria.
Proof. From (27), then the Jacobian matrix of system (8) The corresponding characteristic polynomial of the model (8) at the interior equilibrium point E 1 is given by When (θ, m) ∈ Λ 1 or Λ 2 or C 2 , one can obtain from Theorem 3.1 that the positive equilibrium E 1 is unique for system (8). For E 1 , by Theorem (2.1), the derivative of f (P ) with respect to P 1 satisfies f ′ (P 1 ) = (2aP 1 + b)Z 1 > 0, which is easy to be obtained by Fig. 1, then (2aP 1 + b)(a 2 + m + P 1 )P 1 Z 1 > 0. By using the Routh-Hurwitz criterion, we have the real part of two roots of (30) are negative if h(P 1 ) > 0. Therefore, E 1 is asymptotically stable if h(P 1 ) > 0.
When (θ, m) ∈ Λ 3 , system (8) has another positive equilibrium E 3 except for E 1 , where P 1 < P 3 . In this case, the discussion about the stability of system (8) at E 1 is the same as the above. For E 3 , the Jacobin of system (8) at E 3 can be given by (29) where P 1 is replaced by P 3 . Note that, however, the derivative of f (P ) with respect to P 3 satisfies f ′ (P 3 ) = (2aP 3 + b)Z 3 < 0, which can be obtained easily by Fig. 1(c), then (2aP 3 + b)(a 2 + m + P 3 )P 3 Z 3 < 0, which implies that E 3 is always unstable.
If h(P 1 ) = 0, then (30) has a pair of pure imaginary complex conjugate roots. To discuss the effect of refuges, we here select the parameter m as a bifurcation parameter. Let λ = v(m)±iω(m) be the complex conjugate roots and there exists a critical valuem satisfying H(P 1 (m)) = 0, then the transversality condition satisfies which means that Hopf bifurcation can occur. The proof is finished.
When (θ, m) ∈ C 1 , we have a similar conclusion for E 2 .

4.3.
Saddle-node bifurcation for (θ, m) ∈ C 3 . According to Theorem 3.1, when (θ, m) ∈ C 3 , two positive equilibria E 1 and E 3 coalesce at an unique positive equilibrium denoted by E * (P * , Z * ) with P When (θ, m) ∈ C 3 , we have from (16) that m satisfies And a straightforward calculation gives that one eigenvalue of Jacobian matrix of system (8) at E * is zero, the other is −η(a2+m+P * ) P * h(P * ). By (28), we have that h(P * ) = 0 is equivalent to In the following, we will consider the dynamical behavior of system (8) at the equilibrium E * in two cases: K = K * and K = K * , respectively.
When K = K * , the relevant positive equilibrium E * is a Lyapunov-type equilibrium. The transformation x = P − P * and y = Z − Z * (36) brings E * to the origin, then system (8) in a neighborhood of the origin becomes where φ(x, y) = Z * ax 2 + ax 2 y.
Let the right hand side of the first equation of (37) be zero, we obtain a function for x in terms of y by implicit function theorem [34], x = u(y) = − (P * ) 2 ηh(P * ) y + · · · with u(0) = 0.
which shows that the equilibrium E * is a stable saddle-node [35]. Therefore, we have the following theorem.
Theorem 4.4. If m = m * and K = K * , the equilibrium E * of system (8) is a stable saddle-node, where m * and K * are given in (34) and (35) respectively.

4.4.
Bogdanov-Takens bifurcation for (θ, m) ∈ C 3 . Based on the discussion in subsection 4.3, we have that system (8) at E * has a zero eigenvalues with multiplicity 2 when K = K * , which suggests that the equilibrium E * may undergo a Bogdanov-Takens bifurcation. Through the further analysis, the following theorem will be given and proved. Proof. Similar to the process which leads to (37), then we have where R i0 (i = 1, 2) is C ∞ in (x, y) and R i0 = O(|(x, y)| 3 ).
Hence, the equilibrium E * of system (8) is a cusp of codimension 2 [36]. The proof is completed.

4.5.
Hopf bifurcation. From Theorem 4.2 and 4.3, the positive equilibrium E 1 (or E 2 ) is a center-type non-hyperbolic equilibrium of system (8) when h(P 1 ) = 0 (or h(P 2 ) = 0). In this case, system (8) may undergo a Hopf bifurcation at the equilibrium E 1 (or E 2 ). Let E(P ,Z) be any one of the candidates E 1 and E 2 .
To determine the stability of the positive equilibrium E(P ,Z) and the direction of the Hopf bifurcation, the relevant Liapunov coefficients of the above mentioned equilibria will be computed.
5. Numerical simulations. In this section, the previous theoretical results will be numerically performed by using MATLAB 7.1 and an assistive software Matcont (a graphical Matlab package). The relevant parameter values are referenced from [15,16](see Table 1 for details). Note that let E(P, Z) be any one positive equilibrium of system (8), then the coordinates of the related positive equilibrium E(P * , Z * ) of system (3) is P * = P + m, Z * = β 1 Z.
We first study the impact of toxin released by phytoplankton on the dynamical behaviors of system (3) through system (8) more intuitively, the Hopf bifurcation diagram for system (8) with θ as a bifurcation parameter is drawn in Fig. 4, where the value of m is fixed as 8. Fig. 4 shows that system (8) undergoes two Hopf bifurcation points (H). The values of θ for H from left to right are 0.194742 and 0.201745, where the corresponding first Liapunov numbers are −5.531406e − 4 and −4.126233e − 4 respectively. This means that a Hopf bifurcation occurs, which generates limit cycles for model (8).
From Fig. 5, we find that the related dynamical behavior of system (3) depends on the toxin level released by phytoplankton. The phytoplankton and zooplankton is able to be in a stable state with θ small enough. With the increasing of θ, system (3) has a limit cycle which is possible to be caused by the heavy growth of phytoplankton and the mass death of zooplankton because of high concentrations of toxin, which means that algal bloom is likely to occur. By comparing the above case with less toxin, if one keeps θ increasing, the system then is able to become stable again where the abundance of phytoplankton will be at a higher level. If θ is large enough, system (3) is still stable, but in this case, zooplankton due to 12     the powerful toxic substances dies out while the size of phytoplankton due to the extinction of zooplankton reaches the environmental carrying capacity. Thus, we can conclude that, toxin produced by phytoplankton has a significant effect on the occurring and terminating algal blooms.
Similarly, in order to understand the influence of the refuge of phytoplankton on the dynamical behavior of system (3) more clearly, our numerical experiments are carried out by varying the value of refuge capacity m with other parameter values fixed. With θ = 0.2, Fig. 6 depicts the the equilibrium curve of system (8)  At the same time, the equilibrium level of phytoplankton due to the refuge becomes large as the value of m grows, which shows that the population of phytoplankton can increase because the refuge effect protects phytoplankton from the predation of zooplankton.
From Fig. 7, we observe that the related dynamical behavior of system (3) also relies on the refuge capacity m. If m is not large, there is a limit cycle for system (3). However, system (3) has a global asymptotically stable equilibrium with the limit cycle disappearing as m becomes progressively larger. Keeping m increasing, the system is able to become oscillatory again where the abundance of phytoplankton will be at a higher level by comparing the case with small m value. With larger m, system (3) becomes stable again. If the value of m is large enough, however, zooplankton due to the lack of touchable phytoplankton becomes extinct. But phytoplankton population due to the excessive protection of the refuge effect can achieve the environmental carrying capacity. We therefore come to a conclusion that the refuge also pays an positive role in the growth of phytoplankton and has an impressible effect on the occurrence and termination of algal blooms. With the purpose of exploring the influences of toxin and refuge on the dynamical behavior of our systems in combination, by choosing m and θ as two bifurcation parameters, the corresponding bifurcation diagram is presented in Fig. 8 respectively. One Bogdanov-Takens point (BT) can be detected. According to the theoretical analysis, there is an equilibrium with a double zero eigenvalues for our system at BT point.
Beside, in order to study the complicated dynamics of system (8) at the critical positive equilibrium E(P * , Z * ) in subsection 4.3 and 4.4, we first fix the value of θ (θ = 0.25) and choose K and m as two bifurcation parameters to verify the existence of BT bifurcation of codimension 2 in system (8) by numerical simulations. In this case, we have that there is an unique equilibrium (P * , Z * ) = (3.1257, 1.0368) for (K * , m * ) = (12.7759, 0.6146). By perturbing (K, m) in a small neighborhood of (K * , m * ), Fig. 9 shows that there is an unique limit cycle for small values (ǫ 1 , ǫ 2 ) = (−0.05, −0.03). Second, we select θ as an additional bifurcation parameter to checkout the existence of BT bifurcation of codimension 3. The unique equilibrium (P * , Z * ) = (3.1257, 1.0368) is calculated when (K * , m * , θ * ) = (12.7759, 0.6146, 0.2681). By perturbing (K, m, θ) in a small neighborhood of (K * , m * , θ * ), the numerical simulations display that two limit cycles can appear for  (ǫ 1 , ǫ 2 , ǫ 3 ) = (0.0191099, 0.000836, −0.0001), see Fig. 10-13. Therefore, one can see that BT bifurcation with different codimension can occur in system (8) when the bifurcation parameters are chosen properly, which is consistent with our theoretical analysis as before.
6. Discussion. In this paper, we have studied a model incorporating the effects of both refuge and toxin of phytoplankton on the phytoplankton-zooplankton interactions. Despite its simplicity, our model still gives fascinating insights into the phenomenons occurring in lakes.
Comparing the study about the phytoplankton-zooplankton model with toxin in [15], the complex dynamical behaviors can be observed for a suitable set of parameter values in our model with the effect of phytoplankton refuge. By selecting the related parameter values of toxin release rate and refuge at some certain levels, our  Figure 10. The phase portraits of system (8) by taking (K, m, θ) = (K * , m * , θ * ) = (12.7759, 0.6146, 0.2681) and perturbing (K, m, θ) in a small neighborhood of (K * , m * , θ * ). The partial enlarged details of S 1 ,S 2 and S 3 marked are shown by following fig.  11-13. analysis indicates that two different positive equilibria can exist. And more complicated bifurcations, for example, Bogdanov-Takens bifurcation of high codimension, can be shown. These complexities are unlikely to be explained by some previous models which only consider the effect of toxic chemicals.
The previous results in [15,16,17,18,19] have demonstrated that the level of toxicity plays a decisive role in the termination of plankton blooms. In order to verify those statements, we fix the value of refuge capacity and vary the value of toxin production rate to explore the dynamical behaviors of our proposed system. With the refuge capacity in a proper range, our study suggests that the phenomenon of periodic oscillations disappears and the ecological system becomes more stable when the concentration of toxin overtakes a threshold, which means that the same result can be obtained in our study.
However, our study also shows more complicated interactions between phytoplankton and zooplankton. Especially, when the refuge capacity is fixed in a certain range, the periodic oscillations can be developed from nonexistence to existence as toxin production rate varies. In biologically, this means that various amounts of toxin released by phytoplankton may act as a biological trigger for algal blooms. In addition, when the system is in a stable state for small toxin production rate, the equilibrium level of phytoplankton can increase as the the toxin production rate increases, which is mainly because the death of zooplankton due to the increasing of toxin increases.
For the sake of studying the effect of refuge, keeping toxin production rate at a fixed level, we discuss the dynamical behaviors shown in our model by changing the value of refuge capacity. Our study indicates that the refuge does make a  contribution to increasing the population of phytoplankton and reducing the opportunity of phytoplankton from extinction since it can lead to the decrease of the predation ability of zooplankton, which is consistent with the previous finding about predation-prey interactions in [21]. Besides, our findings suggest that, when refuge capacity is very small, the phenomenon of periodic oscillations can fade away gradually with the increase of phytoplankton population in refuge. When refuge capacity further increases, however, our system can experience the oscillations of the population sizes of phytoplankton and zooplankton and then become stable again. Therefore, the moderate refuge of phytoplankton may serve as a biological control not only for the termination of algal blooms but also for the occurrence of algal blooms.
From our study, it is clear that the dynamical behaviors of system (3) we propose rely on the refuges of phytoplankton and toxin produced by phytoplankton. And the significant impacts of both phytoplankton refuge and toxin on the occurrence or termination of algal bloom should not be ignored in lake modelling.
A good understanding of two natural behaviors of some phytoplankton, refuge and toxin that impact on the dynamical behaviors of lake ecosystem, is an importance when making decisions regarding prediction and control of algal bloom in lakes. Besides refuge and toxin, however, a lot of factors mentioned in [5,6,7,8,9,10,11,12] and references therein may effect the phytoplankton-zooplankton interactions and the occurrence and termination of algal bloom in lakes. Therefore, the mechanism of algal bloom in lakes is still needed to be studied continually.