Bifurcation analysis of a Singular Nutrient-plankton-fish model with taxation, protected zone and multiple delays

A differential algebraic nutrient-plankton-fish model with taxation, free fishing zone, protected zone and multiple delays is investigated in this paper. First, the conditions of existence and control of singularity induced bifurcation are given by regarding economic interest as bifurcation parameter. Meanwhile, the existence of Hopf bifurcations are investigated when migration rates, taxation and the cost per unit harvest are taken as bifurcation parameters respectively. Next, the local stability of the interior equilibrium, existence and properties of Hopf bifurcation are discussed in the different cases of five delays. Furthermore, the optimal tax policy is obtained by using Pontryagin's maximum principle. Finally, some numerical simulations are presented to demonstrate analytical results.


1.
Introduction. In the aquatic ecosystem, there are abundant micro organisms called plankton, which are composed of phytoplankton and zooplankton. During the recent years, the related results of zooplankton-phytoplankton models have been given by many scholars [3,4,7,19,25,32,35,43].
In fact, nutrients plays an important role in the growth of plankton, therefore, a plankton-nutrient model is worth investigating. He and Ruan [16] proposed plankton-nutrient interaction models, and gave some conditions for the global attractivity of the positive equilibrium. Pardo [30] investigated a phytoplankton nutrient system and analyzed the global behavior of system. Wang et al. [36] put forward a nutrient-plankton model with delayed nutrient cycling, and studied stability and the existence of global Hopf bifurcation at the positive equilibrium.
Considering the coexistence of nutrient, plankton and fish, the following model is proposed based on [36], where N (t) is the concentration of nutrients that can be absorbed. P (t), Z(t) and F (t) denote the density of nutrient, phytoplankton, zooplankton and fish at time t, respectively. The parameters N 0 and D are constant input rate and the washout rate for nutrient, respectively. α, β and β 2 are the capture rates from nutrient to phytoplankton, from phytoplankton to zooplankton and from zooplankton to fish, respectively. α 1 , β 1 and β 4 denote the ratio of biomass conversion (satisfying the obvious restriction 0 < α 1 < α, 0 < β 1 < β and 0 < β 4 < β 2 ). d 1 , d 2 and d 3 denote natural death rate for phytoplankton, zooplankton and fish population, respectively. m 1 and m 2 represent death phytoplankton and zooplankton are decomposed into detritus and converted into nutrients. α 2 is the interspecific competition coefficient of the phytoplankton.
Although a lot of research has been conducted over the past few years, very few results are interactions between nutrient, plankton and fish stocks. Fish, in general, is a good source of high-quality protein and contains many vitamins and minerals. Thus, it is consumed by many species, including humans all over the world. In addition, harvesting plays an important role not only in people's life, but also in the dynamics of biological resources [8,17,42]. As we know that the ecological balance of the oceans is being destroyed by overfishing. Thus, some measures must be implemented through relevant departments. Clearly, taxation and marine reserves are the effective measures to prevent overfishing. Some scholars have investigated preypredator systems with taxation or marine reserves, but a few of people introduce taxation and marine reserves into plankton models. Chakraborty et al. [2] described a prey-predator fishery system with prey free fishing zone, protected zone and taxation, and discussed dynamical behavior of the proposed system. Gupta et al. [12] studied period doubling cascades of prey-predator model with nonlinear harvesting, and obtained the conditions of optimal taxation control. Zhang et al. [41] proposed a prey-predator fishery model with reserves, selective harvesting and taxation, and gave an optimal taxation policy.
As we all know, many fishermen in coastal areas live mainly on fishing. However, in recent years, the marine ecosystem has been threatened due to human overfishing. Therefore, it is of great significance to study the impact of commercial gains and economic benefits on ocean systems.
The economic theory of common property resources is addressed by Gordon [10], The dynamic effects of harvesting on ecosystem were studied from the perspective of economics. In order to discuss the economic interest of commercial harvesting, the following equation is proposed Net Economic Revenue (NER) = Total Revenue (TR) − Total Cost (TC). Based on above equation, taxation is nonnegative and bounded i.e. 0 ≤ T < p 1 . Thus, we have where p 1 , c, v and T are the price per unit harvested biomass, the cost per unit harvest, the net economic revenue and the taxation per unit harvested biomass, respectively.
In addition, time delay is usually introduced into mathematical models [1,9,20,21,22,26,27,28,29,33,37,38,39,40] to reflect population interaction and coexistence mechanisms that depend on past history. Over recent years, more and more scholars addressed some singular models with delays in order to investigate the dynamic effects of commercial harvesting on the economic benefits of prey-predator system [20,28,33]. However, few people consider differential algebraic nutrientplankton-fish model with multiple delays. Santra et al. [33] presented a singular predator-prey model with concept super predator, and discussed the conditions of the existence and control singularity induced bifurcation. Meng and Wu [28] put forward a delayed differential algebraic phytoplankton-zooplankton-fish model with taxation and nonlinear fish harvesting, and gave the conditions of singularity induced bifurcation, Hopf bifurcation and optimal tax policy.
Based upon above analysis, by combining the economic interest, taxation and multiple delays, system (2) can be expressed by differential algebraic system where τ 1 and τ 2 are the discrete time period required for the decomposition of phytoplankton and zooplankton population, respectively. τ 3 is the discrete time period required for absorption of nutrients. τ 4 and τ 5 are the discrete time period required for the maturity of phytoplankton and zooplankton, respectively. The initial conditions of system (3) are

System (3) can be transformed into the following matrix form
where The rest of the paper is organized as follows. In Section 2, we analyze the existence of positive solutions when v = 0 and v > 0 respectively. In absence of time delay, we discuss the existence and control of singularity induced bifurcation by considering economic interest as bifurcation parameter. What's more, by taking the migration rates, taxation and the cost of unit harvesting as bifurcation parameters, we discuss the sufficient conditions of Hopf bifurcations. In Section 3, in the presence of time delay, the local stability, the existence and properties Hopf bifurcation of system are discussed based on some relevant theoretical knowledge. In Section 4, Pontryagin's maximum principle is applied to get optimal tax value that can assure maximize benefit and ecosystem balance. Numerical simulations are used to demonstrate the validity of the mathematical conclusion in Section 5. The discussion and conclusion are contained in the last section.
2. Bifurcation and control of system without delay.
2.1. The existence of the positive equilibrium.
For system (3), there is the phenomenon of bioeconomic equilibrium state when v = 0. First, we suppose Thus, if the assumption (H 1 ) holds, then the interior equilibrium is X 0 = (N 0 , P 0 , , and Z 0 satisfies the following equation here By using Routh-Hurwitz criterion [24], Eq.(5) has at least one positive root when M 3 < 0.

Singular induced bifurcation and control.
In this part, we discuss the stability of system (3) at the interior equilibrium X 0 .
When τ = 0, system (3) takes the following form Thus, we can get the following results.
Theorem 2.1. When the bifurcation parameter v increases through zero and , singular induced bifurcation of system (7) takes place around the interior equilibrium X 0 .
Proof. We can obtain the Jacobin matrix of system (7) around X 0 Further, we have that In addition, we can obtain that According to Lemma 4.1 and Theorem 3 in [34], when v increase through 0, one eigenvalue of system (7) moves from C − to C + along real axis by diverging through ∞. Hence, the Theorem 2.1 is proved.
Next, we will design state feedback controller to remove singularity induced bifurcation of system (7) around X 0 . When v = 0, we can calculate rank (J X0 ,ÂJ X0 , AJ X0 ) = 6 according to the leading matrixĀ of system (3) and J X0 . By using Theorem 2-2.1 in [6], we find that system (7) is locally controllable around X 0 . Therefore, based on Theorem 3-1.2 in [6], the following linear feedback controller can be designed to stabilize system (7) around X 0 , Applying such controller into system (7), we can obtain that a controlled system is For system (8), we obtain the following result.
Proof. At first, based on the Jacobin matrix of system (8) around X 0 , the characteristic equation of system (8) around X 0 is We suppose that the assumption is By using Routh-Hurwitz criterion [24], we implies the equilibrium of system (8) is locally asymptotically stable if the assumption (H 4 ) holds.
Remark 1. According to designing controller, we can stabilize the system at the equilibrium, which shows that both sustainable development of marine ecosystem and ideal economic benefits of fish harvesting can be sustained by adjusting commercial harvest effort on fish.

Hopf bifurcation.
In this subsection, we will discuss the existence of Hopf bifurcation aroundX of system (7) when δ 1 is taken as a bifurcation parameter.
3. Bifurcation and its properties of system with time delay.
As we know that the necessary condition of existence of Hopf bifurcation is critical eigenvalue passing through the imaginary axis with non-zero velocity [13]. Therefore, we will check transversal condition.
(ii) If the assumption (H 8 ) holds, then system (3) undergoes a Hopf bifurcation at the positive equilibriumX when τ = τ 0 . That is, system (3) has a branch of periodic solution bifurcating from the zero solution near τ = τ 0 .

BIFURCATION IN A DELAYED SINGULAR MODEL WITH TAXATION AND RESERVE 19
f (u t , µ) = (f 1 (u t , µ),f 2 (u t , µ),f 3 (u t , µ),f 4 (u t , µ),f 5 (u t , µ)) T , and Then L µ is a parameter family of bounded linear operator in C([−1, 0], R 5 ) → R 5 . By the Reisz representation theorem [18], there exists 5 × 5 matrix-valued function In fact, we can choose where δ(θ) is a Dirac delta function. Next, we define and Thus, system (39) can be rewritten as which is an equation of the form we desired.
In order to maintain the sustainable development of marine ecosystem, we apply the state feedback controllerū =k(P − P * ) (k is a feedback gain) to eliminate Hopf bifurcation. The range of k is evaluated based on Routh-Hurwitz criterion [24]. Of course, we only discuss the control of Hopf bifurcation by numerical simulation due to the complexity of theoretical analysis.

Remark 2.
According to the theory of Hopf bifurcation and stability, the sate feedback method can be used to eliminate those unfavorable phenomena and stabilize proposed marine ecosystems. So that biological population can be controlled at states. Remark 3. We first only take the time delay τ 2 as the bifurcation parameter, then take time delay τ 4 as the bifurcation parameter while the time delay τ 2 is taken as its stability interval. At last, we consider the five delays as the same to obtain the existence of Hopf bifurcation and study the properties of Hopf bifurcation under this case. In order to eliminate Hopf bifurcation, we can design the corresponding state feedback controller. We will only give some numerical results in the Section 5.

4.
Optimal tax policy. In this section, we will discuss the optimal tax policy of system (38). Our aim is to get the optimal tax policy that maximize the benefits of fish population without extinction. We noticed that as the tax rate increases, the harvest decreases. Therefore, the taxation plays a decisive role.
When τ 1 = τ 2 = τ 3 = τ 4 = τ 5 = 0, system (38) becomes system (7). We will consider not only the net economic revenue of the harvested but also the social economic revenue of the regulatory agency based on [4,2]. Hence, we havē The present value J(T ) of a continuous time-stream of revenues is given by where δ is the continuous annual discount rate which is fixed by harvesting agencies.
The control problem over an infinite time horizon is given by maxJ(T ) = ∞ 0 e −δtP dt, subject to system (7). We take advantage of the maximum principle to get the optimal solution of this problem. The associated Hamiltonian function H = H(t, N, P, Z, F 1 , F 2 , E, T ) is whereλ i =λ i (t), i = 1, 2, 3, 4, 5 are adjoint variables corresponding to the variables N, P, Z, F 1 , F 2 , E, respectively. The optimal control T which maximizes H must satisfy the following conditions: Assuming that the optimal solution does not occur at T = T max or T = 0 when ∂H ∂T = 0, we can obtainλ 5 (t) = 0. Based on Pontryagin's Maximum Principle [4], the adjoint variables must satisfy the following adjoint equations Therefore, we get In order to obtain an optimal equilibrium, we shall rewrite Eq.(51) as follows based on [5],λ i (t) = µ i (t)e −δt , i = 1, 2, 3, 4, 5.
From Eq.(52), we can get an equation for T . At the same time, we can get the optimal equilibrium solution X optimal . We only obtain numerical solution due to the complexity of theoretical analysis.
Similarly, we also can consider δ 2 , c, T as bifurcation parameters. When δ 1 = 0.25, some results at the positive equilibrium can be depicted in Figs.4-6. For example, from Fig.5, we can see that the critical value of c are c * l = 1.0, c * r = 1.7. That is, when the bifurcation parameter c is less than the value c * l = 1.0, the system is locally asymptotically stable at the positive equilibrium (see Fig. 5(f)). Hopf bifurcation occurs when c * l < c = 1 < c * r (see Fig.5(g)). Further, the system is locally asymptotically stable again once c = 2.0 is bigger than c * r = 1.7 (see Fig.5(h)). Similarly, some results about δ 2 and T can be seen in Fig.4 and Fig.6 respectively.
The Eq.(52) is too complex to solve the expression of T in theory. We only consider numerical example to find the optimal tax value which makes the the interior equilibrium be stable. When δ = 3, c = 1.5 dolars, δ 1 = 0.25, δ 2 = 0.5, v = 0.008 dolars, we get T = 1.61 dolars is a feasible solution of Eq.(52). Hence T optimal = 1.61 dolars and the corresponding point X optimal = (N ,P ,Z,F 1 ,F 2 ) = (2..49, 2.5, 0.49, 0.51, 1.10) is taken as the optimum equilibrium.  6. Conclusions and discussions. Generally, nutrients plays an important role in the growth of plankton. Therefore, a better plankton-nutrient model is worth investigating. In addition, as we know that the ecological balance of the oceans is being destroyed by overfishing. Thus, taxation and marine reserves are the effective measures to prevent overfishing. Most of the researchers have considered plankton-nutrient model. But few scholars study a singular nutrient-plankton-fish model including multiple delays, free fishing zone, protected zone and taxation as this paper. Without considering time delay, we discussed existence of singularity induced bifurcation by regarding economic interest as bifurcation parameter. In order to eliminate singularity induced bifurcation, we designed a linear state feedback controller u(t) = −1(E(t) − 1.5) (see Fig.1), which showed that both sustainable development of nutrient-plankton-fish ecosystem and ideal economic interest of     commercial harvesting could be adjusted. What's more, by taking migration rate, taxation and the cost per unit harvest as bifurcation parameters, we obtained some results that system (7) occurred Hopf bifurcation when δ 1 = δ * 1 = 2.35. System (7) is stable when 0 < δ 1 < 2.35 and becomes stable again when 5.93 ≤ δ 1 (see Figs. 2 and 3). While we chose δ 2 , c, T , we obtained some results (see Figs. 4-6).
With considering time delay, when τ 2 is fixed in [0, τ * 2 ), we obtained the sufficient conditions for Hopf bifurcation by considering τ 4 as a bifurcation parameter, and designed a state feedback controllerū = 0.1(P − P * ) in order to eliminate Hopf bifurcation (see Figs.7 (c) and (f)). By analyzing the associated characteristic transcendental equation of system with τ 1 = τ 2 = τ 3 = τ 4 = τ 5 = τ , we also obtained the corresponding critical time delay is 0.05. It is obvious that system loses local stability aroundX when time delays crosses its corresponding critical value, and Hopf bifurcation occurs (see Fig.8). And by using normal form theory and center manifold theorem, we knew that Hopf bifurcation was supercritical, the bifurcating periodic solutions existed for τ 4 > τ * 4 = 1.22, and bifurcating periodic solutions were stable. Pontryagin's maximum principle was applied to obtain optimal tax value T optimal = 1.61 dolars and optimal interior equilibrium X optimal = (2.49, 2.5, 0.49, 0.51, 1.10), and to maximize the economic benefit which is sum of the net economic revenue of the harvested and social economic revenue of the regulatory agency as well as conservation of the ecosystem.