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

  • Received: 18 December 2015 Published: 01 April 2017
  • MSC : Primary: 92D25, 92D40; Secondary: 34A47

  • 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 Ⅱ 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.

    Citation: Juan Li, Yongzhong Song, Hui Wan, Huaiping Zhu. Dynamical analysis of a toxin-producing phytoplankton-zooplankton model with refuge[J]. Mathematical Biosciences and Engineering, 2017, 14(2): 529-557. doi: 10.3934/mbe.2017032

    Related Papers:

    [1] Xinyou Meng, Jie Li . Stability and Hopf bifurcation analysis of a delayed phytoplankton-zooplankton model with Allee effect and linear harvesting. Mathematical Biosciences and Engineering, 2020, 17(3): 1973-2002. doi: 10.3934/mbe.2020105
    [2] Saswati Biswas, Pankaj Kumar Tiwari, Yun Kang, Samares Pal . Effects of zooplankton selectivity on phytoplankton in an ecosystem affected by free-viruses and environmental toxins. Mathematical Biosciences and Engineering, 2020, 17(2): 1272-1317. doi: 10.3934/mbe.2020065
    [3] Jean-Jacques Kengwoung-Keumo . Competition between a nonallelopathic phytoplankton and an allelopathic phytoplankton species under predation. Mathematical Biosciences and Engineering, 2016, 13(4): 787-812. doi: 10.3934/mbe.2016018
    [4] Zhichao Jiang, Xiaohua Bi, Tongqian Zhang, B.G. Sampath Aruna Pradeep . Global Hopf bifurcation of a delayed phytoplankton-zooplankton system considering toxin producing effect and delay dependent coefficient. Mathematical Biosciences and Engineering, 2019, 16(5): 3807-3829. doi: 10.3934/mbe.2019188
    [5] He Liu, Chuanjun Dai, Hengguo Yu, Qing Guo, Jianbing Li, Aimin Hao, Jun Kikuchi, Min Zhao . Dynamics induced by environmental stochasticity in a phytoplankton-zooplankton system with toxic phytoplankton. Mathematical Biosciences and Engineering, 2021, 18(4): 4101-4126. doi: 10.3934/mbe.2021206
    [6] Wenjie Yang, Qianqian Zheng, Jianwei Shen, Linan Guan . Bifurcation and pattern dynamics in the nutrient-plankton network. Mathematical Biosciences and Engineering, 2023, 20(12): 21337-21358. doi: 10.3934/mbe.2023944
    [7] Zhenyao Sun, Da Song, Meng Fan . Dynamics of a stoichiometric phytoplankton-zooplankton model with season-driven light intensity. Mathematical Biosciences and Engineering, 2024, 21(8): 6870-6897. doi: 10.3934/mbe.2024301
    [8] Ruiqing Shi, Jianing Ren, Cuihong Wang . Stability analysis and Hopf bifurcation of a fractional order mathematical model with time delay for nutrient-phytoplankton-zooplankton. Mathematical Biosciences and Engineering, 2020, 17(4): 3836-3868. doi: 10.3934/mbe.2020214
    [9] Jean-Jacques Kengwoung-Keumo . Dynamics of two phytoplankton populations under predation. Mathematical Biosciences and Engineering, 2014, 11(6): 1319-1336. doi: 10.3934/mbe.2014.11.1319
    [10] Elvira Barbera, Giancarlo Consolo, Giovanna Valenti . A two or three compartments hyperbolic reaction-diffusion model for the aquatic food chain. Mathematical Biosciences and Engineering, 2015, 12(3): 451-472. doi: 10.3934/mbe.2015.12.451
  • 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 Ⅱ 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 phytoplankton-zooplankton 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 Ⅱ, (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 nutrient-phytoplankton-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 nutrient-phytoplankton 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ˊalez-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ˊalez-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 Ⅱ 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 Ⅲ response function instead of a Holling type Ⅱ, 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],

    {dPdt=rP(1PK)β1ϕ1(P)Z,dZdt=β2ϕ1(P)ZdZ, (1)

    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,

    {dPdt=rP(1PK)β1ϕ1(P)Z,dZdt=β2ϕ1(P)ZdZθψ(P)Z, (2)

    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 Ⅱ function, ψ(P)=Pa2+P, to model the distribution of toxic substances, where a2 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 Pm 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 Ⅱ functional response, i.e., ϕ1(P)=Pma1+(Pm), to describe the predation of zooplankton on unprotected phytoplankton, where a1 is the the half saturation constant.

    From the above discussion, we will extend the previous model (2) for toxin-producing phytoplankton (P) with refuge and zooplankton (Z) as follows:

    {dPdt=rP(1PK)β1(Pm)a1+(Pm)Z,dZdt=β2(Pm)a1+(Pm)ZdZθPa2+PZ, (3)

    where all parameters are nonnegative and their biological interpretations are summarized in Table 1.

    Table 1. The biological interpretations of all parameters in system (3) with default values used for numerical studies.
    Par.DescriptionValueUnitReference.
    rGrowth rate of phytoplankton0.2h10.07-0.28 [16]
    KEnvironmental carrying capacity50l1108 [15]
    mRefuge capacityPar.l1Defaulted
    β1Predation rate of zooplankton1 h10.6-1.4[16]
    β2Growth efficiency of zooplankton0.15h10.2-0.5[16]
    dMortality rate of zooplankton0.003h10.015-0.15[16]
    a1Half saturation constant3l1Defaulted
    a2Half saturation constant5.7 l15.7[15]
    θToxin production ratePar.h1Defaulted
     | Show Table
    DownLoad: CSV

    It should be noted from the second equation of model (3) that if β2d, 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 inftP(t)m for any solutions (P(t),Z(t)) with initial values P(0)>m,Z(0)>0. And the second equation of the system (3) implies that Z=0 is invariant. Furthermore, we have the following lemma.

    Lemma 2.1. Solutions of system (3) with initial conditions P(0)>m,Z(0)>0 are positive and bounded for all t>0.

    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 suptP(t)K, (4)

    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

    dwdt=cdPdt+dZdt=crP(1PK)dZθPγ+PZcrPdZ(dZ+rcP)+2crPmin{r,d}w+2cr(K+ε). (5)

    And then we have

    0<w(P,Z)2cr(K+ε)min{r,d}(1emin{r,d}t)+w(P(0),Z(0))emin{r,d}t, (6)

    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,Z)R2+|0cP+Z2cr(K+ε)min{r,d}} is an invariant set. For practical biological significance, we will focus our discussion about system (3) in Ω.


    3. The existence and number of equilibria of system (3)

    In order to simplify system (3), we introduce new variables and rescale the time as:

    ˆP=Pm,ˆZ=β1Z,τ=1(a2+P)(a1+(Pm))t, (7)

    if we denote η=rK, then system (3) becomes the following form (here we still use symbols P,Z,t instead of ˆP,ˆZ,τ respectively)

    {dPdt=(a2+m+P)(η(P+m)(KmP)(a1+P)PZ),dZdt=β2P(a2+m+P)Zd(a1+P)(a2+m+P)Zθ(P+m)(a1+P)Z. (8)

    Note that Ω1={(P,Z)|KP(t)>m,Z(t)>0 for all t>0} for system (3) equivalently transforms to Ω2={(P,Z)|KmP(t)>0,Z(t)>0 for all t>0}, which means that we will study system (8) in Ω2.

    Obviously, E0(Km,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 E0, we have

    Z=η(P+m)(KmP)(a1+P)P, (9)

    where P satisfies the quadratic equation

    f(P)=aP2+bP+c=0, (10)

    with

    a=β2dθ,b=am+(β2d)a2(d+θ)a1,c=a1(da2+dm+θm)<0. (11)

    Of special note here is that P also satisfies 0<P<Km. 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,Km).

    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

    d0=β2(1a1a2),  d1=K2β2(K+a1)(K+a2)(1a1a2),d2=Kβ2K+a1(1a1a2),  d3=Kβ2K+a1, (12)

    where d1<d2<d3<d0 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.

    Case 1. a<0, i.e., θ>β2d

    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).

    Figure 1. Existence of the positive roots of f(P) with f(c)<0.

    One can verify that b/(2a)=0 defines a curve Γ1,

    Γ1:m1(θ)=(d+θ)a1(β2d)a2β2dθ. (13)

    By c<0, if the positive solution of the equation (10) exists for θ>β2d, then b/(2a)>0 must hold, which implies 0<m<m1(θ). From θ>β2d and (13), we can get that if d<d0 with a2>a1, 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<d0.

    With d<d0, along the curve Γ1, the derivative of m1(θ) with respect to θ is

    m1(θ)=(β2d)a2β2a1(β2dθ)2<0. (14)

    In addition, it is easy to have that m1(θ) with respect to θ is monotone decreasing. Besides, θ=β2d is a vertical asymptote of m1(θ) and m1()a1. By calculating, the curve Γ1 travels through two points ((β2d)a2da1a1,0) and ((β2d)(K+a2)da1K+a1,K) in the region Λ.

    Δ=0 can define two curves. By writing the discriminant Δ=b24ac as a function of θ and m, a straightforward calculation leads to

    Δ(m,θ)=b24ac=(am+(β2d)a2+(d+θ)a1)24β2a1a2θ,

    then solving the equation Δ(m,θ)=0 for m, one can get two curves

    Γ2:m2(θ)=(d+θ)a1(β2d)a22β2a1a2θβ2dθ, (15)
    Γ3:m3(θ)=(d+θ)a1(β2d)a2+2β2a1a2θβ2dθ. (16)

    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>m2(θ) or m<m3(θ) and Δ=0 if and only if m=m2(θ)=m3(θ). Combined with the above discussion on the curve Γ1, by comparing, we have m3(θ)<m1(θ)<m2(θ). Thus, (10) may have two positive roots if 0<m<m3(θ) and an unique positive root if 0<m=m3(θ), 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).

    For Γ3, θ=(β2d)a2da1±2β2da1(a2a1)a1 when m3=0, θ=β2d is a vertical asymptote of m3(θ), and the derivative of m3(θ) with respect to θ is

    m3(θ)=(β2a1a2β2a1a2θa1)(β2dθ)((d+θ)a1+(β2d)a22β2a1a2θ)(β2dθ)2<0. (17)

    Solving Km=b/(2a), we get a curve Γ4,

    Γ4:m4(θ)=2K+(β2d)a2(d+θ)a1β2dθ=2Km1(θ), (18)

    where Km>b/(2a) is equivalent to 0<m<m4(θ).

    From (18), θ=β2d is a vertical asymptote of m4(θ), m4(0)=2K+(β2d)a2da1β2d, m4()2K+a1, and θ=(β2d)(2K+a2)da12K+a1 when m4=0. Furthermore, the derivative of m4(θ) with respect to θ is given by the following form

    m4(θ)=(β2d)a2β2a1(β2dθ)2=m1(θ)>0. (19)

    Finally, f(Km)=0 defines a curve Γ5,

    Γ5:m5(θ)=K(Kd+Kθ+da2)a1(β2d)(K+a2)Kθ, (20)

    which implies that f(Km)<0 if and only if m>m5(θ) and f(Km)>0 if and only if 0<m<m5(θ).

    From (20), we have θ=(β2d)(1+a2K) is a vertical asymptote of m5(θ), m5(0)=Kda1β2d<K, m5()K+a1 and θ=(K+a2)(β2KdKda1)K(K+a1) when m5=0. The corresponding derivative of m5(θ) with respect to θ is

    m5(θ)=Ka1β2(K+a2)((β2d)(K+a2)Kθ)2<0. (21)

    Since θ>β2d, we get from the characters of the curve Γ5 that equation (10) might have positive root(s) in the region Λ when m5(β2d)=K(Kβ2+da2)a1(β2d)a2>0 i.e., d<d2.

    When d<d2, 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<d1. Through further validation, the intersection point also satisfies the equation of Γ3. That is, when d<d1, 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

    (θ,m)=(Ma1a2β2(2Ma1a2β2)2K2,2K+(β2d)a2(d+θ)a1β2dθ), (22)

    where M=2K(β2d)(K+a2)+a1a2β2.

    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)).

    Therefore, when θ>β2d, we can conclude that

    (1a) when d1d<d2, then system (8) has an unique positive equilibrium if 0<mm5(θ).

    (1b) when d<d1, then system (8) has an unique positive equilibrium if 0<m=m3(θ)<m4(θ) or if 0<m<m5(θ) or if 0<m=m5(θ)m4(θ) and system (8) has two positive equilibrium if max{0,m5(θ)}<m<min{m3(θ),m4(θ)}, or else, system (8) does not have any positive equilibrium.

    (1c) when dd2, then system (8) does not have any positive equilibrium.

    Case 2. a=0, i.e., θ=β2d.

    For the equation (10), the positive root exists only when f(Km)>0. See Fig. 1(d).

    Based on the definition of the curve Γ5 in the region Λ, f(Km)>0 is equivalent to 0<m<m5(β2d), where m5(β2d)>0 can be equivalently transformed into d<d2. As a consequence, we reach a conclusion that the positive equilibrium is existent and unique for system (8) only when 0<m<m5(β2d) with d<d2.

    Case 3. a>0, i.e., θ<β2d

    By using c<0, it is easy to see that the equation (10) has an unique positive root if and only if f(Km)>0. See Fig. 1(e).

    According to the analysis of the curve Γ5, f(Km)>0 is also equivalent to 0<m<m5(θ) in the case. By the monotonicity of m5(θ), then m5(0)>0, i.e., d<d3, must be satisfied, which means there is an unique positive equilibrium for system (8) only when 0<m<m5(θ) with d<d3.

    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 C1, C2 and C1 as

    C1={(θ,m)|θ=β2d,0m<m5(θ)};C2={(θ,m)|β2d<θ,0m=m5(θ)<m4(θ)};C3={(θ,m)|β2d<θ,0m=m3(θ)m4(θ)}.

    which divide the whole region Λ into the following subregions where the positive equilibrium points exist:

    Λ1={(θ,m)|0θ<β2d,0m<m5(θ)};Λ2={(θ,m)|β2d<θ,0m<m5(θ)};Λ3={(θ,m)|β2d<θ,m5(θ)<m<m3(θ),m<m4(θ)}.

    therefore, we summarize the above discussion about the existence and number of the positive equilibria in the following theorem.

    Theorem 3.1. For system (8) with m3(θ),m4(θ) and m5(θ) defined as above.

    (I) If d<d1, as shown in Fig. 2(a), then

    Figure 2. The existence and number of equilibrium of system (8) on the plane (θ,m) for different values of d, where the dotted line denotes θ=β2d.

    (a) for (θ,m)Λ1Λ2C2, system (8) has an unique positive equilibrium E1(P1,Z1);

    for (θ,m)C1, system (8) has an unique positive equilibrium E2(P2,Z2);

    for (θ,m)Λ3, system (8) has two positive equilibria E1(P1,Z1) and E3(P3,Z3);

    for (θ,m)C3, E1(P1,Z1) and E3(P3,Z3) coalesce at an unique positive equilibrium of multiplicity 2.

    (b) for (θ,m)Λ/(Λ1Λ2Λ3C1C2C3), system (8) has none positive equilibrium.

    (II) If d1d<d2, as shown in Fig. 2(b), then

    (a) for (θ,m)C1, system (8) has an unique positive equilibrium E2(P2,Z2);

    for (θ,m)Λ1Λ2, system (8) has an unique positive equilibrium E1(P1,Z1);

    (b) for (θ,m)Λ/(Λ1Λ2C1), system (8) has none positive equilibrium.

    (III) If d2d<d3, as shown in Fig. 2(c), then

    (a) for (θ,m)Λ1, system (8) has an unique positive equilibrium E1(P1,Z1);

    (b) for (θ,m)Λ/Λ1, system (8) has none positive equilibrium.

    (IV) If dd3, system (8) has none positive equilibrium.

    Where when exist the corresponding equilibria are:

    E1(P1,Z1)=(b+Δ2a,η(mb+Δ)(Km+bΔ)(a1b+Δ)4a2(b+Δ)), (23)
    E2(P2,Z2)=(cb,η(mc)(Km+c)(a1c)b2c), (24)
    E3(P3,Z3)=(bΔ2a,η(mbΔ)(Km+b+Δ)(a1bΔ)4a2(b+Δ)). (25)

    From Fig. 2, it is easy to see that the curve Γ5 is a threshold line which determines the change of the number of positive equilibria.

    According to the definition of the curve Γ5, by simply rearranging the expression of f(Km), we have

    f(Km)=(a2+Km)(da2+dK+θK)(β2(Km)(a2+K)(a2+Km)(da2+dK+θK)1). (26)

    Let R0=β2(Km)(a2+K)(a1+Km)(da2+dK+θK), the ratio of growth to mortality of zooplankton at the critical situation. Based on the preceding assumption 0<m<K, then f(Km)>0 if and only if R0>1 and f(Km)<0 if and only if R0<1.

    For d<d1 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 Δ=0. Let ˜R0=R0|Δ=0, then ˜R0=1+(2a(Km)+b)24a(a1+Km)(da2+dK+θK) and 0<˜R0<1 due to θ>β2d.

    We summarize the above discussion into the following proposition.

    Proposition 1. Suppose d<d1 and m(0,m), then

    (a) if R01, system (8) has an unique positive equilibrium;

    (b) if ˜R0<R0<1, system (8) has two positive equilibria;

    (c) if R0=˜R0<1, two positive equilibria coalesce at an unique equilibrium of multiplicity 2;

    (d) if R0<˜R0<1, system (8) has none positive equilibria.

    Proposition 1 shows that system (8) will undergo a backward bifurcation which occurs at R0=1 (Fig. 3).

    Figure 3. A backward bifurcation occurs when 0<d<d1 and m(0,m).

    4. Stability and bifurcation analysis

    Let E(P,Z) be any one equilibrium, then the Jacobian matrix of system (8) at E is given by

    JE=(η(a2+m+P)Ph(P)P(a2+m+P)(2aP+b)Zf(P)), (27)

    where

    h(P)=2(P)3+(a1+2mK)(P)2+ma1(Km). (28)

    4.1. The boundary equilibrium E0

    Theorem 4.1. The boundary equilibrium E0 is unstable if R0>1, asymptotically stable if R0<1 and globally asymptotically stable if R0<˜R0<1.

    Proof. One can verify that system (8) at E0 has an eigenvalue, ηK(K+a2)(a1+(Km))<0. The other eigenvalue is f(Km). The results follow from the fact that f(Km)>0 is equivalent to R0>1.

    By using Proposition 1, system (8) has no positive equilibria if R0<˜R0<1. Since E0 is local asymptotically stable, it follows from Poincarˊe-Bendixson theorem [33] that E0 is also globally asymptotically stable. We finish this proof.


    4.2. The positive equilibrium E1,E2 and E3 for (θ,m)Λ1, Λ2, Λ3, C2 and C1

    In this subsection, the aim is to determine the stability and bifurcation of the existing positive equilibriums of system (8).

    Theorem 4.2. Assume (θ,m)Λ1Λ2Λ3C2, then the positive equilibrium E1 of system (8) is asymptotically stable when h(P1)>0 and unstable when h(P1)<0, system (8) will experience a Hopf-bifurcation at E1 when h(P1)=0. Moreover E3 when exists is always unstable.

    Proof.From (27), then the Jacobian matrix of system (8) at E1(P1,Z1) is

    JE1=(η(a2+m+P1)P1h(P1)P1(a2+m+P1)(2aP1+b)Z10). (29)

    The corresponding characteristic polynomial of the model (8) at the interior equilibrium point E1 is given by

    λ2+η(a2+m+P1)P1h(P1)λ+(2aP1+b)(a2+m+P1)P1Z1=0. (30)

    When (θ,m)Λ1 or Λ2 or C2, one can obtain from Theorem 3.1 that the positive equilibrium E1 is unique for system (8). For E1, by Theorem (2.1), the derivative of f(P) with respect to P1 satisfies f(P1)=(2aP1+b)Z1>0, which is easy to be obtained by Fig. 1, then (2aP1+b)(a2+m+P1)P1Z1>0. By using the Routh-Hurwitz criterion, we have the real part of two roots of (30) are negative if h(P1)>0. Therefore, E1 is asymptotically stable if h(P1)>0.

    When (θ,m)Λ3, system (8) has another positive equilibrium E3 except for E1, where P1<P3. In this case, the discussion about the stability of system (8) at E1 is the same as the above. For E3, the Jacobin of system (8) at E3 can be given by (29) where P1 is replaced by P3. Note that, however, the derivative of f(P) with respect to P3 satisfies f(P3)=(2aP3+b)Z3<0, which can be obtained easily by Fig. 1(c), then (2aP3+b)(a2+m+P3)P3Z3<0, which implies that E3 is always unstable.

    If h(P1)<0, then E1 is an unstable equilibrium. Moreover, system (8) has at least one closed orbit in Ω2, which can be deduced by Poincaré-Bendixson Theorem [33].

    If h(P1)=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 value ˆm satisfying H(P1(ˆm))=0, then the transversality condition satisfies

    dvdm|m=ˆm=(dvdP1dP1dm)|m=ˆm0, (31)

    where

    dvdP1=η(a2+m+P1)(3P1+m+2a1K)P1, (32)
    dP1dm=am+(β2d)a2+(d+θ)a1Δ2Δ, (33)

    which means that Hopf bifurcation can occur.

    The proof is finished.

    When (θ,m)C1, we have a similar conclusion for E2.

    Theorem 4.3. Assume (θ,m)C1, then the positive equilibrium E2 of system (8) is asymptotically stable when h(P2)>0; E2 is unstable when h(P2)<0, and Hopf bifurcation occurs if h(P2)=0.


    4.3. Saddle-node bifurcation for (θ,m)C3

    According to Theorem 3.1, when (θ,m)C3, two positive equilibria E1 and E3 coalesce at an unique positive equilibrium denoted by E(P,Z) with P=b2a and Z=η(mb)(Km+b)(a1b)4a2b, where a=β2dθ,b=am+(β2d)a2(d+θ)a1.

    When (θ,m)C3, we have from (16) that m satisfies

    m=m=(d+θ)a1(β2d)a2+2β2a1a2θβ2dθ. (34)

    And a straightforward calculation gives that one eigenvalue of Jacobian matrix of system (8) at E is zero, the other is η(a2+m+P)Ph(P). By (28), we have that h(P)=0 is equivalent to

    K=K=P+m+P(P+m)(a1+P)(P)2ma1, (35)

    In the following, we will consider the dynamical behavior of system (8) at the equilibrium E in two cases: KK and K=K, respectively.

    When KK, the relevant positive equilibrium E is a Lyapunov-type equilibrium. The transformation

    x=PPandy=ZZ (36)

    brings E to the origin, then system (8) in a neighborhood of the origin becomes

    {dxdt=η(a2+m+P)Ph(P)xP(a2+m+P)y+ψ(x,y)dydt=ϕ(x,y), (37)

    where

    ψ(x,y)=η(K2ma13P)(a2+m+P)x2(a2+m+2P)xy+η(K3ma1a24P)x3x2yηx4,ϕ(x,y)=Zax2+ax2y.

    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.

    Substituting x=u(y) into ϕ(x,y), we have

    ϕ(u(y),y)=(P)Zab2h2(P)y2+o(y2), (38)

    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 KK, 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)C3

    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.

    Theorem 4.5. If m=m, K=K and K3P+a1+2m, then the equilibrium E of system (8) is a cusp of codimension 2 (a Bogdanov-Takens bifurcation point), where m and K are given in (34) and (35) respectively.

    Proof. Similar to the process which leads to (37), then we have

    {dxdt=P(a2+m+P)y+η(K2ma13P)(a2+m+P)x2+R10(x,y),dydt=Zax2+R20(x,y), (39)

    where Ri0(i=1,2) is C in (x,y) and Ri0=O(|(x,y)|3).

    Introducing new variables

    x1=xandy1=P(a2+m+P)y, (40)

    system (39) becomes

    {dx1dt=y1+a11x1y1+a20x21+R11(x1,y1),dy1dt=b20x21+R21(x1,y1), (41)

    where Ri1(i=1,2) is C in (x1,y1) and Ri1=O(|(x1,y1)|3), and

    a11=a2+m+2PP(a2+m+P)>0,a20=η(K2ma13P)(a2+m+P),b20=P(a2+m+P)Za>0. (42)

    Making the near-identity transformation

    x2=x1,y2=y1+a11x1y1+a20x21+R11(x1,y1), (43)

    we have

    {dx2dt=y2,dy2dt=b20x22+a11y2y2a20x221+a11x2+2a20x2y2+R22(x2,y2), (44)

    where R22 is C in (x1,y2) and R22=O(|(x2,y2)|3).

    Furthermore, utilizing the near-identity transformation

    dt=(1+a11x2)dτ, (45)

    we obtain

    {dx2dτ=(1+a11x2)y2,dy2dτ=b20x22+2a20x2y2+a11y22+R23(x2,y2), (46)

    where R23 is C in (x2,y2) and R23=O(|(x2,y2)|3).

    By a time reparameterization

    u=x2andv=(1+a11x2)y2, (47)

    we obtain

    {dudτ=v,dvdτ=b20u2+2a20uv+R24(u,v), (48)

    where R24 is C in (u,v) and R24=O(|(u,v)|3).

    According to (47) and P=b2a, it is easy to see that b20>0. And note that K3P+a1+2m, we have a200.

    Hence, the equilibrium E of system (8) is a cusp of codimension 2 [36]. The proof is completed.

    Next, we study system (8) by perturbing the parameters (θ,m) in a small neighborhood of (K,m) when K3P+a1+2m. Thus we let

    {K=K+ϵ1,m=m+ϵ2, (49)

    in system (8) and we discuss the bifurcations of the resulting system

    {dPdt=(a2+m+ϵ2+P)[η(P+m+ϵ2)(K+ϵ1mϵ2P)(a1+P)PZ],dZdt=β2P(a2+m+ϵ2+P)Zd(a1+P)(a2+m+ϵ2+P)Zθ(P+m+ϵ2)(a1+P)Z, (50)

    for sufficiently small (ϵ1,ϵ2).

    System (50) has a cusp point (P,Z) if (ϵ1,ϵ2)=(0,0). Using the substitution

    {ˉP=PP,ˉZ=ZZ, (51)

    and expanding the right side of system (50) into Taylor series about the origin, we get

    {dˉPdt=(p00+p10ˉP+p20ˉP2+p30ˉP3+p40ˉP4)+(p01+p11ˉP+p21ˉP2)ˉZ,dˉZdt=(q00+q10ˉP+q20ˉP2)+(q01+q11ˉP+q21ˉP2)ˉZ, (52)

    where

    p00=η[ϵ2(K2m2Pϵ2)+ϵ1(m+ϵ2+P)](a1+P)(a2+m+ϵ2+P),p10=η{(a2+m+P+ϵ2)[ϵ2(K2m4P2a1ϵ2)+ϵ1(a1+2P+m+ϵ2)]+(a1+P)[ϵ2(K2m2Pϵ2)+ϵ1(m+ϵ2+P)]},p01=P(a2+m+P+ϵ2),p11=(a2+m+2P+ϵ2),p20=η{(K2ma13P2ϵ2)(a2+m+P+ϵ2)+ϵ2(K2m4P2a1ϵ2)+ϵ1(a1+a2+2P+2m+2ϵ2)},p21=1,p30=η{(K3m4Pa1a23ϵ2)+ϵ1},p40=η,q00=ϵ2((β2dθ)P(d+θ)a1)Z,q10=ϵ2(β2dθ)Z,q01=ϵ2((β2dθ)P(d+θ)a1),q11=ϵ2(β2dθ),q20=(β2dθ)Z,q21=β2dθ. (53)

    Following the transformation

    {˜P=ˉP,˜Z=(p00+p10ˉP+p20ˉP2+p30ˉP3+p40ˉP4)+(p01+p11ˉP+p21ˉP2)ˉZ, (54)

    system (52) is rewritten in the following form

    {d˜Pdt=˜Z,d˜Zdt=f00+f10˜P+f01˜Z+f11˜P˜Z+f20˜P2+f02˜Z2+O((˜P,˜Z)3), (55)

    where

    f00=p01q00p00q01,f10=p01q10+p11q00+p00q11p10q01,f01=p10+q01p00p11p01,f11=2p20+q11+p00p211p10p01p112p00p01p21p201,f20=p01q20+p11q10+p21q00p20q01p10q11p00q21+p11(p10q01+p00q11)p01p00(p211p21)q01p201,f02=p11. (56)

    By rescaling the time variable t using

    t=(1f02˜Z)τ, (57)

    system (55) becomes

    {d˜Pdτ=(1f02˜Z)˜Z,d˜Zdτ=(1f02˜Z)(f00+f10˜P+f01˜Z+f11˜P˜Z+f20˜P2+f02˜Z2+O((˜P,˜Z)3)). (58)

    Furthermore, let

    {ˆP=˜P,ˆZ=(1f02˜Z)˜Z, (59)

    then we have

    {dˆPdτ=ˆZ,dˆZdτ=f00+(f102f00f02)ˆP+f01ˆZ+(f20+f00f2022f10f02)ˆP2+(f11f01f02)ˆPˆZ+O((ˆP,ˆZ)3). (60)

    By (53), we can compute

    limϵ10,ϵ20(f11f01f02)=η(K2ma13P)(a2+m+P)0, (61)

    thus setting ˆP=Pf01f11f01f02 and renaming P to be ˆP, then system (60) is equivalent to

    {dˆPdτ=ˆZ,dˆZdτ=g00+g10ˆP+g20ˆP2+g11ˆPˆZ+O((ˆP,ˆZ)3). (62)

    where

    g00=f00f01(f102f00f02)f11f01f02+f201(f20+f00f2022f10f02)(f11f01f02)2,g10=f102f00f022f01(f20+f00f2022f10f02)f11f01f02,g20=f20+f00f2022f10f02,g11=f11f01f02. (63)

    From (53) and (63), we are able to have

    limϵ10,ϵ20g11=η(K2ma13P)(a2+m+P)0, (64)

    and

    limϵ10,ϵ20g20=P(a2+m+P)(β2dθ)Z0, (65)

    therefore, by performing the change of new variables: ξ1=g211g20ˆP,ξ2=g311g220ˆZ,t=g20g11τ, which finally leads that system (3) takes the following generic normal form

    {dξ1dt=ξ2,dξ2dt=g00g411g320+g10g211g220ξ1+ξ21+ξ1ξ2+O((^ξ1,^ξ2)3). (66)

    Corollary 1. If m=m, K=K and K=3P+a1+2m, then the equilibrium E3 of system (8) is a cusp with the codimension at least 3 (a degenerate Bogdanov-Takens bifurcation point), where m and K are given in (34) and (35) respectively.


    4.5. Hopf bifurcation

    From Theorem 4.2 and 4.3, the positive equilibrium E1 (or E2) is a center-type non-hyperbolic equilibrium of system (8) when h(P1)=0 (or h(P2)=0). In this case, system (8) may undergo a Hopf bifurcation at the equilibrium E1 (or E2). Let E(˜P,˜Z) be any one of the candidates E1 and E2.

    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.

    Making the transformation

    x=P˜Pandy=Z˜Z (67)

    then by expanding the right hand side of the corresponding system in a Taylor series about the origin, we get

    {dxdt=a01y+a20x2+a11xy+a30x3+a21x2y+O(x4),dydt=b10x+b11xy+b20x2+b21x2y, (68)

    where

    a01=˜P(a2+m+˜P),a11=(a2+m+2˜P),a20=η(K2ma13˜P)(a2+m+˜P),a21=1,a30=η(K3ma1a24˜P),b10=(2a˜P+b)˜Z,b11=(2a˜P+b),b20=˜Za,b21=a.

    a and b here are given in (11).

    Applying the formula of the Liapunov number σ for the focus at the origin of (68) in [35], we obtain

    σ=3π2Δ32(2a01a20b20a11a20b10+3a01a30b10), (69)

    where it follows from the proof of Theorem 4.2 that Δ=a01b10=(2a˜P+b)(a2+m+˜P)˜P˜Z>0.

    Theorem 4.6. For system (8), fix all parameters except m>0 (or θ>0) and suppose there exists mk>0 (or θk>0) such that h(˜P,mk)=0 (or h(˜P,θk)=0) and h(˜P,m)m|m=mk0 (or h(˜P,θ)θ|θ=θk0) for the equilibrium E(˜P,˜Z). If σ0, then a curve of periodic solutions bifurcation from E(˜P,˜Z) such that

    (i) a supercritical Hopf bifurcation occurs at E(˜P,˜Z) if σ<0;

    (ii) a subcritical Hopf bifurcation occurs at E(˜P,˜Z) if σ>0.


    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=β1Z.

    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.531406e4 and 4.126233e4 respectively. This means that a Hopf bifurcation occurs, which generates limit cycles for model (8).

    Figure 4. The Hopf bifurcation diagram of system (8) with m=8 and θ as a bifurcation parameter.

    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 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.

    Figure 5. The variation of phytoplankton and zooplankton with the increasing time and the phase plane diagram of system (3) with m=8 for different θ values, where the initial value is (20,3).

    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) and limit point bifurcation of cycle as m varies. From this figure, we see that three Hopf bifurcation points for our systems appear when m changes from small to large, where the values of m for H from small to large are 0.003364, 2.485566 and 8.220599, and the corresponding first Liapunov numbers are 1.817476e4, 2.280040e4 and 4.253572e04 respectively.

    Figure 6. The Hopf bifurcation diagram of system (8) with θ=0.2 and m as a bifurcation parameter.

    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.

    Figure 7. The variation of phytoplankton and zooplankton with the increasing time and the phase plane diagram of system (3) with θ=0.2 for different m values, where the initial value is (20,3).

    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.

    Figure 8. Bifurcation diagram of system (8) by choosing m and θ as two parameters, where the values of m and θ for BT are 5.4248527 and 0.20944969 respectively. Note LP marked represents the limit point at which two positive equilibria collide into one positive equilibrium, and BP represents branch point at which positive equilibrium can disappears with the increase of the value of the parameter θ orm.

    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.

    Figure 9. The phase portraits of system (8) by perturbing (K,m) in a small neighborhood of (K,m)=(12.7759,0.6146).
    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 S1, S2 and S3 marked are shown by following fig. 11-13.
    Figure 11. The local phase portraits of system (8) for (ϵ1,ϵ2,ϵ3)=(0.01665,0.0001,0.0001) in Fig.10.
    Figure 12. The local phase portraits of system (8) for (ϵ1,ϵ2,ϵ3)=(0.322,0,0.0001) in Fig. 10.
    Figure 13. The local phase portraits of system (8) for (ϵ1,ϵ2,ϵ3)=(0.0191099,0.000836,0.0001) in Fig. 10.

    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 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.


    [1] [ H. Duan,R. Ma,X. Xu,F. Kong,S. Zhang,W. Kong,J. Hao,L. Shang, Two-decade reconstruction of algal blooms in China-Lake Taihu, Environ. Sci. Technol., 43 (2009): 3522-3528.
    [2] [ M. Wines, Spring Rain, Then Foul Algae in Ailing Lake Erie Report of The New York Times, 2013. Available from: http://www.nytimes.com/2013/03/15/science/earth/algae-blooms-threaten-lake-erie.html?&_r=0.
    [3] [ W. McLean and J. Macdonald, CLAS: Colby Liberal Arts Symposium, Lake Erie Algal Blooms 2014. Available from: http://digitalcommons.colby.edu/clas/2014/program/414/.
    [4] [ C. B. Lopez, E. B. Jewett, Q. Dortch, B. T. Walton and H. K. Hudnell, Scientific Assessment of Freshwater Harmful Algal Blooms Interagency Working Group on Harmful Algal Blooms, Hypoxia, and Human Health of the Joint Subcommittee on Ocean Science and Technology. Washington, DC. 2008. Available from: https://www.whitehouse.gov/sites/default/files/microsites/ostp/frshh2o0708.pdf.
    [5] [ E. Beltrami,T. O. Carroll, Modeling the role of viral disease in recurrent phytoplankton blooms, J. Math. Biol., 32 (1994): 857-863.
    [6] [ J. Norberg,D. DeAngelis, Temperature effects on stocks and stability of a phytoplankton-zooplankton model and the dependence on light and nutrients, Ecol. Model., 95 (1997): 75-86.
    [7] [ B. Mukhopadhyay,R. Bhattacharyya, Modelling phytoplankton allelopathy in a nutrient-plankton model with spatial heterogeneity, Ecol. Model., 198 (2006): 163-173.
    [8] [ S. Pal,S. Chatterjee,J. Chattopadhyay, Role of toxin and nutrient for the occurrence and termination of plankton bloom results drawn from field observations and a mathematical model, Biosystems, 90 (2007): 87-100.
    [9] [ T. Saha,M. Bandyopadhyay, Dynamical analysis of toxin producing phytoplankton-zooplankton interactions, Nonlinear Anal-Real., 10 (2009): 314-332.
    [10] [ Y. Lv,Y. Pei,S. Gao,C. Li, Harvesting of a phytoplankton-zooplankton model, Nonlinear Anal-Real., 11 (2010): 3608-3619.
    [11] [ M. Bengfort,U. Feudel,F. M. Hilker,H. Malchow, Plankton blooms and patchiness generated by heterogeneous physical environments, Ecol. Complex., 20 (2014): 185-194.
    [12] [ S. Rana,S. Samanta,S. Bhattacharya,K. Al-Khaled,A. Goswami,J. Chattopadhyay, The effect of nanoparticles on plankton dynamics: A mathematical model, Biosystems, 127 (2015): 28-41.
    [13] [ T. G. Hallam,C. E. Clark,R. R. Lassiter, Effects of toxicants on populations: A qualitative approach Ⅰ. Equilibrium environmental exposure, Ecol. Model., 18 (1983): 291-304.
    [14] [ J. T. Turner,P. A. Tester, Toxic marine phytoplankton, zooplankton grazers, and pelagic food webs, Limnol. Oceanogr., 42 (1997): 1203-1214.
    [15] [ J. Chattopadhyay,R. R. Sarkar,S. Mandal, Toxin producing plankton may act as a biological control for planktonic blooms-field study and mathematical modelling, J. Theor. Biol., 215 (2002): 333-344.
    [16] [ J. Chattopadhyay,R. R. Sarkar,A. El Abdllaoui, A delay differential equation model on harmful algal blooms in the presence of toxic substances, Math. Med. Biol., 19 (2002): 137-161.
    [17] [ R. Pal,D. Basu,M. Banerjee, Modelling of phytoplankton allelopathy with Monod-Haldane-type functional response-A mathematical study, Biosystems, 95 (2009): 243-253.
    [18] [ S. Chakraborty,S. Chatterjee,E. Venturino,J. Chattopadhyay, Recurring plankton bloom dynamics modeled via toxin-producing phytoplankton, J. Biol. Phys., 3 (2008): 271-290.
    [19] [ M. Banerjee,E. Venturino, A phytoplankton-toxic phytoplankton-zooplankton model, Ecol. Complex., 8 (2011): 239-248.
    [20] [ M. M. Mullin,E. F. Stewart,F. J. Fuglister, Ingestion by planktonic grazers as a function of concentration of food1, Limnol. Oceanogr., 20 (1975): 259-262.
    [21] [ J. B. Collings, Bifurcation and stability analysis of a temperature-dependent mite predator-prey interaction model incorporating a prey refuge, B. Math. Biol., 57 (1995): 63-76.
    [22] [ E. González-Olivares,R. Ramos-Jiliberto, Dynamic consequences of prey refuges in a simple model system: more prey, fewer predators and enhanced stability, Ecol. Model., 166 (2003): 135-146.
    [23] [ L. Chen,F. Chen,L. Chen, Qualitative analysis of a predator-rey model with Holling type Ⅱ functional response incorporating a constant prey refuge, Nonlinear anal-real., 11 (2010): 246-252.
    [24] [ T. K. Kar, Stability analysis of a prey-predator model incorporating a prey refuge, Commun. Nonlinear. Sci., 10 (2005): 681-691.
    [25] [ Y. Huang,F. Chen,L. Zhong, Stability analysis of a prey-predator model with Holling type Ⅲ response function incorporating a prey refuge, Appl. Math. Comput., 182 (2006): 672-683.
    [26] [ D. E. Schindler,M. D. Scheuerell, Habitat coupling in lake ecosystems, Oikos, 98 (2002): 177-189.
    [27] [ P. J. Wiles,L. A. van Duren,C. Häse,J. Larsen,J. H. Simpson, Stratification and mixing in the Limfjorden in relation to mussel culture, J. Marine. Syst., 60 (2006): 129-143.
    [28] [ K. S. Cheng,S. B. Hsu,S. S. Lin, Some results on global stability of a predator-prey system, J. Math. Biol., 12 (1981): 115-126.
    [29] [ L. P. Liou,K. S. Cheng, Global stability of a predator-prey system, J. Math. Biol., 26 (1988): 65-71.
    [30] [ S. B. Hsu,T. W. Huang, Global stability for a class of predator-prey systems, SIAM J, Appl. Math., 55 (1995): 763-783.
    [31] [ J. F. Talling, The annual cycle of stratification and phytoplankton growth in Lake Victoria (East Africa), Int. Revue Ges. Hydrobiol., 51 (1966): 545-621.
    [32] [ V. Lakshmikantham,S. Leela, null, Differential and Integral Inequalities: Theory and Applications, Academic press, 1969.
    [33] [ G. Teschl, Ordinary Differential Equations and Dynamical Systems Am. Math. Soc., 2012.
    [34] [ H. K. Khalil and J. W. Grizzle, Nonlinear Systems Upper Saddle River: Prentice hall, 2000.
    [35] [ P. Lawrence, null, Differential Equations and Dynamical Systems, , Springer-Verlag, New York, 1991.
    [36] [ Y. A. Kuznetsov, null, Elements of Applied Bifurcation Theory, , Spring-Verlag, New York, 1995.
  • This article has been cited by:

    1. 圈利 姬, Spatiotemporal Dynamics of a Planktonic Ecological Model with Inhibitory Effect, 2019, 08, 2324-7991, 2050, 10.12677/AAM.2019.812236
    2. Ting Gao, Xinyou Meng, Stability and Hopf bifurcation of a delayed diffusive phytoplankton-zooplankton-fish model with refuge and two functional responses, 2023, 8, 2473-6988, 8867, 10.3934/math.2023445
    3. Kaushik Dehingia, Anusmita Das, Evren Hinçal, Kamyar Hosseini, The Effect of Time Delay on the Dynamics of a Plankton-Nutrient System with Refuge, 2025, 55, 0103-9733, 10.1007/s13538-024-01670-0
  • Reader Comments
  • © 2017 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(3852) PDF downloads(667) Cited by(3)

Article outline

Figures and Tables

Figures(13)  /  Tables(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog