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$\acute{a}$lez-Olivares et al.[22] studied the dynamic consequences of a continuous models with the role of prey's refuge, they analytically proved that the local stability of equilibria and existence of limit cycle. In terms of the work of Gonz$\acute{a}$lez-Olivares, Chen et al. [23] further gave a complete dynamical analysis for the same model in [22] by proving the global stability of the positive equilibrium and the uniqueness of limit cycle when exists. Kumar [24] dealt with a prey-predator model with Holling type Ⅱ 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],

$\left\{\begin{aligned} \frac{dP}{dt}&= rP(1-\frac{P}{K})-\beta_1\phi_1(P)Z,\\ \frac{dZ}{dt}&=\beta_2\phi_1(P)Z-dZ,\\ \end{aligned} \right.$

(1)

where $P$ and $Z$ denote the population size of phytoplankton and zooplankton respectively; $\phi_1(P)$ represents the predatory response function; $r$ is the intrinsic growth rate of phytoplankton; $K$ represents the environmental carrying capacity of phytoplankton; $\beta_1$ is the specific predation rate of zooplankton; $\beta_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,

$\left\{\begin{aligned} \frac{dP}{dt}&= rP(1-\frac{P}{K})-\beta_1\phi_1(P)Z,\\ \frac{dZ}{dt}&=\beta_2\phi_1(P)Z-dZ-\theta\psi(P)Z,\\ \end{aligned} \right.$

(2)

where $\theta$ denotes the rate of toxin production per phytoplankton species; $\psi(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, $\psi(P)=\frac{P}{a_2+P}$, to model the distribution of toxic substances, where $a_2$ is the half saturation constant.

The refuge is favorable for phytoplankton growth by preventing the grazing of zooplankton. For the effect of refuges, one naturally subdivides the total phytoplankton into two groups, one portion is in the refuge and the other portion is outside the refuge. Usually, one would consider that the population size of phytoplankton in refuge as one variables depends on different growth and predation pressure. However, the effectiveness of refuges for phytoplankton depends on various factors, for example, temperature, radiation and winds mainly affect the process of stratification which offers a refuge for phytoplankton in some ways [31]. Due to the complexity of the biological process, we restrict ourselves to study a special case that the refuge capacity denoted by $m$ remains a constant and is taken as a parameter in the model analysis. Besides, we assume that phytoplankton can be able to take full advantage of refuges. Therefore, we consider that the predation of zooplankton on phytoplankton depends on the population size of phytoplankton outside refuges expressed by $P-m$ and the predation of zooplankton on unprotected phytoplankton increases and can approach its maximum biological limit with the saturation of phytoplankton in water body. Motivated by the study in [22], we employ Holling type Ⅱ functional response, i.e., $\phi_1(P)=\frac{P-m}{a_1+(P-m)}$, to describe the predation of zooplankton on unprotected phytoplankton, where $a_1$ is the the half saturation constant.

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

$\left\{\begin{aligned} \frac{dP}{dt}&= rP(1-\frac{P}{K})-\frac{\beta_1(P-m)}{a_1+(P-m)}Z,\\ \frac{dZ}{dt}&=\frac{\beta_2(P-m)}{a_1+(P-m)}Z-dZ-\frac{\theta P}{a_2 +P}Z,\\ \end{aligned} \right.$

(3)

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

It should be noted from the second equation of model (3) that if $\beta_2\leq d$, the population size of zooplankton will continue to drop or remain unchanged over time, which means that the dynamical behaviors of zooplankton can not be described objectively in the actual situation when the conversion rate of phytoplankton is no more than the mortality rate of zooplankton. Therefore, it makes sense to consider only the case $\beta_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 $\displaystyle \liminf_{t \rightarrow \infty}P(t)\geq 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

$\limsup\limits_{t \rightarrow \infty}P(t)\leq K,$

(4)

hence, for sufficiently small $\varepsilon > 0$, there is a sufficiently large $T > 0$ such that $P(t) <K + \varepsilon$ for all $t>T$.

Define the function $w(P, Z)=cP(t)+Z(t)$, where $c=\frac{\beta_2}{\beta_1}$, then the time derivative of $w(P, Z)$ along the solutions of system (3) is given by

$\begin{align}\label{a2} \begin{split} \frac{dw}{dt}&=c \frac{dP}{dt}+\frac{dZ}{dt} \\ &=crP(1-\frac{P}{K})-dZ-\frac{\theta P}{\gamma +P}Z\\ &\leq crP-dZ\\ &\leq -(dZ+rcP)+2crP\\ &\leq -min\{r,d\}w+2cr(K+\varepsilon). \end{split} \end{align}$

(5)

And then we have

$0<w(P,Z)\leq \frac{2cr(K+\varepsilon)}{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 $t$, $w(t)\leq \frac{2cr(K+\varepsilon)}{min\{r,d\}}$ holds. The proof is finished.

This lemma implies that the set $\Omega=\big\{(P, Z)\in R^2_+|0\leq cP+Z\leq\frac{2cr(K+\varepsilon)}{min\{r,d\}}\big\}$ is an invariant set. For practical biological significance, we will focus our discussion about system (3) in $\Omega$.

---

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:

$\hat{P}=P-m, \hat{Z}=\beta_1 Z, \tau=\frac{1}{(a_2+P)(a_1+(P-m))}t,$

(7)

if we denote $\eta=\frac{r}{K}$, then system (3) becomes the following form (here we still use symbols $P, Z, t$ instead of $\hat{P}, \hat{Z}, \tau$ respectively)

$\left\{\begin{aligned} \frac{dP}{dt}&=(a_2+m+P)\big{(}\eta(P+m)(K-m-P)(a_1+P) -PZ\big{)},\\ \frac{dZ}{dt}&=\beta_2 P(a_2+m+P)Z-d(a_1+P)(a_2+m+P)Z-\theta (P+m)(a_1+P)Z. \end{aligned} \right.$

(8)

Note that $\Omega_1=\big\{(P, Z)|K\geq P(t)>m, Z(t)>0$ for all $t>0\big\}$ for system (3) equivalently transforms to $\Omega_2=\big\{(P, Z)|K-m\geq P(t)>0, Z(t)>0$ for all $t>0\big\}$, which means that we will study system (8) in $\Omega_2$.

Obviously, $E_0(K-m, 0)$ is always a boundary equilibrium of system (8).

Suppose $(P^*,Z^*)$ is a positive equilibrium of the system (8), by putting the right-hand side of system (8) to zero, except for $E_0$, we have

$Z^*=\frac{\eta(P^*+m)(K-m-P^*)(a_1+P^*)}{P^*},$

(9)

where $P^*$ satisfies the quadratic equation

$f(P)=aP^2+bP+c=0,$

(10)

with

$\begin{aligned} a&=\beta_2-d-\theta,\\ b&=am+(\beta_2-d)a_2-(d+\theta)a_1,\\ c&=-a_1(da_2+dm+\theta m)<0. \end{aligned}$

(11)

Of special note here is that $P^*$ also satisfies $0<P^*<K-m$. The interior equilibria $(P^*,Z^*)$ of system (8) in $\Omega_2$ exist if and only if the equation (10) has positive root(s) in the interval $(0, K-m)$.

In order to explore the effect of toxin and refuge on the interactions between phytoplankton and zooplankton, we will only choose $\theta$ 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 $\Lambda=\{(\theta, m)|0<\theta, 0<m<K\}$ of plane $(\theta, m)$.

For the convenience of expression, we first define

$\begin{aligned} &d_0=\beta_2(1-\frac{a_1}{a_2}),~~ d_1=\frac{K^2\beta_2}{(K+a_1)(K+a_2)}(1-\frac{a_1}{a_2}),\\ &d_2=\frac{K\beta_2}{K+a_1}(1-\frac{a_1}{a_2}), ~~ d_3=\frac{K\beta_2}{K+a_1}, \end{aligned}$

(12)

where $d_1<d_2<d_3<d_0$ and the meaning of this definition will become clear in the following statement. Since $f(0)=c<0$, $f(P)$ will have positive root(s) in three cases.

Case 1. $a<0$, i.e., $\theta>\beta_2-d$

In this case, we will discuss the existence conditions of the positive root(s) of the equation (10) in details. It is easy to have that the equation (10) has positive roots only and only if any one of three cases is satisfied in Fig. 1(a)-(c).

One can verify that $-b/(2a)=0$ defines a curve $\Gamma_1$,

$\Gamma_1: m_1(\theta)=\frac{(d+\theta)a_1-(\beta_2-d)a_2}{\beta_2-d-\theta}.$

(13)

By $c<0$, if the positive solution of the equation (10) exists for $\theta>\beta_2-d$, then $-b/(2a)>0$ must hold, which implies $0<m<m_1(\theta)$. From $\theta>\beta_2-d$ and (13), we can get that if $d<d_0$ with $a_2>a_1$, the positive root of (10) might exist. Or else, the equation (10) has no any positive root. Therefore, the existence conditions of the positive root of the equation (10) will be further investigated in the following only when $d<d_0$.

With $d<d_0$, along the curve $\Gamma_1$, the derivative of $m_1(\theta)$ with respect to $\theta$ is

$m_1'(\theta)=-\frac{(\beta_2-d)a_2-\beta_2a_1}{(\beta_2-d-\theta)^2}<0.$

(14)

In addition, it is easy to have that $m_1(\theta)$ with respect to $\theta$ is monotone decreasing. Besides, $\theta=\beta_2-d$ is a vertical asymptote of $m_1(\theta)$ and $m_1(\infty)\rightarrow -a_1$. By calculating, the curve $\Gamma_1$ travels through two points $(\frac{(\beta_2-d)a_2-da_1}{a_1},0)$ and $(\frac{(\beta_2-d)(K+a_2)-da_1}{K+a_1}, K)$ in the region $\Lambda$.

$\Delta=0$ can define two curves. By writing the discriminant $\Delta=b^2-4ac$ as a function of $\theta$ and $m$, a straightforward calculation leads to

$\begin{align*} \Delta(m,\theta)& =b^2-4ac\\ &=(am+(\beta_2-d)a_2+(d+\theta)a_1)^2-4\beta_2a_1a_2\theta, \end{align*}$

then solving the equation $\Delta(m,\theta)=0$ for $m$, one can get two curves

$\Gamma_2: m_2(\theta)=\frac{-(d+\theta)a_1-(\beta_2-d)a_2-2\sqrt{\beta_2a_1a_2\theta}}{\beta_2-d-\theta},$

(15)

$\Gamma_3: m_3(\theta)=\frac{-(d+\theta)a_1-(\beta_2-d)a_2+2\sqrt{\beta_2a_1a_2\theta}}{\beta_2-d-\theta}.$

(16)

In order to ensure the existence of the root of the equation (10), we first need $\Delta\geq0$. From the above discussion, it is not difficult to obtain that $\Delta>0$ if and only if $m>m_2(\theta)$ or $m<m_3(\theta)$ and $\Delta=0$ if and only if $m=m_2(\theta)=m_3(\theta)$. Combined with the above discussion on the curve $\Gamma_1$, by comparing, we have $m_3(\theta)<m_1(\theta)<m_2(\theta)$. Thus, (10) may have two positive roots if $0<m<m_3(\theta)$ and an unique positive root if $0<m=m_3(\theta)$, otherwise, there is no root. With these, one can see that the curve $\Gamma_3$ pays a key role in the existence of the positive root(s) of (10).

For $\Gamma_3$, $\theta=\frac{(\beta_2-d)a_2-da_1\pm 2\sqrt{\beta_2d a_1(a_2-a_1)}}{a_1}$ when $m_3=0$, $\theta=\beta_2-d$ is a vertical asymptote of $m_3(\theta)$, and the derivative of $m_3(\theta)$ with respect to $\theta$ is

$m_3'(\theta)=\frac{(\frac{\beta_2a_1a_2}{\sqrt{\beta_2a_1a_2\theta}}-a_1)(\beta_2-d-\theta) -((d+\theta)a_1+(\beta_2-d)a_2-2\sqrt{\beta_2a_1a_2\theta})}{(\beta_2-d-\theta)^2}<0.$

(17)

Solving $K-m=-b/(2a)$, we get a curve $\Gamma_4$,

$\Gamma_4: m_4(\theta)=2K+\frac{(\beta_2-d)a_2-(d+\theta)a_1}{\beta_2-d-\theta}=2K-m_1(\theta),$

(18)

where $K-m>-b/(2a)$ is equivalent to $0<m<m_4(\theta)$.

From (18), $\theta=\beta_2-d$ is a vertical asymptote of $m_4(\theta)$, $m_4(0)=2K+\frac{(\beta_2-d)a_2-da_1}{\beta_2-d}$, $m_4(\infty)\rightarrow 2K+a_1$, and $\theta=\frac{(\beta_2-d)(2K+a_2)-da_1}{2K+a_1}$ when $m_4=0$. Furthermore, the derivative of $m_4(\theta)$ with respect to $\theta$ is given by the following form

$m_4'(\theta)=\frac{(\beta_2-d)a_2-\beta_2a_1}{(\beta_2-d-\theta)^2}=-m_1'(\theta)>0.$

(19)

Finally, $f(K-m)=0$ defines a curve $\Gamma_5$,

$\Gamma_5: m_5(\theta)=K-\frac{(Kd+K\theta+da_2)a_1}{(\beta_2-d)(K+a_2)-K\theta},$

(20)

which implies that $f(K-m)<0$ if and only if $m>m_5(\theta)$ and $f(K-m)>0$ if and only if $0<m<m_5(\theta)$.

From (20), we have $\theta=(\beta_2-d)(1+\frac{a_2}{K})$ is a vertical asymptote of $m_5(\theta)$, $m_5(0)=K-\frac{da_1}{\beta_2-d}<K$, $m_5(\infty)\rightarrow K+a_1$ and $\theta=\frac{(K+a_2)(\beta_2K-dK-da_1)}{K(K+a_1)}$ when $m_5=0$. The corresponding derivative of $m_5(\theta)$ with respect to $\theta$ is

$m_5'(\theta)=-\frac{Ka_1\beta_2(K+a_2)}{((\beta_2-d)(K+a_2)-K\theta)^2}<0.$

(21)

Since $\theta>\beta_2-d$, we get from the characters of the curve $\Gamma_5$ that equation (10) might have positive root(s) in the region $\Lambda$ when $m_5(\beta_2-d)=K-\frac{(K\beta_2+da_2)a_1}{(\beta_2-d)a_2}>0$ i.e., $d<d_2$.

When $d<d_2$, we have that only three curves $\Gamma_3$, $\Gamma_4$ and $\Gamma_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 $\Gamma_4$ and $\Gamma_5$ uniquely intersect in the region $\Lambda$ when $d<d_1$. Through further validation, the intersection point also satisfies the equation of $\Gamma_3$. That is, when $d<d_1$, three curves have an unique common intersection point $(\theta_*, m_*)$. Otherwise, three curves can not intersect at the same point in the region $\Lambda$. Let $(\theta_*, m_*)$ be the intersection point, by computing, then we get

$(\theta_*, m_*)=(\frac{M-\sqrt{a_1a_2\beta_2(2M-a_1a_2\beta_2)}}{2K^2}, 2K+\frac{(\beta_2-d)a_2-(d+\theta_*)a_1}{\beta_2-d-\theta_*}),$

(22)

where $M=2K(\beta_2-d)(K+a_2)+a_1a_2\beta_2$.

According the above analysis, the relative locations of $\Gamma_3$, $\Gamma_4$ and $\Gamma_5$ can be approximately depicted in the region $\Lambda$ (see Fig. 2(a) and (b)).

Therefore, when $\theta>\beta_2-d$, we can conclude that

(1a) when $d_1\leq d<d_2$, then system (8) has an unique positive equilibrium if $0<m\leq m_5(\theta)$.

(1b) when $d<d_1$, then system (8) has an unique positive equilibrium if $0<m=m_3(\theta)<m_4(\theta)$ or if $0<m<m_5(\theta)$ or if $0<m=m_5(\theta)\leq m_4(\theta)$ and system (8) has two positive equilibrium if $max\{0, m_5(\theta)\}<m<min\{m_3(\theta), m_4(\theta)\}$, or else, system (8) does not have any positive equilibrium.

(1c) when $d\geq d_2$, then system (8) does not have any positive equilibrium.

Case 2. $a=0$, i.e., $\theta=\beta_2-d$.

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

Based on the definition of the curve $\Gamma_5$ in the region $\Lambda$, $f(K-m)>0$ is equivalent to $0<m<m_5(\beta_2-d)$, where $m_5(\beta_2-d)>0$ can be equivalently transformed into $d<d_2$. As a consequence, we reach a conclusion that the positive equilibrium is existent and unique for system (8) only when $0<m<m_5(\beta_2-d)$ with $d<d_2$.

Case 3. $a>0$, i.e., $\theta<\beta_2-d$

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

According to the analysis of the curve $\Gamma_5$, $f(K-m)>0$ is also equivalent to $0<m<m_5(\theta)$ in the case. By the monotonicity of $m_5(\theta)$, then $m_5(0)>0$, i.e., $d<d_3$, must be satisfied, which means there is an unique positive equilibrium for system (8) only when $0<m<m_5(\theta)$ with $d<d_3$.

From the discussion above, three curves $\Gamma_3$, $\Gamma_4$ and $\Gamma_5$ are important to determine the existence and the numbers of the positive equilibrium points of system (8). Here we define $C_1$, $C_2$ and $C_1$ as

$\begin{align*} C_1&=\big\{(\theta, m)| \theta=\beta_2-d, 0\leq m<m_5(\theta)\big\};\\ C_2&=\big\{(\theta, m)| \beta_2-d<\theta, 0\leq m=m_5(\theta)<m_4(\theta)\big\};\\ C_3&=\big\{(\theta, m)|\beta_2-d<\theta, 0\leq m=m_3(\theta)\leq m_4(\theta)\big\}. \end{align*}$

which divide the whole region $\Lambda$ into the following subregions where the positive equilibrium points exist:

$\begin{align*} \Lambda_1&=\big\{(\theta, m)| 0\leq\theta<\beta_2-d, 0\leq m<m_5(\theta)\big\};\\ \Lambda_2&=\big\{(\theta, m)| \beta_2-d<\theta, 0\leq m<m_5(\theta)\big\};\\ \Lambda_3&=\big\{(\theta, m)| \beta_2-d<\theta, m_5(\theta)<m<m_3(\theta), m<m_4(\theta)\big\}. \end{align*}$

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 $m_3(\theta), m_4(\theta)$ and $m_5(\theta)$ defined as above.

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

(a) for $(\theta, m)\in \Lambda_1 \bigcup \Lambda_2 \bigcup C_2$, system (8) has an unique positive equilibrium $E_1(P_1, Z_1)$;

for $(\theta, m)\in C_1$, system (8) has an unique positive equilibrium $E_2(P_2, Z_2)$;

for $(\theta, m)\in \Lambda_3$, system (8) has two positive equilibria $E_1(P_1, Z_1)$ and $E_3(P_3, Z_3)$;

for $(\theta, m)\in C_3$, $E_1(P_1, Z_1)$ and $E_3(P_3, Z_3)$ coalesce at an unique positive equilibrium of multiplicity $2$.

(b) for $(\theta, m)\in \Lambda/(\Lambda_1 \bigcup \Lambda_2 \bigcup \Lambda_3 \bigcup C_1\bigcup C_2 \bigcup C_3)$, system (8) has none positive equilibrium.

(II) If $d_1\leq d<d_2$, as shown in Fig. 2(b), then

(a) for $(\theta, m)\in C_1$, system (8) has an unique positive equilibrium $E_2(P_2, Z_2)$;

for $(\theta, m)\in \Lambda_1 \bigcup \Lambda_2$, system (8) has an unique positive equilibrium $E_1(P_1, Z_1)$;

(b) for $(\theta, m)\in \Lambda/(\Lambda_1 \bigcup \Lambda_2 \bigcup C_1)$, system (8) has none positive equilibrium.

(III) If $d_2\leq d<d_3$, as shown in Fig. 2(c), then

(a) for $(\theta, m)\in \Lambda_1$, system (8) has an unique positive equilibrium $E_1(P_1, Z_1)$;

(b) for $(\theta, m)\in \Lambda/\Lambda_1$, system (8) has none positive equilibrium.

(IV) If $d\geq d_3$, system (8) has none positive equilibrium.

Where when exist the corresponding equilibria are:

$E_1(P_1, Z_1)=(\frac{-b+\sqrt{\Delta}}{2a},\frac{\eta(m-b+\sqrt{\Delta})(K-m+b-\sqrt{\Delta})(a_1-b+\sqrt{\Delta})}{4a^2(-b+\sqrt{\Delta})}),$

(23)

$E_2(P_2, Z_2)=(\frac{-c}{b},-\frac{\eta(m-c)(K-m+c)(a_1-c)}{b^2c}),$

(24)

$E_3(P_3, Z_3)=(\frac{-b-\sqrt{\Delta}}{2a},-\frac{\eta(m-b-\sqrt{\Delta})(K-m+b+\sqrt{\Delta})(a_1-b-\sqrt{\Delta})}{4a^2(b+\sqrt{\Delta})}).$

(25)

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

According to the definition of the curve $\Gamma_5$, by simply rearranging the expression of $f(K-m)$, we have

$f(K-m)=(a_2+K-m)(da_2+dK+\theta K)\big{(}\frac{\beta_2(K-m)(a_2+K)}{(a_2+K-m)(da_2+dK+\theta K)}-1\big{)}.$

(26)

Let $R_0=\frac{\beta_2(K-m)(a_2+K)}{(a_1+K-m)(da_2+dK+\theta K)}$, the ratio of growth to mortality of zooplankton at the critical situation. Based on the preceding assumption $0<m<K$, then $f(K-m)>0$ if and only if $R_0>1$ and $f(K-m)<0$ if and only if $R_0<1$.

For $d<d_1$ and $m\in(0, m_*)$, the curve $\Gamma_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 $\Gamma_3$ is defined by $\Delta=0$. Let $\tilde{R}_0=R_0|_{\Delta=0}$, then $\tilde{R}_0=1+\frac{(2a(K-m)+b)^2}{4a(a_1+K-m)(da_2+dK+\theta K)}$ and $0<\tilde{R}_0<1$ due to $\theta>\beta_2-d$.

We summarize the above discussion into the following proposition.

Proposition 1. Suppose $d<d_1$ and $m\in(0, m_*)$, then

(a) if $R_0\geq1$, system (8) has an unique positive equilibrium;

(b) if $\tilde{R}_0<R_0<1$, system (8) has two positive equilibria;

(c) if $R_0=\tilde{R}_0<1$, two positive equilibria coalesce at an unique equilibrium of multiplicity $2$;

(d) if $R_0<\tilde{R}_0<1$, system (8) has none positive equilibria.

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

---

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

$\begin{matrix} J_{E}= \begin{pmatrix} -\frac{\eta(a_2+m+P^*)}{P^*}h(P^*)&-P^*(a_2+m+P^*)\\ (2aP^*+b)Z^*& f(P^*) \end{pmatrix} \end{matrix},$

(27)

where

$\label{a10} h(P^*)=2(P^*)^3+(a_1+2m-K)(P^*)^2+ma_1(K-m).$

(28)

---

4.1. The boundary equilibrium $E_0$

Theorem 4.1. The boundary equilibrium $E_0$ is unstable if $R_0>1$, asymptotically stable if $R_0<1$ and globally asymptotically stable if $R_0<\tilde{R}_0<1$.

Proof. One can verify that system (8) at $E_0$ has an eigenvalue, $-\eta K(K+a_2)(a_1+(K-m))<0$. The other eigenvalue is $f(K-m)$. The results follow from the fact that $f(K-m)>0$ is equivalent to $R_0>1$.

By using Proposition 1, system (8) has no positive equilibria if $R_0<\tilde{R}_0<1$. Since $E_0$ is local asymptotically stable, it follows from Poincar$\acute{e}$-Bendixson theorem [33] that $E_0$ is also globally asymptotically stable. We finish this proof.

---

4.2. The positive equilibrium $E_1, E_2$ and $E_3$ for $(\theta, m)\in \Lambda_1$, $\Lambda_2$, $\Lambda_3$, $C_2$ and $C_1$

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

Theorem 4.2. Assume $(\theta, m)\in \Lambda_1 \bigcup \Lambda_2 \bigcup \Lambda_3 \bigcup C_2$, then the positive equilibrium $E_1$ of system (8) is asymptotically stable when $h(P_1)>0$ and unstable when $h(P_1)<0$, system (8) will experience a Hopf-bifurcation at $E_1$ when $h(P_1)=0$. Moreover $E_3$ when exists is always unstable.

Proof.From (27), then the Jacobian matrix of system (8) at $E_1(P_1, Z_1)$ is

$\begin{matrix} J_{E_1}= \begin{pmatrix} -\frac{\eta(a_2+m+P_1)}{P_1}h(P_1)& -P_1(a_2+m+P_1)\\ (2aP_1+b)Z_1& 0 \end{pmatrix} \end{matrix}.$

(29)

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

$\lambda^2+\frac{\eta(a_2+m+P_1)}{P_1}h(P_1)\lambda+(2aP_1+b)(a_2+m+P_1)P_1Z_1=0.$

(30)

When $(\theta, m)\in \Lambda_1$ or $\Lambda_2$ or $C_2$, one can obtain from Theorem 3.1 that the positive equilibrium $E_1$ is unique for system (8). For $E_1$, by Theorem (2.1), the derivative of $f(P)$ with respect to $P_1$ satisfies $f'(P_1)=(2aP_1+b)Z_1>0$, which is easy to be obtained by Fig. 1, then $(2aP_1+b)(a_2+m+P_1)P_1Z_1>0$. By using the Routh-Hurwitz criterion, we have the real part of two roots of (30) are negative if $h(P_1)>0$. Therefore, $E_1$ is asymptotically stable if $h(P_1)>0$.

When $(\theta, m)\in \Lambda_3$, system (8) has another positive equilibrium $E_3$ except for $E_1$, where $P_1<P_3$. In this case, the discussion about the stability of system (8) at $E_1$ is the same as the above. For $E_3$, the Jacobin of system (8) at $E_3$ can be given by (29) where $P_1$ is replaced by $P_3$. Note that, however, the derivative of $f(P)$ with respect to $P_3$ satisfies $f'(P_3)=(2aP_3+b)Z_3<0$, which can be obtained easily by Fig. 1(c), then $(2aP_3+b)(a_2+m+P_3)P_3Z_3<0$, which implies that $E_3$ is always unstable.

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

If $h(P_1)=0$, then (30) has a pair of pure imaginary complex conjugate roots. To discuss the effect of refuges, we here select the parameter $m$ as a bifurcation parameter. Let $\lambda=v(m)\pm i\omega(m)$ be the complex conjugate roots and there exists a critical value $\hat{m}$ satisfying $H(P_1(\hat{m}))=0$, then the transversality condition satisfies

$\frac{dv}{dm}|_{m=\hat{m}}=\big{(}\frac{dv}{dP_1}\frac{dP_1}{dm}\big{)}|_{m=\hat{m}}\neq0,$

(31)

where

$\frac{dv}{dP_1}=\frac{-\eta(a_2+m+P_1)(3P_1+m+2a_1-K)}{P_1},$

(32)

$\frac{dP_1}{dm}=\frac{am+(\beta_2-d)a_2+(d+\theta)a_1-\sqrt{\Delta}}{2\sqrt{\Delta}},$

(33)

which means that Hopf bifurcation can occur.

The proof is finished.

When $(\theta, m)\in C_1$, we have a similar conclusion for $E_2$.

Theorem 4.3. Assume $(\theta, m)\in C_1$, then the positive equilibrium $E_2$ of system (8) is asymptotically stable when $h(P_2)>0$; $E_2$ is unstable when $h(P_2)<0$, and Hopf bifurcation occurs if $h(P_2)=0$.

---

4.3. Saddle-node bifurcation for $(\theta, m)\in C_3$

According to Theorem 3.1, when $(\theta, m)\in C_3$, two positive equilibria $E_1$ and $E_3$ coalesce at an unique positive equilibrium denoted by $E^*(P^*,Z^*)$ with $P^*=-\frac{b}{2a}$ and $Z^*=-\frac{\eta(m^*-b)(K-m^*+b)(a_1-b)}{4a^2b}$, where $a=\beta_2-d-\theta, b=am^*+(\beta_2-d)a_2-(d+\theta)a_1.$

When $(\theta, m)\in C_3$, we have from (16) that $m$ satisfies

$m=m^*=\frac{-(d+\theta)a_1-(\beta_2-d)a_2+2\sqrt{\beta_2a_1a_2\theta}}{\beta_2-d-\theta}.$

(34)

And a straightforward calculation gives that one eigenvalue of Jacobian matrix of system (8) at $E^*$ is zero, the other is $\frac{-\eta(a_2+m+P^*)}{P^*}h(P^*)$. By (28), we have that $h(P^*)=0$ is equivalent to

$K=K^*=P^*+m^*+\frac{P^*(P^*+m^*)(a_1+P^*)}{(P^*)^2-m^*a_1},$

(35)

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

When $K\neq K^*$, the relevant positive equilibrium $E^*$ is a Lyapunov-type equilibrium. The transformation

$x=P-P^* \text{and} y=Z-Z^*$

(36)

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

$\left\{\begin{aligned} \frac{dx}{dt}&= -\frac{\eta(a_2+m+P^*)}{P^*}h(P^*)x-P^*(a_2+m+P^*)y+\psi(x,y)\\ \frac{dy}{dt}&=\phi(x,y), \end{aligned} \right.$

(37)

where

$\begin{align} &\psi(x,y)=\eta(K-2m-a_1-3P^*)(a_2+m+P^*)x^2-(a_2+m+2P^*)x y&\\&+\eta(K-3m-a_1-a_2-4P^*)x^3-x^2y -\eta x^4,& \\ &\phi(x,y)=Z^*ax^2+ax^2y.& \end{align}$

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)=-\frac{(P^*)^2}{\eta h(P^*)}y+\cdots$ with $u(0)=0$.

Substituting $x=u(y)$ into $\phi(x,y)$, we have

$\phi(u(y),y)=\frac{(P^*)^*Z^*a}{b^2h^2(P^*)}y^2+o(y^2),$

(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 $K\neq K^*$, the equilibrium $E^*$ of system (8) is a stable saddle-node, where $m^*$ and $K^*$ are given in (34) and (35) respectively.

---

4.4. Bogdanov-Takens bifurcation for $(\theta, m)\in C_3$

Based on the discussion in subsection 4.3, we have that system (8) at $E^*$ has a zero eigenvalues with multiplicity 2 when $K=K^*$, which suggests that the equilibrium $E^*$ may undergo a Bogdanov-Takens bifurcation. Through the further analysis, the following theorem will be given and proved.

Theorem 4.5. If $m=m^*$, $K= K^*$ and $K^*\neq 3P^*+a_1+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

$\left\{\begin{aligned} \frac{dx}{dt}&=-P^*(a_2+m+P^*)y+\eta(K-2m-a_1-3P^*)(a_2+m+P^*)x^2\\ &+R_{10}(x,y),\\ \frac{dy}{dt}&=Z^*ax^2+R_{20}(x,y), \end{aligned} \right.$

(39)

where $R_{i0}(i=1,2)$ is $C^{\infty}$ in $(x, y)$ and $R_{i0}=O(|(x,y)|^3)$.

Introducing new variables

$x_1=x \text{and} y_1=-P^*(a_2+m+P^*)y,$

(40)

system (39) becomes

$\left\{\begin{aligned} \frac{dx_1}{dt}&=y_1+a_{11}x_1y_1+a_{20}x^2_1+R_{11}(x_1,y_1),\\ \frac{dy_1}{dt}&=b_{20}x_1^2+R_{21}(x_1,y_1), \end{aligned} \right.$

(41)

where $R_{i1}(i=1,2)$ is $C^{\infty}$ in $(x_1, y_1)$ and $R_{i1}=O(|(x_1,y_1)|^3)$, and

$\begin{align}\label{a20} &a_{11}=\frac{a_2+m+2P^*}{P^*(a_2+m+P^*)}>0,& \\ &a_{20}=\eta(K-2m-a_1-3P^*)(a_2+m+P^*),& \\ &b_{20}=-P^*(a_2+m+P^*)Z^*a>0. \end{align}$

(42)

Making the near-identity transformation

$x_2=x_1, y_2=y_1+a_{11}x_1y_1+a_{20}x^2_1+R_{11}(x_1,y_1),$

(43)

we have

$\left\{\begin{aligned} \frac{dx_2}{dt}&=y_2,\\ \frac{dy_2}{dt}&=b_{20}x_2^2+a_{11}y_2\frac{y_2-a_{20}x^2_2}{1+a_{11}x_2} +2a_{20}x_2y_2+R_{22}(x_2,y_2), \end{aligned} \right.$

(44)

where $R_{22}$ is $C^{\infty}$ in $(x_1, y_2)$ and $R_{22}=O(|(x_2,y_2)|^3)$.

Furthermore, utilizing the near-identity transformation

$dt=(1+a_{11}x_2)d\tau,$

(45)

we obtain

$\left\{\begin{aligned} \frac{dx_2}{d\tau}&=(1+a_{11}x_2)y_2,\\ \frac{dy_2}{d\tau}&=b_{20}x_2^2+2a_{20}x_2y_2+a_{11}y^2_2+R_{23}(x_2,y_2), \end{aligned} \right.$

(46)

where $R_{23}$ is $C^{\infty}$ in $(x_2, y_2)$ and $R_{23}=O(|(x_2,y_2)|^3)$.

By a time reparameterization

$u=x_2 \text{and} v=(1+a_{11}x_2)y_2,$

(47)

we obtain

$\left\{\begin{aligned} \frac{du}{d\tau}&=v,\\ \frac{dv}{d\tau}&=b_{20}u^2+2a_{20}uv+R_{24}(u,v), \end{aligned} \right.$

(48)

where $R_{24}$ is $C^{\infty}$ in $(u, v)$ and $R_{24}=O(|(u,v)|^3)$.

According to (47) and $P^*=-\frac{b}{2a}$, it is easy to see that $b_{20}>0$. And note that $K^*\neq 3P^*+a_1+2m^*$, we have $a_{20}\neq0$.

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 $(\theta, m)$ in a small neighborhood of $(K^*, m^*)$ when $K^*\neq 3P^*+a_1+2m^*$. Thus we let

$\left\{\begin{aligned} K&=K^*+ \epsilon_1,\\ m&=m^*+\epsilon_2, \end{aligned} \right.$

(49)

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

$\left\{\begin{aligned} \frac{dP}{dt}&=(a_2+m^*+\epsilon_2+P)\big{[}\eta(P+m^*+\epsilon_2)(K+\epsilon_1-m^*-\epsilon_2-P)(a_1+P)\\ &-PZ\big{]},\\ \frac{dZ}{dt}&=\beta_2 P(a_2+m^*+\epsilon_2+P)Z-d(a_1+P)(a_2+m^*+\epsilon_2+P)Z\\ & -\theta(P+m^*+\epsilon_2)(a_1+P)Z, \end{aligned} \right.$

(50)

for sufficiently small $(\epsilon_1, \epsilon_2)$.

System (50) has a cusp point $(P^*, Z^*)$ if $(\epsilon_1, \epsilon_2)=(0, 0)$. Using the substitution

$\left\{\begin{aligned} \bar{P}&=P-P^*,\\ \bar{Z}&=Z-Z^*, \end{aligned} \right.$

(51)

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

$\left\{\begin{aligned} \frac{d\bar{P}}{dt}&=(p_{00}+p_{10}\bar{P}+p_{20}\bar{P}^2+p_{30}\bar{P}^3+p_{40}\bar{P}^4)+(p_{01}+p_{11}\bar{P}+p_{21}\bar{P}^2)\bar{Z},\\ \frac{d\bar{Z}}{dt}&=(q_{00}+q_{10}\bar{P}+q_{20}\bar{P}^2)+(q_{01}+q_{11}\bar{P}+q_{21}\bar{P}^2)\bar{Z}, \end{aligned} \right.$

(52)

where

$\begin{align} &p_{00}=\eta[\epsilon_2(K-2m^*-2P^*-\epsilon_2)+\epsilon_1(m^*+\epsilon_2+P^*)](a_1+P^*)& \\ &(a_2+m^*+\epsilon_2+P^*),& \\ &p_{10}=\eta\{(a_2+m^*+P^*+\epsilon_2)[\epsilon_2(K-2m^*-4P^*-2a_1-\epsilon_2)& \\& +\epsilon_1(a_1+2P^*+m^*+\epsilon_2)]& \\& +(a_1+P^*)[\epsilon_2(K-2m^*-2P^*-\epsilon_2)+\epsilon_1(m^*+\epsilon_2+P^*)]\},& \\ &p_{01}=-P^*(a_2+m^*+P^*+\epsilon_2),& \\ &p_{11}=-(a_2+m^*+2P^*+\epsilon_2),& \\ &p_{20}=\eta\{(K-2m^*-a_1-3P^*-2\epsilon_2)(a_2+m^*+P^*+\epsilon_2)& \\& +\epsilon_2(K-2m^*-4P^*-2a_1-\epsilon_2)+\epsilon_1(a_1+a_2+2P^*+2m^*+2\epsilon_2)\},& \\ &p_{21}=-1,& \\ &p_{30}=\eta \{(K-3m^*-4P^*-a_1-a_2-3\epsilon_2)+\epsilon_1 \},& \\ &p_{40}=-\eta,& \\ &q_{00}=\epsilon_2((\beta_2-d-\theta)P^*-(d+\theta)a_1)Z^*,& \\ &q_{10}=\epsilon_2(\beta_2-d-\theta)Z^*,& \\ &q_{01}=\epsilon_2((\beta_2-d-\theta)P^*-(d+\theta)a_1),& \\ &q_{11}=\epsilon_2(\beta_2-d-\theta),& \\ &q_{20}=(\beta_2-d-\theta)Z^*,& \\ &q_{21}=\beta_2-d-\theta. \end{align}$

(53)

Following the transformation

$\left\{\begin{aligned} \tilde{P}&=\bar{P},\\ \tilde{Z}&=(p_{00}+p_{10}\bar{P}+p_{20}\bar{P}^2+p_{30}\bar{P}^3+p_{40}\bar{P}^4)+(p_{01}+p_{11}\bar{P}+p_{21}\bar{P}^2)\bar{Z}, \end{aligned} \right.$

(54)

system (52) is rewritten in the following form

$\left\{\begin{aligned} \frac{d\tilde{P}}{dt}&=\tilde{Z},\\ \frac{d\tilde{Z}}{dt}&=f_{00}+f_{10}\tilde{P}+f_{01}\tilde{Z}+f_{11}\tilde{P}\tilde{Z}+f_{20}\tilde{P}^2+f_{02}\tilde{Z}^2+O(\|(\tilde{P},\tilde{Z})\|^3), \end{aligned} \right.$

(55)

where

$\begin{align}\label{a28} &f_{00}=p_{01}q_{00}-p_{00}q_{01},&\\ &f_{10}=p_{01}q_{10}+p_{11}q_{00}+p_{00}q_{11}-p_{10}q_{01},& \\ &f_{01}=p_{10}+q_{01}-\frac{p_{00}p_{11}}{p_{01}},& \\ &f_{11}=2p_{20}+q_{11}+\frac{p_{00}p_{11}^2-p_{10}p_{01}p_{11}-2p_{00}p_{01}p_{21}}{p_{01}^2},& \\ &f_{20}=p_{01}q_{20}+p_{11}q_{10}+p_{21}q_{00}-p_{20}q_{01}-p_{10}q_{11}-p_{00}q_{21}+\frac{p_{11}(p_{10}q_{01}+p_{00}q_{11})}{p_{01}}& \\&-\frac{p_{00}(p_{11}^2-p_{21})q_{01}}{p_{01}^2},& \\ &f_{02}=p_{11}. \end{align}$

(56)

By rescaling the time variable $t$ using

$t=(1-f_{02}\tilde{Z})\tau,$

(57)

system (55) becomes

$\left\{\begin{aligned} \frac{d\tilde{P}}{d\tau}&=(1-f_{02}\tilde{Z})\tilde{Z},\\ \frac{d\tilde{Z}}{d\tau}&=(1-f_{02}\tilde{Z})(f_{00}+f_{10}\tilde{P}+f_{01}\tilde{Z}+f_{11}\tilde{P}\tilde{Z}+f_{20}\tilde{P}^2+f_{02}\tilde{Z}^2\\ &+O(\|(\tilde{P},\tilde{Z})\|^3)). \end{aligned} \right.$

(58)

Furthermore, let

$\left\{\begin{aligned} \hat{P}&=\tilde{P},\\ \hat{Z}&=(1-f_{02}\tilde{Z})\tilde{Z}, \end{aligned} \right.$

(59)

then we have

$\left\{\begin{aligned} \frac{d\hat{P}}{d\tau}&=\hat{Z},\\ \frac{d\hat{Z}}{d\tau}&=f_{00}+(f_{10}-2f_{00}f_{02})\hat{P}+f_{01}\hat{Z}+(f_{20}+f_{00}f_{02}^2-2f_{10}f_{02})\hat{P}^2\\ & +(f_{11}-f_{01}f_{02})\hat{P}\hat{Z}+O(\|(\hat{P},\hat{Z})\|^3). \end{aligned} \right.$

(60)

By (53), we can compute

${\displaystyle \lim\limits_{\epsilon_1\rightarrow 0, \epsilon_2\rightarrow 0}(f_{11}-f_{01}f_{02})=\eta(K-2m^*-a_1-3P^*)(a_2+m^*+P^*)}\neq 0,$

(61)

thus setting $\hat{P}=P-\frac{f_{01}}{f_{11}-f_{01}f_{02}}$ and renaming $P$ to be $\hat{P}$, then system (60) is equivalent to

$\left\{\begin{aligned} \frac{d\hat{P}}{d\tau}&=\hat{Z},\\ \frac{d\hat{Z}}{d\tau}&=g_{00}+g_{10}\hat{P}+g_{20}\hat{P}^2+ g_{11}\hat{P}\hat{Z}+O(\|(\hat{P},\hat{Z})\|^3). \end{aligned} \right.$

(62)

where

$\begin{align}\label{a32} &g_{00}=f_{00}-\frac{f_{01}(f_{10}-2f_{00}f_{02})}{f_{11}-f_{01}f_{02}}+\frac{f_{01}^2(f_{20}+f_{00}f_{02}^2-2f_{10}f_{02})}{(f_{11}-f_{01}f_{02})^2},&\\ &g_{10}=f_{10}-2f_{00}f_{02}-\frac{2f_{01}(f_{20}+f_{00}f_{02}^2-2f_{10}f_{02})}{f_{11}-f_{01}f_{02}},& \\ &g_{20}=f_{20}+f_{00}f_{02}^2-2f_{10}f_{02},& \\ &g_{11}=f_{11}-f_{01}f_{02}. \end{align}$

(63)

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

${\displaystyle \lim\limits_{\epsilon_1\rightarrow 0, \epsilon_2\rightarrow 0}g_{11}=\eta(K-2m^*-a_1-3P^*)(a_2+m^*+P^*)}\neq 0,$

(64)

and

${\displaystyle \lim\limits_{\epsilon_1\rightarrow 0, \epsilon_2\rightarrow 0}g_{20}=-P^*(a_2+m^*+P^*)(\beta_2-d-\theta)Z^*}\neq0,$

(65)

therefore, by performing the change of new variables: $\xi_1=\frac{g_{11}^2}{g_{20}}\hat{P}, \xi_2=\frac{g_{11}^3}{g_{20}^2}\hat{Z}, t=\frac{g_{20}}{g_{11}}\tau$, which finally leads that system (3) takes the following generic normal form

$\left\{\begin{aligned} \frac{d\xi_1}{dt}&=\xi_2,\\ \frac{d\xi_2}{dt}&=\frac{g_{00}g_{11}^4}{g_{20}^3}+\frac{g_{10}g_{11}^2}{g_{20}^2}\xi_1+\xi_1^2+\xi_1\xi_2+O(\|(\hat{\xi_1},\hat{\xi_2})\|^3). \end{aligned} \right.$

(66)

Corollary 1. If $m=m^*$, $K= K^*$ and $K^*= 3P^*+a_1+2m^*$, then the equilibrium $E_3$ 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 $E_1$ (or $E_2$) is a center-type non-hyperbolic equilibrium of system (8) when $h(P_1)=0$ (or $h(P_2)=0$). In this case, system (8) may undergo a Hopf bifurcation at the equilibrium $E_1$ (or $E_2$). Let $E(\tilde{P}, \tilde{Z})$ be any one of the candidates $E_1$ and $E_2$.

To determine the stability of the positive equilibrium $E(\tilde{P}, \tilde{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-\tilde{P} \text{and} y=Z-\tilde{Z}$

(67)

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

$\left\{\begin{aligned} \frac{dx}{dt}&= a_{01}y +a_{20}x^2+a_{11}xy+a_{30}x^3+a_{21}x^2y +O(x^4),\\ \frac{dy}{dt}&=b_{10}x+b_{11}xy+b_{20}x^2+b_{21}x^2y,\\ \end{aligned} \right.$

(68)

where

$\begin{align} &a_{01}=-\tilde{P}(a_2+m+\tilde{P}),& \\ &a_{11}=-(a_2+m+2\tilde{P}),& \\ &a_{20}=\eta(K-2m-a_1-3\tilde{P})(a_2+m+\tilde{P}),& \\ &a_{21}=-1,& \\ &a_{30}=\eta(K-3m-a_1-a_2-4\tilde{P}),& \\ &b_{10}=(2a\tilde{P}+b)\tilde{Z},& \\ &b_{11}=(2a\tilde{P}+b),& \\ &b_{20}=\tilde{Z}a,& \\ &b_{21}=a.& \end{align}$

$a$ and $b$ here are given in (11).

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

$\begin{aligned} \sigma&=-\frac{3\pi}{2\Delta^{\frac{3}{2}}}(-2a_{01}a_{20}b_{20}-a_{11}a_{20}b_{10}+3a_{01}a_{30}b_{10}), \end{aligned}$

(69)

where it follows from the proof of Theorem 4.2 that $\Delta=-a_{01}b_{10}=(2a\tilde{P}+b)(a_2+m+\tilde{P})\tilde{P}\tilde{Z}>0$.

Theorem 4.6. For system (8), fix all parameters except $m>0$ (or $\theta>0$) and suppose there exists $m_k>0$ (or $\theta_k>0$) such that $h(\tilde{P}, m_k)$=0 (or $h(\tilde{P}, \theta_k)$=0) and $\frac{\partial h(\tilde{P}, m)}{\partial m}|_{m=m_k}\neq 0$ (or $\frac{\partial h(\tilde{P}, \theta)}{\partial \theta}|_{\theta=\theta_k}\neq 0$) for the equilibrium $E(\tilde{P}, \tilde{Z})$. If $\sigma\neq0$, then a curve of periodic solutions bifurcation from $E(\tilde{P}, \tilde{Z})$ such that

(i) a supercritical Hopf bifurcation occurs at $E(\tilde{P}, \tilde{Z})$ if $\sigma<0$;

(ii) a subcritical Hopf bifurcation occurs at $E(\tilde{P}, \tilde{Z})$ if $\sigma>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^*=\beta_1 Z$.

We first study the impact of toxin released by phytoplankton on the dynamical behaviors of system (3) through system (8) more intuitively, the Hopf bifurcation diagram for system (8) with $\theta$ 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 $\theta$ for $H$ from left to right are $0.194742$ and $0.201745$, where the corresponding first Liapunov numbers are $-5.531406e-4$ and $-4.126233e-4$ respectively. This means that a Hopf bifurcation occurs, which generates limit cycles for model (8).

From Fig. 5, we find that the related dynamical behavior of system (3) depends on the toxin level released by phytoplankton. The phytoplankton and zooplankton is able to be in a stable state with $\theta$ small enough. With the increasing of $\theta$, 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 $\theta$ increasing, the system then is able to become stable again where the abundance of phytoplankton will be at a higher level. If $\theta$ 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.

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 $\theta=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.817476e-4$, $-2.280040e-4$ and $-4.253572e-04$ respectively.

At the same time, the equilibrium level of phytoplankton due to the refuge becomes large as the value of $m$ grows, which shows that the population of phytoplankton can increase because the refuge effect protects phytoplankton from the predation of zooplankton.

From Fig. 7, we observe that the related dynamical behavior of system (3) also relies on the refuge capacity $m$. If $m$ is not large, there is a limit cycle for system (3). However, system (3) has a global asymptotically stable equilibrium with the limit cycle disappearing as $m$ becomes progressively larger. Keeping $m$ increasing, the system is able to become oscillatory again where the abundance of phytoplankton will be at a higher level by comparing the case with small $m$ value. With larger $m$, system (3) becomes stable again. If the value of $m$ is large enough, however, zooplankton due to the lack of touchable phytoplankton becomes extinct. But phytoplankton population due to the excessive protection of the refuge effect can achieve the environmental carrying capacity. We therefore come to a conclusion that the refuge also pays an positive role in the growth of phytoplankton and has an impressible effect on the occurrence and termination of algal blooms.

With the purpose of exploring the influences of toxin and refuge on the dynamical behavior of our systems in combination, by choosing $m$ and $\theta$ as two bifurcation parameters, the corresponding bifurcation diagram is presented in Fig. 8 respectively. One Bogdanov-Takens point (BT) can be detected. According to the theoretical analysis, there is an equilibrium with a double zero eigenvalues for our system at BT point.

Beside, in order to study the complicated dynamics of system (8) at the critical positive equilibrium $E(P^*, Z^*)$ in subsection 4.3 and 4.4, we first fix the value of $\theta$ ($\theta=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 $(\epsilon_1, \epsilon_2)=(-0.05, -0.03)$. Second, we select $\theta$ 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^*, \theta^*)=(12.7759,0.6146, 0.2681)$. By perturbing $(K, m, \theta)$ in a small neighborhood of $(K^*, m^*, \theta^*)$, the numerical simulations display that two limit cycles can appear for $(\epsilon_1,\epsilon_2, \epsilon_3)=(0.0191099,0.000836,-0.0001)$, see Fig. 10-13. Therefore, one can see that BT bifurcation with different codimension can occur in system (8) when the bifurcation parameters are chosen properly, which is consistent with our theoretical analysis as before.

---

6. Discussion

In this paper, we have studied a model incorporating the effects of both refuge and toxin of phytoplankton on the phytoplankton-zooplankton interactions. Despite its simplicity, our model still gives fascinating insights into the phenomenons occurring in lakes.

Comparing the study about the phytoplankton-zooplankton model with toxin in [15], the complex dynamical behaviors can be observed for a suitable set of parameter values in our model with the effect of phytoplankton refuge. By selecting the related parameter values of toxin release rate and refuge at some certain levels, our 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.

---