Stability and Hopf bifurcation analysis of a delayed phytoplankton-zooplankton model with Allee effect and linear harvesting

1.   Introduction

Blooms are considered to be a rapid and significant increase in marine phytoplankton populations. It is generally believed that zooplankton is the primary predator of blooms. The zooplankton dissolves phytoplankton blooms in the bud by feeding a large amount of phytoplankton, thereby achieving the purpose of preventing and controlling algal blooms. It is clearly stated in an international research plan that the dynamic changes of zooplankton in marine ecosystems control the change in the total amount of primary productivity. This indicates that zooplankton plays a very important role in the production of phytoplankton and the control of existing quantities. In addition, this is also one of the significant reasons that the interaction between phytoplankton and zooplankton has been widely concerned by many scholars ^{[1,2,3,4,5,6]} in recent years.

Due to the importance of the interaction between phytoplankton and zooplankton in marine biological systems, many scholars have proposed a large number of models to describe this interaction. Moreover, they have continued to study this topic for the past few decades and have achieved very significant results ^{[7,8,9,10,11,12,13]}. Chakraborty et al. ^[7] proposed a mathematical model incorporating nutrient concentration, toxin producing phytoplankton (TPP), non-toxic phytoplankton (NTP), and toxin concentration and explained that TPP species controls the outbreaks of other NTP species when nutrient-deficient conditions are beneficial for the TPP species to release toxin. Saha and Bandyopadhyay ^[8] analyzed a toxin producing phytoplankton-zooplankton model in which the release of toxins from phytoplankton species follows a discrete time variation and discussed the basic dynamics of the system. Banerjee and Venturino ^[9] proposed a phytoplankton-toxic phytoplankton-zooplankton model and found that toxic phytoplankton population does not cause zooplankton population to become extinct. Javidi and Ahmad ^[10] investigated the dynamics of a time fractional order toxic-phytoplankton-phytoplankton-zooplankton system(TPPZS) and did some numerical simulations to validate the theoretical results. Han and Dai ^[11] studied the spatiotemporal pattern caused by cross-diffusion of a toxic-phytoplankton-zooplankton model with nonmonotonic functional response and discussed the effect of toxin-producing rate of toxic-phytoplankton (TPP) species and natural mortality rate of zooplankton species on pattern selection. In addition, Han and Dai proposed a spatiotemporal pattern formation and selection driven by nonlinear cross-diffusion of a toxic-phytoplankton-zooplankton model with Allee effect in ^[12]. Zheng and Sugie ^[13] studied a three-dimensional system including phytoplankton, zooplankton and fish and gave a sufficient condition to ensure that the equilibrium of this three-dimensional system is globally asymptotically stable. Furthermore, they also proved that the equilibrium is asymptotically stable under relatively weak conditions.

Because the seabed has a variety of sediments that can be used as a refuge for the prey, some phytoplankton populations reduce the risk of being caught by zooplankton through refuge. Therefore, some scholars ^{[14,15,16,17,18]} have considered the prey refuge in the predator-prey model in recent years. Kar ^[14] proposed a prey-predator model incorporating a prey refuge and studied the influence of prey refuge on the prey-predator model. Chen and Chen ^[15] investigated a predator-prey model with Holling type Ⅱ functional response incorporating a constant prey refuge and gave the basic dynamical behaviors of the model. Tripathi et al. ^[16] studied a prey-predator model with reserved area and found that the predator species would exist if the value of prey reserve does not exceed a threshold value, beyond which the predator species would become extinct. Ghosh et al. ^[17] analyzed the effect of additional food for predator on the dynamical behaviors of a prey-predator model with prey refuge, they suggested that the possibility of extinction of predator species in high prey refuge ecosystems can be eliminated by providing additional food to predator species. Samanta et al. ^[18] proposed a fractional-order prey-predator model with prey refuge and derived some sufficient conditions to guarantee the global asymptotic stability of predator-extinction equilibrium and co-existing equilibrium.

Li et al. ^[19] considered the factors that phytoplankton protects themselves from being eaten by zooplankton by releasing toxin and taking refuge, and proposed a toxic phytoplankton-zooplankton model with refuge

where $P(t)\; \text{and}\; Z(t)$ denote population size of phytoplankton and zooplankton at time $t$, respectively; $r$ is the intrinsic growth rate; $K$ is the environmental carrying capacity in the absence of zooplankton; $\beta_{1}$, $\beta_{2}$ $(0 < \beta_{2} < \beta_{1})$ and $d$ describe the predation rate, the ratio of biomass conservation and the natural death rate of zooplankton species, respectively; $a_{1}$ and $a_{2}$ are the half saturation constants; $\theta$ represents the rates of toxin production per phytoplankton species; $m$ is the constant refuge capacity of phytoplankton population, then $P-m$ denotes the number of unprotected phytoplankton is captured by zooplankton species. Their results ultimately showed that phytoplankton refuge and toxin have important effects on the occurrence and termination of algal blooms in freshwater lakes.

Actually, the choice of prey species growth function and predator species function response is considered to be the most important element in the prey-predator model. Generally speaking, many modelers choose Logistic growth form as prey species growth function without considering the predator species. However, we all know that the resources in an ecosystem are limited, such as space, food, basic nutrition and so on. Thus the average growth rate becomes a decreasing function of population size as the population size increases gradually. When the number of population reaches the environmental capacity $K$, the average growth rate decreases to zero; in addition, any population number above the value $K$ will have a negative growth rate. However, there are a lot of evidences that the low population density is the opposite ^{[20,21,22,23,24,25]}. This phenomenon is the so-called Allee effect, which is the positive density dependence of population growth at low density ^[22,26]. The main causes of the Allee effect are the lack of a spouse, the reduced vigilance against predators, the adjustment of the environment, the reduction of defenses against predators, and many other reasons ^[20,22]. In general, the Allee effect is expressed by an equation of the following form ^[27]

where $X(t)$ represents population size at time $t$; $r$ and $K$ are the intrinsic growth rate and the environmental carrying capacity, respectively; $K_{0}$ is the critical level of phytoplankton population. When the population density is below the critical level, the population growth rate will decline and population will tend to become extinct.

In recent years, more and more people consider time delay into population biological model ^{[28,29,30,31,32,33,34,35]} to study its influence on dynamical behavior of system. What's more, it is well known that most countries of the world achieve the economic benefits of natural resource management by restoring and maintaining the ecosystem's health, the productivity and biodiversity and the overall quality of life in a manner that integrates social and economic objectives. Of course, it also satisfies the need for humans to benefit from natural resources. The prey-predator systems with harvesting have been widely studied by a large number of scholars ^{[6,34,35,36,37,38,39,40]}. Meng et al. ^[40] studied a predator-prey system with harvesting prey and disease in prey species and found that the optimal harvesting effort is closely related to the incubation period of the infectious disease, and the maximum value of the optimal harvesting decreases with the increase of the time delay. Furthermore, Meng et al. ^[34] investigated a nutrient-plankton model with Holling type Ⅳ, delay, and harvesting, their results indicated that the nutrient increases first and then stabilizes as the harvesting time increases, the number of phytoplankton and zooplankton decreases and even stabilizes as the harvesting time increases. Zhang and Zhao ^[36] considered a diffusive predator-prey system with delays and interval biological parameters and believed that the overfishing can lead to species extinction.

Motivated by the work mentioned-above, we propose a delayed phytoplankton-zooplank-ton model with Allee effect and linear harvesting. The results of this paper can be seen as a complement to system (1.1). As far as we know, although a large number of scholars have studied the interaction between phytoplankton and zooplankton, system (1.1) including the Allee effect, time delay and linear harvesting has not been studied yet.

The content of this paper is organized as follows. A delayed phytoplankton-zooplankton model with Allee effect and linear harvesting is described in section 2. In this section, we give the boundedness of the model. Section 3 demonstrates the existence and stability of the equilibria of system (2.1). A detailed discussion of the Hopf bifurcation in section 4. We not only give the existence and property of Hopf bifurcation when the system has no time delay, but also give the existence and property of Hopf bifurcation when the system has time delay. In section 5, the optimal policy is derived by using Pontryagin's Maximum Principle. In section 6, some numerical simulations are given for illustrating the theoretical results. The problem ends with a brief concluding remark.

2.   Model and its basic properties

In this paper, based on the work of the reference ^[19], we propose a delayed phytoplankton-zooplankton system with Allee effect and linear harvesting described by

where $K_{0}$ ($0 < K_{0} < < K$) is the critical level of the growth of phytoplankton; $q_{1}$ and $q_{2}$ are catchability coefficients of the two species; $E$ represents the harvesting effort. The delay $\tau$ in system (2.1) can be regarded as the maturation period of zooplankton species. For biological significance, system (2.1) must satisfy the following initial conditions

where $(\vartheta_{1}(\theta), \vartheta_{2}(\theta))\in C([-\tau, 0], R_{+}^{2})$ and $C([-\tau, 0], R_{+}^{2})$ represents the Banach space formed by all continuous functions from $[-\tau, 0]$ to $R_{+}^{2}$, here $R_{+}^{2} = \{(x_{1}, x_{2}):x_{1} > m, x_{2} > 0\}$.

When $\tau = 0$, system (2.1) becomes

and system (2.2) must satisfy the following initial conditions

As we all know, the boundedness of the model can ensure that the model has good dynamical behaviors. So we have the following result when $\tau = 0$.

Lemma 2.1. All solutions of system (2.2) with initial conditions $P(0) > m, Z(0) > 0$ that start in $\Theta = \{(P, Z)\mid (m, +\infty)\times(0, +\infty)\} \subset R_{+}^{2}$ are uniformly bounded for all $t > 0$.

Proof. Let $(P(t), Z(t))$ be any solution of system (2.2) with positive initial conditions $P(0) > m$, $Z(0) > \; 0$. Notice that there is $\Omega_{1} = \{(P, Z)\mid K\geq P(t) > m, Z(t) > 0\text{ for all } t > 0\}$ for system (2.2). So we have $P(t)\leq K$ for all $t > 0$.

We define $W(P, Z) = cP(t)+Z(t)$, here $c = \frac{\beta_{2}}{\beta_{1}}$, thus

where $\alpha = \min\{r, \frac{r(K+K_{0})}{KK_{0}}\}$.

That is

So,

When $t\rightarrow\infty$, we have

Therefore, all solutions of system (2.2) enter into the invariant set $\Omega = \Big\{(P, Z)\in \Theta\Big|0 < W\leq\frac{2cr(K+K_{0})}{\alpha K_{0}}\Big\}$. The Lemma 2.1 is proved.

Next, we use the method in ^[41] to give the positiveness of solution of system (2.1) in the case of $\tau > 0$.

Lemma 2.2. All solutions of system (2.1) with initial conditions that start in $\Theta = \{(P, Z)\mid (m, +\infty)\times(0, +\infty)\} \subset R_{+}^{2}$ are positive invariant.

Proof. We consider $(P, Z)$ a noncontinuable solution of system (2.1), defined on $[-\tau, \Gamma)$, where $\Gamma\in\; (0, \infty]$. We have to prove that for all $t\in[0, \Gamma)$, $P(t) > m$, and $Z(t) > 0$. Suppose that is not true. Then, there exists $0 < T < \Gamma$ such that for all $t\in[0, T)$, $P(t) > m$, and $Z(t) > 0$ and either $P(T) = m$ or $Z(T) = 0$. For all $t\in[0, T)$, under initial conditions we have

As $(P, Z)$ is continuous on $[-\tau, T]$, there exists a $M\geq0$ such that for all $t\in[-\tau, T]$,

Taking the limit, as $t\rightarrow T$, we get

which contradicts the fact that either $P(T) = m$, or $Z(T) = 0$. Thus, for all $t\in[0, \Gamma)$, $P(t) > m$ and $Z(t) > 0$. The Lemma 2.2 is proved.

3.   The existence and stability of equilibria

3.1.   The existence of equilibria

In order to get the conditions for the existence of the equilibria of system (2.2), we analyze the following Eq (3.1) which is given by

Obviously, the equilibria of system (2.2) are the intersections of the two equations of (3.1). In the absence of zooplankton, that is $Z = 0$, the first equation of (3.1) becomes

Let $\Delta_{1} = \frac{r^{2}(K_{0}+K)^{2}-4K_{0}Kr(r+q_{1}E)}{K_{0}^{2}K^{2}}$, then there is the following conclusion.

Theorem 3.1. The boundary equilibria of system (2.2) are as follows.

(i) If $\Delta_{1} = 0$, i.e., $E = \frac{r(K_{0}-K)^{2}}{4K_{0}Kq_{1}}$, then system (2.2) has a unique boundary equilibrium given by $\widetilde{E}_{1}(\widetilde{P}_{1}, 0)$, here $\widetilde{P}_{1} = \frac{K_{0}+K}{2}$.

(ii) If $\Delta_{1} > 0$, i.e., $E < \frac{r(K_{0}-K)^{2}}{4K_{0}Kq_{1}}$, then system (2.2) has two distinct boundary equilibria $\widetilde{E}_{2, 3}(\widetilde{P}_{2, 3}, 0)$, here $\widetilde{P}_{2, 3} = \frac{K_{0}+K}{2}\pm\frac{1}{2}\sqrt{(K_{0}+K)^{2}-\frac{4K_{0}K(r+q_{1}E)}{r}}$.

For positive equilibria, from the second equation of (3.1), we have

where $A_{1} = \beta_{2}-\theta-d-q_{2}E$, $A_{2} = \beta_{2}(a_{2}-m)+\theta(m-a_{1})+(d+q_{2}E)(m-a_{1}-a_{2}), A_{3} =$$a_{2}(m-a_{1})(d+q_{2}E)-\beta_{2}ma_{2}$.

Let $\Delta_{2} = A_{2}^{2}-4A_{1}A_{3}$ and $h(P) = A_{1}P^{2}+A_{2}P+A_{3}$. From the first equation of (3.1), we have

where $B_{1} = -r$, $B_{2} = r(m-a_{1}+K_{0}+K)$, $B_{3} = -[r(K_{0}+K)(m-a_{1})+K_{0}K(r+q_{1}E)]$, $B_{4} = K_{0}K(r+q_{1}E)(m-a_{1})$.

Since the signs of $A_{1}$, $A_{2}$ and $A_{3}$ are uncertain, there are many cases where $h(P) = 0$ has positive root(s). For the sake of discussion, let $A_{2} > 0$, i.e., $\beta_{2} < \frac{\theta(m-a_{1})+(d+q_{2}E)(m-a_{1}-a_{2})}{m-a_{2}}$. Furthermore, as long as

$m-a_{1}-a_{2} > 0$. Thus, when $A_{2} > 0$, $h(P) = 0$ will have possible positive root(s).

Through simple analyses, we obtain the following results about the existence of the positive equilibria.

Theorem 3.2. Under the assumption that $A_{2} > 0$ holds, the possible positive equilibria of system (2.2) are as follows.

(i) System (2.2) has one positive equilibrium when one of the following nine conditions is satisfied

(a)$\; \; A_{3} > 0, A_{1} > 0, \Delta_{2} = 0, h(m) > 0, h(K) > 0$;

(b)$\; \; A_{3} > 0, A_{1} > 0, \Delta_{2} > 0, h(m) < 0, h(K) > 0$;

(c)$\; \; A_{3} > 0, A_{1} > 0, \Delta_{2} > 0, h(m) > 0, h(K) < 0$;

(d)$\; \; A_{3} > 0, A_{1} < 0, \Delta_{2} > 0, h(m) > 0, h(K) < 0$;

(e)$\; \; A_{3} < 0, A_{1} = 0$;

(f)$\; \; A_{3} < 0, A_{1} > 0, \Delta_{2} > 0, h(m) < 0, h(K) > 0$;

(g)$\; \; A_{3} < 0, A_{1} < 0, \Delta_{2} = 0, h(m) < 0, h(K) < 0$;

(h)$\; \; A_{3} < 0, A_{1} < 0, \Delta_{2} > 0, h(m) > 0, h(K) < 0$;

(j)$\; \; A_{3} < 0, A_{1} < 0, \Delta_{2} > 0, h(m) < 0, h(K) > 0$.

(ii) System (2.2) has two distinct positive equilibria when one of the following two conditions is satisfied

(k)$\; \; A_{3} > 0, A_{1} > 0, \Delta_{2} > 0, h(m) > 0, h(K) > 0$;

(l)$\; \; A_{3} < 0, A_{1} < 0, \Delta_{2} > 0, h(m) < 0, h(K) < 0$.

In the case of (e) in Theorem 3.2, that is, when $\beta_{2} = \theta+d+q_{2}E$ and $\frac{(m-a_{1})(d+q_{2}E)}{m} < \beta_{2} < \frac{\theta(m-a_{1})+(d+q_{2}E)(m-a_{1}-a_{2})}{m-a_{2}}$ simultaneously hold, system (2.2) has a positive equilibrium $E^{*}(P^{*}, Z^{*})$, where

Remark 1. Due to the density of phytoplankton in system (2.2) is greater than the refuge constant $m$ at any time, so the existence of other possible positive equilibria of system (2.2) is complex, we will not list them here. Therefore, we will focus our discussion about system (2.2) at positive equilibrium $E^{*}(P^{*}, Z^{*})$.

3.2.   The local stability of equilibria

3.2.1.   The local stability of boundary equilibria

When $\tau = 0$, we only discuss the local stability of the boundary equilibrium $\widetilde{E}_{1}(\widetilde{P}_{1}, 0)$ and can use the similar methods to obtain the local stability of the other boundary equilibria.

Theorem 3.3. The boundary equilibrium $\widetilde{E}_{1}(\widetilde{P}_{1}, 0)$ of system (2.2) is locally asymptotically stable if and only if $(H_{1})$ holds, where

$(H_{1}):\gamma_{1} = K_{0}+K-2m > 0, r(K_{0}+K)^{2} > 4K_{0}K(r+q_{1}E+d+q_{2}E).$

Proof. The Jacobian matrix of system (2.2) at the boundary equilibrium $\widetilde{E}_{1}(\widetilde{P}_{1}, 0)$ is given by

Thus, we can get

where

From the above analysis, it is easy to know that if and only if $(H_{1})$ holds, we have $\text{tr}(J_{\widetilde{E}_{1}}) < 0$ and $\text{det}(J_{\widetilde{E}_{1}}) > 0$. Hence, the boundary equilibrium $\widetilde{E}_{1}(\widetilde{P}_{1}, 0)$ of system (2.2) is locally asymptotically stable if and only if $(H_{1})$ holds. This completes the proof.

3.2.2.   The local stability of the positive equilibrium

Here, we discuss the local stability of system (2.2) at the positive equilibrium $E^{*}(P^{*}, Z^{*})$ when $\tau = 0$.

Theorem 3.4. If $(H_{2})$ and $(H_{3})$ simultaneously hold, then the positive equilibrium $E^{*}(P^{*}, Z^{*})$ of system (2.2) is locally asymptotically stable, where

$(H_{2}):\frac{N_{1}+N_{2}+N_{3}+N_{4}+N_{5}+N_{6}}{L_{3}} < q_{1} < \frac{M_{1}+M_{2}+M_{3}+M_{4}+M_{5}+M_{6}+M_{7}M_{8}}{L_{1}+L_{2}},$

$(H_{3}):\gamma_{2} = 3P^{*}-2(K_{0}+K) > 0.$

Proof. The Jacobian matrix of system (2.2) at the positive equilibrium $E^{*}(P^{*}, Z^{*})$ is given by

where

Then characteristic equation of system (2.2) around $E^{*}(P^{*}, Z^{*})$ is

where

here

Through analysis and calculation, we can obtain $\text{tr}(J_{E^{*}}) < 0$ and $\text{det}(J_{E^{*}}) > 0$ if and only if $(H_{2})$ and $(H_{3})$ hold. Thus, we easily get that the positive equilibrium $E^{*}(P^{*}, Z^{*})$ of system (2.2) is locally asymptotically stable when $(H_{2})$ and $(H_{3})$ hold. This completes the proof.

4.   Hopf bifurcation of the positive equilibrium

4.1.   Hopf bifurcation of system without time delay

In this section, we investigate the existence of Hopf bifurcation around $E^{*}(P^{*}, Z^{*})$ of system (2.2). Taking prey refuge $m$ as the bifurcation parameter, the critical value of Hopf bifurcation is a positive root of $\text{tr}(J_{E^{*}}) = 0$, thus $m = m_{H}$ which satisfies $\text{det}(J_{E^{*}})_{m = m_{H}} > 0$. As we know, when the value of $m$ exceeds its critical value $m = m_{H}$, the stability of the positive equilibrium $E^{*}(P^{*}, Z^{*})$ will be changed (Figure 1).

In addition, we verify the transversality condition under which Hopf bifurcation occurs. By some simple calculations, we have that $\frac{\text{d}}{\text{d}m}\text{tr}(J_{E^{*}})_{m = m_{H}}\neq0$, which implies that the stability of the positive equilibrium $E^{*}(P^{*}, Z^{*})$ changes when the parametric restriction $\text{tr}(J_{E^{*}}) = 0$ and the transversality condition are satisfied simultaneously.

Therefore, we draw the following conclusion about the occurrence of Hopf bifurcation at the positive equilibrium $E^{*}(P^{*}, Z^{*})$.

Theorem 4.1. Under the condition that the positive equilibrium $E^{*}(P^{*}, Z^{*})$ exists, the stability of system (2.2) at the positive equilibrium $E^{*}(P^{*}, Z^{*})$ is changed through critical value $m = m_{H}$.

In order to facilitate the discussion of the direction of Hopf bifurcation, we use the related theory of ^[42] to calculate the first Lyapunov number $l_{1}$ at the positive equilibrium $E^{*}(P^{*}, Z^{*})$ of system (2.2).

We transform the positive equilibrium $E^{*}(P^{*}, Z^{*})$ of system (2.2) to the origin by translation $\bar{P} = P-P^{*}$ and $\bar{Z} = Z-Z^{*}$. Then, system (2.2) in a neighborhood of the origin can be obtained

where $m_{110}$, $m_{101}$, $m_{210}$ and $m_{201}$ are the elements of the Jacobian matrix at the equilibrium $E^{*}(P^{*}, Z^{*})$. So system (2.2) experiences the Hopf bifurcation at the equilibrium $E^{*}(P^{*}, Z^{*})$, we have $m_{110}+\; m_{201} = 0$ and $\Delta = \text{det}(J_{E^{*}}) = m_{110}m_{201}-m_{101}m_{210} > 0$. The coefficients $m_{kij}(k = 1, 2)$ are given by

and $R(\bar{P}, \bar{Z})$ and $S(\bar{P}, \bar{Z})$ are power series in $(\bar{P}, \bar{Z})$ with terms $\bar{P}^{i}$, $\bar{Z}^{j}$ satisfying $i+j\geq4$.

The calculation formula of the first Lyapunov number $l_{1}$ ^[42] determining the stability of the limit cycle in the planar system is described by

Based on the above analyses, we conclude the following conclusion about the direction of the limit cycle.

Theorem 4.2. If $l_{1} < 0$, then Hopf bifurcation is supercritical. If $l_{1} > 0$, then Hopf bifurcation is subcritical.

4.2.   Hopf bifurcation of system with time delay

In this section, we will discuss the effect of maturation period of zooplankton $\tau$ on the dynamical behavior of system (2.1).

4.2.1.   Local stability and the existence of Hopf bifurcation

The local stability of the positive equilibrium $E^{*}(P^{*}, Z^{*})$ and the existence of Hopf bifurcation are studied by considering the maturity delay $\tau$ as bifurcation parameter. First, system (2.1) is linearized at positive equilibrium $E^{*}(P^{*}, Z^{*})$, we obtain

where

It is easy to see that corresponding characteristic equation of system (4.2) at positive equilibrium $E^{*}(P^{*}, Z^{*})$ is

where $A = -a_{11}-a_{22}, B = a_{11}a_{22}-a_{12}a_{21}, C = -a_{12}b_{21}$.

When $\tau = 0$, Eq (4.3) becomes

then Eq (4.4) is equivalent to Eq (3.3). Hence, the positive equilibrium $E^{*}(P^{*}, Z^{*})$ is locally asymptotically stable when $(H_{2})$ and $(H_{3})$ hold.

When $\tau\neq0$, in order to get the existence of Hopf bifurcation, let $\lambda = i\omega_{0}(\omega_{0} > 0)$ be a root of Eq (4.3) and substitute it into Eq (4.3), then we have

By separating the real and imaginary parts, we obtain

From which, we have

Let $u = \omega_{0}^{2}$, then Eq (4.6) deduces to

Record Eq (4.7) as $f(u) = u^{2}-(2B-A^{2})u+B^{2}-C^{2}$. From Eq (4.6), it follows that if

$(H4):2B-A^{2} > 0, \; \; B^{2}-C^{2} > 0, \; \; \text{and}\; \; A^{4}-2A^{2}B+4C^{2} = 0$

holds, then Eq (4.6) has a unique positive root $\omega_{0}^{2}$. Substituting $\omega_{0}^{2}$ into Eq (4.5), we obtain

If

$(H5):2B-A^{2} > 0, \; \; B^{2}-C^{2} > 0, \; \; \text{and}\; \; (2B-A^{2})^{2} > 4(B^{2}-C^{2})$

holds, then Eq (4.6) has two different positive roots $\omega_{+}^{2}$ and $\omega_{-}^{2}$. We substitute $\omega_{\pm}^{2}$ into Eq (4.5), then

Let $\lambda$ be the root of Eq (4.3) satisfying $\text{Re}\lambda(\tau_{2n}) = 0$ (rep.$\text{Re}\lambda(\tau_{2k}^{\pm}) = 0$) and $\text{Im}\lambda(\tau_{2n}) = \omega_{0}$ (rep.$\text{Im}\lambda(\tau_{2k}^{\pm}) = \omega_{\pm}$). Then, when $\tau_{0} = \min\{\tau_{2n}, \tau_{2k}\}, n, k = 0, 1, 2\cdot\cdot\cdot$, we have the following conclusion.

Theorem 4.3. If $f^{'}(u) > 0$, since $\frac{\text{d}\text{Re}(\lambda)}{\text{d}\tau}\Big|_{\tau = \tau_{0}}$ and $f^{'}(u)$ have the same sign, then $\frac{\text{d}\text{Re}(\lambda)}{\text{d}\tau}\Big|_{\tau = \tau_{0}} > 0$ holds.

Proof. From Eq (4.3), we can obtain

Notice that $Ce^{-\lambda\tau} = -(\lambda^{2}+A\lambda+B)$. Thus, we obtain that

Substituting $\lambda = i\omega_{0}$ into (4.8), we have

Therefore,

If $f^{\prime}(\omega_{0}^{2})\neq0$, we have $\frac{\text{d}\text{Re}(\lambda)}{\text{d}\tau}\Big|_{\tau = \tau_{0}}\neq0$, then $\frac{\text{d}\text{Re}(\lambda)}{\text{d}\tau}\Big|_{\tau = \tau_{0}} > 0$. If $\frac{\text{d}\text{Re}(\lambda)}{\text{d}\tau}\Big|_{\tau = \tau_{0}} < 0$, then Eq (4.3) has the roots of the real part when $\tau < \tau_{0}$. This contradicts that the positive equilibrium $E^{*}(P^{*}, Z^{*})$ is locally asymptotically stable when $\tau\in(0, \tau_{0})$.

According to the Hopf bifurcation theorem given in ^[43], we can get the following result on stability and bifurcation of system (2.1).

Theorem 4.4. For system (2.1), suppose $(H4)$ and $(H5)$ are satisfied. Then the positive equilibrium $E^{*}(P^{*}, Z^{*})$ is locally asymptotically stable when $\tau\in(0, \tau_{0})$, and system (2.1) undergoes a Hopf bifurcation at $E^{*}(P^{*}, Z^{*})$ when $\tau$ passes through its critical value $\tau = \tau_{0}$.

4.2.2.   Direction and stability of Hopf bifurcation

In the previous section, we have already obtained sufficient conditions to guarantee that system (2.1) has a periodic solution at positive equilibrium $E^{*}(P^{*}, Z^{*})$ when the critical value $\tau = \tau_{0}$. In this section, we will use the normal form theory and the center manifold theorem presented in Hassard et al. ^[43] to obtain the properties of the bifurcating periodic solution.

Let $\tau = \tau_{0}+\mu$, $\mu\in R$. Then $\mu = 0$ is a Hopf bifurcation value of system (2.1). Without loss of generality, we can choose the phase space as $C = C([-\tau, 0], \Theta)$. Let $u(t) = (u_{1}(t), u_{2}(t))^{T}\in \Theta$, $u_{t}(\theta) = u(t+\theta)\in C$, system (2.1) can be rewritten as

where $L_{\mu}:C\rightarrow \Theta$ and $f:R\times C\rightarrow \Theta$ are defined by

where

here

Based to the Reisz representation theorem, there is a bounded variation function $\eta(\bullet, \mu):[-\tau, 0]\rightarrow\; \Theta^{2\times2}$ such that

Then, we have a choice that

where $\delta$ is the Dirac function. For $\varphi\in C^{1}([-\tau, 0], \Theta^{2})$, we define the operator $A(\mu)$ as

Further, we can define the operator $R(\mu)$ as

Then, we can rewrite Eq (4.9) as the following operator equation

For $\psi\in C^{1}([-\tau, 0], \Theta^{2})$, the adjoint operator $A_{0}^{*}$ of $A_{0}$ can be defined as

where $\eta^{T}$ is the transpose of $\eta$, and the domains of $A_{0}$ and $A_{0}^{*}$ are $C^{1}([-\tau, 0], \Theta^{2})$ and $C^{1}([0, \tau], \Theta^{2})$, respectively.

For $\varphi\in C^{1}([-\tau, 0], \Theta^{2})$ and $\psi\in C^{1}([0, \tau], \Theta^{2})$, a bilinear form is given by

where $\eta(\theta) = \eta(\theta, 0)$. Then $A_{0}$ and $A_{0}^{*}$ are a pair of adjoint operators. By the results obtained in the last section, we know that $\pm i\omega_{0}$ are the eigenvalues of $A_{0}$. Thus, they are also eigenvalues of $A_{0}^{*}$. What's more, we assume that $q(\theta)$ is the the eigenvalues of $A_{0}$ corresponding to $i\omega_{0}$, then

By Eq (4.10), Eq (4.13) becomes

From Eq (4.14) we can get

where $V = (v_{1}, v_{2})^{T}\in \Theta^{2}$ is a constant vector. By virtue of Eq (4.15), we may obtain

which yields

On the other hand, if $-i\omega_{0}$ is the eigenvalue of $A_{0}^{*}$, then we have

For the non-zero vector $q^{*}(\xi)$, $\xi\in[0, \tau_{0}]$, we can obtain

Let $q^{*}(\xi) = DV^{*}e^{i\omega_{0}\xi}$, where $\xi\in[0, \tau_{0}]$, and $V^{*} = (v^{*}_{1}, v^{*}_{2})^{T}$ be a constant vector. Similarly, we get

In fact, by Eq (4.12), we have

Actually, we can choose

From Eq (4.16) we obtain that $\langle q^{*}, q\rangle = 1$. More importantly

thus, $\langle q^{*}, \bar{q}\rangle = 0$ is easily verified.

Based on the above analysis, we can come to the following simple conclusion.

Theorem 4.5. Let $q(\theta) = Ve^{i\omega_{0}\theta}$, where $\theta\in[-\tau, 0]$, is the eigenvector of $A_{0}$ corresponding to $i\omega_{0}$, and let $q^{*}(s) = DV^{*}e^{i\omega_{0}s}$, where $s\in[0, \tau]$ is the eigenvector of $A_{0}^{*}$ corresponding to $-i\omega_{0}$. Then

where

Next, we use the same notations as the previous part to study the stability of bifurcating periodic solution. We shall calculate the coordinates to express the center manifold $C_{0}$ at $\mu = 0$. Let $u_{t}$ be the solution of $\overset{.}u(t) = A_{\mu}u_{t} + R_{\mu}u_{t}$ at $\mu = 0$ and define

On the center manifold $C_{0}$, we have

and

where $z$ and $\bar{z}$ are the local coordinates for the center manifold $C_{0}$ in the direction of $q^{*}$ and $\bar{q^{*}}$. Notice that $W$ is real if $u_{t}$ is real. So here we only consider the real solutions.

For the solution $u_{t}\in C_{0}$ of Eq (4.11), we have

Then Eq (4.19) can be rewritten as

where

According to Eqs (4.11), (4.17) and (4.18), we may obtain

We can rewrite Eq (4.22) as

where

Substituting the corresponding series into Eq (4.22) and comparing the coefficients, we have

According to Eqs (4.19) and (4.20), we know

Comparing the coefficients with Eq (4.21), it follows that

where

Since $W_{20}(\theta)$ and $W_{11}(\theta)$ appear in $g_{21}$, we need to compute them further. From Eqs (4.22) and (4.23), we have that for $\theta\in[-\tau, 0)$,

Comparing the coefficients with Eq (4.24) gives that

From Eqs (4.10) and (4.25), we have

and

Solving for $W_{20}(\theta)$ and $W_{11}(\theta)$, we obtain

and

where $E_{1} = (E^{(1)}_{1}, E^{(2)}_{1})$ and $E_{2} = (E^{(1)}_{2}, E^{(2)}_{2})$ are the following two-dimensional constant vectors, and can be determined by setting $\theta = 0$ in $H$.

In fact, the following formula is true at $\theta = 0$,

we have

and

Substituting Eqs (4.28) and (4.29) into Eqs (4.26) and (4.27), respectively, we get

and

Substituting Eq (4.28) into Eq (4.32), we have

Since $i\omega_{0}$ is the eigenvalue of $A_{0}$ and $q(0)$ is the corresponding eigenvector, we get

that is

and

By Eqs (4.30) and (4.34)–(4.36), we obtain

which leads to

where

Similarly, substituting Eqs (4.29) and (4.31) into Eq (4.33), we have

which leads to

where

Through the above calculation and analysis, it is not difficult to see that each $g_{ij}$ depends on the parameters and time delay of system (2.1). Thus, we can obtain the following expressions

which determine the direction of Hopf bifurcation and stability of bifurcating periodic solutions of system (2.1) on the center manifold at the critical values $\tau_{0}$. From the conclusion of Hassard et al. ^[43], we summarize the following main findings.

Theorem 4.6. The values of the parameters $\mu_{2}$, $\beta$ and $T_{2}$ of (4.37) will determine the properties of Hopf bifurcation of system (2.1).

(i) The sign of $\mu_{2}$ determines the direction of Hopf bifurcation: The Hopf bifurcation is supercritical if $\mu_{2} > 0$ and the Hopf bifurcation is subcritical if $\mu_{2} < 0$.

(ii) The sign of $\beta$ determines the stability of bifurcating periodic solutions: The periodic solutions are stable if $\beta < 0$ and unstable if $\beta > 0$.

(iii) The sign of $T _{2}$ determines the period of bifurcating periodic solutions: The period increases if $T _{2} > 0$ and the period decreases if $T _{2} < 0$.

5.   Optimal harvesting policy

When $\tau = 0$, let $\hat{c}$ be the constant harvesting cost per unit effort and $p_{1}$, $p_{2}$ are the constant price per unit biomass of phytoplankton and zooplankton, respectively. Then the net economic revenue to the society is given by

In order to study the optimal harvesting yield of system (2.2), we will consequently maximize the full return from resource management. Hence we take the harvesting effort $E$ as a control variable and consider a objective function defined by the present value $J$ of a continuous time-stream of revenues

where

where $\delta$ is the instantaneous annual discount rate. We consider that the present value of capital flow over time depends on discount rate $\delta$. Here $X = [P, Z]^{T}$ is the vector of state variables, it can be written as $\dot{X} = f(X, E)$, $X(0) = X_{0}$ by according to the state equations of (2.2).

Our optimal control problem is to maximize Eq (5.1) subject to the state equations of (2.2) and to the control constraint

where $E_{\text{max}}$ is the maximum harvesting effort.

Based on the Pontryagin's Maximum Principle ^[44], we can obtain the optimal solution of this problem. We know that the convexity of objective function with respect to $E(t)$, the linearity of the control differential equations and the compactness of the range values of the state variables to ensure the existence of optimal control.

The present value Hamiltonian function of the optimal problem is formulated by

where $\lambda(t) = [\lambda_{1}(t), \lambda_{2}(t)]^{T}$ is the vector of constant or adjoint variables.

Next, we substitute $\lambda^{T}(t)$ and $f(X, E)$ into above Hamiltonian function. Then, the Hamiltonian function becomes

To make $H$ reach the maximum on the control set $0\leq E(t)\leq E_{\text{max}}$, the condition that the Hamiltonian function $H$ must satisfy is presented by

Pontryagin's Maximum Principle ^[44] states that the optimal state trajectory, optimal control, and corresponding adjoint variable vector must satisfy the following adjoint equations

Obviously, the adjoint equations are

For the positive optimal solutions, which satisfy $\dot{P} = \dot{Z} = 0$ (in other words, $P$, $Z$ are not dependent on $t$) and from two equations of system (2.2), we have

From the above analysis, it is obvious that $E$ is also independent of $t$. Furthermore, we get

and

Differentiating Eq (5.2) and replacing value of $\dot{\lambda_{1}}$, $\dot{\lambda_{2}}$, we get

By Eqs (5.2) and (5.6), we can get

where

Now removing $E$ from Eqs (5.3) and (5.4), we obtain

which is the optimal trajectory of the steady state given by the optimal solutions $P = P_{\delta}$, $Z = Z_{\delta}$. Then, we substitute $\lambda_{1}$ and $\lambda_{2}$ into Eq (5.5) and obtain optimal equilibrium level of effort given by

By solving Eqs (5.9) and (5.10) when assigning a certain value to $\delta$, we can obtain the optimal equilibrium level $(P_{\delta}, Z_{\delta})$.

The optimal harvesting effort at any time is determined by

where $E_{\text{min}}$ is the minimum harvesting effort.

6.   Numerical simulations

In this section, we will do some numerical simulations to verify the theoretical results. The values of all parameters in system (2.2) are sourced from Table 1. And the initial values of the system (2.2) are assumed to be $P(0) = 400$, $Z(0) = 800$.

First, according to the case (e) in Theorem 3.2, there is a positive equilibrium $E^{*} = (343.2, 1004.4)$ of system (2.2). In addition, if assumptions $(H2)$ and $(H3)$ hold, the positive equilibrium $E^{*}$ is locally asymptotically stable (Figure 2). We can clearly see that the solution curve starting from different initial values eventually tends to the black point, that is, the positive equilibrium $E^{*} = (343.2, 1004.4)$ (Figure 2b).

Second, for the parameter values given above, we obtain $m_{H} = 286.6$ and the first Lyapunov number $l_{1} = -0.0187 < 0$ by simple calculation. This indicates that system (2.2) has a stable limit cycle around the positive equilibrium $E^{*}$. This result is shown in Figure 3.

Next, for the given parameters, we get $\omega_{0} = 0.2354$ and $\tau_{0} = 0.011$ when $\tau\neq0$. According to Theorem 4.4, we obtain that the positive equilibrium $E^{*} = (343.2, 1004.4)$ is locally asymptotically stable when $\tau = 0.01 < \tau_{0} = 0.011$ (Figure 4). Then, we choose the value of $\tau$ as $\tau = 0.1 > \tau_{0} = 0.011$ and obtain that $C_{1}(0) = -0.0019 - 0.0152i$, $\mu_{2} = 6.4946\times 10^{-21} > 0$, $\beta = -0.0038 < 0$, $T_{2} = 0.0065 > 0$ by using formula (4.37). On the basis of Theorem 4.6, we know that the system (2.1) experiences a supercritical Hopf bifurcation when $\tau$ passes its critical value $\tau_{0}$. Other than this, system (2.1) has stable bifurcating periodic solutions and the period of the bifurcating periodic solutions is increasing. We clearly see from the Figure 5 that the positive equilibrium $E^{*} = (343.2, 1004.4)$ is destabilized through a Hopf bifurcation.

Finally, we consider the following parameter values: $p_{1} = 1$, $p_{2} = 5$, $\hat{c} = 0.1$, $\delta = 0.2$ and the other parameters remain unchanged. Figure 6 shows the solution curve of the state variables with time. Figure 7a, b show the variation curve of the adjoint variables $\lambda_{1}$ and $\lambda_{2}$, respectively. It is easy to see from the Figure 7 that the adjoint variables $\lambda_{1}$ and $\lambda_{2}$ tend ultimately to $0$ with the increase of time. In addition, the effect of the constant refuge capacity of phytoplankton population $m$, the critical value of the growth of phytoplankton $K_{0}$, the half saturation constant $a_{1}$ and the half saturation constant $a_{2}$ on the optimal harvesting effort in Figure 8. It is not difficult to see that the optimal harvesting effort decreases as $m$ increases (Figure 8a), but the optimal harvesting effort increases as $K_{0}$, $a_{1}$ and $a_{2}$ increase (Figure 8b–d).

7.   Conclusions

In this paper, we studied the dynamics of the phytoplankton-zooplankton system in which the growth of phytoplankton is affected by Allee effect and the growth of zooplankton is affected by maturation delay. For the positive equilibrium $E^{*}$, due to the expression of the trace and determinant of its Jacobian matrix are very complicated, so the stability of the positive equilibrium $E^{*}$ is verified by the combination of theoretical derivation and numerical simulation. When the maturation delay of zooplankton is not considered, the strict mathematical proof of the existence of Hopf bifurcation is given by using the relevant bifurcation theory. In addition, we derive the expression of the first Lyapunov number $l_{1}$ that determines the direction of the Hopf bifurcation. Furthermore, when considering the maturation period of phytoplankton, we obtain some properties of the Hopf bifurcation through the normal form theory and the center manifold theorem. Because plankton has certain economic significance, we also consider linear harvesting for both phytoplankton and zooplankton and obtain the optimal harvesting policy by Pontryagin's Maximum Principle in this paper.

Comparing with the study on the phytoplankton-zooplankton system in ^[19], our model is more realistically by considering Allee effect and maturation delay into growth of phytoplankton and zooplankton, respectively. What's more, both populations are linearly harvested. For phytoplankton, we can predict the stability of system (2.2) at the positive equilibrium $E^{*}$ by selecting the parameter values to determine the critical value $m_{H}$ of the refuge capacity $m$. When we choose $m > m_{H}$, system (2.2) is stable at the positive equilibrium $E^{*}$, but when $m < m_{H}$, system (2.2) experiences a supercritical Hopf bifurcation and loses stability. This indicates that increasing the refuge capacity $m$ is important for balancing ecosystem. For zooplankton, we regard the maturation delay of zooplankton as a parameter to predict the stability of system (2.1). The positive equilibrium $E^{*}$ of the system (2.1) is stable, if the maturation delay of zooplankton $\tau < \tau_{0}$, which implies that the densities of phytoplankton and zooplankton will tend to be stable constants, indicating the ecosystem is balanced. But system (2.1) undergoes a Hopf bifurcation at the positive equilibrium $E^{*}$ if the the maturation delay of the zooplankton $\tau > \tau_{0}$, which shows that the densities of phytoplankton and zooplankton species will oscillate periodically.

Through the above analysis, we know that if the refuge capacity $m < m_{H}$ of phytoplankton or the maturation delay $\tau > \tau_{0}$ of zooplankton, the densities of phytoplankton and zooplankton will change periodically, which indicates the system becomes unstable. If the refuge capacity $m > m_{H}$ of phytoplankton or the maturation delay $\tau < \tau_{0}$ of zooplankton, the densities of phytoplankton and zooplankton will tend to be stable, which suggests a state of ecosystem balance. Thus, we can maintain the ecological balance by adjusting the values of the refuge capacity $m$ of phytoplankton and the maturation delay $\tau$ of zooplankton determined by system parameters.

Acknowledgments

The authors are grateful to the anonymous reviewers for their constructive comments. This work is supported by the National Natural Science Foundation of China (Grant Nos. 11661050 and 11861044), and the HongLiu First-class Disciplines Development Program of Lanzhou University of Technology.

Conflict of interest

All authors declare no conflicts of interest in this paper.