#### Abstract

In this paper, a Holling type IV nutrient-plankton model with time delay and linear plankton harvesting is investigated. The existence and local stability of all equilibria of model without time delay are given. Regarding time delay as bifurcation parameter, such system around the interior equilibrium loses its local stability, and Hopf bifurcation occurs when time delay crosses its critical value. In addition, the properties of the bifurcating periodic solutions are investigated based on normal form theory and center manifold theorem. What is more, the global continuation of the local Hopf bifurcation is discussed by using a global Hopf bifurcation result. Furthermore, the optimal harvesting is obtained by the Pontryagin’s Maximum Principle. Finally, some numerical simulations are given to confirm our theoretical analysis.

#### 1. Introduction

In the aquatic ecosystem, all aquatic food chains are based on plankton, which are composed of phytoplankton and zooplankton. As pointed out in [1], phytoplankton in particular have great contribution for our earth; for example, they provide food for marine life and oxygen for human life, and they also absorb more than half the carbon dioxide which may be contributing to global warming.

In the past decade, the issues of phytoplankton-zooplankton models have been investigated by many experts [2–13]. Chattopadhy et al. [2] showed the effect of the delay forming the period of toxic substances on blooms. Wang et al. [4] considered a toxic phytoplankton-zooplankton model with delay-dependent coefficient and investigated influence of delay and harvesting for the system. Yang et al. [8] discussed the impact of reaction-diffusion on the dynamic behaviours of plankton. Sharma et al. [11] studied the bifurcation behaviour analysis of a plankton model with multiple delays. Nevertheless, these models mainly included interactions of phytoplankton-zooplankton.

However, nutrients play an significant role in the growth of plankton. The concentration of nutrients can better explain their mechanism in the phytoplankton-zooplankton system. As far as we know, the nutrient-plankton systems have been considered by a number of scholars [14–21]. Chakrabortyour et al. [14] demonstrated that the toxin producing phytoplankton (TPP) species control the bloom of other nontoxic phytoplankton species when nutrient-deficient conditions are favorable for TPP species to release toxic chemicals. Zhang and Wang [15] discussed Hopf bifurcation and bistability of nutrient-phytoplankton-zooplankton model. Fan et al. [16] proposed a nutrient model for a water ecosystem. Chakraborty et al. [17] investigated spatial dynamics of a nutrient-phytoplankton system with toxic effect on phytoplankton. Wang et al. [18] analyzed a nutrient-plankton model with delayed nutrient cycling. Jang and Allen [19] studied deterministic and stochastic nutrient-phytoplankton-zooplankton models with periodic toxin producing phytoplankton. Das and Ray [20] debated the effect of delay on nutrient cycling in phytoplankton-zooplankton interactions in estuarine system. Sharma et al. [21] discussed the effect of the discrete time delay on the nutrient-plankton interaction and proposed the model as follows:where , , and denote the concentration of nutrient, phytoplankton, and zooplankton population at time , respectively. is the time delay which is incorporated with the assumption that the liberation of toxin is not instantaneous; rather it is mediated by some time lag. Li and Xiao [22] showed a predator-prey model with Holling type IV: , where the predator is specifically not taking nutrient, and the function form of biomass conversion of herbivore is Holling type IV function. They also assumed that the extra mortality in zooplankton due to toxic substance term is also expressed by Holling type IV function.

Generally speaking, there are some economic benefits for some plankton in aquatic food chain. Some phytoplankton such as phylum Diatom, phylum Pyrrophyta, and phylum Chrysophyta can be harvested as bait for aquaculture, and some zooplankton such as jellyfish, krill, and* Acetes* can be harvested for human food. Thus, these plankton play a crucial role in marine conservation and fishery management. As we all know, harvesting is an essential element in our life. On the one hand, the harvesting is one of the economic sources of human economy. On the other hand, it is a method to control the balance of ecosystem. Thus, zooplankton-phytoplankton system with harvesting has attracted the attention of scholars [23–29]. Chakraborty and Das [23] proposed a two-zooplankton one-phytoplankton system with harvesting and analyzed the impact of harvesting on the coexistence and competitive exclusion of competitive predators. Lv et al. [24] proposed a phytoplankton-zooplankton system with harvesting and pointed out that overharvesting could lead to population extinction and appropriate harvesting should guarantee the sustainability of the population types in an aquatic ecosystem.

Some works have already been made to understand the nutrient-plankton system with time delay and toxic phytoplankton, but the study of nutrient-plankton system which contains Holling type IV function, time delay, toxic phytoplankton, and linear harvesting is rarely done so far. Taking account of the effect of the time delay in maturation of zooplankton into nutrient-plankton model can make the model more realistic. The main aim of the present study is to see whether the time delay in maturation of zooplankton and harvesting of plankton can exert an impact on the nutrient-plankton system.

Based on the works [21, 22] and the above discussions, the following model is proposed by introducing harvesting and time delaywhere , , and denote the number of nutrient, phytoplankton, and zooplankton population at time , respectively. is the constant input of nutrient concentration and is washout rate. and are the nutrient uptake rate and conversion rate of nutrient for growth of the phytoplankton, respectively. and are the mortality rates of the phytoplankton and zooplankton population, respectively. and are the maximal ingestion rate and conversion rate of zooplankton. denotes the rate of toxin liberation. is the nutrient recycle rate after death of phytoplankton population. The parameter is the half saturation constant for a Holling type IV function. denotes the time delay in maturation of zooplankton. The constants and are the catchability coefficients of the phytoplankton and zooplankton population, and is the harvesting effort.

The initial conditions of system (2) take the following formwhere , the Banach space of continuous functions mapping the interval into , where .

The organisation of this paper is as follows: in Section 2, preliminary dynamical results of system without delay are analyzed, such as the existence of equilibria, boundness, and the local stability of equilibria. In Section 3, with the presence of time delay, the local stability and the existence of Hopf bifurcation of system (2) are investigated. Based on normal form theory and center manifold theorem, direction of Hopf bifurcation and stability of the bifurcating periodic solutions are also discussed. What is more, the global existence of bifurcating periodic solution is given by using a global Hopf bifurcation result due to Wu [30]. In Section 4, the optimal control problem is formulated by Pontryagin’s Maximum Principle. In Section 5, numerical simulations are given to illustrate the results.

#### 2. Preliminary Results of System without Delay

##### 2.1. Positivity and Boundedness

When , system (2) isThe initial conditions for system (4) take the following form

Lemma 1. *If , solutions of system (4) with initial conditions are defined on and remain positive for all .*

*Proof. *From the first equation of system (4), we obtainHence, we obtainFrom the second equation of system (4), we haveSince , the third equation of system (4) shows that by a standard comparison argument

Theorem 2. *All solutions of system (4) with initial conditions are bounded for all .*

*Proof. *In order to check the boundedness of solutions of system (4), we define the function As derivative of with respect to system (4), we obtainwhere . Therefore, by the comparison theory [31], we haveThus, , and are bounded.

##### 2.2. Equilibria and Their Local Stability

In this subsection, we will discuss the existence of equilibria of system (4) and their local asymptotical stability.

The given system (4) has three nonnegative equilibria, defined as the plankton extinction equilibrium , the zooplankton free equilibrium , and the coexistence equilibrium .

The zooplankton free equilibrium is , whereSince , if , then there is only one possible equilibrium involving phytoplankton but not zooplankton.

The coexistence equilibrium is , where , , and satisfy the following equationsFrom the third equation of (14), we haveIf , and , then (15) has exactly two positive real rootsAnd if , and , then (15) has two equal positive real rootsAgain from the first and second equations of (14), we haveIf , , then . If , , then . And if , , then .

In order to discuss the existence and the local stability of equilibria of system (4), some assumptions are given.

: , , and − (or − ).

: , and − , / − .

: , and .

: /, /, < < /, /, and .

: , /, and .

and .

From the above analysis, we obtain the following results.

Theorem 3. *(i) The plankton extinction equilibrium always exists.**(ii) The zooplankton free equilibrium exists if .**(iii) The coexistence equilibrium (or ) exists if holds; the coexistence equilibria and exist if holds; the coexistence equilibrium exists if () holds.*

According to the results of Theorem 3 and the above assumptions, we have the results as follows.

Theorem 4. *(1) If , then the plankton extinction equilibrium is locally asymptotically stable.**(2) If and , then the zooplankton free equilibrium is locally asymptotically stable.**(3) The coexistence equilibrium (or ) is locally asymptotically stable if conditions and hold. The coexistence equilibria and are locally asymptotically stable if conditions and hold. The coexistence equilibrium is locally asymptotically stable if conditions and hold.*

*Proof. *(1) The characteristic equation on is givenIt is clear that (19) has negative roots , and . So if , thenThus, the plankton extinction equilibrium is locally asymptotically stable. (2) The characteristic equation about is given byIf , thenFurther, the other roots and satisfy that . If , then . Therefore, if and , the zooplankton free equilibrium is locally asymptotically stable.

(3) The characteristic equation about is givenwhereSince , if () and () hold, then , . Therefore all roots of (23) have negative real parts. Based on Routh-Hurwitz criterion we obtain that the coexistence equilibrium (or ) is locally asymptotically stable.

Similarly, if and hold, the coexistence equilibria and are locally asymptotically stable; if and hold, the coexistence equilibrium is locally asymptotically stable.

*Remark 5. *From biological point of view, phytoplankton and zooplankton cannot survive, and the state is stable if the assumption holds; phytoplankton survive but zooplankton cannot survive, and the state is stable if the assumptions and hold; both phytoplankton and zooplankton survive, and the state is stable if the assumptions hold.

#### 3. Hopf Bifurcation and Analysis of System with Delay

##### 3.1. Hopf Bifurcation

In this part, we will only investigate the effect of time delay on the dynamics of system (2) about the coexistence equilibrium . For convenience, we denote the coexistence equilibrium = . We first translate the coexistence equilibrium to the origin by setting , and and drop the bar for convenience of notation. Then, the linearized system of system (2) can be obtainedwhere

Obviously, the characteristic equation of system (25) around takes the following formwhere

For , according to Corollary 2.4 [32], the sum of orders of the zeros of (27) in the open right half-plane can change only when a zero appears on or crosses the imaginary axis. It is obvious that is not a root of (27) when , and then the imaginary axis cannot be crossed by for any . And due to Corollary 2.4 [32], a stability switch (or a cross of the imaginary axis) occurs with for some real .

Now, we assume that is a root of (27). Separating real and imaginary parts, we haveAdding up the squares of both equations, we obtainBy denoting , (30) becomeswhere .

Let . According to the distribution of the roots of (31) and the results of [33], we can easily get the following lemmas.

Lemma 6. *(i) If , (31) has at least one positive root.**(ii) If and , (31) does not have any positive root.**(iii) If and , (31) has positive roots if and only if and .*

Thus, if the condition holds, then . Equation (31) has at least one positive root. Without loss of generality, we assume that it has three real positive roots, which are defined by , respectively. For convenience, we assume . Then (30) has three positive roots From (29), we have

Then is a pair of purely imaginary roots of (27) with . DefineAnd we have the following transversality condition.

Lemma 7. *Let be the root of (27) satisfying , , , . Then*

*Proof. *Differentiating both sides of (27) with respect to , we obtainThenFrom (30), we haveThus we haveSince , , we obtain , , . This ends the proof.

To summarize the above discussions, we have the following result.

Theorem 8. *Suppose that and hold. System (2) at is locally asymptotically stable whenever and unstable when . System (2) undergoes a Hopf bifurcation for .*

##### 3.2. Properties of Hopf Bifurcation

From Theorem 8, we have obtained the conditions for the existence of Hopf bifurcation when . In this subsection, we shall consider the direction of the Hopf bifurcation and the stability of bifurcating periodic solutions of system (2) by using the center manifold and normal form theory presented by Hassard [34].

Let , , and for . Using Taylor series expansion about , system (2) can be rewritten aswhere andwhereLinearization of the delayed system (41) around the origin is givenDefinewhere .

By the Riesz representation theorem [35], there exists matrix-valued function of bounded variation for such thatIn fact, we can choosewhere is a Dirac delta function.

Next, for , we defineandand then system (43) can be rewritten aswhere and for .

For , the adjoint operator of A is defined asand a bilinear inner productwhere . From the above discussion, we can know that are the eigenvalues of , and thus are also the eigenvalues of .

Let be the eigenvector of corresponding to the eigenvalue , and be the eigenvector of corresponding to the eigenvalue . With the conditions and , we obtainandwhere I is a third order identity matrix. Hence, we can obtain , .

In order to satisfy that , we need to choose the value of in the following part. For convenience, we denote . By virtue of (52) and (53), it derives thatHence, we can choose .

Nextly, we will compute the coordinates to detail center manifold at . Let be solution of (49) when . We define

On the center manifold , it derives thatwhere and are local coordinates for in the direction of and , respectively.

We only consider real solutions as is real. If is real, then we getFor the solution of (49), since , from (49) and (51), we havewhereIt follows from (55) and (56) that