Bifurcation analysis in a diffusive phytoplankton–zooplankton model with harvesting

A diffusive phytoplankton–zooplankton model with nonlinear harvesting is considered in this paper. Firstly, using the harvesting as the parameter, we get the existence and stability of the positive steady state, and also investigate the existence of spatially homogeneous and inhomogeneous periodic solutions. Then, by applying the normal form theory and center manifold theorem, we give the stability and direction of Hopf bifurcation from the positive steady state. In addition, we also prove the existence of the Bogdanov–Takens bifurcation. These results reveal that the harvesting and diffusion really affect the spatiotemporal complexity of the system. Finally, numerical simulations are also given to support our theoretical analysis.


Introduction
In marine ecosystems, the plankton consist of two species, phytoplankton and zooplankton, which are the basis of the aquatic food chain [1]. Plankton is very important to the marine ecosystems, and the accumulation of plankton can cause "red tide" [2]. It is one of the most serious environmental problems faced by the whole world. A red tide is badly endangering the health of marine and human life.
In order to better understand plankton and explain the cause of "red tide", the dynamics of plankton systems have been studied by a number of authors [3][4][5][6][7]. The results of studies suggested that plankton systems could exhibit complex dynamic behaviors, such as Hopf bifurcation, global Hopf-bifurcation, Hopf-transcritical bifurcation, and so on. Other references [8][9][10][11] discussed the persistence, positivity, boundedness, and chaos of phytoplankton-zooplankton models. Most of these models were governed by ordinary or delay differential equations.
The level of plankton species changes not only in time but also in space. Hence, interaction and spatial processes of phytoplankton and zooplankton should be taken into account in mathematical models of plankton population dynamic systems [12][13][14]. Consequently, the construction of phytoplankton and zooplankton (prey and predator) models was commonly used by reaction-diffusion systems. According to a widely accepted approach [15][16][17][18], the functioning of phytoplankton and zooplankton (prey and predator) can be described by the following reaction-diffusion system: where (1) P and Z describe the densities of phytoplankton and zooplankton (prey and predator) at moment T and position X, respectively; (2) D 1 and D 2 are diffusivities; (3) f (P, Z) is the local growth and natural mortality of the prey, g(P, Z) describes the functional response for the grazing of phytoplankton by zooplankton; (4) k is the ratio of biomass conversion. The parameter d is the mortality rate of the predator. The choice of the functional response f (P, Z) and g(P, Z) in (1.1) can pick various combinations, which depend on the type of the prey and predator population. Based on the results of field and laboratory observations for plankton systems [19], we assume that the growth of the prey is logistic and g(P, Z) takes the Holling type II functional response [20,21]. According to the above facts, the model system (1.1) can be expressed in the following form: where a 1 describes the maximum predation rate, n is the self-saturation prey density.
Furthermore, some phytoplankton and zooplankton can be harvested for food. Hence, the stocks of edible plankton play an important role for fishery management. For system (1.2), considering constant harvesting, Chang et al. [25] discussed the existence and stability of Hopf bifurcation from the positive constant steady state and derived the direction and stability of bifurcating periodic solutions, and also considered an optimal control problem.
In several types of harvesting, Michaelis-Menten type harvesting is more realistic from biological and economic points of view [26]. In this paper, we investigate the following model with Michaelis-Menten type prey harvesting: where E measures the harvest effort, q represents the catchability, and m 1 and m 2 are appropriate constants. For simplicity, dimensionless variables are introduced. Let then we can obtain the following reaction-diffusion phytoplankton-zooplankton model: In this paper, we consider the following homogeneous reaction-diffusion phytoplankton-zooplankton system with homogeneous Neumann boundary condition: where is bounded in R with smooth boundary ∂ , denotes the Laplacian operator. The organization of the rest of the paper is as follows: In Sect. 2, we obtain global existence and boundedness of system (1.5) and demonstrate the asymptotic behavior of system (1.5). In Sect. 3, we give bifurcation analysis of system (1.5) under different parameters, including Hopf bifurcation, Bogdanov-Takens bifurcation. In Sect. 4, some numerical simulations are included to test and verify our theoretical analysis. Finally, we end the paper with a brief conclusion in Sect. 5.

Dynamical behavior 2.1 Global existence and boundedness
In this subsection, we prove the global existence of solutions of (1.5) and establish a priori bound of the solution.
By using the Neumann boundary condition and adding (2.2) and (2.3), we can obtain Since lim sup t→+∞ u(x, t) ≤ 1, we can obtain lim sup t→+∞ U(x, t) ≤ | |. Thus, for small ε > 0, there exists T 1 > 0 such that which completes the proof of (b).

Existence and stability of the positive constant steady state solution
The steady state solutions of (1.5) satisfy: (2.6) Then (1.5) has the unique positive constant solution (σ , v σ ) by mathematical calculation, where The positive constant solution exists if and only if Further in this section, we fix α, β, β 1 , c, m, σ and take h as the bifurcation parameter. We have discussed the existence and stability of the positive constant steady state solutions, and the bifurcating periodic solutions were affected by the variance of h. The local stability of the positive constant steady state solution can be summarized as follows.
Based on the biological significance, we abandon this case.

Bifurcation analysis using the harvesting h as the parameter
Next, we study Hopf bifurcations from the constant steady state of (1.5) with the Neumann boundary condition on the spatial domain = (0, π).
In this subsection, we show the existence of the spatial homogeneous and nonhomogeneous periodic solutions of system (1.5), and we also obtain the conditions for direction and stability of the Hopf bifurcation. As is well known, the eigenvalue problem -ψ = μψ, x ∈ (0, π), ψ (0) = ψ (π) = 0 has eigenvalues μ n = k 2 (k = 0, 1, 2, . . .) corresponding to the eigenfunction ψ n (x) = cos kx. From the proof of Theorem 2.2, the eigenvalues λ(k) of J k (k ≥ 0) are given by In the following, we will identify the Hopf bifurcation points h H that satisfy the necessary and sufficient condition (H1) in [20]: (H1) There exists k ∈ N 0 = N ∪ {0} such that And then we discuss the spatially nonhomogeneous Hopf bifurcation for k ≥ 1 (H1). In the following, we assume that the condition of Theorem 2. 2 (a(1), b(1), b(3)) is satisfied. The case of the conditions of Theorem 2. 2 (a(2), b(2), b(4)) is similar, thus they are omitted here. From T k = 0, we have then T j (h H j ) = 0 and T i (h H j ) = 0 for any i = j. These h H j satisfy the following inequality: Apparently, if holds, then for all x ≥ 1, g(x) = -d 2 M * x + d 2 d 2 x 2 is positive. Finally, we verify that the transversality condition holds: By applying the Hopf bifurcation theorem [20] and our analysis, we get the following results. Proof Above we have already discussed the existence of Hopf bifurcation. In the following, we will calculate the direction of Hopf bifurcation and the stability of spatially homogeneous periodic solutions by using the method in [20,29]. Take where ω 0 = √ -B * C * (-B * C * > 0). Define the inner product in X C by Let f (u, v) = u(1u) -βuv α+u -hu c+u and g(u, v) = β 1 uv α+umv, at a positive constant solution (σ , v σ ), by calculation, we obtain Then we get the following formulas by simple calculation: 20 , Therefore, According to the results of [20,29], we know that the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions are stable (unstable) when Re(c 1 (λ H 0 )) < 0 (Re(c 1 (λ H 0 )) > 0). In Sect. 4, we test and verify the direction of Hopf bifurcation and the stability of the spatially homogeneous periodic solutions by numerical simulations.

Bifurcation analysis with k as the parameter
In this section, we investigate bifurcation analysis arising at the positive constant steadystate solution by using k as the parameter. Based on (2.8) and the conclusions in [20,[30][31][32], we have the following results. Proof (1) Based on (2.8), we know that there is k 0 = 0 such that T k 0 = A * , D k 0 = -B * C * > 0. According to the analysis of Sect. 3.1, we know that when h = h H 0 , then (2) (a) There exists k 0 > 0, when k 2 0 = A * d 1 +d 2 , k 4 0 d 2 2 + B * C * < 0, then T k 0 = 0, D k 0 > 0. For j = k 0 , when j > k 0 , Hence, Eq. (2.8) has a pair of purely imaginary roots, and all other roots of Eq. (2.8) have negative real parts. Following [20,31], system (1.5) undergoes a Hopf bifurcation at k 0 , which completes the proof. (b) There exists k 0 > 0, when We know that D j 2 = 0 has two roots j 2 = k 2 0 , j 2 = hence, when j > k 0 , then D j > 0. Hence Eq. (2.8) has a double zero root, and all other roots of Eq. (2.8) have negative real parts. Following [30,32], system (1.5) undergoes a Bogdanov-Takens bifurcation at k 0 , which completes the proof.

Numerical simulations
In this section, we illustrate several conclusions by numerical simulations with Matlab. In model (1.5), we choose the following parameters: Then model (1.5) has a unique positive equilibrium (σ , v σ ) ≈ (0.2598, 1.1585). By direct computation, we have h * = 0.7697,h = 0.4201. By Theorem 2.2, we know that (σ , v σ ) is locally asymptotically stable whenh < h < h * (shown in Fig. 1), the positive solution (σ , v σ ) is unstable for 0 < h <h. By Theorem 3.1, we know that when h crossesh = h H 0 , the positive solution (σ , v σ ) loses its stability and Hopf bifurcation occurs (see Fig. 2).  In model (1.5), we choose the following parameters: Model (1.5) has a unique positive equilibrium (σ , v σ ) ≈ (0.2486, 0.3885). By direct computation, we have h * = 0.6226,h = 0.3918. By Theorem 3.2, we know that the diffusive system can undergo the Bogdanov-Takens bifurcation, this indicates that diffusion plays an important role in leading to complex dynamic behaviors. In this paper, the universal unfolding of the system near the Bogdanov-Takens bifurcation point is still not obtained, but we found that system (1.5) can exhibit some complex dynamic behaviors by numerical simulations such as quasi-periodic solutions (shown in Fig. 3).

Conclusion
In this paper, we investigated the dynamics of a diffusive phytoplankton-zooplankton model with nonlinear harvesting. By choosing different parameters, we obtained the parameter ranges of the existence of bifurcations. We have shown that the harvesting and the diffusion have a combined effect on the dynamic behaviors of the system. According to the analysis of the above, we know that if the harvesting h overruns h * , then zooplankton will die out. It means overfishing can break the coexistence of phytoplankton and zooplankton. We have also shown that the system can undergo Hopf bifurcation, Bogdanov-Takens bifurcation, when the parameters of the harvesting and diffusion cross certain critical values. It turns out that, by adjusting the parameters of the harvesting and diffusion, the planktonic ecological system can develop towards a healthy direction which can avoid "red tide". Furthermore, we also studied the stability of Hopf bifurcation by applying the normal form theory and the center manifold theorem.
According to Sect. 3, we know that the diffusive system can exhibit Bogdanov-Takens bifurcation (see Theorem 3.2). Meanwhile, the system may have a quasi-periodic solution near the Bogdanov-Takens bifurcation, which has been verified by some numerical simulations (see Fig. 3). This indicates that diffusion and harvesting can increase the dynamic complexity of system (1.5). The universal unfolding of system (1.5) near the Bogdanov-