I. INTRODUCTION

Marine phytoplankton and zooplankton are essential components of marine ecosystems and support the regular operation of the entire marine ecosystem. The research of marine phytoplankton and animal ecology is conducive to our comprehensive understanding of the status of an aquatic ecosystem. Marine plankton refers to the aquatic organisms suspended in the water and moving with water flow, mainly including phytoplankton and zooplankton, as well as other organisms such as planktonic viruses, planktonic bacteria, and archaea. Phytoplankton is the primary producer in the sea; it converts solar energy into organic energy through photosynthesis, initiates the material circulation and energy flow in the sea, and is the most basic link in the marine food chain. Zooplankton is an essential consumer in the sea; this part of organic matter is utilized through the food chain and further transferred to the upper trophic level through secondary production processes. Therefore, phytoplankton and zooplankton provide food and energy sources for the upper trophic level organisms through the above primary and secondary production processes, supporting the regular operation of the entire marine ecosystem.

Phytoplankton is not only the bottom but also the most crucial component of the marine ecosystem. It is divided into toxic and non-toxic phytoplankton. At the same time, zooplankton can distinguish different types of phytoplankton. To avoid feeding on toxic phytoplankton, which has a similar synergistic behavior with selective grazing in the predator-prey system [1-5]. In marine plankton ecosystems, the hypothetical mechanisms of selective grazing include prey morphology (size, color, shape, and colony formation), intestinal genetic strains, and poisonous chemicals released by prey [6-12]. Thus, the avoidance effect of zooplankton on toxins from toxic phytoplankton and the harmful effects of toxic compounds released by toxic species on their competitors have been studied [13-20].

In this paper, we consider not only the effect of toxin avoidance on species existence, but also the impact of human beings on the harvest of non-toxic phytoplankton and zooplankton is considered, whereas non-toxic phytoplankton on species existence and the human harvest has been applied in many models [21-27]. Since time delay is widely studied in the phytoplankton-zooplankton model [28-31], another essential purpose of our research is to explore the effect of pregnancy delay and toxin onset delay on the dynamic system. Finally, we find that optimal strategies are applied in many models to constrain overfishing [32-33]. Through the research we know that in fisheries, there is a fishing strategy called specific fishing, that is, fishermen catch almost only one particular type of fish or several species associated with it, such as these three species in our article, so we need a feedback mechanism to control this particular capture. Based on the dual phytoplankton-zooplankton system, we consider the optimal tax policy to constrain this particular fishing.

The organizational structure of this paper is as follows. In Section 2, we establish a mathematical model with double time delays for avoiding toxic species by zooplankton in the presence of non-toxic species. And give a parameter explanation in Table 2. In Section 3, we analyze the boundedness and stability of the boundary equilibrium point and the internal equilibrium point in the delay-free model. And obtain the bistability between the equilibrium points. The results are summarized in Table 1 and Fig 1. In Section 4, by analyzing different situations of this double delay model, we obtain the critical value of time delay when the system undergoes Hopf bifurcation. In Section 5, we study the optimal tax policy without time delay using the principle of Pontryagin's maximum. In addition, we use the parameters and initial values given in Table 2 and (6.1) to simulate several cases of double-delay systems in Matlab to verify all theoretical results in Section 6. Lastly, we end this paper with some conclusions and significance in Section 7.

II. MODEL FORMULATION

Considering the toxin refuge of zooplankton, a nontoxic phytoplankton-toxic zooplankton model was proposed in [14]. They showed that avoidance effects can promote the coexistence of non-toxic phytoplankton, toxic phytoplankton and zooplankton. Which can be shown as(with symbols slightly varied):

where $N$, $T$, and $Z$ represent the biomass of nontoxic phytoplankton, toxic phytoplankton, and zooplankton, respectively. $k_1$ and $k_2$ are the environmental carrying capacities of nontoxic phytoplankton (NTP) and toxin-producing phytoplankton (TPP) species, respectively. $r_1$ and $r_2$ represent the constant intrinsic growth rates of $N$ and $T$, respectively. ${\alpha }_1$ and ${\alpha }_2$ measure the competitive effect of $T$ on $N$, and $N$ on $T$, respectively. $w_1$ and $w_2$ represent the rates at which $N$ and $T$ are consumed by $Z$, respectively. $p_1$ and $p_2$ are half-saturation constants for NTP and TPP, respectively. $\beta$ represents the intensity of avoidance of $T$ by $Z$ in the presence of $N$, and $d$ is the natural mortality of zooplankton. As the research merely focuses on a single time model, moreover overfishing has an important impact on the stability of marine ecosystems, human harvest and time delays should be taken into account. The increment in zooplankton population due to predation does not appear immediately after consuming phytoplankton; it takes some time(say ${\tau }_1$ ), which can be regarded as the gestation period in zooplankton. The decrease of zooplankton population caused by ingestion of toxic phytoplankton does not occur immediately. Still, it requires a certain time(say ${\tau }_2$ ), which can be regarded as the reaction time after zooplankton poisoning. Accordingly the bio-economic model with time delays on the interactions of nontoxic phytoplankton, toxic plankton and zooplankton with toxin avoidance effects, which can be shown as follows:

Notes where $N$, $T$, and $Z$ represent the biomass of nontoxic phytoplankton, toxic phytoplankton and zooplankton, respectively. and represent the maturation gestation delay and the toxin onset delay, respectively. $c_1$ and $c_2$ represent the conversion rate of $N$ to $Z$ and $T$ to $Z$, respectively. Due to the experience of human capture, we assume that humans can distinguish between toxic phytoplankton and non-toxic phytoplankton when capturing zooplankton and phytoplankton. So, we put $q_1$ and $q_2$ to represent the fishing coefficients of nontoxic phytoplankton and zooplankton, respectively. And $E$ is the effort used to harvest the population. To investigate the effect of time delay on the dynamic behavior of the model, we will first study the stability of the equilibrium point of the following model without time delay.

III. DYNAMICAL BEHAVIOR OF NON-DELAYED MODEL

a) Positivity and boundedness of the solution

In this subsection, firstly, we shall show the positivity and boundedness of solutions of the system (2.3), which is vital for the biological understanding of the system and the subsequent analysis.

Lemma 3.1. All the solutions with initial values of system (2.3), which start in $R_{+}^3$, are always positive and bounded.

Proof. Firstly, we rewrite the model (2.3) and take the linear as the following form:

where $X(t)=(N,T,Z{)}^T\in R_{+}^3$ and $F(X)$ is simplified as the following

We want to prove that $(N(t),T(t),Z(t))\in R_{+}^3$ for all $t\in [0,+\mathrm{\infty })$. For system (2.3) with initial value , and , we have

which shows that all the solutions of system (2.3) are always positive for all .

Secondly, we prove the boundedness of the solution. Let $(N(t),T(t),Z(t))$ be the solutions of system (2.3), we define a function

Then, by differentiating (3.2) concerning $t$, we obtain

when , we can obtain

noting $\kappa =\frac{c_1{r}_2{k}_1(r_1+\eta {)}^2+c_2{r}_1{k}_2(r_2+\eta {)}^2}{4r_1{r}_2}$, therefore, applying a theorem on differential inequalities [34], we obtain $0\le W\le \frac{\kappa }{\eta }+\frac{W(N(0),T(0),Z(0))}{e^{\eta t}}$, let $t\to +\mathrm{\infty }$, $W(N,T,Z)\le \frac{\kappa }{\eta }$. So, all solutions of system (2.3) enter the region

This shows that every solution of the system is bounded.

b) Equilibrium points and their stability

System (2.3) possesses six different equilibrium points:

• (i) the plankton-free equilibrium, $E_0=(0,0,0)$, which always exists;
• (ii) TPP and zooplankton-free equilibrium, $E_1=(k_1,0,0)$, which is always feasible;
• (iii) NTP and zooplankton-free equilibrium, $E_2=(0,k_2,0)$, which is always feasible;
• (iv) zooplankton-free equilibrium, $E_3=(\hat{N},\hat{T},0)$, where

(v)TPP-free equilibrium $E_4=(\overline{N},0,\overline{Z})$, where

(vi) the interior equilibrium, $E^{\ast }=(N^{\ast },T^{\ast },Z^{\ast })$, where

and $N^{\ast }$ can be obtained from

Next, we illustrate the existence and stability of six equilibria when human harvest and avoidance factor exist simultaneously by solving Jacobi determinant of different equilibria, and summarize them in Table 1.

Equilibria analysis: Obviously, the equilibria $E_0$, $E_1$ and $E_2$ always exist. The zooplankton-free equilibrium $E_3$ exists, let $\hat{N}$ and $\hat{T}$ both be positive, that is and . The TPP-free equilibrium $E_4$ exists, let $\overline{N}$ and $\overline{Z}$ both be positive, that is and . The interior equilibrium point $E^{\ast }$ exists; let $N^{\ast }$, $T^{\ast }$ and $Z^{\ast }$ all be positive, that is , and Eq.(3.4) has at least one positive root.

In the following, we summarize the eigenvalues and local stability conditions around the feasible equilibrium point of each organism of system (2.3).

• (i) The eigenvalues of the plankton-free equilibrium $E_0=(0,0,0)$ are $r_1,r_2$ and $-d-q_{2}E$. Therefore, it is a saddle point and hence always unstable.
• (ii) The eigenvalues of the TPP and zooplankton-free equilibrium $E_1=(k_1,0,0)$ are $-r_1-q_{1}E$, $r_2(1-\frac{k_1{\alpha }_2}{k_2})$ and $\frac{c_1{w}_1{k}_1}{p_1+k_1}-d-q_{2}E$. When $c_1{\omega }_1-d-q_{2}E\le 0$, and hold, $E_1$ is LAS(locally asymptotically stable). On the contrary, if , and hold, we can also obtain $E_1$ is LAS.
• (iii) The eigenvalues of the NTP and zooplankton-free equilibrium $E_2=(0,k_2,0)$ are $r_2(1-\frac{k_2{\alpha }_1}{k_1})-q_{1}E,-r_2$ and $-\frac{c_2{w}_2{k}_2}{p_2+k_2}-d-q_{2}E$, Therefore, $E_2$ is LAS if .
• (iv) The eigenvalues of the zooplankton-free equilibrium $E_3=(\hat{N},\hat{T},0)$ are $\frac{c_1{w}_{1}N}{p_1+\hat{N}}-\frac{c_2{w}_{2}T}{p_2+\hat{T}+\beta \hat{N}}-d-q_{2}E$, ${\lambda }_1$ and ${\lambda }_2$, where ${\lambda }_1$ and ${\lambda }_2$ are the roots of the equation

Therefore, let , ${\lambda }_1$ and ${\lambda }_2$ with negative real parts, that is , and . If the above conditions are satisfied, $E_3$ is LAS.

(v) The eigenvalues of the TPP-free equilibrium $E_4=(\overline{N},0,\overline{Z})$ are $r_2(1-\frac{{\alpha }_2\overline{N}}{k_2})-\frac{w_2\overline{Z}}{p_2+\beta N}$, ${\overline{\lambda }}_1$ and ${\overline{\lambda }}_2$, where ${\overline{\lambda }}_1$ and ${\overline{\lambda }}_2$ are the roots of the equation

where

Therefore, let , ${\overline{\lambda }}_1$ and ${\overline{\lambda }}_2$ with negative real parts, that is and . If the above conditions are satisfied, $E_4$ is LAS.

(vi) By solving the Jacobi determinant of $E^{\ast }$, we can get its characteristic equation as follows

The interior equilibrium $E^{\ast }=(N^{\ast },T^{\ast },Z^{\ast })$ is LAS if

• (a)
• (b)
• (c)

From the calculation of the eigenvalues, obviously, $\beta$ does not affect the stability of $E_1$ and $E_2$. Still, it has a significant impact on the stability of $E_3$ and $E_4$ (because the eigenvalues of $E_1$ and $E_2$ are independent of $\beta$, but related to human harvest). On the other hand, we not only find that the equilibrium point of system (2.3) is affected by human harvest, but also has a particular impact on its stability (it can be seen from the eigenvalue of each equilibrium point).

Next, the biological explanations of the above different equilibria are discussed below. Since all these interpretations are mainly based on local asymptotic stability conditions, initial abundance of all the populations may also play an essential role for the system's dynamics together with the parameters. Different from the biological explanation in [14], we not only consider the effect of $\beta$ on species coexistence, but also human harvest as an essential factor in species coexistence.

• (i) $E_0$: Extinction of all the populations at a time is impossible.
• (ii) $E_1$: From the analysis of research results, whenever the carrying capacity of the NTP population $(k_1)$ stays within the specific threshold values of , both TPP and zooplankton will eventually become extinct from the system. Now, through the analysis of the $k_1$ threshold range, as the intensification of the harvest for zooplankton, the equilibrium point $E_1$ remains stable for a more extensive range of $k_1$, and we can say that over-fishing of zooplankton $(q_{2}E)$ may accelerate the extinction of TPP and zooplankton.
• (iii) $E_2$: If the carrying capacity of NTP population $(k_1)$ stays below the threshold value $\frac{r_2{\alpha }_1{k}_2}{r_2-q_{1}E}$, both NTP and zooplankton eventually extinct. With the competitive effect of TPP on NTP $({\alpha }_1)$, the environmental carrying capacities of toxin-producing phytoplankton $(k_2)$ and harvesting term for NTP and zooplankton

$(q_{1}E)$ increase, respectively. The equilibrium point $E_2$ remains stable for a larger scale of $k_1$; we can say that the possibility of deracinating NTP and zooplankton at a time increases with the increase in ${\alpha }_1$, $k_2$ and $q_{1}E$.

(iv) $E_3$: When the carrying capacity of NTP population $(k_1)$ remains within two threshold values (it can be obtained by the threshold value $(k_1)$ of $E_1$ and $E_2$ ) together with the competitive effects $({\alpha }_1,{\alpha }_2)$, the harvesting term on NTP $(q_{1}E)$ are present and the values of all three are small, the zooplankton population will go extinct on the condition that , whereas both NTP and TPP persist in the system. The chance of zooplankton extinction increases with the decrease in avoidance of TPP by zooplankton $(\beta )$, TPP consumption rate $(w_1)$, the half-saturation constant for TPP $(p_2)$, the harvesting term on zooplankton $(q_{2}E)$ and the zooplankton mortality $(d)$. For a specific parameter setup , we can find a threshold value of the avoidance of TPP by zooplankton , below which the zooplankton population will become extinct. On the contrary, for , the extinction of zooplankton dose not depend on the intensity of avoidance; it maybe has something relationship with the harvest term on zooplankton $(q_{2}E)$.

• (v) $E_4$: If the carrying capacity of NTP population $(k_1)$ remains within two threshold values , then TPP becomes extinct under the condition , whereas both NTP and zooplankton persist in the system. The possibility of TPP extinction increases with the reduction in the avoidance of TPP by zooplankton $(\beta )$, the half-saturation constant for TPP $(p_2)$, and the growth rate of TPP $(r_2)$, decreases with the rise of the competitive effect of $N$ on $T({\alpha }_2)$ and the TPP consumption rate $(w_2)$. Similarly, for a particular parameter setup , we can find a threshold value of the avoidance of TPP by zooplankton , below which TPP may become extinct. On the contrary, for , TPP extinction dose not depend on the avoidance. Because the biological analysis of $E_4$ found that the harvesting term has little impact on the extinction of TPP compared with other equilibrium points. In conclusion, for , TPP extinction dose not depend on the avoidance of TPP by zooplankton $(\beta )$ and harvest term on zooplankton $(q_{2}E)$.
• (vi) $E^{\ast }=(N^{\ast },T^{\ast },Z^{\ast })$: When the competitive effects $({\alpha }_1)$, the fishing coefficients of nontoxic phytoplankton $(q_1)$, the environmental carrying capacities of nontoxic phytoplankton $(k_1)$, and the effort used to harvest the population $(E)$ remain very small, whereas the constant intrinsic growth rates of $N(r_1)$, there may be a possibility of coexistence of all the three species.

Table 1: Existence and stability conditions of the equilibrium points.
Equilibrium	Existence conditions	Stability conditions
E0=(0,0,0)	Always exist	Always unstable
E1=(k1,0,0)	Always exist	(i)c1w1-d-q2E>0,α2>k2/k1,k1< p1(d+q2E)/c1w1-d-q2E,or (ii)c1w1-d-q2E≤0,α2>k2/k1
E2=(0,k2,0)	Always exist	(i) k1<r2α1k2/r2-q1E
E3=(ˆN,ˆT,0)	(i) α2>k2/k1,(ii) α1=(α1α2-1)q1k1E+k1/k2	(i) c1w1ˆN-d-q2E<c2w2ˆT/p2+ˆT+βˆN,(ii) b1>0,ˆc1>0
E4=(ˆN,0,ˆZ)	(i) w1>d+q2E/c1,(ii) k1=r1N/r1-ˆq1E(p1+E)	(i) r2(1-α2ˆN)/p2+βˆN,(ii) a2+ˆb2<0, a2ˆb2+ˆc2>0
E*=(N*,T*,Z*)	(i) k1=q1k1E/r1+N*+α1T*,(ii) c2w2(p1+N*)>c1w1N*(d+q2E)(p1+N*)>0,(iii) positive root of Eq.(3.4) exists	(i) D1>0,(ii) D3>0,(iii) D1D2-D3>0

c) Bistability analysis of equilibrium point

The existence and stability of these equilibrium points are summarized in Table 1 and Fig 1. When , equilibria $E_2=(0,k_2,0)$, $E_3=(\hat{N},\hat{T},0)$, $E_1=(k_1,0,0)$ and $E_4=(\overline{N},0,\overline{Z})$ keep stable for , , and , respectively(Fig.1(a)). Obviously, for $k_1$ at the different equilibria above, the coexistence of NTP, TPP, and zooplankton requires the three ranges , , and , respectively. Therefore, the system exhibits these three possible types of bistability, where

• (i) $E_1$ and $E_2$
• (ii) $E_2$ and $E_4$
• (iii) $E_3$ and $E_4$

The above three types are locally asymptotically stable for different ranges of $k_1$.

For , we can observe the bistability of $E_1$ and $E_2$ (Fig.1(b)(c)). If conditions and hold simultaneous, we can find the bistability of $E_2$ and $E_4$ (Fig.1(d)(e)). On the contrary, if holds, for either or we'll get the existence of stable $E_2$ together with unstable $E_4$. Identically, for together with , and , we can observe the bistability of $E_3$ and $E_4$ (Fig.1(f)-(i)).

Now, let's discuss the importance of avoiding toxic species by zooplankton $(\beta )$ together with the harvesting term $(q_{1}E,q_{2}E)$ for the survival of the different species groups.

Firstly, let's discuss the effect of $\beta$ on three types of bistability. It can be seen from the previous analysis that the stability of $E_1$ and $E_2$ does not depend on the value of $\beta$. However, for the stability of $E_3$ and $E_4$, it is related to the critical value of $\beta$. When $\beta$ is less than this critical value, $E_3$ and $E_4$ remain stable. Thus, $\beta$ does not affect the bistability of $(E_1,E_2)$; when $\beta$ is below some threshold value, we will observe the bistability of $(E_2,E_4)$ and $(E_3,E_4)$, and as the $\beta$ value increases, the original bistability may disappear. and . From these conditions, we can see the establishment of the above conclusion.)

Secondly, let's discuss the effect of the harvesting term $(q_{1}E,q_{2}E)$ on three types of bistability. From the analysis of the previous data, it can be seen that although the stability of $E_1$ and $E_2$ does not depend on the value of $\beta$, when humans overfish NTP and zooplankton, that is, $q_{1}E$ and $q_{2}E$ are too large, it may affect the bistability of $E_1$ and $E_2$. For $E_3$ and $E_4$, although their stability is directly related to the threshold value of $\beta$, the existence of $q_{1}E$ and $q_{2}E$ will also affect the threshold value of $\beta$, further influencing the stability of $E_3$ and $E_4$. Therefore, $q_{1}E$ and $q_{2}E$ may affect the bistability of $(E_1,E_2)$, $(E_2,E_4)$ and $(E_3,E_4)$; the increase of $q_{1}E$ and $q_{2}E$ may also lead to the disappearance of this bistability.

IV. DYNAMICAL BEHAVIOR OF THE DELAYED MODEL

In this section, we focus on the local stability and Hopf bifurcation of the delayed model; the delayed system (2.2) has the following form

Next, assuming ${\mathrm{\Lambda }}_1(t)=N(t)-N^{\ast }$, ${\mathrm{\Lambda }}_2(t)=T(t)-T^{\ast }$, ${\mathrm{\Lambda }}_3(t)=Z(t)-Z^{\ast }$ at the positive equilibrium point, and linearizing the system (2.2), we can obtain

where

We linearize the system(2.2) about positive equilibrium $E^{\ast }=(N^{\ast },T^{\ast },Z^{\ast })$, and get

Notes

Fig. 1: Stability of different equilibria for different ranges of $k_1$. The dotted arrow indicates the range where bistability occurs, (a) means no bistability, (b) and (c) bistability of $E_1$ and $E_2$, (d) and (e) bistability of $E_2$ and $E_4$, (f)-(i) bistability of $E_3$ and $E_4$.

where $N_1,T_1,Z_1$ are small perturbations around the equilibrium point $E^{\ast }=(N^{\ast },T^{\ast },Z^{\ast })$. We have

The characteristic equation for the linearized system (2.2) is obtained as

where

with

Case (1): ${\tau }_1={\tau }_2=0$

In this case, Section 3 covers the analysis of the system when ${\tau }_1={\tau }_2=0$.

Case (2): .

In this case, the characteristic equation(4.4) becomes

putting in Eq.(4.5), and separating the real and imaginary parts, we have

Squaring and adding the equation(4.6), we obtain

Simplifying Eq.(4.7) and substituting ${\omega }^2=\psi$, the above equation can be written as

where

If (H1) holds, Eq.(4.8) has no positive roots, which implies all the roots of Eq.(4.5) have negative real parts. Therefore, $E^{\ast }$ is asymptotically stable for all when (H1) holds.

If (H2) holds, Eq.(4.8) has exactly one positive root ${\omega }_0$, substituting ${\omega }_0$ in Eq.(4.6), we obtain

For the critical value of ${\tau }_2$, we can obtain

For the transversality condition, differentiating Eq.(4.5) with respect to ${\tau }_2$, we get

Solving ${\left(\frac{d\xi }{d{\tau }_2}\right)}^{-1}$, we obtain

Then at ${\tau }_2={\tau }_{20}$ and $\xi =i{\omega }_0$, we can get

Then

here

From this, we can get

If (H3): $AD+BC\ne 0$ holds, the transversal condition $\mathrm{sgn}[\mathrm{Re}(\frac{d\xi }{d{\tau }_2}{)}_{{\tau }_2={\tau }_{20},\xi =i{\omega }_0}{]}^{-1}\ne 0$. From the above analysis, the following theorem can be drawn

Theorem 4.1. For ${\tau }_1=0$ and , we have the following results:

(i) If (H1) holds, then the equilibrium $E^{\ast }$ is asymptotically stable for all .

(ii) If (H3) holds, and (H2) holds, then the equilibrium $E^{\ast }$ is locally asymptotically stable for all together with unstable for and undergoes Hopf bifurcation at ${\tau }_2={\tau }_{20}$.

Case (3): .

In this case, the characteristic equation(4.4) becomes as follows

putting in Eq.(4.12), and separating the real and imaginary parts, we have

Squaring and adding the equation(4.13), we obtain

Based on the calculation method for case (2), we can simplify (4.14) to the following

• (H4): .
• If (H4) holds, Eq.(4.15) has no positive roots, which implies all the roots of Eq.(4.12) have negative real parts. Therefore, $E^{\ast }$ is asymptotically stable for all when (H4) holds.
• (H5): or or .
• If (H5) holds, Eq.(4.15) has exactly one positive root ${\hat{\omega }}_0$, substituting ${\hat{\omega }}_0$ in Eq.(4.13), we obtain

For the critical value of ${\tau }_1$, we can obtain

For the transversality condition, differentiating Eq.(4.13) with respect to ${\tau }_1$, we get

Solving ${\left(\frac{d\xi }{d{\tau }_1}\right)}^{-1}$, we obtain

Then at ${\tau }_1={\tau }_{10}$ and $\xi =i{\omega }_0$, we can get

Now

From this, we can get

If (H6): $A_{\ast }{D}_{\ast }+B_{\ast }{C}_{\ast }\ne 0$ holds, the transversal condition $[\mathrm{Re}(\frac{d\xi }{d{\tau }_1}{)}_{{\tau }_1={\tau }_{10},\xi =i\hat{\omega }}{]}^{-1}\ne 0$. From the above analysis, the following theorem can be drawn

Theorem 4.2. For ${\tau }_2=0$ and , we have the following results:

(i) If (H4) holds, then the equilibrium $E^{\ast }$ is asymptotically stable for all .

(ii) If (H6) and (H5) hold, then the equilibrium $E^{\ast }$ is locally asymptotically stable for all together with unstable for and undergoes Hopf bifurcation at ${\tau }_1={\tau }_{10}$.

Case (4): ${\tau }_1$ is fixed in $(0,{\tau }_{10}]$ and .

We consider the gestation delay ${\tau }_1$ to be stable in the interval $(0,{\tau }_{10}]$, taking ${\tau }_2$ as a control parameter. Let $\xi =u+i\omega$ be the root of Eq.(4.4). Putting this value in Eq.(4.4), separating real and imaginary parts, we obtain

Putting $u=0$ in Eqs.(4.19) and (4.20), we obtain

Squaring and adding Eqs.(4.21) and (4.22) to eliminate ${\tau }_2$, we get

Noting that Eq.(4.23) is transcendental. Now, Eqs.(4.21) and (4.22) can be written as

where

Without losing generality, the Eq.(4.23) has finite positive roots $\widetilde{{\omega }_1},\widetilde{{\omega }_2},\dots ,\widetilde{{\omega }_k}$, for every fixed $\tilde{\omega }$, there exists a sequence $\{{\tau }_{2i}^j|j=0,1,2\dots \}$, where

let ${\tilde{\tau }}_2=min\{{\tau }_{2i}^{(j)}|i=0,1,2,\dots ,k,j=0,1,2\dots \}$, when ${\tau }_2={\tilde{\tau }}_2$, $\tilde{\omega }={\tilde{\omega }}_i{|}_{{\tau }_2={\tilde{\tau }}_2}$, $i=1,2,3,\dots$, the characteristic equation (4.4) has purely imaginary roots $\pm i\tilde{\omega }$. Then, we will verify the transversality condition, differentiating the characteristic equation (4.4) with respect to ${\tau }_2$, we can obtain

Now

where

From this we can get

If (H7): $EH+FG\ne 0$ holds, the transversal condition $\mathrm{sgn}[\mathrm{Re}(\frac{d\xi }{d{\tau }_2}{)}_{{\tau }_2={\tilde{\tau }}_2,\xi =i\tilde{\omega }}{]}^{-1}\ne 0$. From the above analysis, we have the following theorem.

Theorem 4.3. For system(2.2), assume (H7) holds with ${\tau }_1$ is fixed in $(0,{\tau }_{10}]$ and , then the equilibrium $E^{\ast }$ is locally asymptotically stable for ${\tau }_2\in (0,{\tilde{\tau }}_2)$ whereas system (2.2) undergoes Hopf bifurcation at ${\tau }_2={\tilde{\tau }}_2$.

Case(5): ${\tau }_2$ is fixed in $(0,{\tau }_{20}]$ and , so take ${\tau }_1$ as a control parameter; the analysis is the same as case(4), so we omit it.

V. OPTIMAL TAX POLICY

From previous studies, overfishing may lead to the extinction of populations. However, in the society, the adequate protection of the ecosystem is a common problem we need to face. In the face of the increasingly severe harmful effects of overfishing on ecosystems, people began to find the most suitable methods for fishery control in various areas of sustainable development policies, for example, seasonal fishing, property leasing, taxation, licensing fees, etc. Taxes are generally considered to be better than other regulatory approaches, so that we will view the optimal tax policy for the double phytoplankton - single zooplankton system based on model (2.3). Here, we take $E$ as a time-dependent dynamic variable controlled by equations. Therefore, there is the following equation.

Where $Q(t)$ is the amount of capital invested in fisheries at time $t$, $I(t)$ is the total investment rate (in physical form) at time $t$ and $\gamma$ is the constant depreciation rate of capital. Suppose that the effort $E$ at any time is proportional to the instantaneous amount of investment capital. For example, if $Q(t)$ represents the number of standard fishing vessels that can be used, it is reasonable to assume that $Q(t)$ and $E$ should be proportional. When $\epsilon =1$, it can be considered that the maximum fishing capacity $(E)$ is equal to the number of available vessels at time $t$ ( $Q(t)$ ). When $\epsilon =0$, it means that even though there may be fishing boats, the fishing is not expanded; it also reflects the over-exploitation of fisheries. At this time the fish population has been seriously depleted, so fishing vessels can no longer be used. These are simulated capital levels may be adjusted, thus prove the reasonableness of the equation (5.2). Regulators control the development of fisheries by imposing a tax ( ) on the unit biomass of terrestrial fish. When ( ) can be understood as any subsidy to fishermen. Net income of fishermen ('Net income' for short) is $E[(u_1-v)q_{1}N+(u_2-v)q_{2}N-C]$, where $u_i$, $i=1,2$ is the constant price of unit biomass of nontoxic phytoplankton and zooplankton, respectively. $C$ is the fixed cost per unit of harvesting effort.

We assume the gross profit margin on capital investment is proportional to this 'Net income.' So, we have

For $\phi =1$, Eq.(5.2) shows that the highest investment rate at any time is equal to the net income of the fishermen at that time. $\phi =0$ can only be used when the net income of fishermen is negative; that is, current capital assets cannot be divested. If the fishery is operating at a loss and allows capital to be withdrawn, the only owner of the fishery will benefit by allowing the capital assets to be continuously withdrawn, because negative investment means withdrawal of investment, so it is the case of , . By combining Eqs.(5.1) and (5.2), we can get

Fishermen and regulators are two different parts of society. Therefore, the income they receive is society's income accumulated through fisheries. The net economic income to society is

this is equal to the net economic income of fishermen plus the economic income of regulators. Therefore without considering the time delay, Eq.(2.3) can be rewritten as

Next, we will use the principle of Pontryagin's maximum to get the path of the best tax policy. If the fish population stays along this path, then regulators can ensure that their goals are achieved. The goal of regulatory agencies is to maximize the total net income of society as a result of harvesting activities. Specifically, the goal is to maximize revenue over a continuous time stream $(J)$.

where $\delta$ is the discounting factor. Therefore, our goal is to determine an optimal tax $v=v(t)$ that maximizes compliance with Eq.(5.4) and constrains $v_{\min }\le v(t)\le v_{\max }$ on the control variable $v(t)$. When , it will have the effect of accelerating the rate of fishery expansion. The Hamiltonian of the problem is obtained by

where ${\lambda }_1,{\lambda }_2,{\lambda }_3$ and ${\lambda }_4$ are the adjoint variables. For $v\in [v_{min},v_{max}]$, the Hamiltonian must be maximized. Assuming that the control constraint is not bound, that is, the optimal solution does not appear as $v=v_{min}$ or $v=v_{max}$. We can get by singular control [9]

Now, the adjoint equations are