# Dynamical Behaviour of a Nutrient-Plankton Model with Holling Type IV, Delay, and Harvesting.

1. IntroductionIn 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-phytoplanktonzooplankton 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:

[mathematical expression not reproducible] (1)

where N(t), x(t), and y(t) denote the concentration of nutrient, phytoplankton, and Zooplankton population at time t, respectively. [tau] 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: p(u) = cu/(m + [u.sup.2]), 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 delay

[mathematical expression not reproducible], (2)

where N(t), P(t), and Z(t) denote the number of nutrient, phytoplankton, and zooplankton population at time t, respectively. [N.sub.0] is the constant input of nutrient concentration and [alpha] is washout rate. b and [b.sub.1] ([b.sub.1] [less than or equal to] b) are the nutrient uptake rate and conversion rate of nutrient for growth of the phytoplankton, respectively. k and [mu] are the mortality rates of the phytoplankton and zooplankton population, respectively. [beta] and [beta], ([[beta].sub.1] [less than or equal to] [beta]) are the maximal ingestion rate and conversion rate of zooplankton. [theta] denotes the rate of toxin liberation. k, ([k.sub.1] [less than or equal to] k) is the nutrient recycle rate after death of phytoplankton population. The parameter r is the half saturation constant for a Holling type IV function. [tau] denotes the time delay in maturation of zooplankton. The constants [q.sub.1] and [q.sub.2] are the catchability coefficients of the phytoplankton and zooplankton population, and E is the harvesting effort.

The initial conditions of system (2) take the following form

[mathematical expression not reproducible], (3)

where [[phi].sub.1]([??]), [[phi].sub.2]([??]), [[phi].sub.3]([??]) [member of] C([-[tau], 0], [R.sup.3.sub.+]), the Banach space of continuous functions mapping the interval [-[tau], 0] into [R.sup.3.sub.+], where [R.sup.3.sub.+] = {([x.sub.1],[x.sub.2],[x.sub.3]): [x.sub.1] [greater than or equal to] 0, i = 1,2, 3}.

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 [tau] = 0, system (2) is

[mathematical expression not reproducible], (4)

The initial conditions for system (4) take the following form

N(0) = [[phi].sub.1](0) [greater than or equal to] 0, P(O) = [[phi].sub.2](0) [greater than or equal to] 0, Z(0) = [[phi].sub.3](0) [greater than or equal to] 0. (5)

Lemma 1. If r [greater than or equal to] 1, solutions of system (4) with initial conditions are defined on [0, +[infinity]) and remain positive for all t [greater than or equal to] 0.

Proof. From the first equation of system (4), we obtain

[??](i) [greater than or equal to] -[alpha]N - bNP (6)

Hence, we obtain

N(t) [greater than or equal to] [[phi].sub.1](0) exp [[integral].sup.t.sub.0](-[alpha]bP(s))ds > 0. (7)

From the second equation of system (4), we have

[mathematical expression not reproducible]. (8)

Since [[phi].sub.3](0) [greater than or equal to] 0, the third equation of system (4) shows that by a standard comparison argument

[mathematical expression not reproducible]. (9)

Theorem 2. All solutions of system (4) with initial conditions are bounded for all t [greater than or equal to] 0.

Proof. In order to check the boundedness of solutions of system (4), we define the function

W(t) = N(t)+[b/[b.sub.1]]P(t) + [[beta]b/([[beta].sub.1]-[theta])[b.sub.1]] Z(t). (10)

As derivative of W(t) with respect to system (4), we obtain

[mathematical expression not reproducible] (11)

where m = min{[alpha]a, K, [q.sub.1]E - [b.sub.1]k/b, [mu] + [q.sub.2]E}. Therefore, by the comparison theory [31], we have

W(t) [less than or equal to] (W(0)-[N.sub.0]/m)[e.sup.-mt] + [N.sub.0]/m. (12)

Thus, N(t), P(t), and Z(t) 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 [E.sub.0]([N.sub.0]/ [alpha], 0, 0), the zooplankton free equilibrium [E.sub.1]([N.sub.1],[P.sub.1],0), and the coexistence equilibrium [E.sup.*] ([N.sup.*], [P.sup.*], [Z.sup.*]).

The zooplankton free equilibrium is [E.sub.1]([N.sub.1], [P.sub.1], 0), where

[N.sub.1] = k+[q.sub.1]E/[b.sub.1] [P.sub.1] = [alpha](k+[q.sub.1]E)-[b.sub.1][N.sub.0]/[k.sub.1][b.sub.1]-bk-[bq.sub.1]E. (13)

Since [k.sub.1][b.sub.1]-bk-[bq.sub.1]E < 0, if [b.sub.1][N.sub.0] > [alpha](k+[q.sub.1]E), then there is only one possible equilibrium involving phytoplankton but not Zooplankton.

The coexistence equilibrium is [E.sup.*]([N.sup.*], [P.sup.*], [Z.sup.*]), where [N.sup.*], [P.sup.*], and [Z.sup.*] satisfy the following equations

[N.sub.0] - [alpha]N - bNP + [k.sub.1]P = 0, [b.sub.1]N-(k + [q.sub.1]E) - [beta]Z/[r+[P.sup.2]] = 0 ([[beta].sub.1] - [theta])P/r+[P.sup.2] - [mu]+[q.sub.2]E) = 0. (14)

From the third equation of (14), we have

([mu] + [q.sub.2]E)[P.sup.2] - ([[beta].sub.1] - [theta]) P + r ([mu] + [q.sub.2]E) = 0. (15)

If [[beta].sub.1] - [theta] > 0, and [DELTA] = [([[beta].sub.1] - [theta]).sup.2] - 4r[([mu] + [q.sub.2]E).sup.2] > 0, then (15) has exactly two positive real roots

[P.sup.*.sub.1] = [[beta].sub.1]-[theta]-[square root of [DELTA]]/2([mu]+[q.sub.2]E). [P.sup.*.sub.2] = [[beta].sub.1]-[theta]+[square root of [DELTA]]/1([mu]+[q.sub.2]E). (16)

And if [[beta].sub.1] - [theta] > 0, and [DELTA] = 0, then (15) has two equal positive real roots

[P.sup.*.sub.3] = [[beta].sub.1]-[theta]/2([mu]+[q.sub.2]E). (17)

Again from the first and second equations of (14), we have

[mathematical expression not reproducible] (18)

[mathematical expression not reproducible]

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

[mathematical expression not reproducible].

From the above analysis, we obtain the following results.

Theorem 3. (i) The plankton extinction equilibrium [E.sub.0]([N.sub.0]/[alpha], 0,0) always exists.

(ii) The zooplankton free equilibrium [E.sub.1] exists if [b.sub.1][N.sub.0] > [alpha](k + [q.sub.1]E).

(iii) The coexistence equilibrium [E.sup.*.sub.1]([N.sup.*.sub.1],[P.sup.*.sub.1],[Z.sup.*.sub.1]) (or [E.sup.*.sub.2]([N.sup.*.sub.2],[P.sup.*.sub.2],[Z.sup.*.sub.2])) exists if([H.sub.1]) holds; the coexistence equilibria [E.sup.*.sub.1] and [E.sup.*.sub.2] exist if ([H.sub.2]) holds; the coexistence equilibrium [E.sup.*.sub.3] ([N.sup.*.sub.3], [P.sup.*.sub.3], [Z.sup.*.sub.3]) exists if([H.sub.3]) holds.

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

Theorem 4. (1) If [b.sub.1][N.sub.0] < [alpha](k + [q.sub.1]E), then the plankton extinction equilibrium [E.sub.0] is locally asymptotically stable.

(2) If [b.sub.1][N.sub.0] > [alpha](k + [q.sub.1]E) and [mu] >([[beta].sub.1] - [theta])([alpha]k + [alpha][q.sub.1]E- [b.sub.1][N.sub.0]) [([k.sub.1][b.sub.1] -bk-[bq.sub.1]E).sup.2]/[(r([k.sub.1][b.sub.1]-bk-[bq.sub.1]E).sup.2] + [([alpha]k + [alpha][q.sub.1]E- [b.sub.1][N.sub.0]).sup.2]) - [q.sub.1]E, then the zooplankton free equilibrium [E.sub.1] is locally asymptotically stable.

(3) The coexistence equilibrium [E.sup.*.sub.1] (or [E.sup.*.sub.2]) is locally asymptotically stable if conditions ([H.sub.1]) and ([H'.sub.1]) hold. The coexistence equilibria [E.sup.*.sub.1] and [E.sup.*.sub.2] are locally asymptotically stable if conditions ([H.sub.2]) and ([H'.sub.2]) hold. The coexistence equilibrium [E.sup.*.sub.3] is locally asymptotically stable if conditions ([H.sub.3]) and ([H'.sub.3]) hold.

Proof. (1) The characteristic equation on [E.sub.0]([N.sub.0]/[alpha],0,0) is given

([lambda] + [alpha])([lambda] - [b.sub.1][N.sub.0]/[alpha] + k + [q.sub.1]E)([lambda]+[mu]+[Eq.sub.2]) = 0. (19)

It is clear that (19) has negative roots [[lambda].sub.1] = -[alpha], and [[lambda].sub.2] = -([mu] + [Eq.sub.2]). So if [b.sub.1][N.sub.0] < [alpha](k + [q.sub.1]E), then

[[lambda].sub.3] = [b.sub.1][N.sub.0]/[alpha] - k - [q.sub.1]E < 0. (20)

Thus, the plankton extinction equilibrium [E.sub.0] = ([N.sub.0]/[alpha], 0, 0) is locally asymptotically stable. (2) The characteristic equation about [E.sub.1]([N.sub.1], [P.sub.1],0) is given by

[[lambda] - ([[beta].sub.1]-[theta])[P.sub.1]/r+[([P.sub.1]).sup.2] + [mu] + [q.sub.2]E] x [[[lambda].sup.2] + ([alpha] + b[P.sub.1])[lambda] + [b.sub.1][P.sub.1] (b[N.sub.1] - [k.sub.1])] = 0. (21)

If [mu] > ([[beta].sub.1]-[theta])([alpha]k + [alpha][q.sub.1]E-[b.sub.1][N.sub.0])([k.sub.1][b.sub.1]-bk- [bq.sub.1]E)/[(r([k.sub.1][b.sub.1]-bk - [bq.sub.1]E).sup.2] + [([alpha]k + [alpha][q.sub.1]E - [b.sub.1][N.sub.0]).sup.2]) - [q.sub.1]E, then

[[lambda].sub.1] = ([[beta].sub.1,]-[theta])[P.sub.1]/r+[([P.sub.1]).sup.2] - [mu] [q.sub.2]E < 0. (22)

Further, the other roots [[lambda].sub.2] and [[lambda].sub.3] satisfy that [[lambda].sub.2] + [[lambda].sub.3] = - ([alpha] + [b.sub.1][P.sub.1]) < 0. If [[lambda].sub.2][[lambda].sub.3] = [b.sub.1][P.sub.1](-b[N.sub.1] + [k.sub.1]) > 0, then [[lambda].sub.2] <0, [[lambda].sub.3] < 0. Therefore, if [b.sub.1][N.sub.0] > [alpha](k + [q.sub.1]E) and [mu] > ([[beta].sub.1] - [theta])[([alpha]k + [alpha][q.sub.1]E-[b.sub.1][N.sub.0]).sup.2])([k.sub.1][b.sub.1]-bk-[bq.sub.1]E)/[(r([k.sub.1][b.sub.1]-bk- [bq.sub.1]E).sup.2] + [([alpha]k + [alpha][q.sub.1]E - [b.sub.1][N.sub.0]).sup.2]) - [q.sub.1]E, the zooplankton free equilibrium [E.sub.1] = ([N.sub.1],[P.sub.1],0) is locally asymptotically stable.

(3) The characteristic equation about [E.sup.*.sub.1]([N.sup.*.sub.1],[P.sup.*.sub.1],[Z.sup.*.sub.1]) is given

[[lambda].sup.3] + [M.sub.1][[lambda].sup.2] + [M.sub.2][lambda] + [M.sub.3] = 0, (23)

where

[mathematical expression not reproducible]. (24)

Since [M.sub.1][M.sub.2]-[M.sub.3] = [m.sup.2.sub.1][m.sub.2]+[alpha][m.sub.1][m.sub.3]+[P.sup.*.sub.1][m.sub.1][m.sub.6]+([m.sub.2]+[m.sub.3])[M.sub.2], if ([H.sub.1]) and ([H'.sub.1]) hold, then [M.sub.1] > 0, [M.sub.1][M.sub.2] - [M.sub.3] > 0. Therefore all roots of (23) have negative real parts. Based on Routh-Hurwitz criterion we obtain that the coexistence equilibrium [E.sup.*.sub.1]([N.sup.*.sub.1],[P.sup.*.sub.1],[Z.sup.*.sub.1]) (or [E.sup.*.sub.2]([N.sup.*.sub.2],[P.sup.*.sub.2],[Z.sup.*.sub.2])) is locally asymptotically stable.

Similarly, if ([H.sub.2]) and ([H'.sub.2]) hold, the coexistence equilibria [E.sup.*.sub.1]([N.sup.*.sub.1],[P.sup.*.sub.1],[Z.sup.*.sub.1]) and [E.sup.*.sub.2]([N.sup.*.sub.2],[P.sup.*.sub.2],[Z.sup.*.sub.2]) are locally asymptotically stable; if ([H.sub.3]) and ([H'.sub.3]) hold, the coexistence equilibrium [E.sup.*.sub.3] ([N.sup.*.sub.3], [P.sup.*.sub.3,] [Z.sup.*.sub.3]) is locally asymptotically stable.

Remark 5. From biological point of view, phytoplankton and zooplankton cannot survive, and the state is stable if the assumption [b.sub.1][N.sub.0] < [alpha](k + [q.sub.1]E) holds; phytoplankton survive but zooplankton cannot survive, and the state is stable if the assumptions [b.sub.1][N.sub.0] > [alpha](k + [q.sub.1]E) and [mu] > ([[beta].sub.1]-[theta])([alpha]k + [alpha][q.sub.1]2E-[b.sub.1][N.sub.0])([k.sub.1][b.sub.1] -bk-[bq.sub.1]E)/[(r([k.sub.1][b.sub.1]-bk-[bq.sub.1]E).sup.2] + [([alpha]k + [alpha][q.sub.1]E- [b.sub.1][N.sub.0]).sup.2]) - [q.sub.1]E hold; both phytoplankton and zooplankton survive, and the state is stable if the assumptions ([H.sub.1]) - ([H'.sub.3]) 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 [E.sup.*.sub.1]. For convenience, we denote the coexistence equilibrium [E.sup.*.sub.1] = [E.sup.*]. We first translate the coexistence equilibrium [E.sup.*] to the origin by setting [bar.N] = N - [N.sup.*], [bar.P] = P - [P.sup.*] and [bar.Z] = Z - [Z.sup.*] and drop the bar for convenience of notation. Then, the linearized system of system (2) can be obtained

[mathematical expression not reproducible], (25)

where

[mathematical expression not reproducible] (26)

Obviously, the characteristic equation of system (25) around [E.sup.*] takes the following form

[[lambda].sup.3] + [r.sub.1][[lambda].sup.2] + [r.sub.2][lambda] + [r.sub.5] + ([r.sub.3] + [r.sub.4][lambda]) [e.sup.- [lambda][tau]] = 0, (27)

where

[mathematical expression not reproducible]. (28)

For [tau] [not equal to] 0, 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 [lambda] = 0 is not a root of (27) when [[beta].sub.1] [not equal to] [theta], and then the imaginary axis cannot be crossed by [lambda] = 0 for any [tau] > 0. And due to Corollary 2.4 [32], a stability switch (or a cross of the imaginary axis) occurs with [lambda] = [+ or -]i[omega] for some real [omega] > 0.

Now, we assume that [lambda] = i[omega] ([omega] > 0) is a root of (27). Separating real and imaginary parts, we have

[r.sub.3] cos ([omega][tau]) + [r.sub.4][omega]sin ([omega][tau]) = [r.sub.1][[omega].sup.2]-[r.sub.5], [r.sub.4][omega] cos([omega][tau]) - [r.sub.3]sin ([omega][tau]) = [[omega].sup.3] - [r.sub.2][omega]. (29)

Adding up the squares of both equations, we obtain

[[omega].sup.6] + ([r.sup.2.sub.1]-2[r.sub.2])[[omega].sup.4]+([r.sup..sub.2]-2[r.sub.1][r.sub.5]- [r.sup.2.sub.4][[omega].sup.2]+[r.sup.2.sub.5]-[r.sup.2.sub.3] = 0. (30)

By denoting z = [[omega].sup.2], (30) becomes

[z.sup.3] + [p.sub.1][z.sup.2] + [p.sub.2]z + [p.sub.3] = 0, (31)

where [p.sub.1] = 2[r.sub.2], [p.sub.2] = [r.sup.2.sub.2] - 2[r.sub.1][r.sub.5] - [r.sup.2.sub.4],[p.sub.3] = [r.sup.2.sub.5]-[r.sup.2.sub.3].

Let h(z) = [z.sup.3] + [p.sup.1][z.sup.2] + [p.sub.2]z + [p.sub.3]. 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 [p.sub.3] < 0, (31) has at least one positive root.

(ii) If [p.sub.3] [greater than or equal to] 0 and [p.sup.2.sub.1] - 3[p.sub.2] [less than or equal to] 0 (31) does not have any positive root.

(iii) If [p.sub.3] [greater than or equal to] 0 and [q.sup.2.sub.1] - 3[q.sub.2] > 0 (31) has positive roots if and only if [z*.sub.1] = (l/3)(-[p.sub.1] + [square root of ([p.sup.2.sub.1]-3[p.sub.2])]) > 0 and h([z.sup.*.sub.1]) [less than or equal to] 0.

Thus, if the condition (H,) holds, then h(0) = [p.sub.3] < 0. 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 [z.sub.1],[z.sub.2],[z.sub.3], respectively. For convenience, we assume [z.sub.1] < [z.sub.2] < [z.sub.3]. Then (30) has three positive roots [[omega].sub.k] = [bar.[z.sub.k]] = 1,2,3. From (29), we have

[[tau].sup.j.sub.k] = 1/[[omega].sub.k](arcos [r.sub.4][[omega].sup.4.sub.k]+ ([r.sub.1][r.sub.3]-[r.sub.2][r.sub.4])[[omega].sup.2.sub.k]-[r.sub.3][r.sub.5] / [r.sup.2.sub.4][[omega].sup.2.sub.k]+[r.sup.2.sub.3] + 2j[pi]), k = 1,2,3; j = 0,1,2, ... (32)

Then [+ or -]i[[omega].sub.k] is a pair of purely imaginary roots of (27) with [[tau].sup.j.sub.k]. Define

[mathematical expression not reproducible] (33)

And we have the following transversality condition.

Lemma 7. Let [lambda]([tau]) = [alpha]([tau]) + i[omega]([tau]) be the root of (27) satisfying [lambda]([[tau].sup.j.sub.k]) = 0 [omega]([[tau].sup.j.sub.k]) = [[omega].sub.k], k = 1, 2, 3, j = 0, 1, 2, .... Then

[mathematical expression not reproducible]. (34)

Proof. Differentiating both sides of (27) with respect to [tau], we obtain

d[lambda]/d[tau] = [[lambda]e.sup.-[lambda][tau]]([r.sub.4][lambda] + [r.sub.3])/ 3[[lambda].sup.2] + 2[r.sub.1][lambda] + [r.sub.2] - [tau] ([r.sub.4][lambda] + [r.sub.3])[e.sup.-[lambda][tau]] + [r.sub.4][e.sup.-[lambda][tau]]. (35)

Then

[(d[lambda]/d[tau]).sup.-1] = (3[[lambda].sup.2]+2[r.sub.1][lambda]+[r.sub.2][e.sup.[lambda][tau]]+[r.sub.4]/[lambda]([r.sub.4]+[r.sub.3]) - [tau]/[lambda]. (36)

From (30), we have

[mathematical expression not reproducible] (37)

Thus we have

[mathematical expression not reproducible] (38)

Since h(0) < 0, h([infinity]) = [infinity], we obtain h'([z.sub.1]) > 0, h'([z.sub.2]) < 0, h'([z.sub.3]) > 0. This ends the proof.

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

Theorem 8. Suppose that ([H.sub.1]) and ([H'.sub.1]) hold. System (2) at [E.sup.*] is locally asymptotically stable whenever [tau] [member of] [0, [[tau].sub.0]] and unstable when [tau] > [[tau].sub.0]. System (2) undergoes a Hopf bifurcation for [tau] = [[tau].sub.0].

3.2. Properties of Hopf Bifurcation. From Theorem 8, we have obtained the conditions for the existence of Hopf bifurcation when [tau] = [[tau].sub.0]. 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 x(t) = N(-[tau]t) - [N.sup.*], y(t) = P([tau]t) - [P.sup.*], z(t) = Z([tau]t) - [Z.sup.*], [tau] = [[tau].sub.0] + [sigma], [phi](t) = [(x(t), y(t), z(t)).sup.T], and [[phi].sub.t]([??]) = [phi](t + [??]) for [??] [member of] [-1,0]. Using Taylor series expansion about [E.sup.*] = ([N.sup.*], [P.sup.*], [Z.sup.*]), system (2) can be rewritten as

[??](t) = [tau]([G.sub.1][phi](t)+[G.sub.2][phi](t-1) + f), (39)

where

[mathematical expression not reproducible], (40)

[mathematical expression not reproducible], (41)

where

[mathematical expression not reproducible], (42)

Linearization of the delayed system (41) around the origin is given

[??](t) = [tau]([G.sub.1][phi](t) + [G.sub.2][phi](t-1)). (43)

Define

[L.sub.[sigma]]([phi]) = ([[tau].sub.0] + [sigma]) ([G.sub.1][phi](t) + [G.sub.2][phi](t-1)), (44)

Where [phi] [member of] [C.sup.1]([-1,0],[R.sup.3]).

By the Riesz representation theorem [35], there exists 3 x 3 matrix-valued function [zeta]([??], [sigma]) of bounded variation for [??] [member of] [-1,0] such that

[L.sub.[sigma]([phi])] = [[integral].sup.0.sub.-1]d[zeta]([??],[sigma])[phi]([??]), for [phi] [member of] C ([-1,0],[R.sup.3]). (45)

In fact, we can choose

[zeta]([??], [sigma]) = ([[tau].sub.0] + [sigma])[[G.sub.1][delta]([??]) + [G.sub.2][delta]([??]+1)], (46)

where [delta]([??]) is a Dirac delta function.

Next, for [phi] [member of] [C.sup.1]([-1,0],[R.sup.3]), we define

[mathematical expression not reproducible] (47)

and

[mathematical expression not reproducible], (48)

and then system (43) can be rewritten as

[[??].sub.t] = A([sigma])[u.sub.t] + R([sigma])[u.sub.t], (49)

where u = [(x, y, z).sup.T] and [u.sub.t] = u(t + [??]) for [??] [member of] [-1,0].

For [psi] [member of] [C.sup.1]([-1,0],[R.sup.3]), the adjoint operator [A.sup.*] of A is defined as

[mathematical expression not reproducible], (50)

and a bilinear inner product

[mathematical expression not reproducible], (51)

where [zeta]([??]) = [zeta]([??],0). From the above discussion, we can know that [+ or -]i[[omega].sub.0] are the eigenvalues of A(0), and thus [+ or -]i[[omega].sub.0] are also the eigenvalues of [A.sup.*](0).

Let [mathematical expression not reproducible] be the eigenvector of A(0) corresponding to the eigenvalue i[[omega].sub.0][[tau].sub.0], and [mathematical expression not reproducible] be the eigenvector of [A.sup.*](0) corresponding to the eigenvalue -i[[omega].sub.0][[tau].sub.0]. With the conditions A(0)q(0) = i[[omega].sub.0]q(0) and [A.sup.*] (0) [q.sup.*] (0) = -i[[omega].sub.0] [q.sup.*] (0), we obtain

[mathematical expression not reproducible], (52)

and

[mathematical expression not reproducible]. (53)

where I is a third order identity matrix. Hence, we can obtain [q.sub.2] = (i[[omega].sub.0] - [a.sub.11])/[a.sub.12],[q.sub.3] = ((i[[omega].sub.0] - [a.sub.22])[q.sub.2]- [a.sub.21])/[a.sub.23], [q.sup.*.sub.2] = -(i[[omega].sub.0] + [a.sub.11])/[a.sub.21], [q.sup.*.sub.3] = -[a.sub.23][q.sub.3]/i[[omega].sub.0].

In order to satisfy that <[q.sup.*], q> = 1, we need to choose the value of [D.sub.1] in the following part. For convenience, we denote [GAMMA] = [(1,[q.sub.2],[q.sub.3]).sup.T], [[GAMMA].sup.*] = [(1,[q.sup.*.sub.2],[q.sup.*.sub.3]).sup.T]. By virtue of (52) and (53), it derives that

[mathematical expression not reproducible]. (54)

Hence, we can choose [mathematical expression not reproducible]

Nextly, we will compute the coordinates to detail center manifold [C.sub.0] at [sigma] [not equal to] 0. Let [[mu].sub.t] be solution of (49) when [sigma] = 0. We define

z(t) = <[q.sup.*], [u.sub.t]>, W(t, [??]) = [u.sub.t] - zq - [bar.z][bar.q] = [u.sub.t] - 2 Rez (t)q([??]). (55)

On the center manifold [C.sub.0], it derives that

[mathematical expression not reproducible], (56)

where z and [bar.z] are local coordinates for [C.sub.0] in the direction of [q.sup.*] and [bar.[q.sup.*]], respectively.

We only consider real solutions as W is real. If [[mu].sub.t] is real, then we get

<[q.sup.*], W> = <[q.sup.*], [u.sub.t]-zq-[bar.z][bar.q]> = <[q.sup.*], [u.sub.t]> - <[q.sup.*],q> - <[q.sup.*],[bar.q]>. (57)

For the solution [[mu].sub.t] of (49), since [sigma] = 0, from (49) and (51), we have

[mathematical expression not reproducible], (58)

where

[mathematical expression not reproducible] (59)

follows from (55) and (56) that

[mathematical expression not reproducible]. (60)

From (41) and (59), it derives that

[mathematical expression not reproducible] (61)

Here [mathematical expression not reproducible].

According to comparing its coefficients with (59), it derives that

[mathematical expression not reproducible]. (62)

Since [W.sub.11] and [W.sub.20] are still unknown, we still compute them next. By virtue of (49) and (55), we have

[mathematical expression not reproducible]. (63)

Here,

[mathematical expression not reproducible] (64)

By computing and comparing coefficients, we can obtain

[mathematical expression not reproducible]. (65)

From (63), we know that

[mathematical expression not reproducible]. (66)

Using (59) in (66) and comparing coefficient with (64), it derives that

[H.sub.20]([??]) = -[g.sub.20]q([??])-[[bar.g].sub.02][bar.q]([??]), (67)

[H.sub.11]([??]) = -[g.sub.11]q([??])-[[bar.g].sub.11][bar.q]([??]). (68)

From the definition of A(0), based on (65) and (67), we can obtain

[mathematical expression not reproducible]. (69)

Because of [mathematical expression not reproducible], we have

[mathematical expression not reproducible] (70)

Similarly, it derives that from (65) and (66)

[mathematical expression not reproducible] (71)

Here [E.sub.1] = [([E.sup.(1).sub.1], [E.sup.(2).sub.1], [E.sup.(3).sub.1]).sup.T], [E.sub.2] = [([E.sup.(1).sub.2], [E.sup.(2).sub.2]>, [E.sup.(3).sub.2]).sup.T] are three dimensional constant vectors and can be determined by setting [??] = 0 in H(z, [bar.z], [??]).

By above definitions and conditions, we can compute

[mathematical expression not reproducible] (72)

where

[mathematical expression not reproducible], (73)

where

[mathematical expression not reproducible]. (74)

Based on above analysis, the following values can be given

[mathematical expression not reproducible] (75)

By computing (75), we have the following theorem.

Theorem 9. Properties of bifurcating periodic solution in the center manifold at critical value r0 are determined as follows:

(i) If [[mu].sub.2] > 0([[mu].sub.2] < 0), then Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for [tau] > [[tau].sub.0] ([tau] < [[tau].sub.0]).

(ii) Bifurcating periodic solutions are stable (unstable) if [[beta].sub.2] <0 ([[beta].sub.2] > 0).

(iii) Period of the periodic solution increases (decreases) if [T.sub.2] >0 ([T.sub.2] < 0).

3.3. Global Hopf Bifurcation. In the above parts, we have obtained that system (2) undergoes Hopf bifurcation at [F.sup.*]([N.sup.*],[F.sup.*],[Z.sup.*]) when [tau] = [[tau].sup.j.sub.k]. We will study the global continuation of these local bifurcating periodic solutions by using the global Hopf bifurcation theorem for delay differential equations [30], and state that system (2) admits periodic solutions globally for delay [tau] in a finite interval [[[tau].sup.j.sub.k], +[infinity]], k = 1,2,3; j = 1,2,.... This includes the values that are not necessarily near the local Hopf bifurcation values.

For convenience, we give some notations, [R.sup.3.sub.+] = {(N,P,Z) : N > 0, > 0, Z > 0}, X = C([-[tau], 0], [R.sup.3.sub.+]), and x(t) = [(N(t), P(t), Z(t)).sup.T] We can rewrite system (2) as the following function

[mathematical expression not reproducible] (76)

Where [mathematical expression not reproducible].

According to the work of Wu [30], we need to introduce the following notions

[mathematical expression not reproducible]. (77)

Let C([E.sup.*],[[tau].sup.j.sub.k],2[pi]/[[omega].sub.k]) be the connected component of C([E.sup.*],[[tau].sub.j.sub.k],2[pi]/[[omega].sub.k]) in [SIGMA], and denote [Proj.sub.[tau]]([E.sup.*],[[tau].sup.j.sub.k],2[pi]/[[omega].sub.k]) as its projection on [tau] component. [+ or -]i[[omega].sub.k] is a pair purely imaginary roots of (27) with [tau] = [[tau].sup.i.sub.k]. By Theorem 8, we know that C([E.sup.*],[[tau].sup.j.sub.k],2[pi]/[[omega].sub.k]) is nonempty.

Lemma 10. Whenever the equilibrium [E.sup.*] exists, all the nontrivial solutions of system (2) are uniformly bounded.

Proof. We assume that (N(t), P(t), Z(t)) is an arbitrary positive periodic solution of system (2) subject to initial condition (3). Let [t.sub.1] and [t.sub.2] be real quantities defined by N([t.sub.1]) = min|N(i)|, N([t.sub.2]) = max{N(t)}, P([t.sub.1]) = min|P(i)|, P([t.sub.2]) = max{P(t)}, Z([t.sub.1]) = min|Z(i)}, and Z([t.sub.2]) = max|Z(i)}. From the first and second equations of system (2), we have

[mathematical expression not reproducible], (78)

which clearly implies that

[mathematical expression not reproducible], (79)

Since [k.sub.1][b.sub.1] < kb, [[theta].sub.1] > 0.

Similarly, from the second equation of system (2), we have

[mathematical expression not reproducible]. (80)

Thus, we obtain

[beta]P([t.sub.2])P([t.sub.2])/r+P[([t.sub.2]).sup.2] < [b.sub.1]N([t.sub.2])P([t.sub.2]) < [b.sub.1][[theta].sub.1]. (81)

We defined the function [W.sub.1](t) = N(t) + (b/[b.sub.1])P(t) + (b[beta]/[b.sub.1][[beta].sub.1])Z(t), for all t > 0. Then it follows from system (2) that

[mathematical expression not reproducible], (82)

where [??] = min{[alpha],k+[q.sub.1][q.sub.2]E}. Therefore, by comparison with theory [31], we have

[W.sub.1](t) [less than or equal to] ([W.sub.1](0)-[N.sub.0]b[[theta].sub.1]/[??])[e.sup.-[??]t] + [N.sub.0]+b[[theta].sub.1]/[??]. (83)

Thus, all the nontrivial solutions of system (2) are uniformly bounded for t [greater than or equal to] 0.

Lemma 11. System (2) has no nonconstant [tau]-periodic solution if the condition ([H.sub.1]) and

([H.sub.4]) : sup {[a.sub.11] + [a.sub.22] + [absolute value of [a.sub.23]], [absolute value of [a.sub.31] + [a.sub.32]] + [a.sub.11]+ [absolute value of [a.sub.12]], [absolute value of [a.sub.21]] + [a.sub.22]} < 0 (84)

are satisfied.

Proof. In order to explain that system (2) has no nonconstant periodic solution of periodic [tau], we only prove that system (4) has no nonconstant periodic solution (see [36], Lemma 5.5). We rewrite system (4) as [??](t) = f (u(t)). Suppose that u(t) is nonconstant periodic solution for system (4); then () = ([u.sub.1],[u.sub.2],[u.sub.3]) is nonconstant periodic solution of the following equations

[mathematical expression not reproducible]. (85)

From the Lemma 10, we know that the periodic orbits of system (85) corresponding to u(t) are included in the following area

G = {u [member of] [R.sup.3]: [bar.m] < [absolute value of [u.sub.j]] < [bar.M], j = 1,2,3}, (86)

where [bar.m] and [bar.M] are the boundedness of the periodic solutions of system (2) in the G.

The Jacobian matrix of system (85) is

[mathematical expression not reproducible]. (87)

By using the method of [37], its second additive compound matrix is

[mathematical expression not reproducible]. (88)

The second compound system is

[[??].sub.1] = ([a.sub.11] + [a.sub.22])[x.sub.1] + [a.sub.23][x.sub.2], [[??].sub.2] = ([a.sub.31] + [a.sub.32])[x.sub.1] + [a.sub.11][x.sub.2] + [a.sub.12][x.sub.3], [[??].sub.3] = [a.sub.21][x.sub.2] + [a.sub.22][x.sub.3]. (89)

We select vector module for [R.sup.3] as follows

[absolute value of (([x.sub.1],[x.sub.2],[x.sub.3]))] = max {[absolute value of [x.sub.1]], [absolute value of [x.sub.2]], [absolute value of [x.sub.3]]}. (90)

Then, the Lozinski [??] measure [??]([J.sup.[2]](u)) of second additive compound matrix corresponding to a vector module above is

[mu]([J.sup.[2]](u)) = max{[a.sub.11]+[a.sub.22]+[absolute value of [a.sub.23]], [absolute value of [a.sub.31]+[a.sub.32]] + [a.sub.11] + [absolute value of [a.sub.12]], [absolute value of [a.sub.21]], [a.sub.22]}. (91)

Thus, if ([H.sub.4]) holds, there exists an [epsilon] such that [mu](t) [less than or equal to] -[epsilon] < 0 for all u [member of] G.

This shows that the conclusion of Lemma 11 follows from Corollary 3.5 in [38].

Theorem 12. Assume ([H.sub.1]) and ([H.sub.4]) are fulfilled; then there exists [tau] [member of] I, such that system (2) has at least one positive periodic solution when [tau] varies among [[tau].sup.j.sub.1], [[tau].sup.j.sub.2], and [[tau].sup.j.sub.3], j [greater than or equal to] 1, where

I = {[[tau].sup.j.sub.k] > 0 : [[tau].sup.j.sub.k] are the Hopf bifurcation points of system (2)}. (92)

Proof. The characteristic matrix of system (2) at an equilibrium [bar.E] = ([bar.N], [bar.P], [bar.Z]) is given

[mathematical expression not reproducible], (93)

where

[mathematical expression not reproducible] (94)

([bar.E],[bar.[tau]],[bar.T]) is called a center if [??]([bar.E],[bar.[tau]],[bar.T]) = 0 and det([DELTA]([bar.E], [bar.[tau]],[bar.T])(t))(2[pi]i/T) = 0. A center is said to be isolated if it is the only center in some neighbourhood of ([bar.E],[bar.[tau]],[bar.T]). It is obvious that (76) has an equilibrium [x.sup.*] = [E.sup.*] in [R.sup.3.sub.+] under the condition ([H.sub.1]). It follows that ([E.sup.*], [[tau].sup.j.sub.k], 2[pi]/[[omega].sub.k]) is an isolated center from the discussion about the local Hopf bifurcation theorem in Section 3.1. By Lemma 7 and Theorem 8, there exist constants [epsilon] > 0 and [delta] > 0, as well as a smooth curve [lambda] : ([[tau].sup.j.sub.k] - [delta], [[tau].sup.j.sub.k] + [delta]) [right arrow] C such that

det [DELTA]([bar.E],[bar.[tau]],[bar.T])([lambda]([tau])) = 0, [absolute value of ([lambda] - i[[omega].sub.k])] < [epsilon] (95)

when [tau] [member of] [[[tau].sup.j.sub.k]-[delta], [[tau].sup.j.sub.k] + [delta]] holds, and

[mathematical expression not reproducible] (96)

Denotes [T.sub.k] = 2[pi]/[[omega].sub.k], and

[[OMEGA].sub.[epsilon],T] = {(u, T) : 0 < u < [epsilon], [absolute value of (T - [T.sub.k])] < [epsilon]}. (97)

Apparently, if [absolute value of ([tau] - [[tau].sup.j.sub.k])] [less than or equal to] [delta], (u, T) [member of] [partial derivative][[OMEGA].sub.[epsilon],T] such that det [mathematical expression not reproducible] if and only if u = 0, r = [[tau].sup.j.sub.k]. Hence all the assumptions ([A.sub.1])- ([A.sub.4]) of Theorem 3.2 in [30] are satisfied for m = 1. Moreover, if we define

[mathematical expression not reproducible], (98)

then we have the crossing number

[mathematical expression not reproducible] (99)

By Theorem 3.2 of Wu [30], we conclude that the connected component C([E.sup.*], [[tau].sup.j.sub.k],2[pi]/[[omega].sub.k]) through ([E.sup.*],[[tau].sup.j.sub.k],2[pi]/[[omega].sub.k]) in [SIGMA] is nonempty. Meanwhile, we have that

(i) C([E.sup.*],[[tau].sup.j.sub.k], 2[pi]/[[omega].sub.k]) is unbounded or

(ii) C([E.sup.*],[[tau].sup.j.sub.k],2[pi]/[[omega].sub.k]) is bounded, C([E.sup.*],[[tau].sup.j.sub.k],2[pi]/[[omega].sub.k]) [intersection] N(F) is finite, and [mathematical expression not reproducible].

By expression (32), we all know that there exists an integer N > 0, such that 2[pi] < [[tau].sup.j.sub.k][[omega].sub.k] < 2(n + 1)[pi], [[tau].sup.j.sub.k] [member of] J, which implies that [[tau].sup.j.sub.k]/(n + 1) < 2[pi]/[[omega].sub.k] < [[tau].sup.j.sub.k].

Therefore, we have that [tau]/(n + 1) < T < [tau] if ([bar.E], [tau], T) [member of] C([E.sup.*],[[tau].sup.j.sub.k], 2[pi]/[[omega].sub.k]). This fact and Lemma 11 show that the projection of C([E.sup.*],[[tau].sup.j.sub.k],2[pi]/[[omega].sub.k]) onto the T-space is bounded. On the flip side, Lemma 10 implies that the projection of C([E.sup.*],[[tau].sup.j.sub.k],2[pi]/[[omega].sub.k]) onto the x-space is uniformly bounded for [tau] [member of] [Proj.sub.[tau]]([E.sup.*],[[tau].sup.j.sub.k],2[pi]/[[omega].sub.k]). Otherwise, C([E.sup.*],[[tau].sup.j.sub.k],2[pi]/[[omega].sub.k]) is bounded and [Proj.sub.[tau]]([E.sup.*], [[tau].sup.j.sub.k], 2[pi]/[[omega].sub.k]) [subset] [[[tau].sub.k],[infinity]], which contradicts the fact that nonconstant periodic solution and T are both uniformly bounded for [tau] [member of] [Proj.sub.[tau]]([E.sup.*], [[tau].sup.j.sub.k], 2[pi]/[[omega].sub.k]). The proof is completed.

4. Optimal Harvesting

The emphasis of this section is on the profit making aspect of plankton; optimal control addresses the problem of finding control law for a given system such that a certain optimality criterion is satisfied [25]. In the economic perspective, the commercial development of recyclable resources is a critical issue and it is necessary to determine the optimal tradeoff between present and future harvests.

The harvesting problem covers cost of harvesting comprising a set of differential equations that describe the paths of the control variables that minimize the cost functional. It is assumed that price is a function that decreases with increasing biomass. Thus, to maximize the total discounted net revenues from the plankton, the optimal control problem can be formulated as

[mathematical expression not reproducible], (100)

where c is the constant harvesting cost per unit effort, [p.sub.1] and [p.sub.2] are the constant price per unit biomass of the phytoplankton and zooplankton, and [v.sub.1], [v.sub.2] are the economic constants. The convexity of the objective function with respect to E and the compactness of range values of the state variables can be combined to give the existence of the optimal control.

This problem is subject to system (2) and the control constrains 0 < E(t) < [E.sub.max]; here, [E.sub.max] is the maximum effect capacity of the harvesting industry. We want to find the optimal harvesting effect [E.sub.*] such that

J([E.sub.*]) = max {J(E) | E [member of] U}, (101)

where U is the control set defined by

U = {E | E is measurable and 0 [less than or equal to] E [less than or equal to] [E.sub.max] for all t}. (102)

To find the maximum value of J(E), we set X(t) = (N(t),P(t),Z(t)) and [X.sub.[tau]](t) = ([N.sub.[tau]](t),[P.sub.[tau]](t),[Z.sub.[tau]](t)), where [N.sub.[tau]](t) = N(t - [tau]), [P.sub.[tau]] = P(t - [tau]), [Z.sub.[tau]](t) = Z(t - [tau]). Now, the Hamiltonian function of this optimal control problem is given

[mathematical expression not reproducible] (103)

where [[lambda].sub.i] = [[lambda].sub.i](t), i = 1,2,3, are adjoint variables. Transversality conditions are [[lambda].sub.i]([t.sub.f]) = 0, i = 1,2,3.

Pontryagin's Maximum Principle with delay given in [28] is used to derive the necessary conditions for the optimal control. If we consider X(t) and [X.sub.[tau]](t) defined above, then there exists a continuous function [lambda](t) on [0,[t.sub.f]] satisfying the following three functions. The state equation is

X'(t) = [H.sub.[lambda]](X(t), [X.sub.[tau]](t), E(t), [lambda](t)), (104)

the optimality condition is

0 = [H.sub.E](X(t),[X.sub.[tau]](t), E(t), [lambda](t)), (105)

and the adjoint equation is

[mathematical expression not reproducible], (106)

where [H.sub.[lambda]], [H.sub.E], [H.sub.x], and [mathematical expression not reproducible] denote the derivative with respect to [lambda], E, X, [X.sub.t], respectively. Now we use the necessary conditions to Hamiltonian H.

Then there exist adjoint equations

[mathematical expression not reproducible]. (107)

The condition for the optimal harvesting can be obtained from (105), which is

[mathematical expression not reproducible]. (108)

Based on the above analysis, we have the following theorem.

Theorem 13. There exist an optimal control E(t) and corresponding solutions to system (2) that maximin J(E) over U. Furthermore, there exist adjoint variables [[lambda].sub.1](t), [[lambda].sub.2](t), and [[lambda].sub.3](t) which satisfy (107) with the transversality conditions [[lambda].sub.i]([t.sub.f]) = 0, i = 1,2,3. Moreover the optimal control is given by expression (108).

5. Numerical Simulation

In this section, numerical simulations are carried out to verify the theoretical findings by the use of MATLAB. Most values of parameters in system (2) are taken from [21]. [N.sub.0] = 1, [alpha] = 0.485, b = 0.388, [k.sub.1] = 0.0555, [b.sub.1] = 0.35, k = 0.07, [beta] = 0.76, [q.sub.1] = 0.3, E = 0.5, [[beta].sub.1] = 0.66, [mu] = 0.03, [theta] = 0.085, [q.sub.2] = 0.1, r=5.

When [tau] =0, we observe that ([H.sub.1]) holds and the coexisting equilibrium is [E.sup.*] = (1.37,0.78,1.90). Condition ([H'.sub.1]) is also valid. Therefore, the equilibrium [E.sup.*] is locally asymptotically stable (see Figure 1). Then, when [tau] [not equal to] 0, we obtained h(0) = [p.sub.3] = -0.0002 < 0. Hence, (30) has at least one positive root. We obtain [[omega].sub.0] = 0.13, and the critical value [[tau].sub.0] = 8.43. According to Theorem 8, system (2) remains stable when 0 < [tau] < [[tau].sub.0], but system (2) is unstable when [tau] > [[tau].sub.0]. The stability behaviour of system (2) around [E.sup.*] for [tau] = [[tau].sub.0] 8.3 is shown in Figure 2. Hopf bifurcation occurs when [tau] = 9 > [[tau].sub.0], shown in Figure 3.

Again, through the above parameter values, the properties of bifurcating periodic solutions at [[tau].sub.0] are given

[C.sub.1](0) = -307.56 - 114.2, [[mu].sub.2] = 53958, [[beta].sub.2] = -615.12, [T.sub.2] = 19131. (109)

Therefore, Hopf bifurcation is supercritical, bifurcating periodic solutions are stable (see Figure 4), and period of the periodic solution increases according to Theorem 9.

Nextly, we consider a set of parametric values to illustrate optimal control problem as follows: [N.sub.0] = 1, [alpha] = 0.485, b = 0.388, [k.sub.1] = 0.0555, [b.sub.1] = 0.35, k = 0.07, [beta] = 0.76, [q.sub.1] = 0.3, [[beta].sub.1] = 0.66, [mu] = 0.03, [theta] = 0.085, [q.sub.2] = 0.1, r = 5, c = 0.1, a = 0.2, [delta] = 0.2, [v.sub.1] = 0.4, [v.sub.2] = 0.7, [p.sub.1] = 1, [p.sub.2] = 5. We first use forward forth order Runge-Kutta method for the state variables of (4) with nonnegative initial conditions. Then we solve the adjoint variables from (107) by using backward forth order Runge-Kutta method. It is observed that as the constant input of nutrient concentration or the half saturation constant increases, the optimal harvesting effort increases (see (a,b) of Figure 5). However, as the rate of toxin liberation increases, the optimal harvesting effort decreases (see (c) of Figure 5). We observe that the optimal harvesting effort tends to be stable. We see that the nutrient increases firstly and then tends to be stable as time varies in the presence of harvesting (see (a) of Figure 6), and the phytoplankton and Zooplankton population decrease and even tend to be stable as time changes in the presence of harvesting (see (b,c) of Figure 6). All the final values of adjoint variables are zero (see Figure 7). Finally, using the same method in subject (100), we research optimal harvesting policy for model (2). Optimal harvesting effort is mainly dependent on gestation delay r of zooplankton displayed in Figure 8. We can observe that the optimal harvesting effort increases monotonically.

6. Conclusion and Discussion

Sharma et al. [21] discussed the effect of the time delay about the liberation of toxin on nutrient-plankton interaction. In this paper, the dynamic of nutrient-plankton interaction was considered with the gestation delay and linear harvesting. It is assumed that the mortality caused by toxin to zooplankton population is described as Holling type IV function. We obtained that the system lost local stability around E* and Hopf bifurcation occurred when time delay crossed its corresponding critical value. By using normal form theory and center manifold theorem, we obtained the properties of Hopf bifurcation. Next, as time delay varies in some regions, the global existence result of the periodic solution of Hopf bifurcation was given by a global Hopf bifurcation result due to Wu [30]. Pontryagin's Maximum Principle was used to obtain the optimal harvesting policy. We analyzed the effect of different parameter variables on the harvesting effort (see of Figure 5) and obtained the impact of time delay r on the harvesting effort (see of Figure 8). In the sense of biology, optimal harvesting effort increased with increasing of the constant input of nutrient concentration; optimal harvesting effort stayed the same as the half saturation constant changes; optimal harvesting effort decreased with the increasing toxin liberation of plankton. And optimal harvesting effort increased with the longer time delay in maturation of zooplankton.

In the future research, in order to make the model more practical, harvesting might be delayed because zooplankton avoid harvesting, so this factor can be considered into our nutrient-plankton system, which is

[mathematical expression not reproducible], (110)

where [[tau].sub.1] and [[tau].sub.2] are the time delays in maturation of zooplankton and harvesting because zooplankton avoid harvesting, respectively.

https://doi.org/10.1155/2018/9232590

Data Availability

The data used to support the findings of this study are included within the article. The reader can find the corresponding values of some parameters in the beginning of "Numerical Simulation". Most parameter values and data utilized in the numerical simulation section of this paper are taken from the numerical simulation section in [21]. Of course, the reader can freely access the parameter values and data supporting the conclusions of the study in the corresponding article [21].

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Grant No. 11661050, 11461041), the Development Program for HongLiu Outstanding Young Teachers in Lanzhou University of Technology (Q201308), and the Development Program for HongLiu Distinguished Young Scholars in Lanzhou University of Technology.

References

[1] J. Duinker and G. Wefer, "Das C[O.sub.2]-problem und die Rolle des Ozeans," Naturwissenschaften, vol. 81, no. 6, pp. 237-242, 1994.

[2] J. Chattopadhyay, "A delay differential equation model on harmful algal blooms in the presence of toxic substances," Mathematical Medicine and Biology, vol. 19, no. 2, pp. 137-161, 2002.

[3] T. Saha and M. Bandyopadhyay, "Dynamical analysis of toxin producingphytoplankton-zooplankton interactions," Nonlinear Analysis: Real World Applications, vol. 10, no. 1, pp. 314-332, 2009.

[4] Y. Wang, W. Jiang, and H. Wang, "Stability and global Hopf bifurcation in toxic phytoplankton-zooplankton model with delay and selective harvesting," Nonlinear Dynamics, vol. 73, no. 1-2, pp. 881-896, 2013.

[5] J. Zhao and J. Wei, "Stability and bifurcation in a two harmful phytoplankton-zooplankton system," Chaos, Solitons & Fractals, vol. 39, no. 3, pp. 1395-1409, 2009.

[6] P. Wang, M. Zhao, H. Yu, C. Dai, N. Wang, and B. Wang, "Nonlinear dynamics of a marine phytoplankton-zooplankton system," Advances in Difference Equations, vol. 2016, no. 1, 2016.

[7] T. Liao, H. Yu, and M. Zhao, "Dynamics of a delayed phytoplankton-zooplankton system with Crowley-Martin functional response," Advances in Difference Equations, vol. 1, pp. 5-35, 2017.

[8] R. Yang, M. Liu, and C. Zhang, "A diffusive toxin producing phytoplankton model with maturation delay and threedimensional patch," Computers & Mathematics with Applications. An International Journal, vol. 73, no. 5, pp. 824-837, 2017.

[9] Z. Zhang and M. Rehim, "Global qualitative analysis of a phytoplankton-zooplankton model in the presence of toxicity," International Journal of Dynamics and Control, vol. 5, no. 3, pp. 799-810, 2017.

[10] M. Banerjee and E. Venturino, "A phytoplankton-toxic phytoplankton-zooplankton model," Ecological Complexity, vol. 8, no. 3, pp. 239-248, 2011.

[11] A. K. Sharma, A. Sharma, and K. Agnihotri, "Bifurcation behaviors analysis of a plankton model with multiple delays," International Journal of Biomathematics, vol. 9, no. 6, pp. 1-25,2016.

[12] H. Xiang, Y. Wang, and H. Huo, "Analysis of the binge drinking models with demographics and nonlinear infectivity on networks," Journal of Applied Analysis and Computation, vol. 8, no. 5, pp. 1535-1554, 2018.

[13] X.-Y. Meng and Y.-Q. Wu, "Bifurcation and control in a singular phytoplankton-zooplankton-fish model with nonlinear fish harvesting and taxation," International Journal of Bifurcation and Chaos, vol. 28, no. 3,1850042, 24 pages, 2018.

[14] S. Chakraborty, S. Roy, and J. Chattopadhyay, "Nutrient-limited toxin production and the dynamics of two phytoplankton in culture media: a mathematical model," Ecological Modelling, vol. 213, no. 2, pp. 191-201, 2008.

[15] T. Zhang and W. Wang, "Hopf bifurcation and bistability of a nutrient-phytoplankton-zooplankton model," Applied Mathematical Modelling, vol. 36, no. 12, pp. 6225-6235, 2012.

[16] A. Fan, P. Han, and K. Wang, "Global dynamics of a nutrient-plankton system in the water ecosystem," Applied Mathematics and Computation, vol. 219, no. 15, pp. 8269-8276, 2013.

[17] S. Chakraborty, P. K. Tiwari, A. K. Misra, and J. Chattopadhyay, "Spatial dynamics of a nutrient-phytoplankton system with toxic effect on phytoplankton," Mathematical Biosciences, vol. 264, pp. 94-100, 2015.

[18] Y. Wang, H. Wang, and W. Jiang, "Stability switches and global Hopf bifurcation in a nutrient-plankton model," Nonlinear Dynamics, vol. 78, no. 2, pp. 981-994, 2014.

[19] S. R. Jang and E. J. Allen, "Deterministic and stochastic nutrient-phytoplankton-zooplankton models with periodic toxin producing phytoplankton," Applied Mathematics and Computation, vol. 271, pp. 52-67, 2015.

[20] K. Das and S. Ray, "Effect of delay on nutrient cycling in phytoplankton-zooplankton interactions in estuarine system," Ecological Modelling, vol. 215, no. 1-3, pp. 69-76, 2008.

[21] A. Sharma, A. K. Sharma, and K. Agnihotri, "The dynamic of plankton-nutrient interaction with delay," Applied Mathematics and Computation, vol. 231, pp. 503-515, 2014.

[22] Y. Li and D. Xiao, "Bifurcations of a predator-prey system of Holling and Leslie types," Chaos, Solitons & Fractals, vol. 34, no. 2, pp. 606-620, 2007.

[23] K. Chakraborty and K. Das, "Modeling and analysis of a two-zooplankton one-phytoplankton system in the presence of toxicity," Applied Mathematical Modelling: Simulation and Computation for Engineering and Environmental Systems, vol. 39, no. 3-4, pp. 1241-1265, 2015.

[24] Y. F. Lv, Y. Z. Pei, S. J. Gao, and C. G. Li, "Harvesting of a phytoplankton-zooplankton model," Nonlinear Analysis: Real World Applications, vol. 11, no. 5, pp. 3608-3619, 2010.

[25] M. R. Garvie and C. Trenchea, "Optimal control of a nutrient-phytoplankton-zooplankton-fish system," SIAM Journal on Control and Optimization, vol. 46, no. 3, pp. 775-791,2007.

[26] R. P. Gupta and P. Chandra, "Bifurcation analysis of modified Leslie-Gower predator-prey model with Michaelis-Menten type prey harvesting," Journal of Mathematical Analysis and Applications, vol. 398, no. 1, pp. 278-295, 2013.

[27] S. V. Krishna, P. D. N. Srinivasu, and B. Kaymakcalan, "Conservation of an ecosystem through optimal taxation," Bulletin of Mathematical Biology, vol. 60, no. 3, pp. 569-584,1998.

[28] S. Biswas, S. K. Sasmal, S. Samanta, M. Saifuddin, N. Pal, and J. Chattopadhyay, "Optimal harvesting and complex dynamics in a delayed eco-epidemiological model with weak Allee effects," Nonlinear Dynamics, vol. 87, no. 3, pp. 1553-1573, 2017.

[29] T. Das, R. N. Mukherjee, and K. S. Chaudhuri, "Harvesting of a prey-predator fishery in the presence of toxicity," Applied Mathematical Modelling, vol. 33, no. 5, pp. 2282-2292, 2009.

[30] J. H. Wu, "Symmetric functional-differential equations and neural networks with memory," Transactions of the American Mathematical Society, vol. 350, no. 12, pp. 4799-4838, 1998.

[31] V. Lakshmikantham and S. Leela, Differential and Integal Inequalities (Theory and Application): Ordinary Differential Equations I, Academic Press, New York, NY, USA, 1969.

[32] S. G. Ruan and J. J. Wei, "On the zeros of transcendental functions with applications to stability of delay differential equations with two delays," Dynamics of Continuous, Discrete & Impulsive Systems A: Mathematical Analysis, vol. 10, no. 6, pp. 863-874, 2003.

[33] Y. Song, M. Han, and J. Wei, "Stability and Hopf bifurcation analysis on a simplified BAM neural network with delays," Physica D: Nonlinear Phenomena, vol. 200, no. 3-4, pp. 185-204, 2005.

[34] B. D. Hassard, N. D. Kazarinoff, and Y.-H. Wan, Theory and Applications of Hopf Bifurcation, Cambridge University Press, Cambridge, UK, 1981.

[35] Y. Kuang, Delay Differential Equations with Applications in Population Dynamics, Academic Press, New York, NY, USA, 1993.

[36] Y. Qu, J. Wei, and S. Ruan, "Stability and bifurcation analysis in hematopoietic stem cell dynamics with multiple delays," Physica D: Nonlinear Phenomena, vol. 239, no. 20-22, pp. 2011-2024, 2010.

[37] J. S. Muldowney, "Compound matrices and ordinary differential equations," Rocky Mountain Journal of Mathematics, vol. 20, no. 4, pp. 857-872, 1990.

[38] Y. Li and J. S. Muldowney, "On Bendixson's criterion," Journal of Differential Equations, vol. 106, no. 1, pp. 27-39, 1993.

Xin-You Meng (iD), Jiao-Guo Wang, and Hai-Feng Huo (iD)

School of Science, Lanzhou University of Technology, Lanzhou, Gansu 730050, China

Correspondence should be addressed to Xin-You Meng; xymeng@lut.cn

Received 15 August 2018; Accepted 19 November 2018; Published 9 December 2018

Academic Editor: Ewa Pawluszewicz

Caption: Figure 1: Waveform plot and phase of system (4) at the equilibrium (1.37,0.78,1.90). (a) Waveform plot, (b) phase diagram.

Caption: Figure 2: The equilibrium [E.sup.*] of system (2) is locally asymptotically stable with initial value (1.2,1,1.7) when [tau] =8.3 < [[tau].sub.0]. (a) Waveform plot, (b) phase diagram.

Caption: Figure 3: The form of limit cycle of system (2) around [E.sup.*] with initial value (1.2,1,1.6) when [tau] = 9 > [[tau].sub.0]. (a) Waveform plot, (b) phase diagram.

Caption: Figure 4: A stable periodic orbits at [tau] = 13. (a) Regular oscillations of system (2), (b) phase diagram.

Caption: Figure 5: Optimal plankton harvesting with respect to different parameters: (a) the constant input of nutrient concentration [N.sub.0], (b) the half saturation constant r, (c) the rate of toxin liberation [theta].

Caption: Figure 6: Diagrams for the state variables with control for system (4): (a) the nutrient, (b) the phytoplankton population, (c) the zooplankton population.

Caption: Figure 7: Diagrams for the adjoint variables of system (4): (a) [[lambda].sub.1], (b) [[lambda].sub.2], (c) [[lambda].sub.3].

Caption: Figure 8: Variation of the optimal effort E with respect to time with different [tau] in system (2).

Printer friendly Cite/link Email Feedback | |

Title Annotation: | Research Article |
---|---|

Author: | Meng, Xin-You; Wang, Jiao-Guo; Huo, Hai-Feng |

Publication: | Discrete Dynamics in Nature and Society |

Date: | Jan 1, 2018 |

Words: | 9996 |

Previous Article: | On Convergence of Infinite Matrix Products with Alternating Factors from Two Sets of Matrices. |

Next Article: | Complexity Analysis of Dynamic Cooperative Game Models for Supply Chain with the Remanufactured Products. |