Stability and Hopf bifurcation of a delayed diffusive phytoplankton-zooplankton-fish model with refuge and two functional responses

1.   Introduction

Plankton is divided into two groups: zooplankton and phytoplankton. Phytoplankton are the primary producers in aquatic ecological models, as well as the main supplier of dissolved oxygen to phytoplankton blooms. Phytoplankton opens up the food web of aquatic ecosystems. Zooplankton, as economic aquatic animals, constitute an important feed for fish and other economic animals in the middle and upper waters, which is of great significance to the development of fishery^[1,2,3,4].

Plankton has been studied extensively in many ways. In References^[1,2,3,4], the researchers found that zooplankton eat plankton and zooplankton smaller than themselves, or feed on algae, bacteria, copepods and other food scraps. Therefore, there is a predator-prey relationship between zooplankton and phytoplankton. In 1939, Fleming^[5] presented the first mathematical model of plankton. Since then, researchers have done a lot of work on plankton^[6,7,8,9,10], focusing on factors such as nutrients, temperature, light, viral diseases and harvest to understand the bloom and disappearance of algal blooms. However, the researchers found that toxins released by toxic-producing phytoplankton have an effect on the termination of plankton blooms, which means that toxic chemicals can act as biological control for other plankton populations^[11].

When the predator captures the prey, the prey seeks shelter because of a survival instinct. Nature provides shelter to the prey, and this behavior keeps the balance of the predator-prey model^[12,13,14]. For lake ecosystems, prey refuge can stabilize plankton biomass by preventing phytoplankton from being temporarily eaten by zooplankton. Scholars have found that phytoplankton shelters can be obtained through benthic sediments, which can allow phytoplankton to temporarily escape from zooplankton predation. At the same time, the water layer can also form a temporary shelter for phytoplankton, and the shelter can prevent the extinction of the prey population^[15,16,17]. Li et al.^[18] proposed a model with the refuge as follows:

where $P$ and $Z$ represent the number of phytoplankton and zooplankton, respectively. $r$ represents the intrinsic growth rate of phytoplankton, $K$ represents the environmental carrying capacity of phytoplankton, $\beta_{1}$ represents the predation rate, $\beta_{2}$ represents the conversion rate, $d$ represents the mortality rate of zooplankton, $\theta$ represents the toxin release rate, $a_{1}$ and $a_{2}$ represent half-full constants, $m$ represents the number of protected phytoplankton when phytoplankton have the ability to shelter and $P-m$ denotes the number of unprotected phytoplankton that can be preyed upon by zooplankton. They^[18] studied the effect of refuge of phytoplankton on the phytoplankton - zooplankton model.

In nature, the population diffuse from one area to another in order to survive. The predator-prey model with diffusion can generate complex spatial patterns ^{[19,20,21,22]}. There are two types of diffusion: self-diffusion and cross - diffusion. The former refers to the diffusion of one species from the higher - density area to the lower - density area in order to survive, while the latter refers to the diffusion of one species as influenced by other species. Meanwhile, time delay in nature is an important factor affecting the predator - prey model. A dynamical model without time delay can only be an approximation ^[23]. Generally speaking, there are many kinds of delayed factors in the growth process, such as the digestion of food ^[24], the maturation of cells ^[25], pregnancy^[26,27]. These processes are not instantaneous and take time to complete. Later, it is found that the existence of time delay would make the positive steady state of a predator - prey model lose stability, resulting in bifurcation or periodic oscillation^[21,28]. Zhao et al.^[19] proposed a reaction - diffusion model with mature delay:

where $d_{1}$ and $d_{2}$ represent the diffusion coefficients. $\tau$ is the time required for phytoplankton to form a mature cell to release toxins. They^[19] analyzed the stability of the equilibrium, the existence and properties of Hopf bifurcation. They^[19] also found that time delay has an effect on the model (1.2). Recently, Hopf bifurcation has also continued to be investigated in fractional-order dynamical systems ^{[29,30,31,32]} and integer-order differential systems ^{[33,34,35,36,37]}.

In lake ecosystems, fish eat plankton for survival, while zooplankton eat phytoplankton. Thus, phytoplankton, zooplankton and fish form a food chain. There is a lot of work that has done by many researchers^[20,38,39].

The functional response reflects the predator's predation on the prey, which is an important part of the predator-prey model. It can be divided into prey-dependent (Holling I-IV type^[40], Ivlev type ^[41], Rosenzweig type ^[22]) and predator-dependent functional responses (Beddington-DeAngelis type ^[42,43], Crowley-Martin (C-M) type ^[44], Hassell-Varley type ^[45], ratio-dependent type^[46]). A predator-prey model that takes into account interactions between predators is more realistic. In ^[44], Crowley and Martin first proposed the C-M functional response:

where $F(H, P)$ represents the predation rate per predator, $H$ represents the density of the prey per unit of area, $P$ represents the density of the predator per unit of area, $a$ represents the attack coefficient, $b$ is the handling time and $c$ is the interference coefficient. This functional response indicates that there is interference between predators as they feed on and handle prey. In ^[47], the authors proposed a model with the C-M functional response and showed that the system has complex dynamical behavior. In ^[48], the interaction between mature prey and predator is assumed to be the C-M functional response; the authors analyzed the positivity, boundedness and existence of equilibrium points. They not only analyzed the stability behavior of the delayed and non-delayed system, but also discussed the properties of Hopf bifurcation by choosing delay as the bifurcation parameter. In recent years, the researchers have done a great deal of work on this predator-prey models with a C-M functional response^[49,50,51]. Therefore, in this paper, on the basis of Eq (1.3), we will consider that fish predation on zooplankton follows the C-M functional response:

where $Z$ and $F$ represent the number of zooplankton and fish, respectively. $a$ represents interference between zooplankton, $b$ is the interference between fish and $\gamma$ is the maximum rate of fish predation on zooplankton. For different values of $a$ and $b$, we have that (i) Eq (1.4) has turned into a Holling II functional response if $a > 0$ and $b = 0$; (ii) Eq (1.4) expresses a saturation response with respect to a predator if $a = 0$ and $b > 0$; (iii) Eq (1.4) becomes a linear mass-action functional response if $a = 0$ and $b = 0$.

In this paper, we consider the self-diffusion of populations. Let $P(x, t)$, $Z(x, t)$ and $F(x, t)$ be the densities of phytoplankton, zooplankton and fish at the location $x$ and the time $t$, respectively. We also give the following assumptions.

(1) The phytoplankton population follows the logistic growth under the condition of no zooplankton.

(2) Zooplankton preyed upon by fish is a C-M functional response.

(3) Phytoplankton preyed upon by zooplankton is a Holling II functional response. Nature provides shelter to the prey, so a constant proportion $m\in(0, 1)$ of the phytoplankton take refuge, leaving $(1-m)P$ of the unprotected phytoplankton available for zooplankton grazing, following the Holling II functional response, i.e., $f_{1}(P) = \frac{(1-m)P}{a_{1}+(1-m)P}$.

(4) The progress of release of toxins takes the Holling II functional response and considers the mature delay of toxins in the cell, i.e., $f_{2}(P) = \frac{P(t-\tau)}{a_{2}+P(t-\tau)}$.

Based on the above assumptions, a schematic diagram that expresses the interactions of phytoplankton, zooplankton and fish is depicted in Figure 1. The corresponding model is

where $d_{1}$, $d_{2}$ and $d_{3}$ represent the diffusion coefficients of each population, respectively; $\Delta$ is the Laplace operator, $\beta_{1}$ represents the maximum predation rate of phytoplankton by zooplankton, $\beta_{2}$ is the conversion rate, $K$ represents the environmental carrying capacity, $\gamma_{1}$ represents the predation rate, $\gamma_{2}$ represents the conversion rate, $m$ is the refuge proportion of phytoplankton, $\alpha_{1}$ and $\alpha_{2}$ are half saturation constants, $\delta$ represents the release rate of phytoplankton toxins and $g_{1}$ and $g_{2}$ represent the natural mortality rates of zooplankton and fish, respectively. $a$ represents the degree of interference between zooplankton, $b$ represents the degree of interference between fish and $\tau$ represents the mature time delay of phytoplankton toxin release. The Neumann boundary condition indicates that the area is closed and no individuals can move across this area.

We rescale the model (1.5) by

For the sake of convenience, omitting the breaking line, the model (1.5) becomes

where $0 < P < 1,$ $0 < m < 1$ and all parameters are positive.

Our paper is organized as follows. The existence and stability of equilibrium of the model (2.1) are discussed in Section 2. Meanwhile, the occurrence of Hopf bifurcation is given by choosing $m$ as a bifurcation parameter. The existence and properties of Hopf bifurcation of the model (3.1) are discussed in Section 3. In Section 4, Hopf bifurcation of the reaction-diffusion model (1.6) at the positive equilibrium and its properties are analyzed. In Section 5, a numerical simulation is demonstrated to prove the previous theoretical results by using Matlab software. In the last section, we have a brief discussion.

2.   Existence and stability of equilibria of ODE model

In this section, we will investigate the dynamical behavior of the model (1.6) with no delay and no diffusion. That is, the model is

According to the existence theorem of the solution of ordinary differential equations, we can know that the solution of the model (2.1) exists. By Lemma 2.1 in Reference^[52], we give the following lemma to explain the positivity of the solution of the model (2.1).

Lemma 2.1. All solutions of the model (2.1) that start positive remain positive.

Proof. Let $(P(t), Z(t), F(t))$ be any solution of the model (2.1). Assuming that the initial time is $t_{0}$ and one solution of the model (2.1) is at least not positive, then we have the following three cases:

(i) there exists time $t_{1}$ such that $P(t_{0}) > 0, \; P(t_{1}) = 0, \; P'(t_{1})\leq 0, \; Z(t) > 0, \; F(t) > 0, \; t_{0} < t < t_{1}$;

(ii) there exists time $t_{2}$ such that $Z(t_{0}) > 0, \; Z(t_{2}) = 0, \; Z'(t_{2})\leq 0, \; P(t) > 0, \; F(t) > 0, \; t_{0} < t < t_{2}$;

(iii) there exists time $t_{3}$ such that $F(t_{0}) > 0, \; F(t_{3}) = 0, \; F'(t_{3})\leq 0, \; P(t) > 0, \; Z(t) > 0, \; t_{0} < t < t_{3}$.

If the first case is true, then we get $P'(t_{1}) = 0$. This contradicts with $P'(t_{1})\leq 0$. Similarly, we have that $Z'(t_{2}) = 0$, which contradicts with $Z'(t_{2})\leq 0$. And, $F'(t_{3}) = 0$, which contradicts with $F'(t_{3})\leq 0$.

Because of the arbitrariness of $P(t), Z(t)$ and $F(t)$, all solutions of the model (2.1) remain positive for all $t > t_{0}$. Thus, all solutions of the model (2.1) that start positive remain positive.

2.1.   Existence of all equilibria

Obviously, the model (2.1) has three boundary equilibria: one trivial equilibrium $E_{0}(0, 0, 0)$ and two boundary equilibria $E_{1}(1, 0, 0)$ and $E_{2}(P_{2}, Z_{2}, 0)$ if $(H_{1}):\beta_{2}-\delta-g_{1} > 0$ holds, where

here, $B = (\beta_{2}-g_{1})(1-m)\alpha_{2}-(\delta+g_{1})\alpha$.

Assume that the model (2.1) has the coexistence equilibrium $E_{*}(P_{*}, Z_{*}, F_{*})$, where $P_{*}$, $Z_{*}$ and $F_{*}$ satisfy

By the first equation of (2.2), we can obtain

Substituting Eq (2.3) into the third equation of (2.2), we have

If $(H_{2}):(\gamma_{2}-bg_{2})(1-P_{*})[\alpha_{1}+(1-m)P_{*}]-bg_{2}a\beta_{1}(1-m) > 0$ holds, then we have that $F_{*} > 0$.

Substituting Eqs (2.3) and (2.4) into the second equation of (2.2), $P_{*}$ is the positive root of the equation

where

Under the condition $(H_{1})$, we have that $M_{5} > 0$. By the Descarte's rule of signs^[53], Eq (2.5) has only one positive root $P_{*}$ if and only if one of the following terms is satisfied:

Here, we give two figures as examples of the first two cases to verify the conclusion that there is only one positive root (see Figure 2). We first give the assumption $(H_{3})$: one of the conditions (1)–(5) is true. Therefore, the following conclusion can be obtained.

Theorem 2.1. Under the conditions $(H_{1})$, $(H_{2})$ and $(H_{3})$, the model (2.1) has only one positive equilibrium $E_{*}(P_{*}, Z_{*}, F_{*})$, which is determined by Eqs (2.3)–(2.5).

2.2.   Stability of all equilibria

The stability of all equilibria will be analyzed in this part.

The Jacobian matrix of the model (2.1) is

where

The characteristic equation of the model (2.1) is

According to Eq (2.7), we have the following conclusions.

Theorem 2.2. The boundary equilibrium $E_{0}(0, 0, 0)$ of the model (2.1) is always unstable.

Proof. The characteristic equation of the model (2.1) at $E_{0}$ is

It has three roots:

Thus, the boundary equilibrium $E_{0}$ is unstable.

Theorem 2.3. If $(H_{4}):\frac{\beta_{2}(1-m)}{\alpha_{1}+1-m}-\frac{\delta }{\alpha_{2}+1}-g_{1} < 0$ holds, then the boundary equilibrium $E_{1}(1, 0, 0)$ of the model (2.1) is locally asymptotically stable.

Proof. The characteristic equation of the model (2.1) at $E_{1}$ is

It has three roots:

When $(H_{4})$ is true, then $\lambda_{2} < 0$. Therefore, the boundary equilibrium $E_{1}$ is locally asymptotically stable under the condition $(H_{4})$.

The characteristic equation of the model (2.1) at $E_{2}(P_{2}, Z_{2}, 0)$ is

where

Assume that the condition $(H_{5}):s_{1} < 0,$ $s_{2}^2-4s_{3} > 0,$ $s_{2} < 0$ and $s_{3} > 0$ are true; we can get the following result.

Theorem 2.4. Under the assumptions $(H_{1})$ and $(H_{5})$, the boundary equilibrium $E_{2}(P_{2}, Z_{2}, 0)$ of the model (2.1) is locally asymptotically stable.

The characteristic equation of the model (2.1) at $E_{*}(P_{*}, Z_{*}, F_{*})$ is

where

From the Routh-Hurwitz criterion^[54], if $M_{6} > 0,$ $M_{7} > 0$ and $M_{6}M_{7}-M_{8} > 0$, then all solutions of Eq (2.8) have negative real parts. When $M_{6} > 0,$ $M_{7} > 0$ and $M_{6}M_{7}-M_{8} < 0$, Eq (2.8) has one negative root and a pair of complex roots with a positive real part. Assume that $(H_{6}):$ $M_{6} > 0$, $M_{7} > 0$ and $M_{6}M_{7}-M_{8} > 0$; then, the stability of the positive equilibrium $E_{*}$ will be obtained.

Theorem 2.5. Suppose that the conditions $(H_{1})$–$(H_{3})$ are true. If the condition $(H_{6})$ holds, then $E_{*}$ is locally asymptotically stable. Further, $E_{*}$ loses stability when $M_{6}M_{7}-M_{8}$ passes through 0; in other words, the model (2.1) undergoes a Hopf bifurcation at $E_{*}$ when $M_{6}M_{7}-M_{8} = 0$.

Next, we will choose $m$ as the bifurcation parameter to study the occurrence of Hopf bifurcation of the model (2.1) at $E_{*}$. By using the results in Reference^[55], the following result can be obtained.

Theorem 2.6. If the characteristic equation of the model (2.1) at $E_{*}(P_{*}, Z_{*}, F_{*})$ is

where $M_{6}(m)$, $T(m) = M_{6}(m)M_{7}(m)-M_{8}(m)$ and $M_{8}(m)$ are the smooth functions of $m$ and there exists a positive number $m^*$ that satisfies

(1) $M_{6}(m^*) > 0,$ $T(m^*) = 0$ and $M_{8}(m^*) > 0$;

(2) $\frac{{\rm d}T}{{\rm d}m}|_{m = m^*}\neq 0,$

then Hopf bifurcation occurs at $E_{*}(P_{*}, Z_{*}, F_{*})$ when $m = m^*$.

We used Matlab software for numerical simulations to obtain this result in Section 5. Meanwhile, we can also choose $\delta$ as a bifurcation parameter to study the occurrence of Hopf bifurcation of the model (2.1) at $E_{*}(P_{*}, Z_{*}, F_{*})$. Since the discussion process is similar, we will only give the bifurcation diagram in Section 5.

Under some conditions, Hopf bifurcation may not take place. Thus, we will discuss the global asymptotical stability of the positive equilibrium $E_{*}$ as follows.

Theorem 2.7. Suppose that the conditions $(H_{1})$–$(H_{3})$ are true. The positive equilibrium $E_{*}$ of the model (2.1) is globally asymptotically stable if $1-\frac{\beta_{1}(1-m)^2Z_{*}}{\alpha_{1}[\alpha_{1}+(1-m)P_{*}]}-\frac{\delta^2} {[\alpha_{2}+P_{*}]^2} > 0$ and $\frac{[\alpha_{1}+(1-m)P_{*}]\beta_{1}P_{*}\delta}{\beta_{2}\alpha_{1}(\alpha_{2}+P_{*})} -2-\frac{\beta_{1}^2(1-m)P_{*}\gamma_{2}}{\alpha_{1}g_{1}ab} > 0$.

Proof. Let $(P, Z, F)$ be any positive solution of the model (2.1). Define a Lyapunov function

Calculating the derivative of $V(t)$ along the solution of the model (2.1), then we have

Here, it is obvious that $\frac{\gamma_{2}Z_{*}}{(b+F)(a+Z_{*})(b+F_{*})}+\frac{\gamma_{2}^2a^2}{(a+Z)^2(b+F)^2(a+Z_{*})^2} > 0$. If $1-\frac{\beta_{1}(1-m)^2Z_{*}}{\alpha_{1}[\alpha_{1}+(1-m)P_{*}]}-\frac{\delta^2}{[\alpha_{2}+P_{*}]^2} > 0$ and $-2-\frac{\beta_{1}^2(1-m)P_{*}\gamma_{2}}{\alpha_{1}g_{1}ab}+\frac{[\alpha_{1}+(1-m)P_{*}]\beta_{1}P_{*}\delta} {\beta_{2}\alpha_{1}(\alpha_{2}+P_{*})} > 0$ hold, then we have that the coefficients of $(P-P_{*})^2$, $(Z-Z_{*})^2$ and $(F-F_{*})^2$ are always negative. Thus, $\frac{{\rm d}V(t)}{{\rm d}t}$ is negative. This completes the proof.

3.   Hopf bifurcation of DDE model

3.1.   Existence of Hopf bifurcation

Here, under the conditions $(H_{1})$–$(H_{3})$, we choose delay $\tau$ as the bifurcation parameter and study its influence on the stability of the positive equilibrium $E_{*}(P_{*}, Z_{*}, F_{*})$. The model (2.1) is

Let $u_{1}(t) = P(t)-P_{*},$ $u_{2}(t) = Z(t)-Z_{*},$ $u_{3}(t) = F(t)-F_{*}$ and $u(t) = (u_{1}(t), u_{2}(t), u_{3}(t))^T$ $\in{R^3}$. The linearized system of the model (3.1) at $E_{*}$ is

where

Then, the model (3.2) can also be given by

where

Thus, O(0, 0, 0) is the zero equilibrium of the model (3.2). The characteristic equation of the model (3.2) at O is

where

The roots of Eq (3.5) have been discussed above when $\tau = 0$. Next, we will study the effect of delay $\tau(\tau > 0)$ on the model (3.2).

Suppose that $\lambda = i\omega(\omega > 0)$ is a pair of pure imaginary roots of Eq (3.5). Substituting it into Eq (3.5), we can obtain

Separating the real and imaginary parts of Eq (3.6), we can obtain

From Eq (3.7), we can get

Let $z = \omega^2$, $D_{5} = D_{2}^2-2D_{1}$, $D_{6} = D_{1}^2-2D_{0}D_{1}-D_{3}^2$ and $D_{7} = D_{0}^2-D_{4}^2.$ Then, Eq (3.8) can be rewritten as

and we have

If Eq (3.9) has at least one positive root, then Hopf bifurcation takes place. Assume that $(H_{7}):\Delta_{1} = D_{5}^2-3D_{6}\leq0$ and $(H_{8}):\Delta_{1} = D_{5}^2-3D_{6} > 0.$ By Lemmas 2.2 and 4.2 in Reference ^[56], we can get the following results.

Since ${\lim_{z\rightarrow +\infty}f(z)} = +\infty,$ Eq (3.9) has at least one positive root when $D_{7} < 0$.

If $(H_{7})$ holds, then $f(z)$ is monotonically increasing for $z\in[0, +\infty)$; so, when $D_{7}\geq0$ and $(H_{7})$ hold, Eq (3.9) has no positive root for $z\in[0, +\infty)$.

When $D_{7}\geq0$ and $(H_{8})$ hold, Eq (3.10) has two roots, that is, $z_{1}^{*}$ and $z_{2}^{*}$, where

Furthermore, we have

Therefore, we can obtain $z_{1}^{*}$ and $z_{2}^{*}$ as the local minimum and the local maximum of $f(z)$, respectively.

Theorem 3.1. For Eq (3.9), the following conclusions are true.

(1) If $D_{7} < 0$, then Eq (3.9) has at least one positive root.

(2) If the condition $(H_{7})$ holds and $D_{7}\geq0$, then Eq (3.9) has no positive root.

(3) If the condition $(H_{8})$ holds and $D_{7}\geq0$, then Eq (3.9) has two positive roots when $z_{1}^{*} > 0$ and $f(z_{1}^{*})\leq0$.

Without loss of generality, suppose that Eq (3.9) has three positive roots defined by $z_{1},$ $z_{2}$ and $z_{3}$, respectively. Thus, Eq (3.8) has three positive roots $\omega_{1} = \sqrt{z_{1}}$, $\omega_{2} = \sqrt{z_{2}}$ and $\omega_{3} = \sqrt{z_{3}}$. From Eq (3.7), we can get

thus, $\pm i\omega_{k}$ is a pair of purely imaginary roots of Eq (3.5) when $\tau = \tau_{k}^{j}$.

Define

Let $\lambda(\tau) = \alpha(\tau)+i\beta(\tau)$ be the root of Eq (3.5) near $\tau = \tau_{k}^{j}$, and let it satisfy $\alpha(\tau_{k}^{j}) = 0$ and $\beta(\tau_{k}^{j}) = \omega_{k},$ $j = 0, 1, 2, \cdots,$ $k = 1, 2, 3.$

Theorem 3.2. Assume that $z_{k} = \omega_{k}^2$ and $f'(z_{k})\neq0$; then, $\frac{{\rm d}Re\lambda(\tau)}{{\rm d}\tau}|_{\lambda = i\omega_{0}}\neq0$.

Proof. Differentiating Eq (3.5) for $\tau$, we have

Substituting $\lambda(\tau) = i\omega_{0}$ into Eq (3.12), we can get

Therefore, we can get

Thus, we have

This ends the proof of Theorem 3.2.

Theorem 3.3. Under the conditions $(H_{1})$–$(H_{3})$, the dynamics of the model (3.1) at the positive equilibrium $E_{*}$ can be obtained.

(i) If $D_{7}\geq0$ and the condition $(H_{7})$ is true, then $E_{*}$ is locally asymptotically stable for all $\tau > 0$.

(ii) If $D_{7} < 0$ or the condition $(H_{8})$, $D_{7} > 0$, $z_{1}^{*} > 0$ and $f(z_{1}^{*})\leq0$ are true, then $E_{*}$ is locally asymptotically stable when $\tau\in[0, \tau_{0})$; but, $E_{*}$ is unstable when $\tau > \tau_{0}$.

(iii) If the conditions in (ii) are all satisfied and $f'(z_{k})\neq0$ holds, then a Hopf bifurcation occurs at $E_{*}$ when $\tau = \tau_{0}$.

3.2.   Properties of Hopf bifurcation

For $\tau = \tau_{0}$, the existence of Hopf bifurcation has been discussed. The properties of the Hopf bifurcation will be discussed in consideration of the work of Hassard et al. ^[57].

Let $\tau = \tau_{0}+\mu,$ $\mu\in{R}$; we can obtain that $\mu = 0$ is a Hopf bifurcation value of the model (3.1). Let the phase space be $\mathcal{C} = C([-1, 0], R^3)$ and $t\rightarrow \frac{t}{\tau}$; then, the model (3.1) can be expressed as a functional differential equation (FDE) in $\mathcal{C}$ as

where $L_{\mu}:\mathcal{C}\rightarrow {R^3}$ and $F:R\times\mathcal{C}\rightarrow {R^3}.$

Define $\varphi(\theta) = (\varphi_{1}(\theta), \varphi_{2}(\theta), \varphi_{3}(\theta))^{T}\in{R^3}, \theta\in[-1, 0]$ such that

and

where $L_{1}$ and $L_{2}$ are defined in Eq (3.4), and

From the Riesz representation theorem^[58], there is a matrix function whose elements are functions of bounded variation $\eta(\theta, \mu)\in C([-1, 0], \; R^3)$; also,

Let

where $\delta(\theta)$ is a Dirac function, as follows:

For $\varphi\in\mathcal{C} = C([-1, 0], R^3)$, we define that

Therefore, Eq (3.13) becomes

For $\psi\in\mathcal{C}^{*} = C([-1, 0], R^3),$ let

and establish the bilinear inner product

where $\eta(\theta) = \eta(\theta, 0)$. Then, $A(0)$ and $A^{*}(0)$ are adjoint operators. According to Theorem 3.3, we can get that $\pm i\omega_{0}\tau_{0}$ are eigenvalues of $A(0)$. Thus, they are also eigenvalues of $A^{*}(0)$. Assume that $q(\theta) = (1, q_{1}, q_{2})^{T}e^{i\omega_{0}\tau_{0}\theta}$ is the eigenvector of $A(0)$ corresponding to $i\omega_{0}\tau_{0}$; then, $A(0)q(\theta) = i\omega_{0}\tau_{0}q(\theta)$. By the definition above, we have

that is,

Solving it, we can get

Similarly, we assume that $q^{*}(s) = D(1, q_{1}^{*}, q_{2}^{*})e^{-i\omega_{0}\tau_{0}s}$ is the eigenvector of $A^{*}(0)$ corresponding to $-i\omega_{0}\tau_{0}$; then, $A^{*}(0)q^{*}(s) = -i\omega_{0}\tau_{0}q^*(s)$. By the definition above, we can get

By Eq (3.2), $\langle q^*, q\rangle$ can be expressed as

Thus, we can choose

then, $\langle q^*, q\rangle = 1.$

Next, we adopt the ideas of Hassard et al.^[57] to compute the coordinates describing the center manifold $C_{0}$ at $\mu = 0$. We assume that $u_{t}$ is the solution of Eq (3.13) when $\mu = 0$. Let

On the center manifold $C_{0}$, we can get

where

$z$ and $\bar{z}$ express the local coordinates for the center manifold $C_{0}$ in the direction of $q^*$ and $\bar{q}^*$. We can get that $W$ is real when $u_{t}$ is real. For the real solution $u_{t}\in C_{0}$ of Eq (3.17), when $\mu = 0$, we have

Eq (3.20) can also be written as

where

By Eq (3.18), we can get

where

From Eqs (3.22) and (3.23), we obtain

By Eq (3.18), we have

Then, we have

where

Comparing the coefficients of Eqs (3.21) and (3.26), we can get

Since the expression of $g_{21}$ contains $W_{20}(\theta)$ and $W_{11}(\theta)$, we must compute $W_{20}(\theta)$ and $W_{11}(\theta)$. According to Eq (3.22), when $\theta\in[-1, 0)$, we get

Comparing the coefficients of Eqs (3.23) and (3.28), we can receive

By Eq (3.24), we have

Solving Eq (3.30), we can obtain

where $E_{i} = (E_{i}^{(1)}, E_{i}^{(2)}, E_{i}^{(3)})\in R^3$ $(i = 1, 2)$ is a constant vector.

For $\theta = 0$, from Eq (3.22), we have

From Eq (3.23), we can get

According to the meaning of $A(0)$ and Eq (3.24), we have

where $\eta(\theta) = \eta(\theta, 0)$.

From Eqs (3.14), (3.16) and (3.33), we can get

Substituting Eqs (3.31) and (3.32) into Eq (3.34), we have

Thus, all expressions of $g_{ij}$ can be represented in full. Also, we have

which determine the direction of Hopf bifurcation and the stability of bifurcating periodic solutions on the center manifold at $\tau = \tau_{0}$.

Theorem 3.4. From Eq (3.35), we have the following conclusions.

(i) $\mu_{1}$ determines the direction of Hopf bifurcation: if $\mu_{1} > 0$, then the Hopf bifurcation is supercritical; if $\mu_{1} < 0$, then the Hopf bifurcation is subcritical and the bifurcating periodic solutions exist when $\tau > \tau_{0};$

(ii) $\mu_{2}$ determines the stability of the bifurcating periodic solutions: if $\mu_{2} < 0$, then the bifurcating periodic solutions are stable; if $\mu_{2} > 0$, then the bifurcating periodic solutions are unstable;

(iii) $T_{1}$ determines the period of the bifurcating periodic solutions: if $T_{1} > 0$, then the period increases; if $T_{1} < 0$, then the period decreases.

4.   Hopf bifurcation of PDE model

4.1.   Existence of Hopf bifurcation

Next, we will analyze the existence and properties of Hopf bifurcation of the model (1.6). The model (1.6) is linearized at the positive equilibrium $E_{*}$ in the phase space $\mathcal{C} = C([-\tau, 0], R^3)$:

where $D = diag\{d_{1}, d_{2}, d_{3}\}$ and $L_{1}$ and $L_{2}$ are defined in Eq (3.4).

We know that $\Delta$ has the eigenvalues $-(\frac{n}{l})^2$ and $n\in N_{0}$ under Neumann boundary conditions in $[0, l\pi]$. Then, the characteristic equation of the model (4.1) is

where

When $\tau = 0$, Eq (4.2) can be rewritten as

Furthermore, we have

Assume that the following condition holds true: $(H_{9})$: $q_{2n} > 0,$ $q_{2n}(q_{1n}+q_{3n})-(q_{0n}+q_{4n}) > 0,$ $q_{0n}+q_{4n} > 0,$ $n\in N_{0}$. Using the Routh-Hurwitz criterion^[54], we can obtain the next conclusion.

Theorem 4.1. If the conditions $(H_{1})$–$(H_{3})$ and $(H_{9})$ hold, then all roots of Eq (4.3) have negative real parts, that is, the positive equilibrium $E_{*}$ of the model (1.6) is locally asymptotically stable when $\tau = 0$.

When $\tau\neq0$, the time delay may have some effect on the model (1.6). Therefore, we will analyze the effect of delay $\tau$ on the positive equilibrium $E_{*}$. Let $\lambda = i\omega_{2}(\omega_{2} > 0)$ be the solution of Eq (4.2); we can get

Separating the real and imaginary parts of Eq (4.4), we can receive

which follows that

Let $p_{2n} = q_{2n}^2-2q_{1n}$, $p_{1n} = q_{1n}^2-2q_{0n}q_{2n}-q_{3n}^2,$ $p_{0n} = q_{0n}^2-q_{4n}^2$ and $\omega_{2}^2 = y$; we can obtain

Further, if $f(y) = y^3+p_{2n}y^2+p_{1n}y+p_{0n}$, then $f'(y) = 3y^2+2p_{2n}y+p_{1n}$. Assume that $(H_{10}):\Delta_{2} = p_{2n}^2-3p_{1n}\leq0$ and $(H_{11}):\Delta_{2} = p_{2n}^2-3p_{1n} > 0$ are true, here $y_{1}^* = \frac{-p_{1n}+\sqrt{\Delta_{2}}}{3}$ and $y_{2}^* = \frac{-p_{1n}-\sqrt{\Delta_{2}}}{3}$ are the local minimum and the local maximum of Eq (4.7), respectively. Similarly, we have the following conclusion by Lemmas 2.2 and 4.2 in^[56].

Theorem 4.2. For Eq (4.7), the following results are true.

(1) If $p_{0n} < 0$, then Eq (4.7) has at least one positive root.

(2) If $p_{0n}\geq0$ and the condition $(H_{10})$ holds, then Eq (4.7) has no positive root.

(3) If $p_{0n}\geq0$ and the condition $(H_{11})$ holds, then Eq (4.7) has positive roots when $y_{1}^{*} > 0$ and $f(y_{1}^{*})\leq0$.

Without loss of generality, we suppose that it has three positive roots defined by $y_{1n}$, $y_{2}$ and $y_{3n}$. So, Eq (4.6) has three positive roots:

Substituting $\omega_{2}^k(k = 1, 2, 3)$ into Eq (4.5), we can get

Thus, when $\tau = \tau_{kn}^j,$ we get that $\lambda = \pm i\omega_{2}^k$ is a pair of purely imaginary roots of Eq (4.2). Define

Assume that $\lambda(\tau) = \varepsilon_{1}(\tau)+i\varepsilon_{2}(\tau)$ is the root of Eq (4.2) near $\tau = \tau_{kn}^j$, and that it {satisfies} $\varepsilon_{1}(\tau_{kn}^j) = 0$ and $\varepsilon_{2}(\tau_{kn}^j) = \omega_{2}^k$, $k = 1, 2, 3,$ $j\in N_{0},$ $n\in N_{0}$.

Theorem 4.3. If $y_{kn} = (\omega_{2}^k)^2$ and $f'(y_{kn})\neq0,$ then we have that $\frac{{\rm d}Re\lambda(\tau)}{{\rm d}\tau}|_{\lambda = i\omega_{n0}}\neq0.$

Proof. Differentiating Eq (4.2) for $\tau$, we can get

Substituting $\lambda = i\omega_{n0}$ into Eq (4.9), we have

then,

Thus, we can get

The proof of Theorem 4.3 is completed.

By computing as shown above, we have the following conclusion.

Theorem 4.4. Under the conditions $(H_{1})$–$(H_{3})$, for the the positive equilibrium $E_{*}$ of the model (1.6), we have the following:

(i) if the condition $(H_{10})$ is true and $p_{0n}\geq0$, then $E_{*}$ is locally asymptotically stable for all $\tau > 0;$

(ii) if $p_{0n} < 0$ or $p_{0n}\geq0$, $y_{1}^{*} > 0$, $f(y_{1}^*)\leq0$ and the condition $(H_{11})$ holds, then $E_{*}$ is locally asymptotically stable when $\tau\in[0, \tau_{n0})$; but, $E_{*}$ is unstable when $\tau > \tau_{n0};$

(iii) if the conditions in (ii) are all satisfied and $f'(y_{kn})\neq0$, then the spatially homogeneous Hopf bifurcation occurs at $E_{*}$ when $\tau = \tau_{0}$ and $n = 0$; and, the spatially inhomogeneous Hopf bifurcation occurs at $E_{*}$ when $\tau = \tau_{n0}$ and $n > 0$.

4.2.   Direction and stability of Hopf bifurcation

Let $\tau_{n} = \tau_{kn}^j,$ $\omega_{n} = \omega_{2n}^k$ and $\tau = \tau_{n}+\mu_{n},$ $\mu_{n}\in R$. Thus, $\mu_{n} = 0$ is a Hopf bifurcation value of the model (1.6). Let $t\rightarrow \frac{t}{\tau}$; then, the model (1.6) can be expressed as an FDE in $\mathcal{C} = C([-1, 0], R^3)$, as follows:

where $L(\theta):\mathcal{C}\rightarrow X, F_{n}(\mu_{n}, u_{t}):\mathcal{C}\rightarrow X$ satisfies

and

where $F(\mu_{n}, \varphi)$ is defined in Eq (3.15), $L_{1}$ and $L_{2}$ are defined in Eq (3.4) and $\varphi = (\varphi_{1}, \varphi_{2}, \varphi_{3})^T\in \mathcal{C}.$

The linear equation of Eq (4.10) at O(0, 0, 0) is

Let $\Lambda = \{i\omega_{n}\tau_{n}, -i\omega_{n}\tau_{n}\}$ and $z_{t}(\theta)\in C = C([-1, 0], R^3)$; we consider the following FDE:

On the basis of the Riesz representation theorem, there exists a $3\times 3$ matrix function $\eta_{n}(\theta, \mu)$ $(-1\leq\theta\leq0)\in C([-1, 0], R^3)$, and it satisfies

Let

where $\delta(\theta)$ is the Dirac delta function.

We set $C^* = C([0, 1], R^{3*}),$ and $R^{3*}$ is the three - dimensional vector space of row vectors. The bilinear inner product is

$A_{n}(\tau_{n})$ describes the infinitesimal generator of the semigroup induced by the solutions of Eq (4.13), and $A_{n}^*(\tau_{n})$ denotes the formal adjoint generator of $A_{n}(\tau_{n})$ satisfying Eq (4.14). Let $V$ and $V^*$ denote the center spaces of the generators $A_{n}(\tau_{n})$ and $A_{n}^*(\tau_{n})$ corresponding to $\Lambda$, respectively. Therefore, $V^*$ is the adjoint space of $V$, ${\rm dim}V = {\rm dim}V^*$.

Lemma 4.1. Let

then, $p_{1}(\theta)=(1, V_{1}, V_{2})^Te^{i\omega_{n}\tau_{n}\theta}$ and $p_{2}(\theta)=\overline{p_{1}(\theta)}$, $-1\leq\theta\leq0$ form the basis of $V$ associated with $\Lambda$; $p_{1}^*(s)=(1, V_{1}^*, V_{2}^*)e^{-i\omega_{n}\tau_{n}s}$ and $p_{2}^*(s)=\overline{p_{1}^*(s)}$, $0\leq s\leq1$ form the basis of $V^*$ associated with $\Lambda$.

Denote $\Phi = (\Phi_{1}, \Phi_{2})$ and $\Psi^* = (\Psi_{1}^*, \Psi_{2}^*)^T$, where

Suppose that $(\Psi^*, \Phi) = (\Psi_{i}^*, \Phi_{j})(i, j = 1, 2.)$ is the basis $\Psi$ of $V^*$, which satisfies

Thus, we have that $(\Psi, \Phi) = I_{2\times2}$.

Let $f_{n} = (\xi_{n}^1, \xi_{n}^2, \xi_{n}^3),$ where $\xi_{n}^1 = (cos\frac{n}{l}x, 0, 0)^T,$ $\xi_{n}^2 = (0, cos\frac{n}{l}x, 0)^T$ and $\xi_{n}^3 = (0, 0, cos\frac{n}{l}x)^T$. $\xi_{n}^{j}$ $(j = 1, 2, 3)$ denotes the eigenfunctions on $R^3$ of the eigenvalues $-(\frac{n}{l})^2$, $n = 0, 1, 2\cdots$. Define $c_{n}\cdot f_{n} = c_{1}\xi_{n}^1+c_{2}\xi_{n}^2+c_{3}\xi_{n}^3$, $c_{n} = (c_{1}, c_{2}, c_{3})^T, c_{j}\in R,$ $j = 1, 2, 3,$ and the center space of Eq (4.12) is written as

where $\varphi\in C$, $C = P_{CN}C\oplus P_{s}C$ and $P_{s}C$ expresses the complementary subspace of $P_{CN}C$.

According to ^[57] and ^[59], the center space of the linear model of (4.10) with $\mu_{n} = 0$ is expressed as $P_{CN}C$, where

Thus, the solution of the model (4.10) can be written as

where $Q(z(t), \bar{z}(t))(\theta) = W(\frac{z+\bar{z}}{2}, i\frac{z-\bar{z}}{2}, 0)$, $z = x_{1}-ix_{2}.$

By Wu^[59], $z$ satisfies

where

Let

and

By computing and comparing the coefficients, we have

We know that $\int_{0}^{\pi}cos^3{(\frac{n}{l})x}dx = 0$ and $\tilde{g}_{02} = \overline{\tilde{g}}_{20}$. Therefore, we can get that $\tilde{g}_{20} = \tilde{g}_{11} = \tilde{g}_{02} = 0$ when $n = 1, 2, 3\cdots$. When $n = 0$, we have

Considering the expression of $\tilde{g}_{21}$, it contains $Q_{20}(\theta)$ and $Q_{11}(\theta)$, so we must compute $Q_{20}(\theta)$ and $Q_{11}(\theta)$. Seeking the derivative on both sides of Eq (4.16), we have

From Wu ^[59], $Q(z, \bar{z})$ satisfies

where

and

From Eq (4.15) and Eqs (4.17)–(4.20), we can get

Because $A_{\tau_{n}}$ has only two characteristic roots with a zero real part, i.e., $\pm i\omega_{n}\tau_{n}$, Eq (4.20) has a unique solution $Q_{ij}(i+j = 2)$ in $P_{S}C$, which satisfies

From Eq (4.20), we can get that, for $\theta\in[-1, 0)$,

Comparing the coefficients in Eqs (4.20) and (4.22), when $\theta\in[-1, 0)$, we can receive

When $\theta = 0$, we have

Using Eqs (4.24) and (4.25), we can get $Q_{20}(0),$ $Q_{11}(0),$ $Q_{20}(-1)$ and $Q_{11}(-1)$. Because $p_{1}(\theta) = p_{1}(0)e^{i\omega_{n}\tau_{n}\theta},$ $\theta\in[-1, 0),$ from Eqs (4.21)–(4.25), we have

where $E_{3} = (E_{3}^{(1)}, E_{3}^{(2)}, E_{3}^{(3)})\in R^3$ and $E_{4} = (E_{4}^{(1)}, E_{4}^{(2)}, E_{4}^{(3)})\in R^3$ satisfy

Therefore, we can get the following values:

which determine the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions on the center manifold at $\tau = \tau_{n}$.

Theorem 4.5. According to Eq (4.26), we have the following conclusions.

(i) $\mu_{3}$ determines the direction of Hopf bifurcation: if $\mu_{3} > 0$, then the Hopf bifurcation is supercritical; if $\mu_{3} < 0$, then the Hopf bifurcation is subcritical and the bifurcating periodic solutions exist when $\tau > \tau_{n};$

(ii) $\mu_{4}$ determines the stability of the bifurcating periodic solutions: if $\mu_{4} > 0$, then the bifurcating periodic solutions are unstable; if $\mu_{4} < 0$, then the bifurcating periodic solutions are stable;

(iii) $T_{2}$ determines the period of the bifurcating periodic solutions: if $T_{2} > 0$, then the bifurcating periodic solutions increase; if $T_{2} < 0$, then the bifurcating periodic solutions decrease.

5.   Numerical simulation

With the help of Matlab software, the stability of the positive equilibrium $E_{*}(P_{*}, Z_{*}, F_{*})$ was simulated with the given values of all parameters in order to confirm the previous theoretical results.

First, we assume that $P(0) = 0.8,$ $Z(0) = 40$ and $F(0) = 0.1$ for the model (2.1). And, we take the values of all other parameters as follows: $\beta_{1} = 0.016,$ $\beta_{2} = 0.7,$ $\gamma_{1} = 0.0875,$ $\gamma_{2} = 0.075,$ $\alpha_{1} = 0.1,$ $\alpha_{2} = 0.2,$ $a = 0.5,$ $b = 0.25,$ $g_{1} = 0.1,$ $g_{2} = 0.2,$ $m = 0.8$ and $\delta = 0.35$. According to Theorem 2.1, we can know that the model (2.1) has one unique positive equilibrium $E_{*}(0.5043, 31.1137, 0.1191)$. From Theorem 2.5, $E_{*}$ is locally asymptotically stable (see Figure 3).

In the second section, we obtained Theorem 2.6 by referring to ^[55]. Now, we will verify some conclusions by taking $m$ as the bifurcation parameter. When $\delta = 0.2$, the critical value is $m^* = 0.82$, which satisfies $M_{6}(m^*) > 0$, $T(m^*) = 0$, $M_{8}(m^*) > 0$ and $\frac{{\rm d}T}{{\rm d}s}|_{m = m^*}\neq0$. That is, all conditions in Theorem 2.6 are satisfied. Therefore, Hopf bifurcation occurs at $E_{*}$ when $m = m^*$. Meanwhile, we can obtain that $P(t)$ reaches a maximum value and $Z(t)$ and $F(t)$ are always $0$ when $m\geq 0.94$. Therefore, we can get the bifurcation diagram as $m$ changes (see Figure 4). If we choose $\delta$ as the bifurcation parameter when $m = 0.75$, we can get the critical value $\delta^{*} = 0.33$. Therefore, Hopf bifurcation occurs at $E_{*}$ when $\delta = \delta^*$. We can obtain that $P(t)$ reaches a maximum value of 1 and $Z(t)$ and $F(t)$ are always $0$ when $\delta\geq 0.49$. Therefore, we can get the bifurcation diagram as $\delta$ changes (see Figure 5). From Figure 5, when other parameter values are fixed, the density of zooplankton and fish will decrease to 0 whether the refuge capacity of phytoplankton or the probability of toxin release of phytoplankton-produced toxic substances increase to some certain value. Properly increasing the shelter capacity of phytoplankton and the rate of toxin release of by phytoplankton can stabilize the population and reach a stable state. Then, the plankton and fish populations will always exist.

For the model (3.1), we assume that $P(0) = 0.5,$ $Z(0) = 30$ and $F(0) = 0.115$. When $m = 0.8$ and $\delta = 0.25$, we have that $\tau_{0} = 4.9397$, and the model (2.1) has one unique positive equilibrium $E_{*}(0.2671, 35.1379, 0.1197)$ according to Theorem 2.1. From Theorem 3.3, $E_{*}$ is locally asymptotically stable when $\tau\in[0, 4.9397]$, but Hopf bifurcation occurs when $\tau\in[4.9397, +\infty)$. From Eq (3.35), we can know that $c_{1}(0) = -512.86-540.47i < 0,$ $\mu_{1} = 1457 > 0,$ $\mu_{2} = -1025.7 < 0$ and $T_{1} = 153.5416 > 0$. Thus, the Hopf bifurcation is supercritical, the bifurcating periodic solution is stable and the period of the bifurcating periodic solutions is increasing, which can be seen in Figure 6 ($\tau = 1$) and Figure 7 ($\tau = 10$). Here, we give the delay bifurcation diagram (see Figure 8). This means that, if the mature delay exceeds the critical value, the model transitions to unstable from stable. At this moment, the model has a Hopf bifurcation near the equilibrium and unstable behavior occurs among populations. In other words, the presence of the mature delay can destabilize the plankton-fish population.

For the model (1.6), we choose $l = 1$, that is, $x\in(0, \pi)$. The values of other parameters as follows: $d_{1} = 0.2,$ $d_{2} = 0.5,$ $d_{3} = 0.1,$ $m = 0.84$ and $\delta = 0.25$. The positive equilibrium is $E_{*}(0.3794, 38.9575, 0.1202)$. From Eq (4.8), we have that $\tau_{n0} = 6.2832$. Based on Theorem 4.1, $E_{*}$ is also locally asymptotically stable when $\tau = 0$ (see Figure 9), and $E_{*}$ is locally asymptotically stable when $\tau = 3.5 < \tau_{n0} = 6.2832$ (see Figure 10). But, $E_{*}$ is unstable when $\tau = 25 > \tau_{n0} = 6.2832$ (see Figure 11). And, we can compute that $c_{2}(0) = -1141.6- 2543.9i < 0,$ $\mu_{3} = 7257.8 > 0,$ $\mu_{4} = -2283.3 < 0$ and $T_{2} = 646.8790 > 0$. From Theorem 4.5, the Hopf bifurcation is supercritical, the bifurcating periodic solution is stable and the period of the bifurcating periodic solutions is increasing. For the reaction-diffusion model, we can know that the model will transitions to unstable from stable if the mature delay exceeds the critical value. At this moment, the model has a spatially homogeneous Hopf bifurcation or spatially inhomogeneous Hopf bifurcation near the equilibrium and unstable behavior occurs between the populations. At this time, the presence of the mature delay can destabilize the plankton-fish population.

By the theoretical conclusions and numerical simulation, we not only find that the existence of delay will deteriorate the system stability under some conditions, but also that the refuge of the prey and the release of toxins will cause the stability of system be damaged in a reaction-diffusion model with delay, even causing Hopf bifurcation to occur at the positive equilibrium $E_{*}$.

6.   Conclusions

In our paper, we establish a phytoplankton-zooplankton-fish model with mature delay and population diffusion by considering the refuge of phytoplankton, C-M functional response and Holling II functional response. In ^[18], the authors found that the refuge affects the stability of the positive equilibrium. However, in our paper, we not only analyzed the effect of refuge, but also studied the effects of diffusion and delay on the model; we obtained that the stability of the system may be destroyed due to the existence of delay. In ^[19], the authors analyzed Hopf bifurcation caused by delay. However, in our paper, we not only obtained the properties of Hopf bifurcation induced by delay, but also the influence of prey refuge on the population. We determined that the existence of prey refuge can also lead to Hopf bifurcation.

After the parameters in the model were selected, the existence and stability of the equilibrium were analyzed. First, we chose $m$ as a bifurcation parameter to study the dynamical behavior as $m$ changes in the model (2.1). Meanwhile, we consider the effect of the parameter $\delta$ on the positive equilibrium in the model (2.1). Through analysis, it could be obtained that the model undergoes a Hopf bifurcation when $m = m^*$ or $\delta = \delta^*$. We found that, when other parameter values are fixed, the densities of zooplankton and fish will decrease to 0 regardless whether the refuge capacity of phytoplankton or the probability of toxin release of phytoplankton-produced toxic substances increase to a certain value. And, we chose the time delay $\tau$ as the bifurcation parameter and discussed the dynamical behavior of the model without diffusion, or with diffusion, respectively. We give the direction of the Hopf bifurcation and the stability of the bifurcating periodic solution by the center manifold theorem and normal form theory. We found that the model transitions to unstable from stable when the mature delay exceeds the critical value. At this moment, the model has a spatially homogeneous Hopf bifurcation or spatially inhomogeneous Hopf bifurcation near the positive equilibrium and unstable behavior occurs between the populations. In a word, the existence of time delay has a great influence on such a model. Meanwhile, we used Matlab software for numerical simulation to prove our theoretical results.

In this paper, we have discussed the influence of factors such as prey refuge, the disturbance between predators, time delay and diffusion on the model. However, in nature, there are external factors to influence the model, such as changing temperature, environmental pollution, human activities and noise. We did not take these influencing factors into account. Therefore, in the future work, we will introduce the influence of environmental pollution on the model and analyze the dynamical behavior of the phytoplankton-zooplankton model under the influence of environmental pollution. The model is

where $m_{1}$ and $m_{2}$ are the effects of environmental pollution on phytoplankton and zooplankton, respectively. We leave this work for the future.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant Nos. 12161054, 11661050 and 11861044), the National Natural Science Foundation of Gansu Province (Grant No. 20JR10RA156), the Doctoral Foundation of Lanzhou University of Technology and the HongLiu First-Class Disciplines Development Program of Lanzhou University of Technology.

Conflict of interest

The authors declare that they have no conflict of interest.