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
[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 |
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
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.
Based on the prey-predation relationship between phytoplankton-zooplankton, the following general model has been introduced and studied in [28,29,30],
{dPdt=rP(1−PK)−β1ϕ1(P)Z,dZdt=β2ϕ1(P)Z−dZ, | (1) |
where
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(1−PK)−β1ϕ1(P)Z,dZdt=β2ϕ1(P)Z−dZ−θψ(P)Z, | (2) |
where
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,
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
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(1−PK)−β1(P−m)a1+(P−m)Z,dZdt=β2(P−m)a1+(P−m)Z−dZ−θPa2+PZ, | (3) |
where all parameters are nonnegative and their biological interpretations are summarized in Table 1.
Par. | Description | Value | Unit | Reference. |
Growth rate of phytoplankton | 0.2 | 0.07-0.28 [16] | ||
| Environmental carrying capacity | 50 | 108 [15] | |
Refuge capacity | Par. | Defaulted | ||
Predation rate of zooplankton | 1 | | 0.6-1.4[16] | |
Growth efficiency of zooplankton | 0.15 | 0.2-0.5[16] | ||
Mortality rate of zooplankton | 0.003 | 0.015-0.15[16] | ||
Half saturation constant | 3 | Defaulted | ||
| Half saturation constant | 5.7 | | 5.7[15] |
Toxin production rate | Par. | Defaulted |
It should be noted from the second equation of model (3) that if
From the first equation of the system (3), it is easy to derive that
Lemma 2.1. Solutions of system (3) with initial conditions
Proof. Let
lim supt→∞P(t)≤K, | (4) |
hence, for sufficiently small
Define the function
dwdt=cdPdt+dZdt=crP(1−PK)−dZ−θPγ+PZ≤crP−dZ≤−(dZ+rcP)+2crP≤−min{r,d}w+2cr(K+ε). | (5) |
And then we have
0<w(P,Z)≤2cr(K+ε)min{r,d}(1−e−min{r,d}t)+w(P(0),Z(0))e−min{r,d}t, | (6) |
which follows from the Gronwall inequality [32].
Therefore, for sufficiently large
This lemma implies that the set
In order to simplify system (3), we introduce new variables and rescale the time as:
ˆP=P−m,ˆZ=β1Z,τ=1(a2+P)(a1+(P−m))t, | (7) |
if we denote
{dPdt=(a2+m+P)(η(P+m)(K−m−P)(a1+P)−PZ),dZdt=β2P(a2+m+P)Z−d(a1+P)(a2+m+P)Z−θ(P+m)(a1+P)Z. | (8) |
Note that
Obviously,
Suppose
Z∗=η(P∗+m)(K−m−P∗)(a1+P∗)P∗, | (9) |
where
f(P)=aP2+bP+c=0, | (10) |
with
a=β2−d−θ,b=am+(β2−d)a2−(d+θ)a1,c=−a1(da2+dm+θm)<0. | (11) |
Of special note here is that
In order to explore the effect of toxin and refuge on the interactions between phytoplankton and zooplankton, we will only choose
For the convenience of expression, we first define
d0=β2(1−a1a2), d1=K2β2(K+a1)(K+a2)(1−a1a2),d2=Kβ2K+a1(1−a1a2), d3=Kβ2K+a1, | (12) |
where
Case 1.
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
Γ1:m1(θ)=(d+θ)a1−(β2−d)a2β2−d−θ. | (13) |
By
With
m′1(θ)=−(β2−d)a2−β2a1(β2−d−θ)2<0. | (14) |
In addition, it is easy to have that
Δ(m,θ)=b2−4ac=(am+(β2−d)a2+(d+θ)a1)2−4β2a1a2θ, |
then solving the equation
Γ2:m2(θ)=−(d+θ)a1−(β2−d)a2−2√β2a1a2θβ2−d−θ, | (15) |
Γ3:m3(θ)=−(d+θ)a1−(β2−d)a2+2√β2a1a2θβ2−d−θ. | (16) |
In order to ensure the existence of the root of the equation (10), we first need
For
m′3(θ)=(β2a1a2√β2a1a2θ−a1)(β2−d−θ)−((d+θ)a1+(β2−d)a2−2√β2a1a2θ)(β2−d−θ)2<0. | (17) |
Solving
Γ4:m4(θ)=2K+(β2−d)a2−(d+θ)a1β2−d−θ=2K−m1(θ), | (18) |
where
From (18),
m′4(θ)=(β2−d)a2−β2a1(β2−d−θ)2=−m′1(θ)>0. | (19) |
Finally,
Γ5:m5(θ)=K−(Kd+Kθ+da2)a1(β2−d)(K+a2)−Kθ, | (20) |
which implies that
From (20), we have
m′5(θ)=−Ka1β2(K+a2)((β2−d)(K+a2)−Kθ)2<0. | (21) |
Since
When
By solving two equations (18) and (20) simultaneously, we can have that two curves
(θ∗,m∗)=(M−√a1a2β2(2M−a1a2β2)2K2,2K+(β2−d)a2−(d+θ∗)a1β2−d−θ∗), | (22) |
where
According the above analysis, the relative locations of
Therefore, when
(1a) when
(1b) when
(1c) when
Case 2.
For the equation (10), the positive root exists only when
Based on the definition of the curve
Case 3.
By using
According to the analysis of the curve
From the discussion above, three curves
C1={(θ,m)|θ=β2−d,0≤m<m5(θ)};C2={(θ,m)|β2−d<θ,0≤m=m5(θ)<m4(θ)};C3={(θ,m)|β2−d<θ,0≤m=m3(θ)≤m4(θ)}. |
which divide the whole region
Λ1={(θ,m)|0≤θ<β2−d,0≤m<m5(θ)};Λ2={(θ,m)|β2−d<θ,0≤m<m5(θ)};Λ3={(θ,m)|β2−d<θ,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
(I) If
(a) for
for
for
for
(b) for
(II) If
(a) for
for
(b) for
(III) If
(a) for
(b) for
(IV) If
Where when exist the corresponding equilibria are:
E1(P1,Z1)=(−b+√Δ2a,η(m−b+√Δ)(K−m+b−√Δ)(a1−b+√Δ)4a2(−b+√Δ)), | (23) |
E2(P2,Z2)=(−cb,−η(m−c)(K−m+c)(a1−c)b2c), | (24) |
E3(P3,Z3)=(−b−√Δ2a,−η(m−b−√Δ)(K−m+b+√Δ)(a1−b−√Δ)4a2(b+√Δ)). | (25) |
From Fig. 2, it is easy to see that the curve
According to the definition of the curve
f(K−m)=(a2+K−m)(da2+dK+θK)(β2(K−m)(a2+K)(a2+K−m)(da2+dK+θK)−1). | (26) |
Let
For
We summarize the above discussion into the following proposition.
Proposition 1. Suppose
(a) if
(b) if
(c) if
(d) if
Proposition 1 shows that system (8) will undergo a backward bifurcation which occurs at
Let
JE=(−η(a2+m+P∗)P∗h(P∗)−P∗(a2+m+P∗)(2aP∗+b)Z∗f(P∗)), | (27) |
where
h(P∗)=2(P∗)3+(a1+2m−K)(P∗)2+ma1(K−m). | (28) |
Theorem 4.1. The boundary equilibrium
Proof. One can verify that system (8) at
By using Proposition 1, system (8) has no positive equilibria if
In this subsection, the aim is to determine the stability and bifurcation of the existing positive equilibriums of system (8).
Theorem 4.2. Assume
Proof.From (27), then the Jacobian matrix of system (8) at
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
λ2+η(a2+m+P1)P1h(P1)λ+(2aP1+b)(a2+m+P1)P1Z1=0. | (30) |
When
When
If
If
dvdm|m=ˆm=(dvdP1dP1dm)|m=ˆm≠0, | (31) |
where
dvdP1=−η(a2+m+P1)(3P1+m+2a1−K)P1, | (32) |
dP1dm=am+(β2−d)a2+(d+θ)a1−√Δ2√Δ, | (33) |
which means that Hopf bifurcation can occur.
The proof is finished.
When
Theorem 4.3. Assume
According to Theorem 3.1, when
When
m=m∗=−(d+θ)a1−(β2−d)a2+2√β2a1a2θβ2−d−θ. | (34) |
And a straightforward calculation gives that one eigenvalue of Jacobian matrix of system (8) at
K=K∗=P∗+m∗+P∗(P∗+m∗)(a1+P∗)(P∗)2−m∗a1, | (35) |
In the following, we will consider the dynamical behavior of system (8) at the equilibrium
When
x=P−P∗andy=Z−Z∗ | (36) |
brings
{dxdt=−η(a2+m+P∗)P∗h(P∗)x−P∗(a2+m+P∗)y+ψ(x,y)dydt=ϕ(x,y), | (37) |
where
ψ(x,y)=η(K−2m−a1−3P∗)(a2+m+P∗)x2−(a2+m+2P∗)xy+η(K−3m−a1−a2−4P∗)x3−x2y−ηx4,ϕ(x,y)=Z∗ax2+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],
Substituting
ϕ(u(y),y)=(P∗)∗Z∗ab2h2(P∗)y2+o(y2), | (38) |
which shows that the equilibrium
Therefore, we have the following theorem.
Theorem 4.4. If
Based on the discussion in subsection 4.3, we have that system (8) at
Theorem 4.5. If
Proof. Similar to the process which leads to (37), then we have
{dxdt=−P∗(a2+m+P∗)y+η(K−2m−a1−3P∗)(a2+m+P∗)x2+R10(x,y),dydt=Z∗ax2+R20(x,y), | (39) |
where
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
a11=a2+m+2P∗P∗(a2+m+P∗)>0,a20=η(K−2m−a1−3P∗)(a2+m+P∗),b20=−P∗(a2+m+P∗)Z∗a>0. | (42) |
Making the near-identity transformation
x2=x1,y2=y1+a11x1y1+a20x21+R11(x1,y1), | (43) |
we have
{dx2dt=y2,dy2dt=b20x22+a11y2y2−a20x221+a11x2+2a20x2y2+R22(x2,y2), | (44) |
where
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
By a time reparameterization
u=x2andv=(1+a11x2)y2, | (47) |
we obtain
{dudτ=v,dvdτ=b20u2+2a20uv+R24(u,v), | (48) |
where
According to (47) and
Hence, the equilibrium
Next, we study system (8) by perturbing the parameters
{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+ϵ1−m∗−ϵ2−P)(a1+P)−PZ],dZdt=β2P(a2+m∗+ϵ2+P)Z−d(a1+P)(a2+m∗+ϵ2+P)Z−θ(P+m∗+ϵ2)(a1+P)Z, | (50) |
for sufficiently small
System (50) has a cusp point
{ˉP=P−P∗,ˉZ=Z−Z∗, | (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(K−2m∗−2P∗−ϵ2)+ϵ1(m∗+ϵ2+P∗)](a1+P∗)(a2+m∗+ϵ2+P∗),p10=η{(a2+m∗+P∗+ϵ2)[ϵ2(K−2m∗−4P∗−2a1−ϵ2)+ϵ1(a1+2P∗+m∗+ϵ2)]+(a1+P∗)[ϵ2(K−2m∗−2P∗−ϵ2)+ϵ1(m∗+ϵ2+P∗)]},p01=−P∗(a2+m∗+P∗+ϵ2),p11=−(a2+m∗+2P∗+ϵ2),p20=η{(K−2m∗−a1−3P∗−2ϵ2)(a2+m∗+P∗+ϵ2)+ϵ2(K−2m∗−4P∗−2a1−ϵ2)+ϵ1(a1+a2+2P∗+2m∗+2ϵ2)},p21=−1,p30=η{(K−3m∗−4P∗−a1−a2−3ϵ2)+ϵ1},p40=−η,q00=ϵ2((β2−d−θ)P∗−(d+θ)a1)Z∗,q10=ϵ2(β2−d−θ)Z∗,q01=ϵ2((β2−d−θ)P∗−(d+θ)a1),q11=ϵ2(β2−d−θ),q20=(β2−d−θ)Z∗,q21=β2−d−θ. | (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=p01q00−p00q01,f10=p01q10+p11q00+p00q11−p10q01,f01=p10+q01−p00p11p01,f11=2p20+q11+p00p211−p10p01p11−2p00p01p21p201,f20=p01q20+p11q10+p21q00−p20q01−p10q11−p00q21+p11(p10q01+p00q11)p01−p00(p211−p21)q01p201,f02=p11. | (56) |
By rescaling the time variable
t=(1−f02˜Z)τ, | (57) |
system (55) becomes
{d˜Pdτ=(1−f02˜Z)˜Z,d˜Zdτ=(1−f02˜Z)(f00+f10˜P+f01˜Z+f11˜P˜Z+f20˜P2+f02˜Z2+O(‖(˜P,˜Z)‖3)). | (58) |
Furthermore, let
{ˆP=˜P,ˆZ=(1−f02˜Z)˜Z, | (59) |
then we have
{dˆPdτ=ˆZ,dˆZdτ=f00+(f10−2f00f02)ˆP+f01ˆZ+(f20+f00f202−2f10f02)ˆP2+(f11−f01f02)ˆPˆZ+O(‖(ˆP,ˆZ)‖3). | (60) |
By (53), we can compute
limϵ1→0,ϵ2→0(f11−f01f02)=η(K−2m∗−a1−3P∗)(a2+m∗+P∗)≠0, | (61) |
thus setting
{dˆPdτ=ˆZ,dˆZdτ=g00+g10ˆP+g20ˆP2+g11ˆPˆZ+O(‖(ˆP,ˆZ)‖3). | (62) |
where
g00=f00−f01(f10−2f00f02)f11−f01f02+f201(f20+f00f202−2f10f02)(f11−f01f02)2,g10=f10−2f00f02−2f01(f20+f00f202−2f10f02)f11−f01f02,g20=f20+f00f202−2f10f02,g11=f11−f01f02. | (63) |
From (53) and (63), we are able to have
limϵ1→0,ϵ2→0g11=η(K−2m∗−a1−3P∗)(a2+m∗+P∗)≠0, | (64) |
and
limϵ1→0,ϵ2→0g20=−P∗(a2+m∗+P∗)(β2−d−θ)Z∗≠0, | (65) |
therefore, by performing the change of new variables:
{dξ1dt=ξ2,dξ2dt=g00g411g320+g10g211g220ξ1+ξ21+ξ1ξ2+O(‖(^ξ1,^ξ2)‖3). | (66) |
Corollary 1. If
From Theorem 4.2 and 4.3, the positive equilibrium
To determine the stability of the positive equilibrium
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=η(K−2m−a1−3˜P)(a2+m+˜P),a21=−1,a30=η(K−3m−a1−a2−4˜P),b10=(2a˜P+b)˜Z,b11=(2a˜P+b),b20=˜Za,b21=a. |
Applying the formula of the Liapunov number
σ=−3π2Δ32(−2a01a20b20−a11a20b10+3a01a30b10), | (69) |
where it follows from the proof of Theorem 4.2 that
Theorem 4.6. For system (8), fix all parameters except
(i) a supercritical Hopf bifurcation occurs at
(ii) a subcritical Hopf bifurcation occurs at
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
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
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
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
At the same time, the equilibrium level of phytoplankton due to the refuge becomes large as the value of
From Fig. 7, we observe that the related dynamical behavior of system (3) also relies on the refuge capacity
With the purpose of exploring the influences of toxin and refuge on the dynamical behavior of our systems in combination, by choosing
Beside, in order to study the complicated dynamics of system (8) at the critical positive equilibrium
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. |
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 |
Par. | Description | Value | Unit | Reference. |
Growth rate of phytoplankton | 0.2 | 0.07-0.28 [16] | ||
| Environmental carrying capacity | 50 | 108 [15] | |
Refuge capacity | Par. | Defaulted | ||
Predation rate of zooplankton | 1 | | 0.6-1.4[16] | |
Growth efficiency of zooplankton | 0.15 | 0.2-0.5[16] | ||
Mortality rate of zooplankton | 0.003 | 0.015-0.15[16] | ||
Half saturation constant | 3 | Defaulted | ||
| Half saturation constant | 5.7 | | 5.7[15] |
Toxin production rate | Par. | Defaulted |
Par. | Description | Value | Unit | Reference. |
Growth rate of phytoplankton | 0.2 | 0.07-0.28 [16] | ||
| Environmental carrying capacity | 50 | 108 [15] | |
Refuge capacity | Par. | Defaulted | ||
Predation rate of zooplankton | 1 | | 0.6-1.4[16] | |
Growth efficiency of zooplankton | 0.15 | 0.2-0.5[16] | ||
Mortality rate of zooplankton | 0.003 | 0.015-0.15[16] | ||
Half saturation constant | 3 | Defaulted | ||
| Half saturation constant | 5.7 | | 5.7[15] |
Toxin production rate | Par. | Defaulted |