## 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): $$ \left\{ \begin{array}{l} \frac {d N}{d t} = r _ {1} N \left(1 - \frac {N + \alpha_ {1} T}{k _ {1}}\right) - \frac {w _ {1} N Z}{p _ {1} + N}, \\ \frac {d T}{d t} = r _ {2} T \left(1 - \frac {T + \alpha_ {2} N}{k _ {2}}\right) - \frac {w _ {2} T Z}{p _ {2} + T + \beta N}, \\ \frac {d Z}{d t} = \frac {w _ {1} N Z}{p _ {1} + N} - \frac {w _ {2} T Z}{p _ {2} + T + \beta N} - d Z, \\ N (0) \geq 0, T (0) \geq 0, Z (0) \geq 0 \end{array} \right. \tag {2.1} $$ 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: $$ \left\{ \begin{array}{l} \frac {d N}{d t} = r _ {1} N \left(1 - \frac {N + \alpha_ {1} T}{k _ {1}}\right) - \frac {w _ {1} N Z}{p _ {1} + N} - q _ {1} E N, \\ \frac {d T}{d t} = r _ {2} T \left(1 - \frac {T + \alpha_ {2} N}{k _ {2}}\right) - \frac {w _ {2} T Z}{p _ {2} + T + \beta N}, \\ \frac {d Z}{d t} = \frac {c _ {1} w _ {1} N \left(t - \tau_ {1}\right) Z \left(t - \tau_ {1}\right)}{p _ {1} + N \left(t - \tau_ {1}\right)} - \frac {c _ {2} w _ {2} T \left(t - \tau_ {2}\right) Z \left(t - \tau_ {2}\right)}{p _ {2} + T \left(t - \tau_ {2}\right) + \beta N \left(t - \tau_ {2}\right)} - d Z - q _ {2} E Z, \\ N (0) \geq 0, T (0) \geq 0, Z (0) \geq 0 \end{array} \right. \tag {2.2} $$ Notes where $N$, $T$, and $Z$ represent the biomass of nontoxic phytoplankton, toxic phytoplankton and zooplankton, respectively. $\tau_{1}(\tau_{1} > 0)$ and $\tau_{2}(\tau_{2} > 0)$ 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. $$ \left\{ \begin{array}{l} \frac {d N}{d t} = r _ {1} N \left(1 - \frac {N + \alpha_ {1} T}{k _ {1}}\right) - \frac {w _ {1} N Z}{p _ {1} + N} - q _ {1} E N, \\ \frac {d T}{d t} = r _ {2} T \left(1 - \frac {T + \alpha_ {2} N}{k _ {2}}\right) - \frac {w _ {2} T Z}{p _ {2} + T + \beta N}, \\ \frac {d Z}{d t} = \frac {c _ {1} w _ {1} N Z}{p _ {1} + N} - \frac {c _ {2} w _ {2} T Z}{p _ {2} + T + \beta N} - d Z - q _ {2} E Z, \\ N (0) \geq 0, T (0) \geq 0, Z (0) \geq 0. \end{array} \right. \tag {2.3} $$ ## 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: $$ \frac {d X}{d t} = F (X), \tag {3.1} $$ where $X(t) = (N,T,Z)^T\in R_+^3$ and $F(X)$ is simplified as the following $$ F (X) = \left[ \begin{array}{c} F _ {1} (X) \\F _ {2} (X) \\F _ {3} (X) \end{array} \right] = \left[ \begin{array}{c} r _ {1} N (1 - \frac {N + \alpha_ {1} T}{k _ {1}}) - \frac {w _ {1} N Z}{p _ {1} + N} - q _ {1} E N \\r _ {2} T (1 - \frac {T + \alpha_ {2} N}{k _ {2}}) - \frac {w _ {2} T Z}{p _ {2} + T + \beta N} \\\frac {c _ {1} w _ {1} N Z}{p _ {1} + N} - \frac {c _ {2} w _ {2} T Z}{p _ {2} + T + \beta N} - d Z - q _ {2} E Z \end{array} \right]. $$ We want to prove that $(N(t), T(t), Z(t)) \in R_{+}^{3}$ for all $t \in [0, +\infty)$. For system (2.3) with initial value $N(0) > 0$, $T(0) > 0$ and $Z(0) > 0$, we have $$ N (t) = N (0) \exp \bigl \{\int_ {0} ^ {t} \bigl [ r _ {1} (1 - \frac {N (s) + \alpha_ {1} T (s)}{k _ {1}}) - \frac {w _ {1} Z (s)}{p _ {1} + N (s)} - q _ {1} E ] d s \bigr \}, $$ $$ T (t) = T (0) \exp \{\int_ {0} ^ {t} [ r _ {2} (1 - \frac {T (s) + \alpha_ {1} N (s)}{k _ {2}}) - \frac {w _ {2} Z (s)}{p _ {2} + T (s) + \beta N (s)} ] d s \}, $$ $$ Z (t) = Z (0) \exp \bigl \{\int_ {0} ^ {t} \bigl [ \frac {c _ {1} w _ {1} N (s)}{p _ {1} + N (s)} - \frac {c _ {2} w _ {2} T (s)}{p _ {2} + T (s) + \beta N (s)} - d - q _ {2} E \bigr ] d s \bigr \}, $$ which shows that all the solutions of system (2.3) are always positive for all $t > 0$. 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 $$ W (t) = c _ {1} N (t) + c _ {2} T (t) + Z (t). \tag {3.2} $$ Then, by differentiating (3.2) concerning $t$, we obtain $$ \begin{array}{l} \frac {d W}{d t} + \eta W = c _ {1} r _ {1} N (1 - \frac {N + \alpha_ {1} T}{k _ {1}}) + c _ {2} r _ {2} T (1 - \frac {T + \alpha_ {1} N}{k _ {2}}) - \frac {2 c _ {2} w _ {2} T Z}{p _ {2} + T + \beta N} - d Z - q _ {2} E Z - c _ {1} q _ {1} E N \\+ c _ {1} \eta N + c _ {2} \eta T + \eta Z, \\\leq c _ {1} r _ {1} N \left(1 - \frac {N}{k _ {1}}\right) + c _ {2} r _ {2} T \left(1 - \frac {T}{k _ {2}}\right) - d Z + c _ {1} \eta N + c _ {2} \eta T + \eta Z, \\= - \frac {c _ {1} r _ {1} N ^ {2}}{k _ {1}} + (r _ {1} + \eta) c _ {1} N - \frac {c _ {2} r _ {2} T ^ {2}}{k _ {2}} + (r _ {2} + \eta) c _ {2} T + (\eta - d) Z, \\\leq \frac {c _ {1} k _ {1} (r _ {1} + \eta) ^ {2}}{4 r _ {1}} + \frac {c _ {2} k _ {2} (r _ {2} + \eta) ^ {2}}{4 r _ {2}} + (\eta - d) Z, \\\leq \frac {c _ {1} r _ {2} k _ {1} (r _ {1} + \eta) ^ {2} + c _ {2} r _ {1} k _ {2} (r _ {2} + \eta) ^ {2}}{4 r _ {1} r _ {2}} + (\eta - d) Z, \\\end{array} $$ when $\eta - d < 0$, we can obtain $$ \frac {d W}{d t} + \eta W \leq \frac {c _ {1} r _ {2} k _ {1} (r _ {1} + \eta) ^ {2} + c _ {2} r _ {1} k _ {2} (r _ {2} + \eta) ^ {2}}{4 r _ {1} r _ {2}}, $$ noting $\kappa = \frac{c_1r_2k_1(r_1 + \eta)^2 + c_2r_1k_2(r_2 + \eta)^2}{4r_1r_2}$, therefore, applying a theorem on differential inequalities [34], we obtain $0\leq W\leq \frac{\kappa}{\eta} +\frac{W(N(0),T(0),Z(0))}{e^{\eta t}}$, let $t\to +\infty$, $W(N,T,Z)\leq \frac{\kappa}{\eta}$. So, all solutions of system (2.3) enter the region $$ D = \left\{(N, T, Z) \in \mathbb {R} _ {+} ^ {3}: 0 \leq W (N, T, Z) \leq \frac {\kappa}{\eta} \right\}. \tag {3.3} $$ 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} = (\widehat{N},\widehat{T},0)$, where $$ \widehat {N} = \frac {\alpha_ {1} k _ {2} - k _ {1}}{\alpha_ {1} \alpha_ {2} - 1} - \frac {q _ {1} k _ {1} E}{r _ {1}}, \quad \widehat {T} = \frac {\alpha_ {2} k _ {1} - k _ {2}}{\alpha_ {1} \alpha_ {2} - 1}; $$ (v)TPP-free equilibrium $E_4 = (\bar{N},0,\bar{Z})$, where $$ \bar {N} = \frac {\left(q _ {2} E + d\right) p _ {1}}{c _ {1} w _ {1} - d - q _ {2} E}, \quad \bar {Z} = \frac {r _ {1} \left(k _ {1} - \bar {N}\right) - q _ {1} k _ {1} E \left(p _ {1} + E\right)}{k _ {1} w _ {1}}; $$ (vi) the interior equilibrium, $E^{*} = (N^{*}, T^{*}, Z^{*})$, where $$ T^{*} = \frac{c_{1} w_{1} N^{*} - (d + q_{2} E)(p_{1} + N^{*})(p_{2} + \beta N^{*})}{(c_{2} w_{ 2} + d + q_{2} E)(p_{1} + N^{*}) - c_{1} w_{1} N^{*}}, \quad Z^{*} = \frac{(p_{1} + N^{*}) r_{1} (k_{1} - N^{*} - \alpha_{1} T^{*}) - q_{1} k_{1} E}{k_{1} w_{1}}; $$ and $N^{*}$ can be obtained from $$ r _ {2} \left(p _ {2} + T ^ {*} + \beta N ^ {*}\right) \left(k _ {2} - T ^ {*} - \alpha_ {2} N ^ {*}\right) - w _ {2} k _ {2} Z ^ {*} = 0. \tag {3.4} $$ 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 $\alpha_2 > \frac{k_2}{k_1}$ and $\alpha_1 > \frac{(\alpha_1\alpha_2 - 1)q_1k_1E}{r_1k_1} + \frac{k_1}{k_2}$. The TPP-free equilibrium $E_4$ exists, let $\bar{N}$ and $\bar{Z}$ both be positive, that is $w_1 > \frac{d + q_2E}{c_1}$ and $k_{1} > \frac{r_{1}N}{r_{1} - q_{1}E(p_{1} + E)}$. The interior equilibrium point $E^{*}$ exists; let $N^{*}$, $T^{*}$ and $Z^{*}$ all be positive, that is $k_{1} > \frac{q_{1}k_{1}E}{r_{1}} + N^{*} + \alpha_{1}T^{*}$, $c_{2}w_{2}(p_{1} + N^{*}) > c_{1}w_{1}N^{*} - (d + q_{2}E)(p_{1} + N^{*}) > 0$ 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_2E$. 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_1w_1k_1}{p_1 + k_1} -d - q_2E$. When $c_{1}\omega_{1} - d - q_{2}E\leq 0$, and $\alpha_{2} > \frac{k_{2}}{k_{1}}$ hold, $E_{1}$ is LAS(locally asymptotically stable). On the contrary, if $c_{1}\omega_{1} - d - q_{2}E > 0$, $\alpha_{2} > \frac{k_{2}}{k_{1}}$ and $k_{1} < \frac{p_{1}(d + q_{2}E)}{c_{1}w_{1} - d - q_{2}E}$ 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\big(1 - \frac{k_2\alpha_1}{k_1}\big) - q_1E, -r_2$ and $-\frac{c_2w_2k_2}{p_2 + k_2} -d - q_2E$, Therefore, $E_{2}$ is LAS if $\mathrm{k}_1 < \frac{r_2\alpha_1k_2}{r_2 - q_1E}$. - (iv) The eigenvalues of the zooplankton-free equilibrium $E_{3} = (\hat{N},\hat{T},0)$ are $\frac{c_1w_1N}{p_1 + \hat{N}} -\frac{c_2w_2T}{p_2 + \hat{T} + \beta\hat{N}} -d - q_2E$, $\lambda_{1}$ and $\lambda_{2}$, where $\lambda_{1}$ and $\lambda_{2}$ are the roots of the equation $$ \lambda^ {2} + \bar {b} _ {1} \lambda + \bar {c} _ {1} = 0, \tag {1} $$ $$ \begin{array}{l} \bar {b} _ {1} = - [ r _ {2} - r _ {1} + \frac {r _ {1} k _ {2} (2 \hat {N} + \alpha_ {1} \hat {T}) - r _ {2} k _ {1} (2 \hat {T} + \alpha_ {2} \hat {N})}{k _ {1} k _ {2}} ], \\\bar {c} _ {1} = r _ {1} r _ {2} [ 1 - (2 \hat {T} + \alpha_ {2} \hat {N}) (2 \hat {N} + \alpha_ {1} \hat {T}) ] [ \frac {1}{(2 \hat {N} + \alpha_ {1} \hat {T}) k _ {2}} + \frac {1}{(2 \hat {T} + \alpha_ {2} \hat {N}) k _ {1}} - \frac {1}{k _ {1} k _ {2}} ] \\+ q _ {1} r _ {2} E \left(\frac {k _ {1} (2 \hat {T} + \alpha_ {2} \hat {N}) - r _ {1} \alpha_ {1} ^ {2} \hat {N} \hat {T}}{k _ {1} k _ {2}} - 1\right). \\\end{array} $$ Therefore, let $\frac{c_1w_1\hat{N}}{p_1 + \hat{N}} - \frac{c_2w_2\hat{T}}{p_2 + \hat{T} + \beta\hat{N}} - d - q_2E < 0$, $\lambda_1$ and $\lambda_2$ with negative real parts, that is $\frac{c_1w_1\hat{N}}{p_1 + \hat{N}} - d - q_2E < \frac{c_2w_2\hat{T}}{p_2 + \hat{T} + \beta\hat{N}}$, $\bar{b}_1 > 0$ and $\bar{c}_1 > 0$. If the above conditions are satisfied, $E_3$ is LAS. (v) The eigenvalues of the TPP-free equilibrium $E_4 = (\bar{N}, 0, \bar{Z})$ are $r_2(1 - \frac{\alpha_2 \bar{N}}{k_2}) - \frac{w_2 \bar{Z}}{p_2 + \beta N}$, $\bar{\lambda}_1$ and $\bar{\lambda}_2$, where $\bar{\lambda}_1$ and $\bar{\lambda}_2$ are the roots of the equation $$ \lambda^ {2} - \left(\tilde {a} _ {2} + \tilde {b} _ {2}\right) \lambda + \tilde {a} _ {2} \tilde {b} _ {2} + \tilde {c} _ {2} = 0, \tag {2} $$ where $$ \tilde {a} _ {2} = (r _ {1} (1 - \frac {2 \bar {N}}{k _ {1}}) - \frac {w _ {1} p _ {1} \bar {Z}}{(p _ {1} + \bar {N}) ^ {2}} - q _ {1} E), $$ $$ \tilde {b} _ {2} = (\frac {c _ {1} w _ {1} \bar {N}}{p _ {1} + \bar {N}} - d - q _ {2} E), \quad \tilde {c} _ {2} = \frac {c _ {1} w _ {1} ^ {2} p _ {1} \bar {N} \bar {Z}}{(p _ {1} + \bar {N}) ^ {3}}. $$ Therefore, let $r_2(1 - \frac{\alpha_2\bar{N}}{k_2}) - \frac{w_2\bar{Z}}{p_2 + \beta\bar{N}} < 0$, $\bar{\lambda}_1$ and $\bar{\lambda}_2$ with negative real parts, that is $(\tilde{a}_2 + \tilde{b}_2) < 0$ and $\tilde{a}_2\tilde{b}_2 + \tilde{c}_2 > 0$. If the above conditions are satisfied, $E_4$ is LAS. (vi) By solving the Jacobi determinant of $E^{*}$, we can get its characteristic equation as follows $$ \lambda^ {3} + D _ {1} \lambda^ {2} + D _ {2} \lambda + D _ {3} = 0. \tag {3} $$ The interior equilibrium $E^{*} = (N^{*}, T^{*}, Z^{*})$ is LAS if - (a) $D_{1} > 0$ - (b) $D_{3} > 0$ - (c) $D_{1}D_{2} - D_{3} > 0$ $$ \begin{array}{l} D _ {1} = - \{r _ {2} [ 1 - \frac {(2 T ^ {*} + \alpha_ {2} N ^ {*})}{k _ {1}} ] - \frac {w _ {2} Z ^ {*} (p _ {2} + \beta N ^ {*})}{(p _ {2} + T ^ {*} + \beta N ^ {*}) ^ {2}} + r _ {1} [ 1 - \frac {(2 N ^ {*} + \alpha_ {1} T ^ {*})}{k _ {1}} ] - \frac {w _ {2} p _ {1} Z ^ {*}}{(p _ {1} + N ^ {*}) ^ {2}} - q _ {1} E \} \\- \left(\frac {c _ {1} w _ {1} N ^ {*}}{p _ {1} + N ^ {*}} - \frac {c _ {2} w _ {2} T ^ {*}}{p _ {2} + T ^ {*} + \beta N ^ {*}} - d - q _ {2} E\right), \\\end{array} $$ $$ D _ {2} = \{\frac {c _ {1} w _ {1} ^ {2} p _ {1} N ^ {*} Z ^ {*}}{(p _ {1} + N ^ {*}) ^ {3}} + \frac {c _ {2} w _ {1} w _ {2} \beta N ^ {*} T ^ {*} Z ^ {*}}{(p _ {2} + T ^ {*} + \beta N ^ {*}) ^ {2} (p _ {1} + N ^ {*})} - \frac {c _ {2} w _ {2} ^ {2} T ^ {*} Z ^ {*} (p _ {2} + \beta N ^ {*})}{(p _ {2} + T ^ {*} + \beta N ^ {*}) ^ {3}} \} $$ $$ \begin{array}{l} + \{r _ {1} [ 1 - \frac {(2 N ^ {*} + \alpha_ {1} T ^ {*})}{k _ {1}} ] - \frac {w _ {2} p _ {1} Z ^ {*}}{(p _ {1} + N ^ {*}) ^ {2}} - q _ {1} E \} \times \{r _ {2} [ 1 - \frac {(2 T ^ {*} + \alpha_ {2} N ^ {*})}{k _ {1}} ] - \frac {w _ {2} Z ^ {*} (p _ {2} + \beta N ^ {*})}{(p _ {2} + T ^ {*} + \beta N ^ {*}) ^ {2}} \} \\+ \frac {r _ {1} \alpha_ {1} N ^ {*}}{k _ {1}} (\frac {r _ {1} \alpha_ {1} T ^ {*}}{k _ {2}} + \frac {w _ {2} \beta T ^ {*} Z ^ {*})}{(p _ {2} + T ^ {*} + \beta N ^ {*}) ^ {2}}) + \{r _ {2} [ 1 - \frac {(2 T ^ {*} + \alpha_ {2} N ^ {*})}{k _ {1}} ] - \frac {w _ {2} Z ^ {*} (p _ {2} + \beta N ^ {*})}{(p _ {2} + T ^ {*} + \beta N ^ {*}) ^ {2}} \\+ r _ {1} [ 1 - \frac {(2 N ^ {*} + \alpha_ {1} T ^ {*})}{k _ {1}} ] - \frac {w _ {2} p _ {1} Z ^ {*}}{(p _ {1} + N ^ {*}) ^ {2}} - q _ {1} E \} \times \{\frac {c _ {1} w _ {1} N ^ {*}}{p _ {1} + N ^ {*}} - \frac {c _ {2} w _ {2} T ^ {*}}{p _ {2} + T ^ {*} + \beta N ^ {*}} - d - q _ {2} E \}, \\\end{array} $$ $$ \begin{array}{l} D _ {3} = - \{\frac {c _ {1} w _ {1} p _ {1} Z ^ {*}}{(p _ {1} + N ^ {*}) ^ {2}} - \frac {c _ {2} w _ {2} \beta T ^ {*} Z ^ {*}}{(p _ {2} + T ^ {*} + \beta N ^ {*}) ^ {2}} \} \times \{- \frac {r _ {1} \alpha_ {1} w _ {2} T ^ {*}}{k _ {1} (p _ {2} + T ^ {*} + \beta N ^ {*})} + \frac {w _ {1} N ^ {*}}{p _ {1} + N ^ {*}} \times (r _ {2} (1 - \frac {(2 T ^ {*} + \alpha_ {2} N ^ {*})}{k _ {2}}) \\- \frac {w _ {2} Z ^ {*} (p _ {2} + \beta N ^ {*})}{(p _ {2} + T ^ {*} + \beta N ^ {*}) ^ {2}}) - (\frac {c _ {2} w _ {2} Z ^ {*} (p _ {2} + \beta N ^ {*})}{(p _ {2} + T ^ {*} + \beta N ^ {*}) ^ {2}}) \times (\frac {w _ {2} T ^ {*}}{p _ {2} + T ^ {*} + \beta N ^ {*}}) \times [ - r _ {1} (1 - \frac {(2 N ^ {*} + \alpha_ {1} T ^ {*})}{k _ {1}}) + \frac {w _ {1} p _ {1} Z ^ {*}}{(p _ {1} + N ^ {*}) ^ {2}} + q _ {1} E ] \} \\+ \frac {w _ {1} N ^ {*}}{p _ {1} + N ^ {*}} \times (\frac {r _ {1} \alpha_ {1} T ^ {*}}{k _ {2}} + \frac {w _ {2} \beta T ^ {*} Z ^ {*}}{(p _ {2} + T ^ {*} + \beta N ^ {*}) ^ {2}}) - (\frac {c _ {2} w _ {2} Z ^ {*} (p _ {2} + \beta N ^ {*})}{(p _ {2} + T ^ {*} + \beta N ^ {*}) ^ {2}}) \times \{- \frac {r _ {1} w _ {2} T ^ {*}}{p _ {2} + T ^ {*} + \beta N ^ {*}} + \frac {r _ {1} w _ {2} (2 N ^ {*} + \alpha_ {1} T ^ {*}) T ^ {*}}{k _ {1} (p _ {2} + T ^ {*} + \beta N ^ {*})} \\+ \frac {w _ {1} w _ {2} p _ {1} T ^ {*} Z ^ {*}}{(p _ {2} + T ^ {*} + \beta N ^ {*}) (p _ {1} + N ^ {*}) ^ {2}} + \frac {w _ {2} q _ {1} E T ^ {*}}{p _ {2} + T ^ {*} + \beta N ^ {*}} + \frac {r _ {1} \alpha_ {1} w _ {1} N ^ {*} T ^ {*}}{k _ {2} (p _ {1} + N ^ {*})} + \frac {w _ {1} w _ {2} \beta N ^ {*} T ^ {*} Z ^ {*}}{(p _ {2} + T ^ {*} + \beta N ^ {*}) ^ {2} (p _ {1} + N ^ {*})} \bigg \} \\+ \left\{r _ {1} \left(1 - \frac {\left(2 N ^ {*} + \alpha_ {1} T ^ {*}\right)}{k _ {1}}\right) - \frac {w _ {2} p _ {1} Z ^ {*}}{\left(p _ {1} + N ^ {*}\right) ^ {2}} - q _ {1} E \right\} \times \left\{r _ {2} \left(1 - \frac {\left(2 T ^ {*} + \alpha_ {2} N ^ {*}\right)}{k _ {1}}\right) - \frac {w _ {2} Z ^ {*} \left(p _ {2} + \beta N ^ {*}\right)}{\left(p _ {2} + T ^ {*} + \beta N ^ {*}\right) ^ {2}} \right\} \\+ \frac {r _ {1} \alpha_ {1} N ^ {*}}{k _ {1}} \times (\frac {r _ {1} \alpha_ {1} T ^ {*}}{k _ {2}} + \frac {w _ {2} \beta T ^ {*} Z ^ {*}}{(p _ {2} + T ^ {*} + \beta N ^ {*}) ^ {2}}). \\\end{array} $$ 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 $\frac{k_2}{\alpha_2} < k_1 < \frac{p_1(d + q_2E)}{c_1w_1 - d - q_2E}$, 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_1k_2}{r_2 - q_1E}$, 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 $\frac{r_2\alpha_1k_2}{r_2 - q_1E} < k_1 < \frac{k_2}{\alpha_2}$ (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 $\frac{c_1w_1\hat{N}}{p_1 + \hat{N}} -d - q_2E < \frac{c_2w_2\hat{T}}{p_2 + \hat{T} + \beta\hat{N}}$, 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 $\left(\frac{c_1w_1\hat{N}}{p_1 + \hat{N}} -(d + q_2E) > 0\right)$, we can find a threshold value of the avoidance of TPP by zooplankton $(\beta < \frac{(c_2w_2\hat{T})(p_1 + \hat{N})}{(\hat{N})(c_1w_1\hat{N} - (d + q_2E)(p_1 + \hat{N}))} -\frac{p_2 + \hat{T}}{\hat{N}})$, below which the zooplankton population will become extinct. On the contrary, for $\frac{c_1w_1\hat{N}}{p_1 + \hat{N}} -(d + q_2E) < 0$, 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 $\left(\frac{(d + q_2E)p_1}{c_1w_1 - d - q_2E} < k_1 < \frac{(d + q_2E)(p_1) + c_1w_1p_1}{c_1w_1 - d - q_2E}\right)$, then TPP becomes extinct under the condition $\left(\frac{r_2(k_2 - \alpha_2\bar{N})}{k_2} < \frac{w_2\bar{Z}}{p_2 + \beta\bar{N}}\right)$, 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 $(k_{2} - \alpha_{2}\bar{N} >0)$, we can find a threshold value of the avoidance of TPP by zooplankton $(\beta < \frac{k_{2}w_{2}\bar{Z}}{\bar{N}r_{2}(k_{2} - \alpha_{2}\bar{N})} -\frac{p_{2}}{N})$, below which TPP may become extinct. On the contrary, for $k_{2} - \alpha_{2}\bar{N} < 0$, 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 $k_{2} - \alpha_{2}\bar{N} < 0$, TPP extinction dose not depend on the avoidance of TPP by zooplankton $(\beta)$ and harvest term on zooplankton $(q_{2}E)$. - (vi) $E^{*} = (N^{*}, T^{*}, Z^{*})$: 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><tr><td colspan="3">Table 1: Existence and stability conditions of the equilibrium points.</td></tr><tr><td>Equilibrium</td><td>Existence conditions</td><td>Stability conditions</td></tr><tr><td>E0=(0,0,0)</td><td>Always exist</td><td>Always unstable</td></tr><tr><td>E1=(k1,0,0)</td><td>Always exist</td><td>(i)c1w1-d-q2E&gt;0,α2&gt;k2/k1,k1&lt; p1(d+q2E)/c1w1-d-q2E,or (ii)c1w1-d-q2E≤0,α2&gt;k2/k1</td></tr><tr><td>E2=(0,k2,0)</td><td>Always exist</td><td>(i) k1&lt;r2α1k2/r2-q1E</td></tr><tr><td>E3=(ˆN,ˆT,0)</td><td>(i) α2&gt;k2/k1,(ii) α1=(α1α2-1)q1k1E+k1/k2</td><td>(i) c1w1ˆN-d-q2E&lt;c2w2ˆT/p2+ˆT+βˆN,(ii) b1&gt;0,ˆc1&gt;0</td></tr><tr><td>E4=(ˆN,0,ˆZ)</td><td>(i) w1&gt;d+q2E/c1,(ii) k1=r1N/r1-ˆq1E(p1+E)</td><td>(i) r2(1-α2ˆN)/p2+βˆN,(ii) a2+ˆb2&lt;0, a2ˆb2+ˆc2&gt;0</td></tr><tr><td>E*=(N*,T*,Z*)</td><td>(i) k1=q1k1E/r1+N*+α1T*,(ii) c2w2(p1+N*)&gt;c1w1N*(d+q2E)(p1+N*)&gt;0,(iii) positive root of Eq.(3.4) exists</td><td>(i) D1&gt;0,(ii) D3&gt;0,(iii) D1D2-D3&gt;0</td></tr></table> ### c) Bistability analysis of equilibrium point The existence and stability of these equilibrium points are summarized in Table 1 and Fig 1. When $c_{1}w_{1} - d - q_{2}E > 0$, equilibria $E_{2} = (0,k_{2},0)$, $E_{3} = (\widehat{N},\widehat{T},0)$, $E_{1} = (k_{1},0,0)$ and $E_{4} = (\bar{N},0,\bar{Z})$ keep stable for $(0 < k_{1} < \frac{r_{2}\alpha_{1}k_{2}}{r_{2} - q_{1}E})$, $(\frac{r_{2}\alpha_{1}k_{2}}{r_{2} - q_{1}E} < k_{1} < \frac{k_{2}}{\alpha_{2}})$, $(\frac{k_{2}}{\alpha_{2}} < k_{1} < \frac{p_{1}(d + q_{2}E)}{c_{1}w_{1} - d - q_{2}E})$ and $(\frac{(d + q_2E)p_1}{c_1w_1 - d - q_2E} < k_1 < \frac{(d + q_2E)(p_1) + c_1w_1p_1}{c_1w_1 - d - q_2E})$, respectively(Fig.1(a)). Obviously, for $k_{1}$ at the different equilibria above, the coexistence of NTP, TPP, and zooplankton requires the three ranges $(k_{1} > \frac{r_{2}\alpha_{1}k_{2}}{r_{2} - q_{1}E})$, $(k_{1} < \frac{k_{2}}{\alpha_{2}})$, and $(k_{1} > \frac{(d + q_{2}E)p_{1}}{c_{1}w_{1} - d - q_{2}E})$, 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 $\frac{k_2}{\alpha_2} < k_1 < \min \left\{\frac{r_2\alpha_1k_2}{r_2 - q_1E},\frac{(d + q_2E)p_1}{c_1w_1 - d - q_2E}\right\}$, we can observe the bistability of $E_{1}$ and $E_{2}$ (Fig.1(b)(c)). If conditions $\frac{(d + q_2E)p_1}{c_1w_1 - d - q_2E} < k_1 < \min \left\{\frac{r_2\alpha_1k_2}{r_2 - q_1E},\frac{(d + q_2E)p_1 + c_1w_1p_1}{c_1w_1 - d - q_2E}\right\}$ and $\left(\frac{r_2(k_2 - \alpha_2\bar{N})}{k_2} < \frac{w_2\bar{Z}}{p_2 + \beta N}\right)$ hold simultaneous, we can find the bistability of $E_{2}$ and $E_{4}$ (Fig.1(d)(e)). On the contrary, if $\frac{(d + q_2E)p_1}{c_1w_1 - d - q_2E} < k_1 < \frac{r_2\alpha_1k_2}{r_2 - q_1E}$ holds, for either $k_{1} > \frac{(d + q_{2}E)(p_{1}) + c_{1}w_{1}p_{1}}{c_{1}w_{1} - d - q_{2}E}$ or $\frac{r_2(k_2 - \alpha_2\bar{N})}{k_2} >\frac{w_2\bar{Z}}{p_2 + \beta\bar{N}}$ we'll get the existence of stable $E_{2}$ together with unstable $E_{4}$. Identically, for $\max \left\{\frac{r_2\alpha_1k_2}{r_2 - q_1E},\frac{(d + q_2E)p_1}{c_1w_1 - d - q_2E}\right\} < k_1 < \min \left\{\frac{k_2}{\alpha_2},\frac{(d + q_2E)p_1 + c_1w_1p_1}{c_1w_1 - d - q_2E}\right\}$ together with $\alpha_{1}\alpha_{2} < 1$, $\frac{c_{1}w_{1}\hat{N}}{p_{1} + \hat{N}} -d - q_{2}E < \frac{c_{2}w_{2}\hat{T}}{p_{2} + \hat{T} + \beta\hat{N}}$ and $\frac{r_2(k_2 - \alpha_2\bar{N})}{k_2} < \frac{w_2\bar{Z}}{p_2 + \beta N}$, 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. $(\frac{r_{2}(k_{2} - \alpha_{2}\hat{N})}{k_{2}} > \frac{w_{2}\bar{Z}}{p_{2} + \beta\hat{N}}, \frac{c_{1}w_{1}\hat{N}}{p_{1} + \hat{N}} - d - q_{2}E < \frac{c_{2}w_{2}\hat{T}}{p_{2} + \hat{T} + \beta\hat{N}}$ and $\frac{r_{2}(k_{2} - \alpha_{2}\hat{N})}{k_{2}} < \frac{w_{2}\bar{Z}}{p_{2} + \beta\hat{N}}$. 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 $$ \frac{\mathrm{d}U(t)}{\mathrm{d}t} = F(U(t),U(t-\tau_1),U(t-\tau_2)), $$ $$ U (t) = [ N (t), T (t), Z (t) ], \quad U (t - \tau_ {1}) = [ N (t - \tau_ {1}), T (t - \tau_ {1}), Z (t - \tau_ {1}) ], $$ $$ U (t - \tau_ {2}) = [ N (t - \tau_ {2}), T (t - \tau_ {2}), Z (t - \tau_ {2}) ]. $$ Next, assuming $\Lambda_1(t) = N(t) - N^*$, $\Lambda_2(t) = T(t) - T^*$, $\Lambda_3(t) = Z(t) - Z^*$ at the positive equilibrium point, and linearizing the system (2.2), we can obtain $$ \frac {\mathrm {d}}{\mathrm {d} t} \left( \begin{array}{l} \Lambda_ {1} (t) \\\Lambda_ {2} (t) \\\Lambda_ {3} (t) \end{array} \right) = L \left( \begin{array}{l} N (t) \\T (t) \\Z (t) \end{array} \right) + M \left( \begin{array}{l} N (t - \tau_ {1}) \\T (t - \tau_ {1}) \\Z (t - \tau_ {1}) \end{array} \right) + S \left( \begin{array}{l} N (t - \tau_ {2}) \\T (t - \tau_ {2}) \\Z (t - \tau_ {2}) \end{array} \right), \tag {4.2} $$ where $$ L = \left(\frac {\partial F}{\partial U (t)}\right) _ {E ^ {*}}, \quad M = \left(\frac {\partial F}{\partial U (t - \tau_ {1})}\right) _ {E ^ {*}}, \quad S = \left(\frac {\partial F}{\partial U (t - \tau_ {2})}\right) _ {E ^ {*}}. $$ We linearize the system(2.2) about positive equilibrium $E^{*} = (N^{*},T^{*},Z^{*})$, and get $$ \frac {\mathrm {d} U (t)}{\mathrm {d} t} = L U (t) + M U \left(t - \tau_ {1}\right) + S U \left(t - \tau_ {2}\right), \tag {4.3} $$ 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}$. $$ L = \left( \begin{array}{c c c} l _ {1 1} & l _ {1 2} & l _ {1 3} \\l _ {2 1} & l _ {2 2} & l _ {2 3} \\0 & 0 & l _ {3 3} \end{array} \right), \quad M = \left( \begin{array}{c c c} 0 & 0 & 0 \\0 & 0 & 0 \\m _ {3 1} & 0 & m _ {3 3} \end{array} \right), \quad S = \left( \begin{array}{c c c} 0 & 0 & 0 \\0 & 0 & 0 \\s _ {3 1} & s _ {3 2} & s _ {3 3} \end{array} \right), \quad U = \left( \begin{array}{c} N _ {1} (\cdot) \\T _ {1} (\cdot) \\Z _ {1} (\cdot) \end{array} \right), $$ where $N_{1}, T_{1}, Z_{1}$ are small perturbations around the equilibrium point $E^{*} = (N^{*},T^{*},Z^{*})$. We have $$ l _ {1 1} = \frac {- r N}{k _ {1}} + \frac {w _ {1} Z N}{(p _ {1} + N) ^ {2}} - q _ {1} E, l _ {1 2} = \frac {r _ {1} \alpha_ {1} N}{k _ {1}}, l _ {1 3} = - \frac {w _ {1} N}{p _ {1} + N}, $$ $$ l _ {2 1} = \frac {r _ {2} \alpha_ {2} T}{k _ {1}} + \frac {w _ {2} \beta T Z}{(p _ {2} + T + \beta N) ^ {2}}, l _ {2 2} = r _ {2} - \frac {(2 r _ {2} T + r _ {2} \alpha_ {1} N)}{k _ {2}}, $$ $$ l _ {2 3} = - \frac {w _ {2} T}{(p _ {2} + T + \beta N)}, l _ {3 3} = - d - q _ {2} E, m _ {3 1} = \frac {c _ {1} w _ {1} p _ {1} Z}{(p _ {1} + N) ^ {2}}, m _ {3 3} = \frac {c _ {1} w _ {1} N}{(p _ {1} + N)}, $$ $$ s _ {3 1} = \frac {c _ {2} w _ {2} \beta T Z}{\left(p _ {2} + T + \beta N\right) ^ {2}}, s _ {3 2} = \frac {c _ {2} w _ {2} Z \left(p _ {2} + \beta N\right)}{\left(p _ {2} + T + \beta N\right) ^ {2}}, s _ {3 3} = \frac {c _ {2} w _ {2} T}{\left(p _ {2} + T + \beta N\right)}. $$ The characteristic equation for the linearized system (2.2) is obtained as $$ D (\xi , \tau_ {1}, \tau_ {2}) \equiv P (\xi) + Q (\xi) \mathrm {e} ^ {- \xi \tau_ {1}} + R (\xi) \mathrm {e} ^ {- \xi \tau_ {2}} = 0, \tag {4.4} $$ where $$ P(\xi) = \xi^{3} + A_{2}\xi^{2} + A_{1}\xi + A_{0},\quad Q(\xi) = B_{2}\xi^{2} + B_{1}\xi + B_{0},\quad R(\xi) = C_{2}\xi^{2} + C_{1}\xi + C_{0}, $$ with $$ A _ {2} = - \left(l _ {3 3} + l _ {2 2} - l _ {1 1}\right), \quad A _ {1} = l _ {1 1} l _ {2 2} + l _ {1 1} l _ {3 3} + l _ {2 2} l _ {3 3} - l _ {1 2} l _ {2 1}, \quad A _ {0} = - l _ {1 1} l _ {2 2} l _ {3 3} + l _ {1 2} l _ {2 1} l _ {3 3} $$ $$ B _ {2} = - m _ {3 3}, \quad B _ {1} = - l _ {1 1} m _ {3 3} - l _ {2 2} m _ {3 3} - l _ {1 3} m _ {3 1}, \quad B _ {0} = + l _ {1 3} l _ {2 2} m _ {3 1} + l _ {1 1} l _ {2 2} m _ {3 3} + l _ {1 2} l _ {2 1} m _ {3 3} - l _ {1 2} l _ {2 3} m _ {3 1}, $$ $$ C _ {2} = - s _ {3 3}, \quad C _ {1} = - l _ {1 3} s _ {3 1} + l _ {1 1} s _ {3 3} - l _ {2 2} l _ {2 3} s _ {3 2} - l _ {2 2} s _ {3 3}, $$ $$ C _ {0} = l _ {1 1} s _ {3 3} + l _ {1 1} l _ {2 3} s _ {3 2} + l _ {1 2} l _ {2 1} s _ {3 3} + l _ {1 3} l _ {2 2} s _ {3 1} - l _ {1 2} l _ {2 3} s _ {3 1} - l _ {1 3} l _ {2 1} s _ {3 2}. $$ 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): $\tau_{1} = 0, \tau_{2} > 0$. In this case, the characteristic equation(4.4) becomes $$ \begin{array}{l} D (\xi , \tau_ {2}) \equiv P (\xi) + Q (\xi) + R (\xi) \mathrm {e} ^ {- \xi \tau_ {2}} \\\equiv \xi^ {3} + A _ {2} \xi^ {2} + A _ {1} \xi + A _ {0} + B _ {2} \xi^ {2} + B _ {1} \xi + B _ {0} + \left(C _ {2} \xi^ {2} + C _ {1} \xi + C _ {0}\right) \mathrm {e} ^ {- \xi \tau_ {2}} = 0, \tag {4.5} \\\end{array} $$ putting $\xi = i\omega (\omega >0)$ in Eq.(4.5), and separating the real and imaginary parts, we have $$ \begin{array}{l} - \left(A _ {2} + B _ {2}\right) \omega^ {2} + \left(A _ {0} + B _ {0}\right) = \left(C _ {2} \omega^ {2} - C _ {0}\right) \cos (\omega \tau_ {2}) - C _ {1} \omega \sin (\omega \tau_ {2}), \\- \omega^ {3} + \left(A _ {1} + B _ {1}\right) \omega = \left(C _ {0} - C _ {2} \omega^ {2}\right) \sin \left(\omega \tau_ {2}\right) - C _ {1} \omega \cos \left(\omega \tau_ {2}\right). \tag {4.6} \\\end{array} $$ Squaring and adding the equation(4.6), we obtain $$ \left[ - \left(A _ {2} + B _ {2}\right) \omega^ {2} + \left(A _ {0} + B _ {0}\right) \right] ^ {2} + \left[ - \omega^ {3} + \left(A _ {1} + B _ {1}\right) \omega \right] ^ {2} = \left(C _ {2} \omega^ {2} - C _ {0}\right) ^ {2} + \left(C _ {1} \omega\right) ^ {2}. \tag {4.7} $$ Simplifying Eq.(4.7) and substituting $\omega^2 = \psi$, the above equation can be written as $$ \Psi (\psi) \equiv \psi^ {3} + a _ {2} \psi^ {2} + a _ {1} \psi + a _ {0} = 0, \tag {4.8} $$ where $$ a _ {2} = - (A _ {2} + B _ {2}) ^ {2} - 2 (A _ {1} + B _ {1}) - C _ {2} ^ {2}, a _ {1} = (A _ {1} + B _ {1}) ^ {2} - 2 (A _ {0} + B _ {0}) (A _ {2} + B _ {2}) - 2 C _ {0} C _ {2} - C _ {1} ^ {2}, a _ {0} = - C _ {0} ^ {2}. $$ $$ \begin{array}{l l} \text {(H 1)}: & a _ {2} > 0, a _ {0} > 0, a _ {2} a _ {1} - a _ {0} > 0. \end{array} $$ 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^{*}$ is asymptotically stable for all $\tau_{2} > 0$ when (H1) holds. $$ \begin{array}{l} \text {(H 2)}: a _ {2} < 0, a _ {1} < 0, a _ {0} < 0 \text {o r} a _ {2} > 0, a _ {1} < 0, a _ {0} < 0 \text {o r} a _ {2} > 0, a _ {1} > 0, a _ {0} < 0. \end{array} $$ If (H2) holds, Eq.(4.8) has exactly one positive root $\omega_0$, substituting $\omega_0$ in Eq.(4.6), we obtain $$ \begin{array}{l} - \left(A _ {2} + B _ {2}\right) \omega_ {0} ^ {2} + \left(A _ {0} + B _ {0}\right) = \left(C _ {2} \omega_ {0} ^ {2} - C _ {0}\right) \cos \left(\omega_ {0} \tau_ {2}\right) - C _ {1} \omega_ {0} \sin \left(\omega_ {0} \tau_ {2}\right), \\- \omega_ {0} ^ {3} + (A _ {1} + B _ {1}) \omega_ {0} = \left(C _ {0} - C _ {2} \omega_ {0} ^ {2}\right) \sin \left(\omega_ {0} \tau_ {2}\right) - C _ {1} \omega_ {0} \cos \left(\omega_ {0} \tau_ {2}\right). \tag {4.9} \\\end{array} $$ For the critical value of $\tau_{2}$, we can obtain $$ \tau_{2j} = \frac{1}{\omega_0} \arccos \left\{ \frac{[C_1 + C_2(A_2 + B_2)]\omega_0^4 + [C_1(A_1 + B_1) - C_0(A_2 + B_2) - C_2(A_0 + B_0)]\omega_0^2 + C_0(A_0 + B_0)}{-(C_0 - C_2\omega_0^2)^2 - (C_1\omega_0)^2} \right\} + \frac{2j\pi}{\omega_0} $$ $$ j = 0, 1, 2 \dots . $$ For the transversality condition, differentiating Eq.(4.5) with respect to $\tau_{2}$, we get $$ \frac {d \xi}{d \tau_ {2}} = \frac {\xi (C _ {2} \xi^ {2} + C _ {1} \xi + C _ {0}) e ^ {- \xi \tau_ {2}}}{3 \xi^ {2} + 2 A _ {2} \xi + A _ {1} + (2 B _ {2} \xi + B _ {1}) + (2 C _ {2} \xi + C _ {1}) e ^ {- \xi \tau_ {2}}}. $$ Solving $\left(\frac{d\xi}{d\tau_2}\right)^{-1}$, we obtain $$ (\frac {d \xi}{d \tau_ {2}}) ^ {- 1} = \frac {3 \xi^ {2} + 2 A _ {2} \xi + A _ {1} + (2 B _ {2} \xi + B _ {1}) + (2 C _ {2} \xi + C _ {1}) e ^ {- \xi \tau_ {2}}}{\xi (C _ {2} \xi^ {2} + C _ {1} \xi + C _ {0}) e ^ {- \xi \tau_ {2}}}. $$ Then at $\tau_{2} = \tau_{20}$ and $\xi = i\omega_0$, we can get $$ \begin{array}{l} [ \mathrm {R e} (\frac {d \xi}{d \tau_ {2}}) _ {\tau_ {2} = \tau_ {2 0}, \xi = i \omega_ {0}} ] ^ {- 1} = \mathrm {R e} [ \frac {3 (i \omega_ {0}) ^ {2} + (2 A _ {2} + B _ {2}) (i \omega_ {0}) + A _ {1} + B _ {1}}{(i \omega_ {0}) (C _ {2} (i \omega_ {0}) ^ {2} + C _ {1} (i \omega_ {0}) + C _ {0}) (\cos (\omega_ {0} \tau_ {2 0}) - i \sin (\omega_ {0} \tau_ {2 0}))} ] \\+ \mathrm {R e} [ \frac {2 C _ {2} (i \omega_ {0}) + C _ {1}}{(i \omega_ {0}) (C _ {2} (i \omega_ {0}) ^ {2} + C _ {1} (i \omega_ {0}) + C _ {0})} ]. \\\end{array} $$ $$ \left[ \mathrm{Re} \left( \frac{d \xi}{d \tau_2} \right) _{\tau_2 = \tau_{20}, \xi = i \omega_0} \right] ^{-1} = \mathrm{Re} \left[ \frac{M_R + M_I i}{N_R + N_I i} \right] + \mathrm{Re} \left[ \frac{Q_R + Q_I i}{P_R + P_I i} \right] = \frac{M_R N_R + M_I N_I}{N_R^2 + N_I^2} + \frac{Q_R P_R + Q_I P_I}{P_R^2 + P_I^2}, $$ $$ M _ {R} = - 3 \omega_ {0} ^ {2} + A _ {1} + B _ {1}, \quad M _ {I} = 2 (A _ {2} + B _ {2}) \omega_ {0}, \quad N _ {R} = \left(C _ {0} \omega_ {0} - C _ {2} \omega_ {0} ^ {3}\right) \sin (\omega_ {0} \tau_ {2 0}) - C _ {1} \omega_ {0} ^ {2} \cos (\omega_ {0} \tau_ {2 0}), $$ $$ N _ {I} = \left(C _ {0} \omega_ {0} - C _ {2} \omega_ {0} ^ {3}\right) \cos \left(\omega_ {0} \tau_ {2 0}\right) + C _ {1} \omega_ {0} ^ {2} \sin \left(\omega_ {0} \tau_ {2 0}\right), \quad Q _ {R} = C _ {1}, \quad Q _ {I} = 2 C _ {2} \omega_ {0}, $$ $$ P _ {R} = - C _ {1} \omega_ {0} ^ {2}, \quad P _ {I} = C _ {0} \omega_ {0} - C _ {2} \omega_ {0} ^ {3}. $$ Then $$ \left[ \operatorname {R e} \left(\frac {d \xi}{d \tau_ {2}}\right) _ {\tau_ {2} = \tau_ {2 0}, \xi = i \omega_ {0}} \right] ^ {- 1} = \frac {A}{B} + \frac {C}{D} = \frac {A D + B C}{B D}, \tag {4.11} $$ here $$ A = M _ {R} N _ {R} + M _ {I} N _ {I}, \quad B = N _ {R} ^ {2} + N _ {I} ^ {2}, $$ $$ C = Q _ {R} P _ {R} + Q _ {I} P _ {I}, \quad D = P _ {R} ^ {2} + P _ {I} ^ {2}. $$ From this, we can get $$ \mathrm {s g n} [ \mathrm {R e} (\frac {d \xi}{d \tau_ {2}}) _ {\tau_ {2} = \tau_ {2 0}, \xi = i \omega_ {0}} ] ^ {- 1} = \mathrm {s g n} [ A D + B C ]. $$ If (H3): $AD + BC \neq 0$ holds, the transversal condition $\mathrm{sgn}[\mathrm{Re}(\frac{d\xi}{d\tau_2})_{\tau_2 = \tau_{20},\xi = i\omega_0}]^{-1} \neq 0$. From the above analysis, the following theorem can be drawn Theorem 4.1. For $\tau_{1} = 0$ and $\tau_{2} > 0$, we have the following results: (i) If (H1) holds, then the equilibrium $E^{*}$ is asymptotically stable for all $\tau_{2} > 0$. (ii) If (H3) holds, and (H2) holds, then the equilibrium $E^{*}$ is locally asymptotically stable for all $\tau_{2} < \tau_{20}$ together with unstable for $\tau_{2} > \tau_{20}$ and undergoes Hopf bifurcation at $\tau_{2} = \tau_{20}$. Case (3): $\tau_{1} > 0, \tau_{2} = 0$. In this case, the characteristic equation(4.4) becomes as follows $$ \begin{array}{l} D (\xi , \tau_ {1}) \equiv P (\xi) + R (\xi) + Q (\xi) e ^ {- \xi \tau_ {1}} \\\equiv \xi^ {3} + A _ {2} \xi^ {2} + A _ {1} \xi + A _ {0} + \left(B _ {2} \xi^ {2} + \left(C _ {2} \xi^ {2} + C _ {1} \xi + C _ {0}\right) + B _ {1} \xi + B _ {0}\right) e ^ {- \xi \tau_ {1}} = 0. \tag {4.12} \\\end{array} $$ putting $\xi = i\omega (\omega >0)$ in Eq.(4.12), and separating the real and imaginary parts, we have $$ \begin{array}{l} - \left(A _ {2} + C _ {2}\right) \omega^ {2} + \left(A _ {0} + C _ {0}\right) = \left(B _ {2} \omega^ {2} - B _ {0}\right) \cos (\omega \tau_ {1}) - B _ {1} \omega \sin (\omega \tau_ {1}), \\- \omega^ {3} + \left(A _ {1} + C _ {1}\right) \omega = \left(B _ {0} - B _ {2} \omega^ {2}\right) \sin \left(\omega \tau_ {1}\right) - B _ {1} \omega \cos \left(\omega \tau_ {1}\right). \tag {4.13} \\\end{array} $$ Squaring and adding the equation(4.13), we obtain $$ [ - (A _ {2} + C _ {2}) \omega^ {2} + (A _ {0} + C _ {0}) ] ^ {2} + [ - \omega^ {3} + (A _ {1} + C _ {1}) \omega ] ^ {2} = (B _ {2} \omega^ {2} - B _ {0}) ^ {2} + (B _ {1} \omega) ^ {2}. \tag {4.14} $$ Based on the calculation method for case (2), we can simplify (4.14) to the following $$ \Psi (\quad) \equiv^ {3} + b _ {2} ^ {2} + b _ {1} + b _ {0} = 0, \tag {4.15} $$ $$ b _ {2} = - (A _ {2} + C _ {2}) ^ {2} - 2 (A _ {1} + C _ {1}) - B _ {2} ^ {2}, \quad b _ {1} = (A _ {1} + C _ {1}) ^ {2} - 2 (A _ {0} + C _ {0}) (A _ {2} + C _ {2}) - 2 B _ {0} B _ {2} - B _ {1} ^ {2}, \quad b _ {0} = - B _ {0} ^ {2}. $$ - (H4): $b_{2} > 0, b_{0} > 0, b_{2}b_{1} - b_{0} > 0$. - 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^{*}$ is asymptotically stable for all $\tau_{1} > 0$ when (H4) holds. - (H5): $b_{2} < 0, b_{1} < 0, b_{0} < 0$ or $b_{2} > 0, b_{1} < 0, b_{0} < 0$ or $b_{2} > 0, b_{1} > 0, b_{0} < 0$. - If (H5) holds, Eq.(4.15) has exactly one positive root $\hat{\omega}_0$, substituting $\hat{\omega}_0$ in Eq.(4.13), we obtain $$ \begin{array}{l} - \left(A _ {2} + C _ {2}\right) \hat {\omega_ {0}} ^ {2} + \left(A _ {0} + C _ {0}\right) = \left(B _ {2} \hat {\omega_ {0}} ^ {2} - B _ {0}\right) \cos \left(\hat {\omega_ {0}} \tau_ {1}\right) - B _ {1} \hat {\omega_ {0}} \sin \left(\hat {\omega_ {0}} \tau_ {1}\right), \\- \hat {\omega} _ {0} ^ {3} + \left(A _ {1} + C _ {1}\right) \hat {\omega} _ {0} = \left(B _ {0} - B _ {2} \hat {\omega} _ {0} ^ {2}\right) \sin \left(\hat {\omega} _ {0} \tau_ {1}\right) - B _ {1} \hat {\omega} _ {0} \cos \left(\hat {\omega} _ {0} \tau_ {1}\right). \tag {4.16} \\\end{array} $$ For the critical value of $\tau_{1}$, we can obtain $$ \tau_ {\mathrm {1 j}} = \frac {1}{\hat {\omega} _ {0}} \arccos \{\frac {[ B _ {1} + B _ {2} (A _ {2} + C _ {2}) ] \hat {\omega} _ {0} ^ {4} + [ B _ {1} (A _ {1} + C _ {1}) - C _ {0} (A _ {2} + C _ {2}) - B _ {2} (A _ {0} + C _ {0}) ] \hat {\omega} _ {0} ^ {2} + B _ {0} (A _ {0} + C _ {0})}{- (B _ {0} - B _ {2} \hat {\omega} _ {0} ^ {2}) ^ {2} - (B _ {1} \hat {\omega} _ {0}) ^ {2}} \} + \frac {2 j \pi}{\hat {\omega} _ {0}}, $$ $$ j = 0, 1, 2 \dots . \tag {4.17} $$ For the transversality condition, differentiating Eq.(4.13) with respect to $\tau_{1}$, we get $$ \frac {d \xi}{d \tau_ {1}} = \frac {\xi (B _ {2} \xi^ {2} + B _ {1} \xi + B _ {0}) e ^ {- \xi \tau_ {1}}}{3 \xi^ {2} + 2 A _ {2} \xi + A _ {1} + (2 C _ {2} \xi + C _ {1}) + (2 B _ {2} \xi + B _ {1}) e ^ {- \xi \tau_ {1}}}. $$ Solving $\left(\frac{d\xi}{d\tau_1}\right)^{-1}$, we obtain $$ \left(\frac {d \xi}{d \tau_ {1}}\right) ^ {- 1} = \frac {3 \xi^ {2} + 2 A _ {2} \xi + A _ {1} + (2 C _ {2} \xi + C _ {1}) + (2 B _ {2} \xi + B _ {1}) e ^ {- \xi \tau_ {1}}}{\xi \left(B _ {2} \xi^ {2} + B _ {1} \xi + B _ {0}\right) e ^ {- \xi \tau_ {1}}}. $$ Then at $\tau_{1} = \tau_{10}$ and $\xi = i\omega_0$, we can get $$ \begin{array}{l} [ \mathrm {R e} (\frac {d \xi}{d \tau_ {1}}) _ {\tau_ {1} = \tau_ {1 0}, \xi = i \hat {\omega_ {0}}} ] ^ {- 1} = \mathrm {R e} [ \frac {3 (i \hat {\omega_ {0}}) ^ {2} + (2 A _ {2} + C _ {2}) (i \hat {\omega_ {0}}) + A _ {1} + C _ {1}}{(i \hat {\omega_ {0}}) (B _ {2} (i \hat {\omega_ {0}}) ^ {2} + B _ {1} (i \hat {\omega_ {0}}) + B _ {0}) (\cos (\hat {\omega_ {0}} \tau_ {1 0}) - i \sin (\hat {\omega_ {0}} \tau_ {1 0}))} ] \\+ \operatorname {R e} \left[ \frac {2 B _ {2} (i \hat {\omega_ {0}}) + B _ {1}}{(i \hat {\omega_ {0}}) (B _ {2} (i \hat {\omega_ {0}}) ^ {2} + B _ {1} (i \hat {\omega_ {0}}) + B _ {0})} \right]. \\\end{array} $$ Now $$ [ \mathrm {R e} (\frac {d \xi}{d \tau_ {1}}) _ {\tau_ {1} = \tau_ {1 0}, \xi = i \omega_ {0}} ] ^ {- 1} = \mathrm {R e} [ \frac {\widehat {M _ {R}} + \widehat {M _ {I}} i}{\widehat {N _ {R}} + \widehat {N _ {I}} i} ] + \mathrm {R e} [ \frac {\widehat {Q _ {R}} + \widehat {Q _ {I}} i}{\widehat {P _ {R}} + \widehat {P _ {I}} i} ] = \frac {\widehat {M _ {R}} \widehat {N _ {R}} + \widehat {M _ {I}} \widehat {N _ {I}}}{\widehat {N _ {R}} ^ {2} + \widehat {N _ {I}} ^ {2}} + \frac {\widehat {Q _ {R}} \widehat {P _ {R}} + \widehat {Q _ {I}} \widehat {P _ {I}}}{\widehat {P _ {R}} ^ {2} + \widehat {P _ {I}} ^ {2}}, $$ $$ \begin{array}{l} \widehat {M _ {R}} = - 3 \hat {\omega_ {0}} ^ {2} + A _ {1} + C _ {1}, \widehat {M _ {I}} = 2 (A _ {2} + C _ {2}) \hat {\omega_ {0}}, \widehat {N _ {R}} = (B _ {0} \hat {\omega_ {0}} - B _ {2} \hat {\omega_ {0}} ^ {3}) \sin (\hat {\omega_ {0}} \tau_ {1 0}) - C _ {1} \hat {\omega_ {0}} ^ {2} \cos (\hat {\omega_ {0}} \tau_ {1 0}), \\\widehat {N _ {I}} = \left(B _ {0} \widehat {\omega_ {0}} - B _ {2} \widehat {\omega_ {0}} ^ {3}\right) \cos (\widehat {\omega_ {0}} \tau_ {1 0}) + B _ {1} \widehat {\omega_ {0}} ^ {2} \sin (\widehat {\omega_ {0}} \tau_ {1 0}), \quad \widehat {Q _ {R}} = B _ {1}, \quad \widehat {Q _ {I}} = 2 B _ {2} \widehat {\omega_ {0}}, \\\widehat {P _ {R}} = - B _ {1} \hat {\omega_ {0}} ^ {2}, \quad \widehat {P _ {I}} = B _ {0} \hat {\omega_ {0}} - B _ {2} \hat {\omega_ {0}} ^ {3}. \\\end{array} $$ $$ \left[ \operatorname {R e} \left(\frac {d \xi}{d \tau_ {1}}\right) _ {\tau_ {1} = \tau_ {1 0}, \xi = i \hat {\omega}} \right] ^ {- 1} = \frac {A _ {*}}{B _ {*}} + \frac {C _ {*}}{D _ {*}} = \frac {A _ {*} D _ {*} + B _ {*} C _ {*}}{B _ {*} D _ {*}}, \tag {4.18} $$ $$ A _ {*} = \widehat {M _ {R} N _ {R}} + \widehat {M _ {I}} \widehat {N _ {I}}, \quad B _ {*} = \widehat {N _ {R}} ^ {2} + \widehat {N _ {I}} ^ {2}, $$ $$ C _ {*} = \widehat {Q _ {R}} \widehat {P _ {R}} + \widehat {Q _ {I}} \widehat {P _ {I}}, D _ {*} = \widehat {P _ {R}} ^ {2} + \widehat {P _ {I}} ^ {2}. $$ From this, we can get $$ \left[ \mathrm {R e} (\frac {d \xi}{d \tau_ {1}}) _ {\tau_ {1} = \tau_ {1 0}, \xi = i \dot {\omega}} \right] ^ {- 1} = \mathrm {s g n} [ A _ {*} D _ {*} + B _ {*} C _ {*} ]. $$ If (H6): $A_{*}D_{*} + B_{*}C_{*} \neq 0$ holds, the transversal condition $[\mathrm{Re}(\frac{d\xi}{d\tau_1})_{\tau_1 = \tau_{10},\xi = i\hat{\omega}}]^{-1} \neq 0$. From the above analysis, the following theorem can be drawn Theorem 4.2. For $\tau_{2} = 0$ and $\tau_{1} > 0$, we have the following results: (i) If (H4) holds, then the equilibrium $E^{*}$ is asymptotically stable for all $\tau_{1} > 0$. (ii) If (H6) and (H5) hold, then the equilibrium $E^{*}$ is locally asymptotically stable for all $\tau_{1} < \tau_{10}$ together with unstable for $\tau_{1} > \tau_{10}$ and undergoes Hopf bifurcation at $\tau_{1} = \tau_{10}$. Case (4): $\tau_{1}$ is fixed in $(0, \tau_{10}]$ and $\tau_{2} > 0$. 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 $$ \begin{array}{l} u ^ {3} - 3 u \omega^ {2} + A _ {2} (u ^ {2} - \omega^ {2}) + A _ {1} u + A _ {0} + (B _ {2} u ^ {2} - B _ {2} \omega^ {2} + B _ {1} u + B _ {0}) e ^ {- u \tau_ {1}} \cos (\omega \tau_ {1}) \\+ \left(2 B _ {2} u \omega + B _ {1} \omega\right) e ^ {- u \tau_ {1}} \sin \left(\omega \tau_ {1}\right) + \left(C _ {2} u ^ {2} - C _ {2} \omega^ {2} + C _ {1} u + C _ {0}\right) e ^ {- u \tau_ {1}} \tag {4.19} \\\end{array} $$ $$ \cos (\omega \tau_ {2}) + (2 C _ {2} u \omega + C _ {1} \omega) \sin (\omega \tau_ {2}) = 0. $$ $$ \begin{array}{l} 3 u ^ {2} \omega - \omega^ {3} + 2 A _ {2} u \omega + A _ {1} \omega - \left(B _ {2} u ^ {2} - B _ {2} \omega^ {2} + B _ {1} u + B _ {0}\right) \sin (\omega \tau_ {1}) + (2 B _ {2} u \omega \\+ B _ {1} \omega) e ^ {- u \tau_ {1}} \cos (\omega \tau_ {1}) - \left(C _ {2} u ^ {2} - C _ {2} \omega^ {2} + C _ {1} u + C _ {0}\right) \sin (\omega \tau_ {2}) + (2 C _ {2} u \omega \\+ C _ {1} \omega) e ^ {- u \tau_ {2}} \cos (\omega \tau_ {2}) = 0. \\\end{array} $$ Putting $u = 0$ in Eqs.(4.19) and (4.20), we obtain $$ A _ {2} \omega^ {2} - A _ {0} = (- B _ {2} \omega^ {2} + B _ {0}) \cos (\omega \tau_ {1}) + (C _ {0} - C _ {2} \omega^ {2}) \cos (\omega \tau_ {2}) + B _ {1} \omega \sin (\omega \tau_ {1}) + C _ {1} \omega \sin (\omega \tau_ {2}). \qquad (4. 2 1) $$ $$ \omega^ {3} - A _ {1} \omega = - \left(B _ {0} - B _ {2} \omega^ {2}\right) \sin \left(\omega \tau_ {1}\right) + B _ {1} \omega \cos \left(\omega \tau_ {1}\right) - \left(C _ {0} - C _ {2} \omega^ {2}\right) \sin \left(\omega \tau_ {2}\right) + C _ {1} \omega \cos \left(\omega \tau_ {2}\right). \tag {4.22} $$ Squaring and adding Eqs.(4.21) and (4.22) to eliminate $\tau_{2}$, we get $$ \omega^ {6} + \tilde {a} _ {4} \omega^ {4} + \tilde {a} _ {3} \omega^ {3} + \tilde {a} _ {2} \omega^ {2} + \tilde {a} _ {0} = 0, \tag {4.23} $$ $$ \tilde {a} _ {4} = - \left(B _ {2} ^ {2} + C _ {2} ^ {2} - A _ {2} ^ {2}\right), \quad \tilde {a} _ {3} = - 2 \left(B _ {2} C _ {1} - B _ {1} C _ {2}\right) \sin (\omega \tau_ {1} - \omega \tau_ {2}), $$ $$ \tilde {a} _ {2} = - \left(\left(B _ {1} ^ {2} - 2 B _ {0} B _ {2} + C _ {1} ^ {2} - 2 C _ {0} C _ {2}\right) + 2 \left(B _ {1} C _ {1} - 2 A _ {0} A _ {2} - A _ {1} ^ {2} - B _ {2}\right)\right) \cos (\omega \tau_ {1} - \omega \tau_ {2}), $$ $$ \tilde {a} _ {0} = - \left(B _ {0} ^ {2} + C _ {0} ^ {2} - A _ {0} ^ {2}\right). $$ Noting that Eq.(4.23) is transcendental. Now, Eqs.(4.21) and (4.22) can be written as $$ \delta_ {1} \cos (\omega \tau_ {2}) + \delta_ {2} \sin (\omega \tau_ {2}) = \delta_ {3} + \delta_ {4} \cos (\omega \tau_ {1}) + \delta_ {5} \sin (\omega \tau_ {1}), \tag {4.24} $$ $$ - \delta_ {2} \cos (\omega \tau_ {2}) + \delta_ {1} \sin (\omega \tau_ {2}) = \delta_ {6} - \delta_ {5} \cos (\omega \tau_ {1}) + \delta_ {4} \sin (\omega \tau_ {1}), \tag {4.25} $$ where $$ \delta_ {1} = C _ {2} \omega^ {2} - C _ {0}, \quad \delta_ {2} = - C _ {1} \omega , $$ $$ \delta_ {3} = A _ {0} - A _ {2} \omega^ {2}, \quad \delta_ {4} = B _ {0} - B _ {2} \omega^ {2}, $$ $$ \delta_ {5} = B _ {1} \omega , \quad \delta_ {6} = \omega^ {3} - A _ {1} \omega . $$ Without losing generality, the Eq.(4.23) has finite positive roots $\widetilde{\omega_1},\widetilde{\omega_2},\dots,\widetilde{\omega_k}$, for every fixed $\widetilde{\omega}$, there exists a sequence $\{\tau_{2i}^j |j = 0,1,2\dots \}$, where $$ \tau_ {2 i} ^ {(j)} = \frac {1}{\tilde {\omega} _ {i}} \tan^ {- 1} \left[ \frac {\left(\delta_ {1} \delta_ {4} + \delta_ {2} \delta_ {4}\right) \sin \left(\tilde {\omega} _ {i} \tau_ {1}\right) - \left(\delta_ {1} \delta_ {5} - \delta_ {2} \delta_ {4}\right) \cos \left(\tilde {\omega} _ {i} \tau_ {1}\right) + \delta_ {1} \delta_ {6} + \delta_ {2} \delta_ {3}}{\left(\delta_ {1} \delta_ {5} - \delta_ {2} \delta_ {4}\right) \sin \left(\tilde {\omega} _ {i} \tau_ {1}\right) + \left(\delta_ {2} \delta_ {5} + \delta_ {1} \delta_ {4}\right) \cos \left(\tilde {\omega} _ {i} \tau_ {1}\right) + \delta_ {1} \delta_ {3} - \delta_ {2} \delta_ {4}} + \frac {k \pi}{\tilde {\omega} _ {i}} \right. $$ $$ j = 0, 1, 2, \dots $$ let $\widetilde{\tau}_2 = \min\{\tau_{2i}^{(j)} | i = 0,1,2,\ldots,k, j = 0,1,2\dots\}$, when $\tau_2 = \widetilde{\tau}_2$, $\widetilde{\omega} = \widetilde{\omega}_i|_{\tau_2 = \widetilde{\tau}_2}$, $i = 1,2,3,\ldots$, the characteristic equation (4.4) has purely imaginary roots $\pm i\widetilde{\omega}$. Then, we will verify the transversality condition, differentiating the characteristic equation (4.4) with respect to $\tau_2$, we can obtain $$ \begin{array}{l} [ \mathrm {R e} (\frac {d \xi}{d \tau_ {2}}) _ {\tau_ {2} = \widetilde {\tau} _ {2}, \xi = i \widetilde {\omega}} ] ^ {- 1} = \mathrm {R e} [ \frac {3 (i \widetilde {\omega}) ^ {2} + 2 A _ {2} (i \widetilde {\omega}) + A _ {1}}{(i \widetilde {\omega}) (C _ {2} (i \widetilde {\omega}) ^ {2} + C _ {1} (i \widetilde {\omega}) + C _ {0}) (\cos (\widetilde {\omega} \widetilde {\tau} _ {2}) - i \sin (\widetilde {\omega} \widetilde {\tau} _ {2}))} ] \\+ \operatorname {R e} \left[ \frac {2 C _ {2} (i \widetilde {\omega}) + C _ {1}}{(i \widetilde {\omega}) \left(C _ {2} (i \widetilde {\omega}) ^ {2} + C _ {1} (i \widetilde {\omega}) + C _ {0}\right)} \right]. \\\end{array} $$ Now $$ \left[ \operatorname {R e} \left(\frac {d \xi}{d \tau_ {2}}\right) _ {\tau_ {2} = \tilde {\tau} _ {2}, \xi = i \tilde {\omega}} \right] ^ {- 1} = \operatorname {R e} \left[ \frac {M _ {R} + M _ {I} i}{N _ {R} + N _ {I} i} \right] + \operatorname {R e} \left[ \frac {Q _ {R} + Q _ {I} i}{P _ {R} + P _ {I} i} \right] = \frac {M _ {R} N _ {R} + M _ {I} N _ {I}}{N _ {R} ^ {2} + N _ {I} ^ {2}} + \frac {Q _ {R} P _ {R} + Q _ {I} P _ {I}}{P _ {R} ^ {2} + P _ {I} ^ {2}}, $$ where $$ M _ {R} = - 3 \widetilde {\omega} ^ {2} + A _ {1}, \quad M _ {I} = 2 A _ {2} \widetilde {\omega}, \quad N _ {R} = \left(C _ {0} \widetilde {\omega} - C _ {1} \widetilde {\omega} ^ {2} - C _ {2} \widetilde {\omega} ^ {3}\right) \sin (\widetilde {\omega} \bar {\tau} _ {2}) $$ $$ N _ {I} = \left(C _ {0} \widetilde {\omega} - C _ {2} \widetilde {\omega} ^ {3}\right) \cos (\widetilde {\omega} \bar {\tau} _ {2}) + C _ {1} \widetilde {\omega} ^ {2} \sin (\widetilde {\omega} \bar {\tau} _ {2}), \quad Q _ {R} = C _ {1}, \quad Q _ {I} = 2 C _ {2} \widetilde {\omega}, $$ $$ P _ {R} = - C _ {1} \widetilde {\omega} ^ {2}, \quad P _ {I} = C _ {0} \widetilde {\omega} - C _ {2} \widetilde {\omega} ^ {3}. $$ $$ \left[ \operatorname {R e} \left(\frac {d \xi}{d \tau_ {2}}\right) _ {\tau_ {2} = \tilde {\tau} _ {2}, \xi = i \tilde {\omega}} \right] ^ {- 1} = \frac {E}{F} + \frac {G}{H} = \frac {E H + F G}{F H}, \tag {4.27} $$ $$ E = M _ {R} N _ {R} + M _ {I} N _ {I}, \quad F = N _ {R} ^ {2} + N _ {I} ^ {2}, $$ $$ G = Q _ {R} P _ {R} + Q _ {I} P _ {I}, \quad H = P _ {R} ^ {2} + P _ {I} ^ {2}. $$ From this we can get $$ \operatorname {s g n} \left[ \operatorname {R e} \left(\frac {d \xi}{d \tau_ {2}}\right) _ {\tau_ {2} = \tilde {\tau} _ {2}, \xi = i \tilde {\omega}} \right] ^ {- 1} = \operatorname {s g n} \left[ E H + F G \right]. $$ If (H7): $EH + FG \neq 0$ holds, the transversal condition $\mathrm{sgn}[\mathrm{Re}(\frac{d\xi}{d\tau_2})_{\tau_2 = \widetilde{\tau}_2,\xi = i\widetilde{\omega}}]^{-1} \neq 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 $\tau_{2} > 0$, then the equilibrium $E^{*}$ is locally asymptotically stable for $\tau_{2} \in (0, \widetilde{\tau}_{2})$ whereas system (2.2) undergoes Hopf bifurcation at $\tau_{2} = \widetilde{\tau}_{2}$. Case(5): $\tau_{2}$ is fixed in $(0, \tau_{20}]$ and $\tau_{1} > 0$, 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. $$ E (t) = \varepsilon Q (t), \quad 0 \leq \varepsilon \leq 1, \quad \frac {d Q}{d t} = I (t) - \gamma Q (t), \quad Q (0) = Q _ {0}. \tag {5.1} $$ 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 $\varepsilon = 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 $\varepsilon = 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 ( $v > 0$ ) on the unit biomass of terrestrial fish. When ( $v < 0$ ) can be understood as any subsidy to fishermen. Net income of fishermen ('Net income' for short) is $E[(u_1 - v)q_1N + (u_2 - v)q_2N - 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 $$ I = E \varphi \left[ \left(u _ {1} - v\right) q _ {1} N + \left(u _ {2} - v\right) q _ {2} Z - C \right], \quad 0 \leq \varphi < 1. \tag {5.2} $$ For $\varphi = 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. $\varphi = 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 $I < 0$, $\varphi > 0$. By combining Eqs.(5.1) and (5.2), we can get $$ \frac {d E}{d t} = E \left\{\varepsilon \varphi \left[ \left(u _ {1} - v\right) q _ {1} N + \left(u _ {2} - v\right) q _ {2} Z - C \right] - \gamma \right\}. \tag {5.3} $$ 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 $$ M E = E \left[ \left(u _ {1} - v\right) q _ {1} N + \left(u _ {2} - v\right) q _ {2} Z - C \right] + E \left[ v \left(q _ {1} N\right) + v \left(q _ {2} N\right) \right], $$ 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 $$ \left\{ \begin{array}{l} \frac {d N}{d t} = r _ {1} N \left(1 - \frac {N + \alpha_ {1} T}{k _ {1}}\right) - \frac {w _ {1} N Z}{p _ {1} + N} - q _ {1} E N, \\ \frac {d T}{d t} = r _ {2} N \left(1 - \frac {T + \alpha_ {2} N}{k _ {2}}\right) - \frac {w _ {2} T Z}{p _ {2} + T + \beta N}, \\ \frac {d Z}{d t} = \frac {c _ {1} w _ {1} N Z}{p _ {1} + N} - \frac {c _ {2} w _ {2} T Z}{p _ {2} + T + \beta N} - d Z - q _ {2} E Z, \\ \frac {d E}{d t} = E \{\varepsilon \varphi [ (u _ {1} - v) q _ {1} N + (u _ {2} - v) q _ {2} Z - C ] - \gamma \}. \end{array} \right. \tag {5.4} $$ 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)$. $$ J = \int_ {0} ^ {+ \infty} E (t) e ^ {- \delta t} [ u _ {1} q _ {1} N + u _ {2} q _ {2} Z - C ] d t, \tag {5.5} $$ 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_{\mathrm{min}} \leq v(t) \leq v_{\mathrm{max}}$ on the control variable $v(t)$. When $v_{\mathrm{min}} < 0$, it will have the effect of accelerating the rate of fishery expansion. The Hamiltonian of the problem is obtained by $$ \begin{array}{l} H = (u _ {1} q _ {1} N + u _ {2} q _ {2} Z - C) E e ^ {- \delta t} + \lambda_ {1} N [ r _ {1} (1 - \frac {N + \alpha_ {1} T}{k _ {1}}) - \frac {w _ {1} Z}{p _ {1} + N} - q _ {1} E ] \\+ \lambda_ {2} \left[ r _ {2} T \left(1 - \frac {T + \alpha_ {1} N}{k _ {2}}\right) - \frac {w _ {2} T Z}{p _ {2} + T + \beta N} \right] + \lambda_ {3} \left[ \frac {c _ {1} w _ {1} N Z}{p _ {1} + N} - \frac {c _ {2} w _ {2} T Z}{p _ {2} + T + \beta N} - d Z - q _ {2} E Z \right] \tag {5.6} \\+ \lambda_ {4} E \{\varepsilon \varphi [ (u _ {1} - v) q _ {1} N + (u _ {2} - v) q _ {2} Z - C ] - \gamma \}, \\\end{array} $$ 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] $$ \frac {\partial H}{\partial v} = - \lambda_ {4} E \varepsilon \varphi \left(q _ {1} N + q _ {2} Z\right) = 0 \Rightarrow \lambda_ {4} = 0. \tag {5.7} $$ Now, the adjoint equations are $$ \begin{array}{l} \frac {d \lambda_ {1}}{d t} = - \frac {\partial H}{\partial N} = - [ u _ {1} q _ {1} E e ^ {- \delta t} + \lambda_ {1} (r _ {1} - \frac {2 r _ {1} N + r _ {1} \alpha_ {1} T}{k _ {1}} - \frac {w _ {1} p _ {1} Z}{(p _ {1} + N) ^ {2}} - q _ {1} E) + \lambda_ {2} [ \frac {w _ {2} \beta T Z}{(p _ {2} + T + \beta N) ^ {2}} - \frac {r _ {2} \alpha_ {2} T}{k _ {2}} ] ]. \\+ \lambda_ {3} \left(\frac {c _ {1} w _ {1} p _ {1} Z}{\left(p _ {1} + N\right) ^ {2}} + \frac {c _ {2} w _ {2} \beta T Z}{\left(p _ {2} + T + \beta N\right) ^ {2}}\right), \\\end{array} $$ $$ \begin{array}{l} \frac {d \lambda_ {2}}{d t} = - \frac {\partial H}{\partial T} = - [ \lambda_ {1} (\frac {r _ {1} \alpha_ {1} N}{k _ {1}}) + \lambda_ {2} [ r _ {2} (1 - \frac {2 T + \alpha_ {2} N}{k _ {2}}) - \frac {w _ {2} Z (p _ {2} + \beta N)}{(p _ {2} + T + \beta N) ^ {2}} ] - \lambda_ {3} (\frac {c _ {2} w _ {2} Z (p _ {2} + \beta N)}{(p _ {2} + T + \beta N) ^ {2}}), \\\frac {d \lambda_ {3}}{d t} = - \frac {\partial H}{\partial Z} = - [ u _ {2} q _ {2} E e ^ {- \delta t} - \lambda_ {1} (\frac {w _ {1} N}{p _ {1} + N}) - \lambda_ {2} (\frac {w _ {2} T}{p _ {2} + T + \beta N}) + \lambda_ {3} (\frac {c _ {1} w _ {1} N}{p _ {1} + N} - \frac {c _ {2} w _ {2} T}{p _ {2} + T + \beta N} - d - q _ {2} E) ], \\\frac {d \lambda_ {4}}{d t} = - \frac {\partial H}{\partial E} = - [ (u _ {1} q _ {1} N + u _ {2} q _ {2} Z - C) e ^ {- \delta t} - \lambda_ {1} q _ {1} N - \lambda_ {3} q _ {2} Z ]. \tag {5.8} \\\end{array} $$ Now start with Eqs.(5.8) and (5.7), using the equilibrium equation we have $$ \frac {d \lambda_ {1}}{d t} = - u _ {1} q _ {1} E e ^ {- \delta t} - \lambda_ {1} [ - \frac {r _ {1} N}{k _ {1}} + \frac {w _ {1} N Z}{(p _ {1} + N) ^ {2}} ] - \lambda_ {2} [ \frac {w _ {2} \beta T Z}{(p _ {2} + T + \beta N) ^ {2}} - \frac {r _ {2} \alpha_ {2} T}{k _ {2}} ] - \lambda_ {3} [ \frac {c _ {1} w _ {1} p _ {1} Z}{(p _ {1} + N) ^ {2}} + \frac {c _ {2} w _ {2} \beta T Z}{(p _ {2} + T + \beta N) ^ {2}} ], $$ $$ \frac {d \lambda_ {2}}{d t} = - \lambda_ {1} [ \frac {r _ {1} \alpha_ {1} N}{k _ {1}} ] - \lambda_ {2} [ \frac {w _ {2} T Z}{(p _ {2} + T + \beta N) ^ {2}} ] - \lambda_ {3} [ \frac {c _ {2} w _ {2} Z (p _ {2} + \beta N)}{(p _ {2} + T + \beta N) ^ {2}} ], $$ $$ \frac {d \lambda_ {3}}{d t} = - u _ {2} q _ {2} E e ^ {- \delta t} + \lambda_ {1} \left(\frac {w _ {1} N}{p _ {1} + N}\right) + \lambda_ {2} \left(\frac {w _ {2} T}{p _ {2} + T + \beta N}\right). \tag {5.9} $$ Using the second and third equations of Equation (5.9) from the fourth equation of Equation (5.8), we can obtain $\frac{d\lambda_1}{dt} = M_1e^{-\delta t} + M_2\lambda_1 + M_3\lambda_2$, where $$ M _ {1} = \frac {(C - u _ {1} q _ {1} N) \delta + u _ {2} q _ {2} Z (q _ {2} E - \delta)}{q _ {1} N}, M _ {2} = - \frac {w _ {1} q _ {2} N Z}{(p _ {1} + N) q _ {1} N}, M _ {3} = - \frac {w _ {2} q _ {2} T Z}{(p _ {2} + T + \beta N) q _ {1} N}. $$ The solution of this linear equation is $$ \lambda_ {1} = N _ {0} e ^ {- M _ {2} t} - \frac {M _ {1} e ^ {- \delta t}}{M _ {2} + \delta} - \frac {M _ {3} \lambda_ {2}}{M _ {2}}. \tag {5.10} $$ Using the same method as above, we can get $$ \lambda_ {3} = I _ {0} e ^ {H _ {2} t} - \frac {H _ {1} e ^ {- \delta t}}{H _ {2} + \delta}, \tag {5.11} $$ where $$ H _ {1} = \left[ \frac {(C - u _ {2} q _ {2} Z) \delta - q _ {1} N (u _ {1} \delta + M _ {1})}{q _ {2} Z} + \frac {M _ {1} M _ {2} q _ {1} N}{(M _ {2} + \delta) q _ {2} Z} \right], \quad H _ {2} = \frac {M _ {2} M _ {3} q _ {1} N}{q _ {2} M _ {2} Z}. $$ Identically $$ \frac {d \lambda_ {2}}{d t} = R _ {1} e ^ {- \delta t} + R _ {2} \lambda_ {2}, \tag {5.12} $$ $$ R _ {1} = \frac {M _ {1}}{M _ {2} + \delta} + \frac {H _ {1}}{H _ {2} + \delta} \left(\frac {c _ {2} w _ {2} Z (p _ {2} + \beta N)}{\left(p _ {2} + T + \beta N\right) ^ {2}}\right), \quad R _ {2} = \frac {M _ {3}}{M _ {2}} \left(\frac {r _ {2} \alpha_ {1} N}{k _ {1}}\right) - \frac {w _ {2} T Z}{\left(p _ {2} + T + \beta N\right) ^ {2}}. $$ So we can get $\lambda_{1}$ $$ \lambda_ {1} = N _ {0} e ^ {M _ {2} t} - \frac {M _ {1} e ^ {- \delta t}}{M _ {2} + \delta} - \frac {M _ {3} (W _ {0} e ^ {R _ {2} t} - \frac {R _ {1} e ^ {- \delta t}}{R _ {2} + \delta})}{M _ {2}}. $$ The shadow price $\lambda_1 e^{-\delta t}$ is bounded as $t \to \infty$, $N_0 = 0$ and $W_0 = 0$, then we can obtain $$ \lambda_ {1} = - \frac {M _ {1} e ^ {- \delta t}}{M _ {2} + \delta} - \frac {M _ {3}}{M _ {2}} \left(e ^ {R _ {2} t} - \frac {R _ {1} e ^ {- \delta t}}{R _ {2} + \delta}\right). \tag {5.13} $$ Now use Eqs.(5.11), (5.12) and (5.13) in the first of Eq.(5.9), we have $$ \begin{array}{l} [ \frac {(C - u _ {1} q _ {1} N ^ {*}) \delta + u _ {2} q _ {2} Z ^ {*} (q _ {2} E ^ {*} - \delta)}{q _ {1} N ^ {*}} ] e ^ {- \delta t} + \frac {w _ {2} q _ {2} N ^ {*} Z ^ {*}}{(p _ {1} + N ^ {*}) q _ {1} N ^ {*}} [ \frac {M _ {1} e ^ {- \delta t}}{M _ {2} + \delta} - \frac {M _ {3}}{M _ {2}} (e ^ {R _ {2} t} - \frac {R _ {1} e ^ {- \delta t}}{R _ {2} + \delta}) ] \\+ \left[ \frac {w _ {2} q _ {2} T ^ {*} Z ^ {*}}{\left(p _ {2} + T ^ {*} + \beta N ^ {*}\right) q _ {1} N ^ {*}} \right] \left[ \frac {R _ {1} e ^ {- \delta t}}{R _ {2} + \delta} \right] + u _ {1} q _ {1} E ^ {*} e ^ {- \delta t} + \left[ \frac {M _ {1} e ^ {- \delta t}}{M _ {2} + \delta} - \frac {M _ {3}}{M _ {2}} \left(e ^ {R _ {2} t} - \frac {R _ {1} e ^ {- \delta t}}{R _ {2} + \delta}\right) \right] \left[ - \frac {r _ {1} N ^ {*}}{k _ {1}} + \frac {w _ {1} N ^ {*} Z ^ {*}}{\left(p _ {1} + N ^ {*}\right) ^ {2}} \right] \tag {5.14} \\= (\frac {R _ {1} e ^ {- \delta t}}{R _ {2} + \delta}) [ \frac {w _ {2} \beta T ^ {*} Z ^ {*}}{(p _ {2} + T ^ {*} + \beta N ^ {*}) ^ {2}} - \frac {r _ {2} \alpha_ {2} T ^ {*}}{k _ {2}} ] + (\frac {H _ {1} e ^ {- \delta t}}{H _ {2} + \delta}) [ \frac {c _ {2} w _ {2} Z ^ {*} (p _ {2} + \beta N ^ {*})}{(p _ {2} + T ^ {*} + \beta N ^ {*}) ^ {2}} ]. \\\end{array} $$ Because of the computational complexity, our optimal equilibrium solution can be expressed as $$ T ^ {*} = \frac {\left[ \left(c _ {1} w _ {1} - \delta\right) N ^ {*} - \delta p _ {1} \right] \left(p _ {2} + \beta N ^ {*}\right)}{\left[ \left(c _ {2} w _ {2} - \delta\right) p _ {1} + \left(c _ {2} w _ {2} - c _ {1} w _ {2} - \delta\right) N ^ {*} \right]}, $$ $$ Z ^ {*} = r _ {1} \left(\frac {p _ {1} + N ^ {*}}{w _ {1} k _ {1}}\right) \left(k _ {1} - N ^ {*} - \alpha_ {1} T ^ {*}\right). \tag {5.15} $$ $N^{*}$ available from the following equation $$ r _ {2} \left(k _ {2} - T ^ {*} - \alpha_ {2} N ^ {*}\right) \left(p _ {2} + T ^ {*} + \beta N ^ {*}\right) - w _ {2} k _ {2} Z ^ {*} = 0. \tag {5.16} $$ $E^{*}$ available from the following equation $$ \frac {r _ {1}}{q _ {1}} \left(1 - \frac {N ^ {*} + \alpha_ {1} T ^ {*}}{k _ {1}}\right) - \frac {w _ {1} Z ^ {*}}{q _ {1} \left(p _ {1} + N ^ {*}\right)} = \frac {c _ {1} w _ {1} N ^ {*}}{q _ {2} \left(p _ {1} + N ^ {*}\right)} - \frac {c _ {2} w _ {2} T ^ {*}}{q _ {2} \left(p _ {2} + T ^ {*} + \beta N ^ {*}\right)} - \frac {d}{q _ {2}}. \tag {5.17} $$ From the complex calculation results, it can be seen that $T^{*}$ and $Z^{*}$ are functions of $v$. Therefore, we can express this function as follows $$ \begin{array}{l} \left[ \frac {\left(C - u _ {1} q _ {1} N ^ {*}\right) \delta + u _ {2} q _ {2} Z ^ {*} \left(q _ {2} E ^ {*} - \delta\right)}{q _ {1} N ^ {*}} \right] e ^ {- \delta t} + \frac {w _ {2} q _ {2} N ^ {*} Z ^ {*}}{\left(p _ {1} + N ^ {*}\right) q _ {1} N ^ {*}} \left[ \frac {M _ {1} e ^ {- \delta t}}{M _ {2} + \delta} - \frac {M _ {3}}{M _ {2}} \left(e ^ {R _ {2} t} - \frac {R _ {1} e ^ {- \delta t}}{R _ {2} + \delta}\right) \right] \\+ \left[ \frac {w _ {2} q _ {2} T ^ {*} Z ^ {*}}{\left(p _ {2} + T ^ {*} + \beta N ^ {*}\right) q _ {1} N ^ {*}} \right] \left[ \frac {R _ {1} e ^ {- \delta t}}{R _ {2} + \delta} \right] + u _ {1} q _ {1} E ^ {*} e ^ {- \delta t} + \left[ \frac {M _ {1} e ^ {- \delta t}}{M _ {2} + \delta} - \frac {M _ {3}}{M _ {2}} \left(e ^ {R _ {2} t} - \frac {R _ {1} e ^ {- \delta t}}{R _ {2} + \delta}\right) \right] \left[ - \frac {r _ {1} N ^ {*}}{k _ {1}} + \frac {w _ {1} N ^ {*} Z ^ {*}}{\left(p _ {1} + N ^ {*}\right) ^ {2}} \right] \tag {5.18} \\- (\frac {R _ {1} e ^ {- \delta t}}{R _ {2} + \delta}) [ \frac {w _ {2} \beta T ^ {*} Z ^ {*}}{(p _ {2} + T ^ {*} + \beta N ^ {*}) ^ {2}} - \frac {r _ {2} \alpha_ {2} T ^ {*}}{k _ {2}} ] - (\frac {H _ {1} e ^ {- \delta t}}{H _ {2} + \delta}) [ \frac {c _ {2} w _ {2} Z ^ {*} (p _ {2} + \beta N ^ {*})}{(p _ {2} + T ^ {*} + \beta N ^ {*}) ^ {2}} ] = f (v). \\\end{array} $$ If $v^{*}$ exists, let $v = v^{*}$ be the solution of $f(v)$. Using the value of $v^{*}$, we can get the optimal solution $(N(v^{*}), T(v^{*}), Z(v^{*}), E(v^{*}))$. Here, we establish the existence of an optimal equilibrium solution satisfying the necessary condition of the maximum principle. As Clark [23] pointed out, it is complicated to find the optimal path composed of explosive control and unbalanced singular control. Because the current model is much more complex than Clark's model, we only consider an optimal equilibrium solution. If we can begin to get $F_{c}^{*} = (N(v^{*}), T(v^{*}), Z(v^{*}), E(v^{*}))$ at any initial state in $[0, S]$ to reach its maximum benefit in a limited time $S_{0}$. The period $[0, S]$ may be a planning cycle, or it may be the shortest cycle closest to $F_{c}^{*}$, so we take $S$ to be the shortest time to reach $F_{c}^{*}$. Let $(N_{0}, T_{0}, Z_{0}, E_{0}) \in R_{+}^{4} / \{0\}$, $(N(v^{*}), T(v^{*}), Z(v^{*}), E(v^{*})) \in R_{+}^{4} / \{0\}$ be the optimal equilibrium. Now, we seek min $S_{0}(v)$ subject to the solution to Eq.(5.5). $$ N (0) = N _ {0}, \quad T (0) = T _ {0}, \quad Z (0) = Z _ {0}, \quad E (0) = E _ {0}, \quad N (S _ {0}) = N (v ^ {*}), \quad T (S _ {0}) = T (v ^ {*}), \quad Z (S _ {0}) = Z (v ^ {*}), $$ $$ E \left(S _ {0}\right) = E \left(v ^ {*}\right), \quad (N, T, Z, E) \in R _ {+} ^ {4} / \{0 \}, \quad t \in [ 0, S _ {0} ]. \tag {5.19} $$ Using the maximum principle, we will get the adjoint variables $\lambda_1$, $\lambda_2$, $\lambda_3$ and $\lambda_4$ as $$ \lambda_ {1} ^ {\prime} = - \frac {\partial H}{\partial N}, \quad \lambda_ {2} ^ {\prime} = - \frac {\partial H}{\partial T}, \quad \lambda_ {3} ^ {\prime} = - \frac {\partial H}{\partial Z}, \quad \lambda_ {4} ^ {\prime} = - \frac {\partial H}{\partial E}. \tag {5.20} $$ The adjoint variables $\lambda_{1},\lambda_{2},\lambda_{3}$ and $\lambda_4$ satisfies the another condition $$ M \{N (t), T (t), Z (t), E (t), \lambda_ {1} (t), \lambda_ {2} (t), \lambda_ {3} (t), \lambda_ {4} (t) \} = 0, \quad t \in [ 0, S _ {0} ], \tag {5.21} $$ where $$ M\{N(t),T(t),Z(t),E(t),\lambda_ {1}(t),\lambda_ {2}(t),\lambda_ {3}(t),\lambda_ {4}(t),t\} = \Bigg{\lceil}\sup_{v_{\min},v_{\max}}\Bigg{\rceil}H(N(t),T(t),Z(t),E(t),\lambda_ {1}(t),\lambda_ {2}(t),\lambda_ {3}(t),\lambda_ {4}(t),v). $$ Eq.(5.19) specifies a set of initial conditions for $\lambda_1, \lambda_2, \lambda_3$ and $\lambda_4$, and Eq.(5.20) uses these initial conditions to determine the unique solution of $\lambda_1, \lambda_2, \lambda_3$ and $\lambda_4$. Therefore, it is easy to obtain the optimal tax policy as follows: $$ \bar{v}(t) = \left\{ \begin{array}{l} v_{\max}, \quad \text{for all } t \in [0, S_0] \text{ if } \frac{\partial H}{\partial v} > 0, \\ v_{\min}, \quad \text{for all } t \in [0, S_0] \text{ if } \frac{\partial H}{\partial v} < 0. \end{array} \right. \tag{5.22} $$ The optimal path in $[0, S]$ is the solution of Eq.(5.5) in its elementary state. We will now combine these two stages to obtain the optimal tax policy and optimal path in an infinite range: $$ v (t) = \bar {v} (t), t \in [ 0, S _ {0} ], v (t) = v ^ {*}, t > S _ {0}, \Gamma (t) = \bar {\Gamma} (t), t \in [ 0, S _ {0} ], \Gamma (t) = F _ {c} = (N _ {b}, T _ {b}, Z _ {b}, E _ {b}), t > S _ {0}. $$ From the above analysis, we can easily observe the following points: - (i) From Eqs.(5.7) and (5.11)-(5.13), we note that $\lambda_{i}e^{-\delta t}$, $(i = 1,2,3,4)$, where $\lambda_{i}$ is an adjoint variable, which remains unchanged in an optimal balance time interval, therefore, they satisfy the transversal condition, that is, they remain bounded to $t\to \infty$ - (ii) Considering the coexistence equilibrium point $F_{c} = (N_{b},T_{b},Z_{b},E_{b})$, The fourth equation of Eq.(5.8) can be written as $$ \left(\lambda_ {1} q _ {1} N _ {b} + \lambda_ {2} q _ {2} Z _ {b}\right) = \left(u _ {1} q _ {1} N _ {b} + u _ {2} q _ {2} Z _ {b} - C\right) e ^ {- \delta t}. $$ This means that the total harvest cost per unit of user's effort is equal to the discount value of the future price under the steady state effort level. (iii) From Eqs.(5.11) and (5.13), we can obtain $$ u _ {1} q _ {1} N _ {b} + u _ {2} q _ {2} Z _ {b} - C = - [ (\frac {M _ {1} e ^ {- \delta t}}{M _ {2} + \delta} - \frac {M _ {3}}{M _ {2}} (\frac {R _ {1} e ^ {- \delta t}}{R _ {2} + \delta})) q _ {1} N + (\frac {H _ {1} e ^ {- \delta t}}{H _ {2} + \delta}) q _ {2} Z ] \rightarrow 0, a s \delta \rightarrow \infty . $$  Fig. 2: The optimal solution of (5.5) for $v = 0.867$. This shows that the unlimited discount rate leads to the complete dissipation of the net economic income to the society, $(u_{1}q_{1}N_{b} + u_{2}q_{2}Z_{b} - C)E = 0$. We also observe that for a zero discount rate, the present value of the continuous time flow reaches its maximum. Due to the complexity of its calculation and to explain our optimal tax policy more intuitively, we continue to study it through numerical simulation. If $r_1 = 6, r_2 = 5, \alpha_1 = 0.2, \alpha_2 = 0.2, k_1 = 100, k_2 = 190, w_1 = 0.3, w_2 = 0.3, p_1 = 50, p_2 = 50, d = 0.2, c_1 = 0.45, c_2 = 0.45, \beta = 0.3, q_1 = 0.2, q_2 = 0.2, \varepsilon = 0.2, \varphi = 0.5, \gamma = 0.2, C = 2, u_1 = 0.1, u_2 = 0.2$, and the discounting factor $\delta = 0.045$ in appropriate units, based on the selection of the above parameter values, we can get the optimal tax is $v = 0.867$. In Fig.2, we get the optimal solution. Therefore, we have used the principle of Pontryagin's maximum to obtain the optimal path of the optimal tax policy, which not only ensures the maximum goal of the regulatory authority, but also the stability of the ecosystem. ## VI. NUMERICAL SIMULATIONS In this section, we will use Matlab to numerically simulate the impact of various parameters on species and the stability of steady state. Therefore, the initial conditions and parameter settings in Table 2 are used for the numerical analysis of the system (2.3). First, we give the time series diagram of $N$, $T$ and $Z$ with short period and long period, then the impact of different $\beta$, $q_{1}E$ and $q_{2}E$ on the survival of species were investigated. Lastly, we study the changes in equilibrium stability with varying delays of time. <table><tr><td colspan="3">Table 2: All the biological descriptions of the parameters are given below:</td></tr><tr><td>Parameter</td><td>Environmental Interpretation</td><td>Value</td></tr><tr><td>(N0,T0,Z0)</td><td>Initial concentrations</td><td>(500,200,1000)</td></tr><tr><td>r1</td><td>Intrinsic growth rate of NTP</td><td>0.56</td></tr><tr><td>r2</td><td>Intrinsic growth rate of TPP</td><td>0.49</td></tr><tr><td>α1</td><td>Competitive effect of TPP on NTP</td><td>0.1</td></tr><tr><td>α2</td><td>Competitive effect of NTP on TPP</td><td>0.1</td></tr><tr><td>k1</td><td>Carrying capacity of NTP</td><td>5600</td></tr><tr><td>k2</td><td>Carrying capacity of TPP</td><td>4900</td></tr><tr><td>w1</td><td>NTP consumption rate</td><td>0.5</td></tr><tr><td>w2</td><td>TPP consumption rate</td><td>0.5</td></tr><tr><td>p1</td><td>Half saturation constants for NTP</td><td>30</td></tr><tr><td>p2</td><td>Half saturation constants for TPP</td><td>30</td></tr><tr><td>c1</td><td>the conversion rate of N to Z</td><td>0.45</td></tr><tr><td>c2</td><td>the conversion rate of T to Z</td><td>0.45</td></tr><tr><td>β</td><td>Intensity of avoidance</td><td>-</td></tr><tr><td>d</td><td>Zooplankton mortality rate</td><td>0.05</td></tr></table> ### a) Time series analysis In Fig.3, we plot the time series of $\beta = 0, \beta = 10, \beta = 1000$ in the first ten days, where the other parameter values and initial conditions are the same as in Table 2. When $q_{1} = q_{2} = 0$ and $\beta = 0$, we can observe that NTP and TPP tend to perish at a fast linear speed. It is obvious that when $\beta$ increases to 10, the concentrate of TPP will first increase to a certain concentration, then decrease and finally tend to extinction, while at this time, NTP still maintains a rapid decline rate until it is extinct(fig.3(a)(b)). On the contrary, when $\beta = 0$, we take $q_{1} = 0.4$, $q_{2} = 1.2$, and $q_{1} = 2$, $q_{2} = 2.5$, respectively. We can observe that with the increase of $q_{1}$ and $q_{2}$, NTP and zooplankton tend to become extinct at a faster rate of decline, while TPP increases more rapidly(fig.3(c)(d)). Based on the values of $q_{1}$ and $q_{2}$ of (fig.3(c)(d)), we increase $\beta$ to 10. Through comparison, we can find that the curves of NTP and zooplankton have almost no change, but the increasing speed of TPP is still accelerated(fig.3(e)(f)). To further explore the influence of $\beta$, we fixed $q_{1}$ and $q_{2}$ as 2 and 2.5, respectively. And increased the value of $\beta$ from 10 to 1000. At this time, We can observe that the concentration of NTP, TPP and zooplankton has almost no change(fig.3(g)(h)). Finally, when $\beta$ exists and is fixed at 10, we increase the concentrations of $q_{1}$ and $q_{2}$ to 6 and 8, respectively. At this time, we can observe that NTP and zooplankton accelerate the decline rate, while TPP has no obvious change(fig.3(i)(j)). In Fig.4, we draw a long-term time series diagram of the system (2.3). We fixed that $q_{1}$ and $q_{2}$ are both 0. In fig.4(a)(b), we can observe the dynamic change of $\beta$ from 0 to 10. First, we take $\beta = 0$, in fig.4(a), we will find the extinction of TPP, while NTP and zooplankton oscillate in the form of limit cycles. Next, we increase $\beta$ to 10, observe the fig.4(b), all species are in a coexistence state, and the system is stabilized to a periodic orbit. These periods show large oscillations of all populations. Secondly, when we fix $\beta = 0$ and increase $q_{1} = q_{2} = 0.1$ to $q_{1} = q_{2} = 0.36$, we can find that when $q_{1}$ and $q_{2}$ are within a certain range, NTP and TPP will coexist, and zooplankton will tend to become extinct(fig.4(c)(d)). Finally, when we fix $\beta = 10$ and increase $q_{1} = q_{2} = 0.36$ to $q_{1} = q_{2} = 0.37$, we will find that the coexistence of NTP and TPP disappears, and then only TPP exists and tends to be stable, while NTP and zooplankton tend to be extinct(fig.4(e)(f)). ### b) Double time delay analysis Now, to explore the influence of pregnancy delay $(\tau_{1})$ and toxin onset delay $(\tau_{2})$ on the stability of equilibrium point in different cases. First, we need to set a set of parameters as follows $$ r _ {1} = 2, r _ {2} = 3, \alpha_ {1} = 0. 3, \alpha_ {2} = 0. 1, k _ {1} = 2 5 0 0, k _ {2} = 3 0 0 0, w _ {1} = w _ {2} = 0. 5, p _ {1} = p _ {2} = 5 0, $$ $$ c _ {1} = c _ {2} = 0. 4 5, d = 0. 0 5, \beta = 0. 5, q _ {1} = 0. 2, q _ {2} = 0. 3, E = 1. \tag {61} $$ With initial values $(N_0, T_0, Z_0) = (400,300,500)$, we perform numerical simulations to verify the theoretical results of the previous delayed system (2.2). For these parameters, we take (6.1) into the delayed system (2.2), the complex dynamical behavior of the system has been observed with time delay. Case i: when $\tau_{1} = 0$, $\tau_{2} > 0$, in this case, $[\mathrm{Re}(\frac{d\xi}{d\tau_2})_{\tau_2 = \tau_{20},\xi = i\omega_0}]^{-1} > 0$, the transversality condition is contented. So when $\tau_{2} < \tau_{20}$ (Fig.5(a)(b)), the positive equilibrium $E^{*}$ is locally asymptotically stable, the positive equilibrium $E^{*}$ is unstable when $\tau_{2} > \tau_{20}$ (Fig.6(a)(b)), when $\tau_{2} = \tau_{20}$, the system undergoes Hopf bifurcation around the positive equilibrium $E^{*}$. (Fig.5(a)(b)) shows the trajectories and phase portrait of system (2.2) for $\tau_{1} = 0$, $\tau_{2} = 1$. It can be clearly seen that the system (2.2) will converge to the positive equilibrium point $E^{*}$. And (Fig.6(a)(b)) shows the trajectories and phase portrait of the system (2.2) for $\tau_{1} = 0$, $\tau_{2} = 1.08$. In this case, the delay system (2.2) has a periodic solution near the positive equilibrium point $(E^{*})$. Case ii: when $\tau_{1} > 0$, $\tau_{2} = 0$, we change the values of $k_{1}$ and $k_{2}$ in (6.1) to $k_{1} = 150$, $k_{2} = 250$, and the others remain unchanged. $[\mathrm{Re}(\frac{d\xi}{dt_1})_{\tau_1 = \tau_{10},\xi = i\hat{\omega}_0}]^{-1} > 0$, the transversality condition is satisfied. (Fig.7(a)(b)) shows the trajectories and phase portrait of the system (2.2) for $\tau_{1} = 0.7$, $\tau_{2} = 0$. It can be seen that although the final equilibrium point tends to be stable, there is no oscillation, indicating that there is no periodic solution in this case. Case iii: when $\tau_{1} = 0.9$ in stable interval $(0, \tau_{10})$, and take $\tau_{2} > 0$ as the parameter, $[\mathrm{Re}(\frac{d\xi}{d\tau_2})_{\tau_2 = \tilde{\tau}_2,\xi = i\tilde{\omega}}]^{-1} \neq 0$, the transversality condition is satisfied. So when $$ \tau_ {2} < \widetilde {\tau} _ {2} $$ (Fig.8(a)(b)), the positive equilibrium $E^{*}$ is locally asymptotically stable, the positive equilibrium $E^{*}$ is unstable when $\tau_{2} > \widetilde{\tau}_{2}$ (Fig.9(a)(b)), when $\tau_{2} = \widetilde{\tau}_{2}$, the system undergoes Hopf bifurcation around the positive equilibrium $E^{*}$. (Fig.8(a)(b)) shows the trajectories and phase portrait of the system (2.2) for $\tau_{1} = 0.9$, $\tau_{2} = 1.06$. It can be clearly seen that the system (2.2) will converge to the positive equilibrium point $E^{*}$. And (Fig.9(a)(b)) shows the trajectories and phase portrait of the system (2.2) for $\tau_{1} = 0.9$, $\tau_{2} = 1.09$; we find the delayed system (2.2) has periodic solutions near the positive equilibrium point $E^{*}$ in this case. Therefore, through the above numerical simulation, we can evidently find the system is stable for small values of the delay, but as the value of delay crosses its critical value, the system loses its stability and undergoes Hopf-bifurcation. Thus the limit cycle exists for $\tau_{1} > \tau_{10}$, $\tau_{2} > \tau_{20}$ and $\tau_{2} > \widetilde{\tau}_{2}$. ## VII. DISCUSSION           Fig. 3: The dynamic changes of the system (1) with different $\beta$, $q_{1}$ and $q_{2}$ in the first 10 days, other parameter values and initial conditions are the same as Table 2. (a)(b): In the case of $q_{1} = q_{2} = 0$, $\beta = 0$ and $\beta = 10$, the TPP concentration will fluctuate and the NTP concentration will barely change. (c)(d): For $\beta = 0$, the concentrations of $q_{1}$ and $q_{2}$ increase, and both NTP and TPP concentrations accelerate towards extinction. (e)(f): Based on (c)(d), for $\beta = 10$, TPP reached a higher flowering concentration, while NTP still maintained a lower concentration. (g)(h): Based on (f), for $\beta = 1000$, NTP and TPP concentrations are almost unchanged. (i)(j): for $\beta = 10$, we increase the concentrations of $q_{1}$ and $q_{2}$ to 6 and 8, respectively. NTP and zooplankton accelerate the decline rate, while TPP has no obvious change.       Fig. 4: The long-term dynamics of the system (2.1), all other parameter values are the same as Table 2. (a): When $q_{1} = q_{2} = 0$, NTP and zooplankton with initial concentrations (500,200,1000) oscillate and TPP populations become extinct. (b): For $\beta = 10$, all populations survive and the system stabilizes to a limit cycle. (c)(d): For $\beta = 0$, $0 \leq q_{1} = q_{2} \leq 0.36$, NTP and TPP can coexist. (e)(f): when we fix $\beta = 10$ and increase $q_{1} = q_{2} = 0.36$ to $q_{1} = q_{2} = 0.37$, we will find that the coexistence of NTP and TPP disappears, and then only TPP exists and tends to be stable, while NTP and zooplankton tend to be extinct.    Fig. 5: The behavior of the system(2.2) for $\tau_{1} = 0, \tau_{2} = 1$ with other parameters chosen in (6.1).   Fig. 6: The behavior of the system(2.2) for $\tau_{1} = 0, \tau_{2} = 1.08$ with other parameters chosen in (6.1).  Fig. 7: The behavior of the system(2.2) for $\tau_{1} = 0.7, \tau_{2} = 0$ with other parameters chosen in (6.1).    Fig. 8: The behavior of the system(2.2) for $\tau_{1} = 0.9, \tau_{2} = 1.06$ with other parameters chosen in (6.1).  Fig. 9: The behavior of the system(2.2) for $\tau_{1} = 0.9, \tau_{2} = 1.09$ with other parameters chosen in (6.1). The predator avoidance effect always attracts ecologists to investigate it. In the aquatic system, zooplankton lives in the environment full of toxic and non-toxic bait (phytoplankton). To make toxic phytoplankton, nontoxic phytoplankton and zooplankton coexist, the avoidance behavior of zooplankton against toxic phytoplankton is an important research topic. In this paper, we consider a biological model with two delays in which zooplankton avoids poisonous phytoplankton in the presence of nontoxic phytoplankton. For this model of poisonous avoidance, due to the avoidance coefficient of zooplankton to toxic phytoplankton, the growth density of zooplankton and toxic phytoplankton is nonlinear. When the poisonous avoidance coefficient is high, the density of poisonous phytoplankton will increase in proportion, and finally tend to be stable. We also consider the impact of human harvest on the coexistence of these three species, the form of avoidance and human harvest have biological significance, which we also analyzed. According to this article, we analyze the positive and boundedness of the system solution without time delay at first. In the bounded area, the densities of nontoxic phytoplankton (NTP), toxic phytoplankton (TPP) and zooplankton (zooplankton) are all non negative. Then we analyze the bistability of the equilibrium points. From fig.1, we can see the bistability of each equilibrium point in different $k_{1}$ ranges. For the dynamic behavior of double time-delay systems, we analyze the local stability and the existence of Hopf bifurcation. Taking the pregnancy delay $\tau_{1}$ and the toxin onset delay $\tau_{2}$ as the bifurcation parameters, the critical value of the time delay for the Hopf bifurcation of the system under different conditions is obtained. We find that the system is stable when the time delay is less than this critical value $(\tau_{1}^{0},\tau_{2}^{0},\tau_{10}^{*}$ and $\tau_{20}^{*}$, respectively), but when we increase the time delay to more than this critical value, the system will become unstable, and then Hopf bifurcation occurs at the critical time. Considering the practical significance of the research, in section 5, we use the principle of Pontryagin's maximum to study the optimal tax policy of the system without time delay, we obtained the optimal path of the optimal tax policy. 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. - [1] A. Mitra, K.J. Flynn, Accounting for variation in prey selectivity by zooplankton, Ecol. Model. 199 (2006) 82-92. - [2] S.B. Linhart, J.D. Roberts, S.A. Shumake, R. Johnson, Avoidance of prey by captive coyotes punished with electric shock, In: Proceedings of the Vertebrate Pest Conference,escholarship, pp. 7 (1976) 302-330. - [3] S. Ghorai, B. Chakraborty, N. Bairagi, Preferential selection of zooplankton and emergence of spatiotemporal patterns in plankton population, Chaos. Soliton. Fract. 153 (2021) 111471. - [4] X.J. Liu, C.H. Tang, C.K. Wong, Microzooplankton grazing and selective feeding during bloom periods in the Tolo Harbour area as revealed by HPLC pigment analysis, J. Sea Res. 90 (2014) 83-94. - [5] Y.L. Zheng, X. Gong, H.W. Gao, Selective grazing of zooplankton on phytoplankton defines rapid algal succession and blooms in oceans, Ecol.1 Modelling. 468 (2022) 109947. - [6] K.G. Porter, Selective grazing and differential digestion of algae by zooplankton, Nature. 244 (1973) 179-180. - [7] B.W. Frost, Effect of size and density of food particles on the feeding behaviour of the marine planktonic copepod Calanuspacificus, Limnol. Oceanogr. 17 (1972) 1752-1765. - [8] Q.Y. Zhao, S.T. Liu, D.D. Tian, Dynamic behavior analysis of phytoplankton zooplankton system with cell size and time delay, Chaos. Solitons. Fractals. 113 (2018) 160-168. - [9] C.J. Zilverberg, J. Angerer, J. Williams, L.J. Metz, K. Harmony, Sensitivity of diet choices and environmental outcomes to a selective grazing algorithm, Ecol. Model. 390 (2018) 10-22. - [10] S. Uye, K. Takamatsu, Feeding interactions between planktonic copepods and red-tide flagellates from Japanese coastal waters, Mar. Ecol. Prog. Ser. 59 (1990) 97-107. - [11] J. Sole, E. Garcia-Ladona, M. Estrada, The role of selective predation in harmful algal blooms, J. Mar. Syst. 62 (2006) 46-54. - [12] K. Agnihotri, H. Kaur, Optimal control of harvesting effort in a phytoplankton zooplankton model with infected zooplankton under the influence of toxicity, Math. Comput. Simul. 19 (2021) 946-964. - [13] J. Chattopadhyay, Effect of toxic substances on a two-species competitive system, Ecol. Model. 84 (13) (1996) 287-289. - [14] S. Chakraborty, S. Bhattacharya, U. Feudel, J. Chattopadhyay, The role of avoidance by zooplankton for survival and dominance of toxic phytoplankton, Ecol. Complexity 11 (2012) 144-153. - [15] S. Khare, O.P. Misra, J. Dhar, Role of toxin producing phytoplankton on a plankton ecosystem, Nonlinear Anal. 4 (3) (2010) 496-502. - [16] N. Turriff, J.A. Runge, A.D. Cembella, Toxin accumulation and feeding behaviour of the planktonic-copepod Calanus jinmarchicus exposed to the red -tide dinoflagellate Alexandrium excavatum, Mar. Biol. 123 (1995) 55-64. - [17] S.J. Jang, J. Baglama, J. Rick, Nutrient-phytoplankton-zooplankton models with a toxin, Math. Comput. Modelling. 43 (12) (2006) 105-118. - [18] B. Dubey, J. Hussain, A model for the allelopathic effect on two competing species, Ecol. Model. 129 (23) (2000) 195-207. - [19] G.P. Samanta, A stochastic two species competition model: Nonequilibrium fluctuation and stability, Int. J. Stoch. Anal. 2011 (2011) 1-7. - [20] S. Roy, J. Chattopadhyay, Toxin-allelepathy among phytoplankton species prevents competitive exclusion, J. Biol. Syst. 15 (01) (2007) 73-93. - [21] Y. Pei, Y. Lv, C. Li, Evolutionary consequences of harvesting for a two-zooplankton one-phytoplankton system, Appl. Math. Model. 36 (4) (2012) 1752-1765. - [22] K. Chakraborty, S. Das, T. Kar, Optimal control of effort of a stage structured preypredator fishery model with harvesting, Nonlinear Anal. RWA. 12 (6) (2011) 3452-3467. - [23] C.W. Clark, Bioeconomic Modelling and Fisheries Management, New York (USA) Wiley. 1985. - [24] P. Panja, S.K. Mondal, D.K. Jana, Effects of toxicants on Phytoplankton-Zooplankton-Fish dynamics and harvesting, Chaos. Solitons. Fractals. 104 (2017) 389-399. - [25] F.F. Zhang, J.M. Sun, W. Tian, Spatiotemporal pattern selection in a nontoxic-phytoplankton-toxic-phytoplankton - zooplankton model with toxin avoidance effects, App. Math. Comput. 423 (2022) 127007. - [26] D. Xiao, L. Jennings, Bifurcations of a ratio-dependent predatorprey system with constant rate harvesting, Siam. Appl. Math. 65 (2005) 737-753. - [27] Y. Lv, Y. Pei, S. Gao, C. Li, Harvesting of a phytoplankton zooplankton model, Nonlinear Anal: Real World Appl. 11 (2010) 3608-3619. - [28] A.K. Sharma, A. Sharma, K. Agnihotri, Bifurcation behaviors analysis of a plankton model with multiple delays, Int. J. Biomath. 9 (06) (2016) 1650086. - [29] K. Agnihotri, H. Kaur, The dynamics of viral infection in toxin producing phytoplankton and zooplankton system with time delay, Chaos Solitons. Fractals. 118 (2019) 122-133. - [30] A. Mondal, A.K. Pal, G.P. Samanta, Rich dynamics of non-toxic phytoplankton, toxic phytoplankton and zooplankton system with multiple gestation delays, Int. J. Dyn. Control. 8 (2020) 112-131. - [31] P. Panday, S. Samanta, N. Pal, Delay induced multiple stability switch and chaos in a predator-prey model with fear effect, Math. Comput. Simulation. 172 (2020) 134-158. - [32] T. Kar, Conservation of a fishery through optimal taxation: a dynamic reaction model, Commun. Nonlinear Sci. Numer. Simul. 10 (2005) 121-131. - [33] S. Sarkar, Optimal fishery harvesting rules under uncertainty, Resour. Energy Econ. 31 (2009) 272-286. - [34] J. Sotomayor, Generic bifurcations of dynamical systems, In Dynamical systems, Academic Press. 1973.

Generating HTML Viewer...