\documentclass[11pt,twoside]{article}\makeatletter \IfFileExists{xcolor.sty}% {\RequirePackage{xcolor}}% {\RequirePackage{color}} \usepackage{colortbl} \usepackage{wrapfig} \usepackage{ifxetex} \ifxetex \usepackage{fontspec} \usepackage{xunicode} \catcode`⃥=\active \def⃥{\textbackslash} \catcode`❴=\active \def❴{\{} \catcode`❵=\active \def❵{\}} \def\textJapanese{\fontspec{Noto Sans CJK JP}} \def\textChinese{\fontspec{Noto Sans CJK SC}} \def\textKorean{\fontspec{Noto Sans CJK KR}} \setmonofont{DejaVu Sans Mono} \else \IfFileExists{utf8x.def}% {\usepackage[utf8x]{inputenc} \PrerenderUnicode{–} }% {\usepackage[utf8]{inputenc}} \usepackage[english]{babel} \usepackage[T1]{fontenc} \usepackage{float} \usepackage[]{ucs} \uc@dclc{8421}{default}{\textbackslash } \uc@dclc{10100}{default}{\{} \uc@dclc{10101}{default}{\}} \uc@dclc{8491}{default}{\AA{}} \uc@dclc{8239}{default}{\,} \uc@dclc{20154}{default}{ } \uc@dclc{10148}{default}{>} \def\textschwa{\rotatebox{-90}{e}} \def\textJapanese{} \def\textChinese{} \IfFileExists{tipa.sty}{\usepackage{tipa}}{} \fi \def\exampleFont{\ttfamily\small} \DeclareTextSymbol{\textpi}{OML}{25} \usepackage{relsize} \RequirePackage{array} \def\@testpach{\@chclass \ifnum \@lastchclass=6 \@ne \@chnum \@ne \else \ifnum \@lastchclass=7 5 \else \ifnum \@lastchclass=8 \tw@ \else \ifnum \@lastchclass=9 \thr@@ \else \z@ \ifnum \@lastchclass = 10 \else \edef\@nextchar{\expandafter\string\@nextchar}% \@chnum \if \@nextchar c\z@ \else \if \@nextchar l\@ne \else \if \@nextchar r\tw@ \else \z@ \@chclass \if\@nextchar |\@ne \else \if \@nextchar !6 \else \if \@nextchar @7 \else \if \@nextchar (8 \else \if \@nextchar )9 \else 10 \@chnum \if \@nextchar m\thr@@\else \if \@nextchar p4 \else \if \@nextchar b5 \else \z@ \@chclass \z@ \@preamerr \z@ \fi \fi \fi \fi \fi \fi \fi \fi \fi \fi \fi \fi \fi \fi \fi \fi} \gdef\arraybackslash{\let\\=\@arraycr} \def\@textsubscript#1{{\m@th\ensuremath{_{\mbox{\fontsize\sf@size\z@#1}}}}} \def\Panel#1#2#3#4{\multicolumn{#3}{){\columncolor{#2}}#4}{#1}} \def\abbr{} \def\corr{} \def\expan{} \def\gap{} \def\orig{} \def\reg{} \def\ref{} \def\sic{} \def\persName{}\def\name{} \def\placeName{} \def\orgName{} \def\textcal#1{{\fontspec{Lucida Calligraphy}#1}} \def\textgothic#1{{\fontspec{Lucida Blackletter}#1}} \def\textlarge#1{{\large #1}} \def\textoverbar#1{\ensuremath{\overline{#1}}} \def\textquoted#1{‘#1’} \def\textsmall#1{{\small #1}} \def\textsubscript#1{\@textsubscript{\selectfont#1}} \def\textxi{\ensuremath{\xi}} \def\titlem{\itshape} \newenvironment{biblfree}{}{\ifvmode\par\fi } \newenvironment{bibl}{}{} \newenvironment{byline}{\vskip6pt\itshape\fontsize{16pt}{18pt}\selectfont}{\par } \newenvironment{citbibl}{}{\ifvmode\par\fi } \newenvironment{docAuthor}{\ifvmode\vskip4pt\fontsize{16pt}{18pt}\selectfont\fi\itshape}{\ifvmode\par\fi } \newenvironment{docDate}{}{\ifvmode\par\fi } \newenvironment{docImprint}{\vskip 6pt}{\ifvmode\par\fi } \newenvironment{docTitle}{\vskip6pt\bfseries\fontsize{22pt}{25pt}\selectfont}{\par } \newenvironment{msHead}{\vskip 6pt}{\par} \newenvironment{msItem}{\vskip 6pt}{\par} \newenvironment{rubric}{}{} \newenvironment{titlePart}{}{\par } \newcolumntype{L}[1]{){\raggedright\arraybackslash}p{#1}} \newcolumntype{C}[1]{){\centering\arraybackslash}p{#1}} \newcolumntype{R}[1]{){\raggedleft\arraybackslash}p{#1}} \newcolumntype{P}[1]{){\arraybackslash}p{#1}} \newcolumntype{B}[1]{){\arraybackslash}b{#1}} \newcolumntype{M}[1]{){\arraybackslash}m{#1}} \definecolor{label}{gray}{0.75} \def\unusedattribute#1{\sout{\textcolor{label}{#1}}} \DeclareRobustCommand*{\xref}{\hyper@normalise\xref@} \def\xref@#1#2{\hyper@linkurl{#2}{#1}} \begingroup \catcode`\_=\active \gdef_#1{\ensuremath{\sb{\mathrm{#1}}}} \endgroup \mathcode`\_=\string"8000 \catcode`\_=12\relax \usepackage[a4paper,twoside,lmargin=1in,rmargin=1in,tmargin=1in,bmargin=1in,marginparwidth=0.75in]{geometry} \usepackage{framed} \definecolor{shadecolor}{gray}{0.95} \usepackage{longtable} \usepackage[normalem]{ulem} \usepackage{fancyvrb} \usepackage{fancyhdr} \usepackage{graphicx} \usepackage{marginnote} \renewcommand{\@cite}[1]{#1} \renewcommand*{\marginfont}{\itshape\footnotesize} \def\Gin@extensions{.pdf,.png,.jpg,.mps,.tif} \pagestyle{fancy} \usepackage[pdftitle={Balancing Coexistence: Ecological Dynamics and Optimal Tax Policies in a Dual Phytoplankton-Zooplankton System Influenced by Toxin Avoidance and Harvesting}, pdfauthor={}]{hyperref} \hyperbaseurl{} \paperwidth210mm \paperheight297mm \def\@pnumwidth{1.55em} \def\@tocrmarg {2.55em} \def\@dotsep{4.5} \setcounter{tocdepth}{3} \clubpenalty=8000 \emergencystretch 3em \hbadness=4000 \hyphenpenalty=400 \pretolerance=750 \tolerance=2000 \vbadness=4000 \widowpenalty=10000 \renewcommand\section{\@startsection {section}{1}{\z@}% {-1.75ex \@plus -0.5ex \@minus -.2ex}% {0.5ex \@plus .2ex}% {\reset@font\Large\bfseries}} \renewcommand\subsection{\@startsection{subsection}{2}{\z@}% {-1.75ex\@plus -0.5ex \@minus- .2ex}% {0.5ex \@plus .2ex}% {\reset@font\Large}} \renewcommand\subsubsection{\@startsection{subsubsection}{3}{\z@}% {-1.5ex\@plus -0.35ex \@minus -.2ex}% {0.5ex \@plus .2ex}% {\reset@font\large}} \renewcommand\paragraph{\@startsection{paragraph}{4}{\z@}% {-1ex \@plus-0.35ex \@minus -0.2ex}% {0.5ex \@plus .2ex}% {\reset@font\normalsize}} \renewcommand\subparagraph{\@startsection{subparagraph}{5}{\parindent}% {1.5ex \@plus1ex \@minus .2ex}% {-1em}% {\reset@font\normalsize\bfseries}} \def\l@section#1#2{\addpenalty{\@secpenalty} \addvspace{1.0em plus 1pt} \@tempdima 1.5em \begingroup \parindent \z@ \rightskip \@pnumwidth \parfillskip -\@pnumwidth \bfseries \leavevmode #1\hfil \hbox to\@pnumwidth{\hss #2}\par \endgroup} \def\l@subsection{\@dottedtocline{2}{1.5em}{2.3em}} \def\l@subsubsection{\@dottedtocline{3}{3.8em}{3.2em}} \def\l@paragraph{\@dottedtocline{4}{7.0em}{4.1em}} \def\l@subparagraph{\@dottedtocline{5}{10em}{5em}} \@ifundefined{c@section}{\newcounter{section}}{} \@ifundefined{c@chapter}{\newcounter{chapter}}{} \newif\if@mainmatter \@mainmattertrue \def\chaptername{Chapter} \def\frontmatter{% \pagenumbering{roman} \def\thechapter{\@roman\c@chapter} \def\theHchapter{\roman{chapter}} \def\thesection{\@roman\c@section} \def\theHsection{\roman{section}} \def\@chapapp{}% } \def\mainmatter{% \cleardoublepage \def\thechapter{\@arabic\c@chapter} \setcounter{chapter}{0} \setcounter{section}{0} \pagenumbering{arabic} \setcounter{secnumdepth}{6} \def\@chapapp{\chaptername}% \def\theHchapter{\arabic{chapter}} \def\thesection{\@arabic\c@section} \def\theHsection{\arabic{section}} } \def\backmatter{% \cleardoublepage \setcounter{chapter}{0} \setcounter{section}{0} \setcounter{secnumdepth}{2} \def\@chapapp{\appendixname}% \def\thechapter{\@Alph\c@chapter} \def\theHchapter{\Alph{chapter}} \appendix } \newenvironment{bibitemlist}[1]{% \list{\@biblabel{\@arabic\c@enumiv}}% {\settowidth\labelwidth{\@biblabel{#1}}% \leftmargin\labelwidth \advance\leftmargin\labelsep \@openbib@code \usecounter{enumiv}% \let\p@enumiv\@empty \renewcommand\theenumiv{\@arabic\c@enumiv}% }% \sloppy \clubpenalty4000 \@clubpenalty \clubpenalty \widowpenalty4000% \sfcode`\.\@m}% {\def\@noitemerr {\@latex@warning{Empty `bibitemlist' environment}}% \endlist} \def\tableofcontents{\section*{\contentsname}\@starttoc{toc}} \parskip0pt \parindent1em \def\Panel#1#2#3#4{\multicolumn{#3}{){\columncolor{#2}}#4}{#1}} \newenvironment{reflist}{% \begin{raggedright}\begin{list}{} {% \setlength{\topsep}{0pt}% \setlength{\rightmargin}{0.25in}% \setlength{\itemsep}{0pt}% \setlength{\itemindent}{0pt}% \setlength{\parskip}{0pt}% \setlength{\parsep}{2pt}% \def\makelabel##1{\itshape ##1}}% } {\end{list}\end{raggedright}} \newenvironment{sansreflist}{% \begin{raggedright}\begin{list}{} {% \setlength{\topsep}{0pt}% \setlength{\rightmargin}{0.25in}% \setlength{\itemindent}{0pt}% \setlength{\parskip}{0pt}% \setlength{\itemsep}{0pt}% \setlength{\parsep}{2pt}% \def\makelabel##1{\upshape ##1}}% } {\end{list}\end{raggedright}} \newenvironment{specHead}[2]% {\vspace{20pt}\hrule\vspace{10pt}% \phantomsection\label{#1}\markright{#2}% \pdfbookmark[2]{#2}{#1}% \hspace{-0.75in}{\bfseries\fontsize{16pt}{18pt}\selectfont#2}% }{} \def\TheFullDate{1970-01-01 (revised: 01 January 1970)} \def\TheID{\makeatother } \def\TheDate{1970-01-01} \title{Balancing Coexistence: Ecological Dynamics and Optimal Tax Policies in a Dual Phytoplankton-Zooplankton System Influenced by Toxin Avoidance and Harvesting} \author{}\makeatletter \makeatletter \newcommand*{\cleartoleftpage}{% \clearpage \if@twoside \ifodd\c@page \hbox{}\newpage \if@twocolumn \hbox{}\newpage \fi \fi \fi } \makeatother \makeatletter \thispagestyle{empty} \markright{\@title}\markboth{\@title}{\@author} \renewcommand\small{\@setfontsize\small{9pt}{11pt}\abovedisplayskip 8.5\p@ plus3\p@ minus4\p@ \belowdisplayskip \abovedisplayskip \abovedisplayshortskip \z@ plus2\p@ \belowdisplayshortskip 4\p@ plus2\p@ minus2\p@ \def\@listi{\leftmargin\leftmargini \topsep 2\p@ plus1\p@ minus1\p@ \parsep 2\p@ plus\p@ minus\p@ \itemsep 1pt} } \makeatother \fvset{frame=single,numberblanklines=false,xleftmargin=5mm,xrightmargin=5mm} \fancyhf{} \setlength{\headheight}{14pt} \fancyhead[LE]{\bfseries\leftmark} \fancyhead[RO]{\bfseries\rightmark} \fancyfoot[RO]{} \fancyfoot[CO]{\thepage} \fancyfoot[LO]{\TheID} \fancyfoot[LE]{} \fancyfoot[CE]{\thepage} \fancyfoot[RE]{\TheID} \hypersetup{citebordercolor=0.75 0.75 0.75,linkbordercolor=0.75 0.75 0.75,urlbordercolor=0.75 0.75 0.75,bookmarksnumbered=true} \fancypagestyle{plain}{\fancyhead{}\renewcommand{\headrulewidth}{0pt}} \date{} \usepackage{authblk} \providecommand{\keywords}[1] { \footnotesize \textbf{\textit{Index terms---}} #1 } \usepackage{graphicx,xcolor} \definecolor{GJBlue}{HTML}{273B81} \definecolor{GJLightBlue}{HTML}{0A9DD9} \definecolor{GJMediumGrey}{HTML}{6D6E70} \definecolor{GJLightGrey}{HTML}{929497} \renewenvironment{abstract}{% \setlength{\parindent}{0pt}\raggedright \textcolor{GJMediumGrey}{\rule{\textwidth}{2pt}} \vskip16pt \textcolor{GJBlue}{\large\bfseries\abstractname\space} }{% \vskip8pt \textcolor{GJMediumGrey}{\rule{\textwidth}{2pt}} \vskip16pt } \usepackage[absolute,overlay]{textpos} \makeatother \usepackage{lineno} \linenumbers \begin{document} \author[1]{ Yuqin} \author[2]{Wensheng Yang} \affil[1]{ Fujian Normal University} \renewcommand\Authands{ and } \date{\small \em Received: 1 January 1970 Accepted: 1 January 1970 Published: 1 January 1970} \maketitle \begin{abstract} In recent years, the impact of toxic phytoplankton on ecological balance has attracted more and more ecologists to study. In this paper, we develop and analyze a model with three interacting species, poisonous and nontoxic phytoplankton, and zooplankton, including zooplankton avoiding toxic phytoplankton in the presence of nontoxic phytoplankton, and the impact of human harvest on the coexistence of these three species. We first introduce the poisonous avoidance coefficient ?? and the human harvest of nontoxicphytoplankton and zooplankton to investigate its impact on species coexistence. We not only find that ?? has a particular effect on the coexistence of these three species. But also that human harvest is an essential factor determining the coexistence of these three species. Secondly, pregnancy delay () and toxin onset delay ( ) are introduced to explore the influence of time delay on the behavior of dynamic systems. When the delay value exceeds its critical value, the system will lose stability and go through Hopf bifurcation. After that, we use the principle of Pontryagin's maximum to study the optimal tax policy without delay. We obtained the optimal path of the optimal tax policy. Finally, we carry out numerical simulations to verify the theoretical results. \end{abstract} \keywords{toxic phytoplankton; human harvest; time delay; optimal tax policy; hopf bifurcation} \begin{textblock*}{18cm}(1cm,1cm) % {block width} (coords) \textcolor{GJBlue}{\LARGE Global Journals \LaTeX\ JournalKaleidoscope\texttrademark} \end{textblock*} \begin{textblock*}{18cm}(1.4cm,1.5cm) % {block width} (coords) \textcolor{GJBlue}{\footnotesize \\ Artificial Intelligence formulated this projection for compatibility purposes from the original article published at Global Journals. However, this technology is currently in beta. \emph{Therefore, kindly ignore odd layouts, missed formulae, text, tables, or figures.}} \end{textblock*} \begin{textblock*}{10cm}(1.05cm,3cm) {{\textit{CrossRef DOI of original article:}} \underline{}} \end{textblock*}\let\tabcellsep& \section[{I. Introduction}]{I. Introduction}\par 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.\par 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 mechanisms of selective grazing include prey morphology (size, color, shape, and colony formation), intestinal genetic strains, and poisonous chemicals released by prey \hyperref[b7]{[6]}\hyperref[b8]{[7]}\hyperref[b9]{[8]} {\ref [9]}\hyperref[b11]{[10]}\hyperref[b12]{[11]}\hyperref[b13]{[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 \hyperref[b14]{[13]}\hyperref[b15]{[14]}\hyperref[b16]{[15]}\hyperref[b17]{[16]}\hyperref[b18]{[17]}\hyperref[b19]{[18]}\hyperref[b20]{[19]}\hyperref[b21]{[20]}.\par 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 \hyperref[b22]{[21]}\hyperref[b23]{[22]}\hyperref[b24]{[23]}\hyperref[b25]{[24]}\hyperref[b26]{[25]}\hyperref[b27]{[26]}\hyperref[b28]{[27]}. Since time delay is widely studied in the phytoplankton-zooplankton model \hyperref[b29]{[28]}\hyperref[b30]{[29]}\hyperref[b31]{[30]}\hyperref[b32]{[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 \hyperref[b33]{[32]}\hyperref[b34]{[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.\par 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 \hyperref[tab_2]{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 {\ref 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 \hyperref[tab_2]{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.\par Considering the toxin refuge of zooplankton, a nontoxic phytoplankton-toxic zooplankton model was proposed in \hyperref[b15]{[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):? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? dN dt = r 1 N (1 - N + ? 1 T k 1 ) - w 1 N Z p 1 + N , dT dt = r 2 T (1 - T + ? 2 N k 2 ) - w 2 T Z p 2 + T + ?N , dZ dt = w 1 N Z p 1 + N - w 2 T Z p 2 + T + ?N -dZ, N (0) ? 0, T (0) ? 0, Z(0) ? 0,\textbf{(2.1)}\par 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 toxinproducing phytoplankton (TPP) species, respectively. r 1 and r 2 represent the constant intrinsic growth rates of N and T , respectively. ? 1 and ? 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. ? 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 ? 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 ? 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:? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? dN dt = r 1 N (1 - N + ? 1 T k 1 ) - w 1 N Z p 1 + N -q 1 EN, dT dt = r 2 T (1 - T + ? 2 N k 2 ) - w 2 T Z p 2 + T + ?N , dZ dt = c 1 w 1 N (t -? 1 )Z(t -? 1 ) p 1 + N (t -? 1 ) - c 2 w 2 T (t -? 2 )Z(t -? 2 ) p 2 + T (t -? 2 ) + ?N (t -? 2 ) -dZ -q 2 EZ, N (0) ? 0, T (0) ? 0, Z(0) ? 0,\textbf{(2.2)}\par where N , T , and Z represent the biomass of nontoxic phytoplankton, toxic phytoplankton and zooplankton, respectively. ? 1 (? 1 > 0) and ? 2 (? 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.? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? dN dt = r 1 N (1 - N + ? 1 T k 1 ) - w 1 N Z p 1 + N -q 1 EN, dT dt = r 2 T (1 - T + ? 2 N k 2 ) - w 2 T Z p 2 + T + ?N , dZ dt = c 1 w 1 N Z p 1 + N - c 2 w 2 T Z p 2 + T + ?N -dZ -q 2 EZ, N (0) ? 0, T (0) ? 0, Z(0) ? 0.\textbf{(2.3)}\par 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.\par All the solutions with initial values of system (2.3), which start in R 3 + , are always positive and bounded.\par Proof. Firstly, we rewrite the model (2.3) and take the linear as the following form:dX dt = F (X),\textbf{(3.1)}\par where X(t) = (N, T, Z) T ? R 3 + and F (X) is simplified as the following F (X) = ? ? F 1 (X) F 2 (X) F 3 (X) ? ? = ? ? ? ? ? ? ? r 1 N (1 -N +?1T k1 ) -w1N Z p1+N -q 1 EN r 2 T (1 -T +?2N k2 ) -w2T Z p2+T +?N c1w1N Z p1+N -c2w2T Z p2+T +?N -dZ -q 2 EZ ? ? ? ? ? ? ? . \section[{Notes}]{Notes}\par We want to prove that (N (t), T (t), Z(t)) ? R 3 + for all t ? [0, +?). For system (2.3) with initial value N (0) > 0, T (0) > 0 and Z(0) > 0, we haveN (t) = N (0) exp\{ t 0 [r 1 (1 -N (s)+?1T (s) k1 ) -w1Z(s) p1+N (s) -q 1 E]ds\}, T (t) = T (0) exp\{ t 0 [r 2 (1 -T (s)+?1N (s) k2 ) - w2Z(s) p2+T (s)+?N (s) ]ds\}, Z(t) = Z(0) exp\{ t 0 [ c1w1N (s) p1+N (s) - c2w2T (s) p2+T (s)+?N (s) -d -q 2 E]ds\},\par which shows that all the solutions of system (2.3) are always positive for all t > 0.\par 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 functionW (t) = c 1 N (t) + c 2 T (t) + Z(t). (3.2)\par Then, by differentiating (3.2) concerning t, we obtaindW dt + ?W = c 1 r 1 N (1 - N + ? 1 T k 1 ) + c 2 r 2 T (1 - T + ? 1 N k 2 ) - 2c 2 w 2 T Z p 2 + T + ?N -dZ -q 2 EZ -c 1 q 1 EN + c 1 ?N + c 2 ?T + ?Z, ? c 1 r 1 N (1 - N k 1 ) + c 2 r 2 T (1 - T k 2 ) -dZ + c 1 ?N + c 2 ?T + ?Z, = - c 1 r 1 N 2 k 1 + (r 1 + ?)c 1 N - c 2 r 2 T 2 k 2 + (r 2 + ?)c 2 T + (? -d)Z, ? c 1 k 1 (r 1 + ?) 2 4r 1 + c 2 k 2 (r 2 + ?) 2 4r 2 + (? -d)Z, ? c 1 r 2 k 1 (r 1 + ?) 2 + c 2 r 1 k 2 (r 2 + ?) 2 4r 1 r 2 + (? -d)Z, when ? -d < 0, we can obtain dW dt + ?W ? c1r2k1(r1+?) 2 +c2r1k2(r2+?) 2 4r1r2 , noting ? = c1r2k1(r1+?) 2 +c2r1k2(r2+?) 2 4r1r2\par , therefore, applying a theorem on differential inequalities \hyperref[b35]{[34]}, we obtain0 ? W ? ? ? + W (N (0),T (0),Z(0)) e ?t\par , let t ? +?, W (N, T, Z) ? ? ? . So, all solutions of system (2.3) enter the regionD = \{(N, T, Z) ? R 3 + : 0 ? W (N, T, Z) ? ? ? \}.\textbf{(3.3)}\par This shows that every solution of the system is bounded. System (2.3) possesses six different equilibrium points:\par (i) the plankton-free equilibrium, E 0 = (0, 0, 0), which always exists;\par (ii) TPP and zooplankton-free equilibrium, E 1 = (k 1 , 0, 0), which is always feasible;\par (iii) NTP and zooplankton-free equilibrium, E 2 = (0, k 2 , 0), which is always feasible;\par (iv) zooplankton-free equilibrium, E 3 = ( N , T , 0), whereN = ? 1 k 2 -k 1 ? 1 ? 2 -1 - q 1 k 1 E r 1 , T = ? 2 k 1 -k 2 ? 1 ? 2 -1 ;\par (v)TPP-free equilibrium E 4 = ( N , 0, Z), whereN = (q 2 E + d)p 1 c 1 w 1 -d -q 2 E , Z = r 1 (k 1 -N ) -q 1 k 1 E(p 1 + E) k 1 w 1 ;\par (vi)the interior equilibrium, E * = (N * , T * , Z * ), whereT * = c 1 w 1 N * -(d + q 2 E)(p 1 + N * )(p 2 + ?N * ) (c 2 w 2 + d + q 2 E)(p 1 + N * ) -c 1 w 1 N * , Z * = (p 1 + N * )r 1 (k 1 -N * -? 1 T * ) -q 1 k 1 E k 1 w 1 ;\par and N * can be obtained fromr 2 (p 2 + T * + ?N * )(k 2 -T * -? 2 N * ) -w 2 k 2 Z * = 0. (3.4)\par 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 {\ref 1}.\par Equilibria analysis: Obviously, the equilibria E 0 , E 1 and E 2 always exist. The zooplankton-free equilibrium E 3 exists, let N and T both be positive, that is ? 2 > k2 k1 and ? 1 > (?1?2-1)q1k1E r1k1 + k1 k2 . The TPP-free equilibrium E 4 exists, let N and Z both be positive, that is w 1 > d+q2E c1 and k 1 > r1N r1-q1E(p1+E) . The interior equilibrium point E * exists; let N * , T * and Z * all be positive, that isk 1 > q1k1E r1 + N * + ? 1 T * , c 2 w 2 (p 1 + N * ) > c 1 w 1 N * -(d + q 2 E)(p 1 + N * )\par > 0 and Eq.(3.4) has at least one positive root.\par In the following, we summarize the eigenvalues and local stability conditions around the feasible equilibrium point of each organism of system (2.3).\par (i) The eigenvalues of the plankton-free equilibrium E 0 = (0, 0, 0) are r 1 , r 2 and -d -q 2 E. Therefore, it is a saddle point and hence always unstable.\par (ii) The eigenvalues of the TPP and zooplankton-free equilibriumE 1 = (k 1 , 0, 0) are -r 1 -q 1 E, r 2 (1 -k1?2 k2 ) and c1w1k1 p1+k1 -d -q 2 E. When c 1 ? 1 -d -q 2 E\par ? 0, and ? 2 > k2 k1 hold, E 1 is LAS(locally asymptotically stable). On the contrary, if c 1 ? 1 -d -q 2 E > 0, ? 2 > k2 k1 and k 1 < p1(d+q2E) c1w1-d-q2E hold, we can also obtain E 1 is LAS. (iii) The eigenvalues of the NTP and zooplankton-free equilibriumE 2 = (0, k 2 , 0) are r 2 (1 -k2?1 k1 ) -q 1 E, -r 2 and -c2w2k2 p2+k2 -d -q 2 E, Therefore, E 2 is LAS if k 1 < r2?1k2 r2-q1E .\par (iv) The eigenvalues of the zooplankton-free equilibriumE 3 = ( N , T , 0) are c1w1 N p1+ N -c2w2 T p2+ T +? N -d -q 2 E\par , ? 1 and ? 2 , where ? 1 and ? 2 are the roots of the equation? 2 + b1 ? + c1 = 0,\textbf{(1)}\par b) Equilibrium points and their stability\par Noteswhere b1 = -[r 2 -r 1 + r 1 k 2 (2 N + ? 1 T ) -r 2 k 1 (2 T + ? 2 N ) k 1 k 2 ], c1 = r 1 r 2 [1 -(2 T + ? 2 N )(2 N + ? 1 T )][ 1 (2 N + ? 1 T )k 2 + 1 (2 T + ? 2 N )k 1 - 1 k 1 k 2 ] + q 1 r 2 E( k 1 (2 T + ? 2 N ) -r 1 ? 1 2 N T k 1 k 2 -1). Therefore, let c1w1 N p1+ N -c2w2 T p2+ T +? N -d -q 2 E < 0, ? 1 and ? 2 with negative real parts, that is c1w1 N p1+ N -d -q 2 E < c2w2 T p2+ T +?\par N , b1 > 0 and c1 > 0. If the above conditions are satisfied, E 3 is LAS.\par (v) The eigenvalues of the TPP-free equilibriumE 4 = ( N , 0, Z) are r 2 (1 -?2 N k2 ) -w2 Z p2+?\par N , ?1 and ?2 , where ?1 and ?2 are the roots of the equation? 2 -(ã 2 + b2 )? + ã2 b2 + c2 = 0,\textbf{(2)}\par whereã2 = (r 1 (1 -2 N k1 ) -w1p1 Z (p1+ N ) 2 -q 1 E), b2 = ( c1w1 N p1+ N -d -q 2 E), c2 = c1w1 2 p1 N Z (p1+ N ) 3 . Therefore, let r 2 (1 -?2 N k2 ) -w2 Z p2+?\par N < 0, ?1 and ?2 with negative real parts, that is (ã 2 + b2 ) < 0 and ã2 b2 + c2 > 0. If the above conditions are satisfied, E 4 is LAS.\par (vi)By solving the Jacobi determinant of E * , we can get its characteristic equation as follows? 3 + D 1 ? 2 + D 2 ? + D 3 = 0.\textbf{(3)}\par The interior equilibriumE * = (N * , T * , Z * ) is LAS if (a) D 1 > 0, (b) D 3 > 0, (c) D 1 D 2 -D 3 > 0,\par whereD 1 = -\{r 2 [1 - (2T * + ? 2 N * ) k 1 ] - w 2 Z * (p 2 + ?N * ) (p 2 + T * + ?N * ) 2 + r 1 [1 - (2N * + ? 1 T * ) k 1 ] - w 2 p 1 Z * (p 1 + N * ) 2 -q 1 E\} -( c 1 w 1 N * p 1 + N * - + r 1 ? 1 N * k 1 ( r 1 ? 1 T * k 2 + w 2 ?T * Z * ) (p 2 + T * + ?N * ) 2 ) + \{r 2 [1 - (2T * + ? 2 N * ) k 1 ] - w 2 Z * (p 2 + ?N * ) (p 2 + T * + ?N * ) 2 + r 1 [1 - (2N * + ? 1 T * ) k 1 ] - w 2 p 1 Z * (p 1 + N * ) 2 -q 1 E\} × \{ c 1 w 1 N * p 1 + N * - c 2 w 2 T * p 2 + T * + ?N * -d -q 2 E\}, D 3 =-\{ c 1 w 1 p 1 Z * (p 1 + N * ) 2 - c 2 w 2 ?T * Z * (p 2 + T * + ?N * ) 2 \} × \{- r 1 ? 1 w 2 T * k 1 (p 2 + T * + ?N * ) + w 1 N * p 1 + N * × (r 2 (1 - (2T * + ? 2 N * ) k 2 ) - w 2 Z * (p 2 + ?N * ) (p 2 + T * + ?N * ) 2 )-( c 2 w 2 Z * (p 2 + ?N * ) (p 2 + T * + ?N * ) 2 )×( w 2 T * p 2 + T * + ?N * )×[-r 1 (1- (2N * + ? 1 T * ) k 1 )+ w 1 p 1 Z * (p 1 + N * ) 2 +q 1 E]\} + w 1 N * p 1 + N * ×( r 1 ? 1 T * k 2 + w 2 ?T * Z * (p 2 + T * + ?N * ) 2 )-( c 2 w 2 Z * (p 2 + ?N * ) (p 2 + T * + ?N * ) 2 )× \{- r 1 w 2 T * p 2 + T * + ?N * + r 1 w 2 (2N * + ? 1 T * )T * k 1 (p 2 + T * + ?N * ) + w 1 w 2 p 1 T * Z * (p 2 + T * + ?N * )(p 1 + N * ) 2 + w 2 q 1 ET * p 2 + T * + ?N * + r 1 ? 1 w 1 N * T * k 2 (p 1 + N * ) + w 1 w 2 ?N * T * Z * (p 2 + T * + ?N * ) 2 (p 1 + N * ) \} + \{r 1 (1 - (2N * + ? 1 T * ) k 1 ) - w 2 p 1 Z * (p 1 + N * ) 2 -q 1 E\} × \{r 2 (1 - (2T * + ? 2 N * ) k 1 ) - w 2 Z * (p 2 + ?N * ) (p 2 + T * + ?N * ) 2 \} + r 1 ? 1 N * k 1 × ( r 1 ? 1 T * k 2 + w 2 ?T * Z * (p 2 + T * + ?N * ) 2 ).\par From the calculation of the eigenvalues, obviously, ? 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 ?, 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).\par 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 \hyperref[b15]{[14]}, we not only consider the effect of ? on species coexistence, but also human harvest as an essential factor in species coexistence.\par (i)E 0 : Extinction of all the populations at a time is impossible.\par (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 k2 ?2 < k 1 < p1(d+q2E) c1w1-d-q2E , 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. + \{r 1 [1 - (2N * + ? 1 T * ) k 1 ] - w 2 p 1 Z * (p 1 + N * ) 2 -q 1 E\} × \{r 2 [1 - (2T * + ? 2 N * ) k 1 ] - w 2 Z * (p 2 + ?N * ) (p 2 + T * + ?N * ) 2 \} (iii)E 2 :\par If the carrying capacity of NTP population (k 1 ) stays below the threshold value r2?1k2 r2-q1E , both NTP and zooplankton eventually extinct. With the competitive effect of TPP on NTP (? 1 ), the environmental carrying capacities of toxin-producing phytoplankton (k 2 ) and harvesting term for NTP and zooplankton \hyperref[b15]{[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. \section[{Ref}]{Ref}\par (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 ? 1 , k 2 and q 1 E.\par (iv)E 3 : When the carrying capacity of NTP population (k 1 ) remains within two threshold values r2?1k2 r2-q1E < k 1 < k2 ?2 (it can be obtained by the threshold value (k 1 ) of E 1 and E 2 ) together with the competitive effects (? 1 , ? 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 c1w1 N p1+ N -d -q 2 E < c2w2 T p2+ T +? 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 (?), 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 ( c1w1 N p1+ N -(d + q 2 E) > 0), we can find a threshold value of the avoidance of TPP by zooplankton (? <(c2w2 T )(p1+ N ) ( N )(c1w1 N -(d+q2E)(p1+ N )) -p2+ T N )\par , below which the zooplankton population will become extinct. On the contrary, for c1w1 N p1+ N -(d + q 2 E) < 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).\par (v)E 4 : If the carrying capacity of NTP population (k 1 ) remains within two threshold values( (d+q2E)p1 c1w1-d-q2E < k 1 < (d+q2E)(p1)+c1w1p1 c1w1-d-q2E\par ), then TPP becomes extinct under the condition ( r2(k2-?2 N ) k2 < w2 Z p2+? N ), 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 (?), 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 (? 2 ) and the TPP consumption rate (w 2 ). Similarly, for a particular parameter setup (k 2 -? 2 N > 0), we can find a threshold value of the avoidance of TPP by zooplankton (?< k2w2 Z N r2(k2-?2 N ) -p2 N )\par , below which TPP may become extinct. On the contrary, for k 2 -? 2 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 -? 2 N < 0, TPP extinction dose not depend on the avoidance of TPP by zooplankton (?) and harvest term on zooplankton (q 2 E).\par (vi)E * = (N * , T * , Z * ): When the competitive effects (? 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.\par Existence and stability conditions of the equilibrium points. \section[{Equilibrium Existence conditions}]{Equilibrium Existence conditions}\par Stability conditions E 0 = (0, 0, 0) Always exist Always unstableE 1 = (k 1 , 0, 0) Always exist (i) c 1 w 1 -d -q 2 E > 0, ? 2 > k 2 k 1 , k 1 < p 1 (d+q 2 E) c 1 w 1 -d-q 2 E , or (ii) c 1 w 1 -d -q 2 E ? 0, ? 2 > k 2 k 1 E 2 = (0, k 2 , 0) Always exist (i) k 1 < r 2 ? 1 k 2 r 2 -q 1 E E 3 = ( N, T , 0) (i) ? 2 > k 2 k 1 , (ii) ? 1 > (? 1 ? 2 -1)q 1 k 1 E r 1 k 1 + k 1 k 2 (i) c 1 w 1 N p 1 + N -d -q 2 E < c 2 w 2 T p 2 + T +? N , (ii) b1 > 0, c1 > 0 E 4 = ( N, 0, Z) (i) w 1 > d+q 2 E c 1 , (ii) k 1 > r 1 N r 1 -q 1 E(p 1 +E) (i) r 2 (1 -? 2 N k 2 ) < w 2 Z p 2 +? N , (ii) ã2 + b2 < 0, ã2 b2 + c2 > 0 E * = (N * , T * , Z * ) (i) k 1 > q 1 k 1 E r 1 + N * + ? 1 T * , (ii) c 2 w 2 (p 1 + N * ) > c 1 w 1 N * -(d + q 2 E)(p 1 + N * ) > 0, (iii) positive root of Eq.(3.4) exists (i) D 1 > 0 , (ii) D 3 > 0, (iii) D 1 D 2 -D 3 > 0 Table 1: \section[{Notes}]{Notes}\par The existence and stability of these equilibrium points are summarized in Table {\ref 1} and {\ref Fig} 1. When c 1 w 1 -d -q 2 E > 0, equilibria E 2 = (0, k 2 , 0), E 3 = ( N , T , 0), E 1 = (k 1 , 0, 0) and E 4 = ( N , 0, Z) keep stable for (0 < k 1 < r2?1k2 r2-q1E ), ( r2?1k2 r2-q1E < k 1 < k2 ?2 ), ( k2 ?2 < k 1 < p1(d+q2E) c1w1-d-q2E ) and ( (d+q2E)p1 c1w1-d-q2E < k 1 < (d+q2E)(p1)+c1w1p1 c1w1-d-q2E\par ), respectively(Fig. \hyperref[fig_0]{1}(a)). Obviously, for k 1 at the different equilibria above, the coexistence of NTP, TPP, and zooplankton requires the three ranges (k 1 > r2?1k2 r2-q1E ), (k 1 < k2 ?2 ), and (k 1 > (d+q2E)p1 c1w1-d-q2E ), 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 .\par The above three types are locally asymptotically stable for different ranges of k 1 .For k2 ?2 < k 1 < min\{ r2?1k2 r2-q1E , (d+q2E)p1 c1w1-d-q2E \}, we can observe the bistability of E 1 and E 2 (Fig.1(b)(c)). If conditions (d+q2E)p1 c1w1-d-q2E < k 1 < min\{ r2?1k2 r2-q1E , (d+q2E)p1+c1w1p1 c1w1-d-q2E\par \} and ( r2(k2-?2 N )k2 < w2 Z p2+?\par N ) hold simultaneous, we can find the bistability of E 2 and E 4 (Fig. \hyperref[fig_0]{1(d)(e)}). On the contrary, if(d+q2E)p1 c1w1-d-q2E < k 1 < r2?1k2 r2-q1E holds, for either k 1 > (d+q2E)(p1)+c1w1p1 c1w1-d-q2E or r2(k2-?2 N ) k2 > w2 Z p2+? N , we'll get the existence of stable E 2 together with unstable E 4 . Identically, for max\{ r2?1k2 r2-q1E , (d+q2E)p1 c1w1-d-q2E \} < k 1 < min\{ k2 ?2 , (d+q2E)p1+c1w1p1 c1w1-d-q2E\par \} together with? 1 ? 2 < 1, c1w1 N p1+ N -d -q 2 E < c2w2 T p2+ T +? N and r2(k2-?2 N ) k2 < w2 Z p2+?\par N , we can observe the bistability of E 3 and E 4 (Fig. \hyperref[fig_0]{1(f)-(i)}). Now, let's discuss the importance of avoiding toxic species by zooplankton (?) together with the harvesting term (q 1 E, q 2 E) for the survival of the different species groups.\par Firstly, let's discuss the effect of ? 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 ?. However, for the stability of E 3 and E 4 , it is related to the critical value of ?. When ? is less than this critical value, E 3 and E 4 remain stable. Thus, ? does not affect the bistability of (E 1 , E 2 ); when ? is below some threshold value, we will observe the bistability of (E 2 , E 4 ) and (E 3 , E 4 ), and as the ? value increases, the original bistability may disappear.( r2(k2-?2 N ) k2 > w2 Z p2+? N , c1w1 N p1+ N -d -q 2 E < c2w2 T p2+ T +? N and r2(k2-?2 N ) k2 < w2 Z p2+?\par N . From these conditions, we can see the establishment of the above conclusion.)\par 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 ?, 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 ?, the existence of q 1 E and q 2 E will also affect the threshold value of ?, 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.\par In this section, we focus on the local stability and Hopf bifurcation of the delayed model; the delayed system (2.2) has the following formdU (t) dt = F (U (t), U (t -? 1 ), U (t -? 2 )),\textbf{(4.1)}\par where U (t) = [N (t), T (t), Z(t)], U (t -? 1 ) = [N (t -? 1 ), T (t -? 1 ), Z(t -? 1 )], U (t -? 2 ) = [N (t -? 2 ), T (t -? 2 ), Z(t -? 2 )]. \section[{Notes}]{Notes}\par Next, assuming ?1 (t) = N (t) -N * , ? 2 (t) = T (t) -T * , ? 3 (t) = Z(t)\par -Z * at the positive equilibrium point, and linearizing the system (2.2), we can obtaind dt ? ? ? 1 (t) ? 2 (t) ? 3 (t) ? ? = L ? ? N (t) T (t) Z(t) ? ? + M ? ? N (t -? 1 ) T (t -? 1 ) Z(t -? 1 ) ? ? + S ? ? N (t -? 2 ) T (t -? 2 ) Z(t -? 2 ) ? ? ,\textbf{(4.2)}\par whereL = ?F ?U (t) E * , M = ?F ?U (t -? 1 ) E * , S = ?F ?U (t -? 2 ) E * .\par We linearize the system(2.2) about positive equilibrium E * = (N * , T * , Z * ), and getdU (t) dt = LU (t) + M U (t -? 1 ) + SU (t -? 2 ),\textbf{(4.3)}\par Fig. \hyperref[fig_0]{1}:Notes where L = ? ? l 11 l 12 l 13 l 21 l 22 l 23 0 0 l 33 ? ? , M = ? ? 0 0 0 0 0 0 m 31 0 m 33 ? ? , S = ? ? 0 0 0 0 0 0 s 31 s 32 s 33 ? ? , U = ? ? ? ? N 1 (?) T 1 (?) Z 1 (?) ? ? ? ? ,\par where N 1 , T 1 , Z 1 are small perturbations around the equilibrium point E * = (N * , T * , Z * ). We havel 11 = -rN k 1 + w 1 ZN (p 1 + N ) 2 -q 1 E, l 12 = r 1 ? 1 N k 1 , l 13 = - w 1 N p 1 + N , l 21 = r 2 ? 2 T k 1 + w 2 ?T Z (p 2 + T + ?N ) 2 , l 22 = r 2 - (2r 2 T + r 2 ? 1 N ) k 2 , l 23 = - w 2 T (p 2 + T + ?N ) , l 33 = -d -q 2 E, m 31 = c 1 w 1 p 1 Z (p 1 + N ) 2 , m 33 = c 1 w 1 N (p 1 + N ) , s 31 = c 2 w 2 ?T Z (p 2 + T + ?N ) 2 , s 32 = c 2 w 2 Z(p 2 + ?N ) (p 2 + T + ?N ) 2 , s 33 = c 2 w 2 T (p 2 + T + ?N ) .\par The characteristic equation for the linearized system (2.2) is obtained asD(?, ? 1 , ? 2 ) ? P (?) + Q(?)e -??1 + R(?)e -??2 = 0,\textbf{(4.4)}\par where Case (1):P (?) = ? 3 + A 2 ? 2 + A 1 ? + A 0 , Q(?) = B 2 ? 2 + B 1 ? + B 0 , R(?) = C 2 ? 2 + C 1 ? + C 0 ,? 1 = ? 2 = 0.\par In this case, Section 3 covers the analysis of the system when ? 1 = ? 2 = 0.\par Case (2):? 1 = 0, ? 2 > 0.\par In this case, the characteristic equation(4.4) becomesD(?, ? 2 ) ? P (?) + Q(?) + R(?)e -??2 ? ? 3 + A 2 ? 2 + A 1 ? + A 0 + B 2 ? 2 + B 1 ? + B 0 + (C 2 ? 2 + C 1 ? + C 0 )e -??2 = 0,\textbf{(4.5)}\par putting ? = i?(? > 0) in Eq.(4.5), and separating the real and imaginary parts, we have-(A 2 + B 2 )? 2 + (A 0 + B 0 ) = (C 2 ? 2 -C 0 ) cos(?? 2 ) -C 1 ? sin(?? 2 ), -? 3 + (A 1 + B 1 )? = (C 0 -C 2 ? 2 ) sin(?? 2 ) -C 1 ? cos(?? 2 ). (4.6)\par Squaring and adding the equation(4.6), we obtain[-(A 2 + B 2 )? 2 + (A 0 + B 0 )] 2 + [-? 3 + (A 1 + B 1 )?] 2 = (C 2 ? 2 -C 0 ) 2 + (C 1 ?) 2 . (4.7)\par Simplifying Eq.(4.7) and substituting ? 2 = , the above equation can be written as?( ) ? 3 + a 2 2 + a 1 + a 0 = 0,\textbf{(4.8}\par ) wherea 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 ) -2C 0 C 2 -C 2 1 , a 0 = -C 2 0 .\par (H1):a 2 > 0, a 0 > 0, a 2 a 1 -a 0 > 0.\par If (H1) holds, Eq.(4.8) has no positive roots, which implies all the roots of Eq.(4.5) have negative real parts.\par Therefore, E * is asymptotically stable for all ? 2 > 0 when (H1) holds.\par (H2): a 2 < 0, a 1 < 0, a 0 < 0 or a 2 > 0, a 1 < 0, a 0 < 0 or a 2 > 0, a 1 > 0, a 0 < 0. If (H2) holds, Eq.(4.8) has exactly one positive root ? 0 , substituting ? 0 in Eq.(4.6), we obtain-(A 2 + B 2 )? 0 2 + (A 0 + B 0 ) = (C 2 ? 0 2 -C 0 ) cos(? 0 ? 2 ) -C 1 ? 0 sin(? 0 ? 2 ), -? 0 3 + (A 1 + B 1 )? 0 = (C 0 -C 2 ? 0 2 ) sin(? 0 ? 2 ) -C 1 ? 0 cos(? 0 ? 2 ). (4.9)\par For the critical value of ? 2 , we can obtain? 2j = 1 ? 0 arccos \{ [C 1 + C 2 (A 2 + B 2 )]? 0 4 + [C 1 (A 1 + B 1 )-C 0 (A 2 + B 2 )-C 2 (A 0 + B 0 )]? 0 2 + C 0 (A 0 + B 0 ) -(C 0 -C 2 ? 0 2 ) 2 -(C 1 ? 0 ) 2 \}+ 2j? ? 0 , j = 0, 1, 2 ? ? ? . (4.10)\par For the transversality condition, differentiating Eq.(4.5) with respect to ? 2 , we getd? d? 2 = ?(C 2 ? 2 + C 1 ? + C 0 )e -??2 3? 2 + 2A 2 ? + A 1 + (2B 2 ? + B 1 ) + (2C 2 ? + C 1 )e -??2 .\par Solving ( d? d?2 ) -1 , we obtain( d? d? 2 ) -1 = 3? 2 + 2A 2 ? + A 1 + (2B 2 ? + B 1 ) + (2C 2 ? + C 1 )e -??2 ?(C 2 ? 2 + C 1 ? + C 0 )e -??2 .\par Then at ? 2 = ? 20 and ? = i? 0 , we can get [Re( d? d? 2 ) ?2=?20,?=i?0 ] -1 = Re[ 3(i? 0 ) 2 + (2A 2 + B 2 )(i? 0 ) + A 1 + B 1 (i? 0 )(C 2 (i? 0 ) 2 + C 1 (i? 0 ) + C 0 )(cos(? 0 ? 20 ) -i sin(? 0 ? 20 )) ] + Re[ 2C 2 (i? 0 ) + C 1 (i? 0 )(C 2 (i? 0 ) 2 + C 1 (i? 0 ) + C 0 ) ]. Now [Re( d? d? 2 ) ?2=?20,?=i?0 ] -1 = Re[ M R + M I i N R + N I i ] + Re[ Q R + Q I i P R + P I i ] = M R N R + M I N I N R 2 + N I 2 + Q R P R + Q I P I P R 2 + P IM R = -3? 0 2 + A 1 + B 1 , M I = 2(A 2 + B 2 )? 0 , N R = (C 0 ? 0 -C 2 ? 0 3 ) sin(? 0 ? 20 ) -C 1 ? 0 2 cos(? 0 ? 20 ), N I = (C 0 ? 0 -C 2 ? 0 3 ) cos(? 0 ? 20 ) + C 1 ? 0 2 sin(? 0 ? 20 ), Q R = C 1 , Q I = 2C 2 ? 0 , P R = -C 1 ? 0 2 , P I = C 0 ? 0 -C 2 ? 0 3 . Then [Re( d? d? 2 ) ?2=?20,?=i?0 ] -1 = A B + C D = AD + BC BD ,\textbf{(4.11)}\par hereA = M R N R + M I N I , B = N R 2 + N I 2 , C = Q R P R + Q I P I , D = P R 2 + P I 2 .\par From this, we can getsgn[Re( d? d? 2 ) ?2=?20,?=i?0 ] -1 = sgn[AD + BC].\par If (H3): AD + BC = 0 holds, the transversal condition sgn[Re( d? d?2 ) ?2=?20,?=i?0 ] -1 = 0. From the above analysis, the following theorem can be drawn For ? 1 = 0 and ? 2 > 0, we have the following results: (i)If (H1) holds, then the equilibrium E * is asymptotically stable for all ? 2 > 0.\par (ii)If (H3) holds, and (H2) holds, then the equilibrium E * is locally asymptotically stable for all ? 2 < ? 20 together with unstable for ? 2 > ? 20 and undergoes Hopf bifurcation at ? 2 = ? 20 . \section[{Case (3):}]{Case (3):}? 1 > 0, ? 2 = 0.\par In this case, the characteristic equation(4.4) becomes as followsD(?, ? 1 ) ? P (?) + R(?) + Q(?)e -??1 ? ? 3 + A 2 ? 2 + A 1 ? + A 0 + (B 2 ? 2 + (C 2 ? 2 + C 1 ? + C 0 ) + B 1 ? + B 0 )e -??1 = 0. (\textbf{4.12)}\par putting ? = i?(? > 0) in Eq.(4.12), and separating the real and imaginary parts, we have-(A 2 + C 2 )? 2 + (A 0 + C 0 ) = (B 2 ? 2 -B 0 ) cos(?? 1 ) -B 1 ? sin(?? 1 ), -? 3 + (A 1 + C 1 )? = (B 0 -B 2 ? 2 ) sin(?? 1 ) -B 1 ? cos(?? 1 ). (4.13)\par Squaring and adding the equation(4.13), we obtain[-(A 2 + C 2 )? 2 + (A 0 + C 0 )] 2 + [-? 3 + (A 1 + C 1 )?] 2 = (B 2 ? 2 -B 0 ) 2 + (B 1 ?) 2 . (4.14)\par Based on the calculation method for case (2), we can simplify (4.14) to the following?( ) ? 3 + b 2 2 + b 1 + b 0 = 0,\textbf{(4.15)}\par whereb 2 = -(A 2 + C 2 ) 2 -2(A 1 + C 1 ) -B 2 2 , b 1 = (A 1 + C 1 ) 2 -2(A 0 + C 0 )(A 2 + C 2 ) -2B 0 B 2 -B 2 1 , b 0 = -B 2 0 .\par Theorem 4.1. \section[{Notes}]{Notes}(H4): b 2 > 0, b 0 > 0, b 2 b 1 -b 0 > 0.\par 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 ? 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.\par If (H5) holds, Eq.(4.15) has exactly one positive root ?0 , substituting ?0 in Eq.(4.13), we obtain-(A 2 + C 2 ) ?0 2 + (A 0 + C 0 ) = (B 2 ?0 2 -B 0 ) cos( ?0 ? 1 ) -B 1 ?0 sin( ?0 ? 1 ), -?0 3 + (A 1 + C 1 ) ?0 = (B 0 -B 2\textbf{?0}\par 2 ) sin( ?0 ? 1 ) -B 1 ?0 cos( ?0 ? 1 ). (4.16)\par For the critical value of ? 1 , we can obtain? 1j = 1 ?0 arccos\{ [B 1 + B 2 (A 2 + C 2 )] ?0 4 + [B 1 (A 1 + C 1 )-C 0 (A 2 + C 2 )-B 2 (A 0 + C 0 )] ?0 2 +B 0 (A 0 + C 0 ) -(B 0 -B 2 ?0 2 ) 2 -(B 1 ?0 ) 2 \}+ 2j? ?0 , j = 0, 1, 2 ? ? ? . (4.17)\par For the transversality condition, differentiating Eq.(4.13) with respect to ? 1 , we getd? d? 1 = ?(B 2 ? 2 + B 1 ? + B 0 )e -??1 3? 2 + 2A 2 ? + A 1 + (2C 2 ? + C 1 ) + (2B 2 ? + B 1 )e -??1 .\par Solving ( d? d?1 ) -1 , we obtain( d? d? 1 ) -1 = 3? 2 + 2A 2 ? + A 1 + (2C 2 ? + C 1 ) + (2B 2 ? + B 1 )e -??1 ?(B 2 ? 2 + B 1 ? + B 0 )e -??1 .\par Then at ? 1 = ? 10 and ? = i ?0 , we can get[Re( d? d? 1 ) ?1=?10,?=i ?0 ] -1 = Re[ 3(i ?0 ) 2 + (2A 2 + C 2 )(i ?0 ) + A 1 + C 1 (i ?0 )(B 2 (i ?0 ) 2 + B 1 (i ?0 ) + B 0 )(cos( ?0 ? 10 ) -i sin( ?0 ? 10 )) ] + Re[ 2B 2 (i ?0 ) + B 1 (i ?0 )(B 2 (i ?0 ) 2 + B 1 (i ?0 ) + B 0 ) ]. Now [Re( d? d? 1 ) ?1=?10,?=i ?0 ] -1 = Re[ M R + M I i N R + N I i ] + Re[ Q R + Q I i P R + P I i ] = M R N R + M I N I N R 2 + N I 2 + Q R P R + Q I P I P R 2 + P I 2 ,\par whereM R = -3 ?0 2 + A 1 + C 1 , M I = 2(A 2 + C 2 ) ?0 , N R = (B 0 ?0 -B 2\textbf{?0}\par 3 ) sin( ?0 ? 10 ) -C 1 ?0 2 cos( ?0 ? 10 ),N I = (B 0 ?0 -B 2 ?0 3 ) cos( ?0 ? 10 ) + B 1 ?0 2 sin( ?0 ? 10 ), Q R = B 1 , Q I = 2B 2 ?0 , P R = -B 1 ?0 2 , P I = B 0 ?0 -B 2 ?0 3 . Then [Re( d? d? 1 ) ?1=?10,?=i? ] -1 = A * B * + C * D * = A * D * + B * C * B * D * ,\textbf{(4.18)}\par Balancing Coexistence: Ecological Dynamics and Optimal Tax Policies in a Dual Phytoplankton-Zooplankton System Influenced by Toxin Avoidance and Harvesting \section[{Notes}]{Notes}\par hereA * = M R N R + M I N I , B * = N R 2 + N I 2 , C * = Q R P R + Q I P I , D * = P R 2 + P I 2 .\par From this, we can get[Re( d? d? 1 ) ?1=?10,?=i? ] -1 = sgn[A * D * + B * C * ].\par If (H6): A * D * + B * C * = 0 holds, the transversal condition [Re( d? d?1 ) ?1=?10,?=i? ] -1 = 0. From the above analysis, the following theorem can be drawn For ? 2 = 0 and ? 1 > 0, we have the following results: (i)If (H4) holds, then the equilibrium E * is asymptotically stable for all ? 1 > 0.\par (ii)If (H6) and (H5) hold, then the equilibrium E * is locally asymptotically stable for all ? 1 < ? 10 together with unstable for ? 1 > ? 10 and undergoes Hopf bifurcation at ? 1 = ? 10 . ? 1 is fixed in (0, ? 10 ] and ? 2 > 0. We consider the gestation delay ? 1 to be stable in the interval (0, ? 10 ], taking ? 2 as a control parameter. Let ? = u + i? be the root of Eq. {\ref (4.4)}. Putting this value in Eq.(4.4), separating real and imaginary parts, we obtain u 3 -3u? 2 + A 2 (u 2 -? 2 ) + A 1 u + A 0 + (B 2 u 2 -B 2 ? 2 + B 1 u + B 0 )e -u?1 cos(?? 1 +(2B 2 u? + B 1 ?)e -u?1 sin(?? 1 ) + (C 2 u 2 -C 2 ? 2 + C 1 u + C 0 )e -u?1 cos(?? 2 ) + (2C 2 u? + C 1 ?) sin(?? 2 ) = 0. (4.19) 3u 2 ? -? 3 + 2A 2 u? + A 1 ? -(B 2 u 2 -B 2 ? 2 + B 1 u + B 0 ) sin(?? 1 ) + (2B 2 u? +B 1 ?)e -u?1 cos(?? 1 ) -(C 2 u 2 -C 2 ? 2 + C 1 u + C 0 ) sin(?? 2 ) + (2C 2 u? +C 1 ?)e -u?2 cos(?? 2 ) = 0.A 2 ? 2 -A 0 = (-B 2 ? 2 + B 0 ) cos(?? 1 ) + (C 0 -C 2 ? 2 ) cos(?? 2 ) + B 1 ? sin(?? 1 ) + C 1 ? sin(?? 2 ).\par (4.21) ? 3 -A 1 ? = -(B 0 -B 2 ? 2 ) sin(?? 1 ) + B 1 ? cos(?? 1 ) -(C 0 -C 2 ? 2 ) sin(?? 2 ) + C 1 ? cos(?? 2 ). (\textbf{4}? 6 + ã4 ? 4 + ã3 ? 3 + ã2 ? 2 + ã0 = 0, (4.23) where ã4 = -(B 2 2 + C 2 2 -A 2 2 ), ã3 = -2(B 2 C 1 -B 1 C 2 ) sin(?? 1 -?? 2 ), ã2 = -((B 1 2 -2B 0 B 2 + C 1 2 -2C 0 C 2 ) + 2(B 1 C 1 -2A 0 A 2 -A 2 1 -B 2 )) cos(?? 1 -?? 2 ), ã0 = -(B 0 2 + C 0 2 -A 0 2 ). \section[{Notes}]{Notes}\par -? 2 cos(?? 2 ) + ? 1 sin(?? 2 ) = ? 6 -? 5 cos(?? 1 ) + ? 4 sin(?? 1 ), {\ref (4.25)} where? 1 = C 2 ? 2 -C 0 , ? 2 = -C 1 ?, ? 3 = A 0 -A 2 ? 2 , ? 4 = B 0 -B 2 ? 2 , ? 5 = B 1 ?, ? 6 = ? 3 -A 1 ?.\par Without losing generality, the Eq.( \hyperref[formula_37]{4}.23) has finite positive roots ? 1 , ? 2 , ? ? ? , ? k , for every fixed ?, there exists a sequence \{? j 2i |j = 0, 1, 2...\}, where? (j) 2i = 1 ?i tan -1 [ (? 1 ? 4 + ? 2 ? 4 ) sin(? i ? 1 ) -(? 1 ? 5 -? 2 ? 4 ) cos(? i ? 1 ) + ? 1 ? 6 + ? 2 ? 3 (? 1 ? 5 -? 2 ? 4 ) sin(? i ? 1 ) + (? 2 ? 5 + ? 1 ? 4 ) cos(? i ? 1 ) + ? 1 ? 3 -? 2 ? 4 + k? ?i j = 0, 1, 2, ? ? ? (4.26) let ? 2 = min\{? (j) 2i |i = 0, 1, 2, ...k, j = 0, 1, 2...\}, when ? 2 = ? 2 , ? = ? i | ?2= ?2 , i = 1, 2, 3, .\par .., the characteristic equation (4.4) has purely imaginary roots ±i ?. Then, we will verify the transversality condition, differentiating the characteristic equation (4.4) with respect to ? 2 , we can obtain[Re( d? d? 2 ) ?2= ?2,?=i ? ] -1 = Re[ 3(i ?) 2 + 2A 2 (i ?) + A 1 (i ?)(C 2 (i ?) 2 + C 1 (i ?) + C 0 )(cos( ? ? 2 ) -i sin( ? ? 2 )) ] + Re[ 2C 2 (i ?) + C 1 (i ?)(C 2 (i ?) 2 + C 1 (i ?) + C 0 ) ]. Now [Re( d? d? 2 ) ?2= ?2,?=i ? ] -1 = Re[ M R + M I i N R + N I i ] + Re[ Q R + Q I i P R + P I i ] = M R N R + M I N I N R 2 + N I 2 + Q R P R + Q I P I P R 2 + P I 2 ,\par whereM R = -3 ? 2 + A 1 , M I = 2A 2 ?, N R = (C 0 ? -C 1 ? 2 -C 2 ? 3 ) sin( ?? 2 )\par ,N I = (C 0 ? -C 2 ? 3 ) cos( ?? 2 ) + C 1 ? 2 sin( ?? 2 ), Q R = C 1 , Q I = 2C 2 ?, P R = -C 1 ? 2 , P I = C 0 ? -C 2 ? 3 . Then [Re( d? d? 2 ) ?2= ?2,?=i ? ] -1 = E F + G H = EH + F G F H ,\textbf{(4.27)}\par here E = M R N R + M I N I , F = N R 2 + N I 2 , G = Q R P R + Q I P I , H = P R 2 + P I 2 . \section[{Notes}]{Notes}\par For system(2.2), assume (H7) holds with ? 1 is fixed in (0, ? 10 ] and ? 2 > 0, then the equilibrium E * is locally asymptotically stable for ? 2 ? (0, ? 2 ) whereas system (2.2) undergoes Hopf bifurcation at ? 2 = ? 2 .\par Case(5): ? 2 is fixed in (0, ? 20 ] and ? 1 > 0, so take ? 1 as a control parameter; the analysis is the same as case(4), so we omit it.\par 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 {\ref (2.3)}. Here, we take E as a time-dependent dynamic variable controlled by equations. Therefore, there is the following equation.E(t) = ?Q(t), 0 ? ? ? 1, dQ dt = I(t) -?Q(t), Q(0) = Q 0 .\textbf{(5.1)}\par 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 ? 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 ? = 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 ? = 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 1 N + (u 2 -v)q 2 N -C],\par 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.\par We assume the gross profit margin on capital investment is proportional to this 'Net income.' So, we haveI = E?[(u 1 -v)q 1 N + (u 2 -v)q 2 Z -C], 0 ? ? < 1.\par (5.2)\par For ? = 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. ? = 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, ? > 0. By combining Eqs.(5.1) and (5.2), we can getdE dt = E\{??[(u 1 -v)q 1 N + (u 2 -v)q 2 Z -C] -?\}.\textbf{(5.3)}\par From this we can getsgn[Re( d? d? 2 ) ?2= ?2,?=i ? ] -1 = sgn[EH + F G].\par If (H7): EH + F G = 0 holds, the transversal condition sgn[Re( d? d?2 ) ?2= ?2,?=i ? ] -1 = 0. From the above analysis, we have the following theorem.\par Theorem 4.3. \section[{V.}]{V.}\par Optimal Tax Policy \section[{Notes}]{Notes}\par 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 isM E = E[(u 1 -v)q 1 N + (u 2 -v)q 2 Z -C] + E[v(q 1 N ) + v(q 2 N )],\par this is equal to the net economic income of fishermen plus the economic income of regulators. Therefore without considering the time delay, Eq.( \hyperref[formula_0]{2}.3) can be rewritten as? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? dN dt = r 1 N (1 - N + ? 1 T k 1 ) - w 1 N Z p 1 + N -q 1 EN, dT dt = r 2 N (1 - T + ? 2 N k 2 ) - w 2 T Z p 2 + T + ?N , dZ dt = c 1 w 1 N Z p 1 + N - c 2 w 2 T Z p 2 + T + ?N -dZ -q 2 EZ, dE dt = E\{??[(u 1 -v)q 1 N + (u 2 -v)q 2 Z -C] -?\}.\par (5.4)\par 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 = +? 0 E(t)e -?t [u 1 q 1 N + u 2 q 2 Z -C]dt,\textbf{(5.5)}\par where ? is the discounting factor. Therefore, our goal is to determine an optimal tax v = v(t) that maximizes compliance with Eq.(5.4) and constrains v min ? v(t) ? v max on the control variable v(t). When v min < 0, it will have the effect of accelerating the rate of fishery expansion. The Hamiltonian of the problem is obtained byH = (u 1 q 1 N + u 2 q 2 Z -C)Ee -?t + ? 1 N [r 1 (1 - N + ? 1 T k 1 ) - w 1 Z p 1 + N -q 1 E] +? 2 [r 2 T (1 -T +?1N k2 ) -w2T Z p2+T +?N ] + ? 3 [ c1w1N Z p1+N -c2w2T Z p2+T +?N -dZ -q 2 EZ] +? 4 E\{??[(u 1 -v)q 1 N + (u 2 -v)q 2 Z -C] -?\},\textbf{(5.6)}\par where ? 1 , ? 2 , ? 3 and ? 4 are the adjoint variables. For v ? [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 {\ref [9]} ?H ?v = -? 4 E??(q 1 N + q 2 Z) = 0 ? ? 4 = 0.\par (5.7)\par Now, the adjoint equations ared? 1 dt = - ?H ?N = -[u 1 q 1 Ee -?t +? 1 (r 1 - 2r 1 N + r 1 ? 1 T k 1 - w 1 p 1 Z (p 1 + N ) 2 -q 1 E) + ? 2 [ w 2 ?T Z (p 2 + T + ?N ) 2 - r 2 ? 2 T k 2 ] + ? 3 ( c 1 w 1 p 1 Z (p 1 + N ) 2 + c 2 w 2 ?T Z (p 2 + T + ?N ) 2 )\par ,d? 2 dt = - ?H ?T = -[? 1 ( r 1 ? 1 N k 1 ) + ? 2 [r 2 (1 - 2T + ? 2 N k 2 ) - w 2 Z(p 2 + ?N ) (p 2 + T + ?N ) 2 ] -? 3 ( c 2 w 2 Z(p 2 + ?N ) (p 2 + T + ?N ) 2 ), d? 3 dt = - ?H ?Z = -[u 2 q 2 Ee -?t -? 1 ( w 1 N p 1 + N ) -? 2 ( w 2 T p 2 + T + ?N ) + ? 3 ( c 1 w 1 N p 1 + N - c 2 w 2 T p 2 + T + ?N -d -q 2 E)], d? 4 dt = - ?H ?E = -[(u 1 q 1 N + u 2 q 2 Z -C)e -?t -? 1 q 1 N -? 3 q 2 Z].\textbf{(5.8)}\par Now start with Eqs.(5.8) and (5.7), using the equilibrium equation we haved? 1 dt =-u 1 q 1 Ee -?t -? 1 [- r 1 N k 1 + w 1 N Z (p 1 + N ) 2 ] -? 2 [ w 2 ?T Z (p 2 + T + ?N ) 2 - r 2 ? 2 T k 2 ] -? 3 [ c 1 w 1 p 1 Z (p 1 + N ) 2 + c 2 w 2 ?T Z (p 2 + T + ?N ) 2 ]\par ,d? 2 dt = -? 1 [ r 1 ? 1 N k 1 ] -? 2 [ w 2 T Z (p 2 + T + ?N ) 2 ] -? 3 [ c 2 w 2 Z(p 2 + ?N ) (p 2 + T + ?N ) 2 ]\par ,d? 3 dt = -u 2 q 2 Ee -?t + ? 1 ( w 1 N p 1 + N ) + ? 2 ( w 2 T p 2 + T + ?N\par ).\par (5.9)\par Using the second and third equations of Equation (5.9) from the fourth equation of Equation (5.8), we can obtaind?1 dt = M 1 e -?t + M 2 ? 1 + M 3 ? 2 ,\par whereM 1 = (C -u 1 q 1 N )? + u 2 q 2 Z(q 2 E -?) q 1 N , M 2 = - w 1 q 2 N Z (p 1 + N )q 1 N , M 3 = - w 2 q 2 T Z (p 2 + T + ?N )q 1 N .\par The solution of this linear equation is? 1 = N 0 e -M2t - M 1 e -?t M 2 + ? - M 3 ? 2 M 2 .\par (5.10)\par Using the same method as above, we can get? 3 = I 0 e H2t - H 1 e -?t H 2 + ? ,\textbf{(5.11)}\par whereH 1 = [ (C -u 2 q 2 Z)? -q 1 N (u 1 ? + M 1 ) q 2 Z + M 1 M 2 q 1 N (M 2 + ?)q 2 Z ], H 2 = M 2 M 3 q 1 N q 2 M 2 Z .\par Identicallyd? 2 dt = R 1 e -?t + R 2 ? 2 ,\textbf{(5.12)}\par whereR 1 = M 1 M 2 + ? + H 1 H 2 + ? ( c 2 w 2 Z(p 2 + ?N ) (p 2 + T + ?N ) 2 ), R 2 = M 3 M 2 ( r 2 ? 1 N k 1 ) - w 2 T Z (p 2 + T + ?N ) 2 .\par So we can get ? 1? 1 = N 0 e M2t - M 1 e -?t M 2 + ? - M 3 (W 0 e R2t -R1e -?t R2+? ) M 2 .\par The shadow price ? 1 e -?t is bounded as t ? ?, N 0 = 0 and W 0 = 0, then we can obtain? 1 = - M 1 e -?t M 2 + ? - M 3 M 2 (e R2t - R 1 e -?t R 2 + ? ).\par (5.13) Now use Eqs.(5.11), (5.12) and (5.13) in the first of Eq.(5.9), we have[ (C -u 1 q 1 N * )? + u 2 q 2 Z * (q 2 E * -?) q 1 N * ]e -?t + w 2 q 2 N * Z * (p 1 + N * )q 1 N * [ M 1 e -?t M 2 + ? - M 3 M 2 (e R2t - R 1 e -?t R 2 + ? )] +[ w2q2T * Z * (p2+T * +?N * )q1N * ][ R1e -?t R2+? ] + u 1 q 1 E * e -?t + [ M1e -?t M2+? -M3 M2 (e R2t -R1e -?t R2+? )][-r1N * k1 + w1N * Z * (p1+N * ) 2 ] = ( R1e -?t R2+? )[ w2?T * Z * (p2+T * +?N * ) 2 -r2?2T * k2 ] + ( H1e -?t H2+? )[ c2w2Z * (p2+?N * ) (p2+T * +?N * ) 2 ].\par (5.14)\par Because of the computational complexity, our optimal equilibrium solution can be expressed asT * = [(c 1 w 1 -?)N * -?p 1 ](p 2 + ?N * ) [(c 2 w 2 -?)p 1 + (c 2 w 2 -c 1 w 2 -?)N * ] , Z * = r 1 ( p1+N * w1k1 )(k 1 -N * -? 1 T * ).\textbf{(5.15)}\par N * available from the following equationr 2 (k 2 -T * -? 2 N * )(p 2 + T * + ?N * ) -w 2 k 2 Z * = 0.\par (5.16)\par E * available from the following equationr 1 q 1 (1 - N * + ? 1 T * k 1 ) - w 1 Z * q 1 (p 1 + N * ) = c 1 w 1 N * q 2 (p 1 + N * ) - c 2 w 2 T * q 2 (p 2 + T * + ?N * ) - d q 2 .\par (5.17)\par 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[ (C -u 1 q 1 N * )? + u 2 q 2 Z * (q 2 E * -?) q 1 N * ]e -?t + w 2 q 2 N * Z * (p 1 + N * )q 1 N * [ M 1 e -?t M 2 + ? - M 3 M 2 (e R2t - R 1 e -?t R 2 + ? )] +[ w2q2T * Z * (p2+T * +?N * )q1N * ][ R1e -?t R2+? ] + u 1 q 1 E * e -?t + [ M1e -?t M2+? -M3 M2 (e R2t -R1e -?t R2+? )][-r1N * k1 + w1N * Z * (p1+N * ) 2 ] -( R1e -?t R2+? )[ w2?T * Z * (p2+T * +?N * ) 2 -r2?2T * k2 ] -( H1e -?t H2+? )[ c2w2Z * (p2+?N * ) (p2+T * +?N * ) 2 ] = f (v).\par (5.18)\par 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 \hyperref[b24]{[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 \section[{VI. Numerical Simulations}]{VI. Numerical Simulations} \section[{Notes}]{Notes}\par In Fig. \hyperref[fig_9]{3}, we plot the time series of ? = 0, ? = 10, ? = 1000 in the first ten days, where the other parameter values and initial conditions are the same as in Table \hyperref[tab_2]{2}. When q 1 = q 2 = 0 and ? = 0, we can observe that NTP and TPP tend to perish at a fast linear speed. It is obvious that when ? 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. \hyperref[fig_9]{3(a)(b)}). On the contrary, when ? = 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. \hyperref[fig_9]{3(c)(d)}). Based on the values of q 1 and q 2 of (fig. \hyperref[fig_9]{3(c)(d})), we increase ? 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. \hyperref[fig_9]{3}(e)(f)). To further explore the influence of ?, we fixed q 1 and q 2 as 2 and 2.5, respectively. And increased the value of ? from 10 to 1000. At this time, We can observe that the concentration of NTP, TPP and zooplankton has almost no change(fig. \hyperref[fig_9]{3(g)(h})). Finally, when ? 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. \hyperref[fig_9]{3}(i)(j)).\par In Fig. \hyperref[fig_11]{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. \hyperref[fig_11]{4}(a)(b), we can observe the dynamic change of ? from 0 to 10. First, we take ? = 0, in fig. \hyperref[fig_11]{4}(a), we will find the extinction of TPP, while NTP and zooplankton oscillate in the form of limit cycles. Next, we increase ? to 10, observe the fig. \hyperref[fig_11]{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 ? = 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. \hyperref[fig_11]{4(c)(d)}). Finally, when we fix ? = 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. \hyperref[fig_11]{4}(e)(f)). Now, to explore the influence of pregnancy delay (? 1 ) and toxin onset delay(? 2 ) on the stability of equilibrium point in different cases. First, we need to set a set of parameters as followsr 1 = 2, r 2 = 3, ? 1 = 0.3, ? 2 = 0.1, k 1 = 2500, k 2 = 3000, w 1 = w 2 = 0.5, p 1 = p 2 = 50, c 1 = c 2 = 0.45, d = 0.05, ? = 0.5, q 1 = 0.2, q 2 = 0.3, E = 1. (\textbf{6 1)}\par 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.\par Case i: when ? 1 = 0, ? 2 > 0, in this case, [Re( d? d?2 ) ?2=?20,?=i?0 ] -1 > 0, the transversality condition is contented. So when ? 2 < ? 20 (Fig. {\ref 5(a)(b)}), the positive equilibrium E * is locally asymptotically stable, the positive equilibrium E * is unstable when ? 2 > ? 20 (Fig. {\ref 6}(a)(b)), when ? 2 = ? 20 , the system undergoes Hopf bifurcation around the positive equilibrium E * . (Fig. {\ref 5(a)(b})) shows the trajectories and phase portrait of system (2.2) for ? 1 = 0, ? 2 = 1. It can be clearly seen that the system (2.2) will converge to the positive equilibrium point E * . And (Fig. {\ref 6(a)(b})) shows the trajectories and phase portrait of the system (2.2) for ? 1 = 0, ? 2 = 1.08. In this case, the delay system (2.2) has a periodic solution near the positive equilibrium point (E * ).\par Case ii : when ? 1 > 0, ? 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. [Re( d? d?1 ) ?1=?10,?=i ?0 ] -1 > 0, the transversality condition is satisfied. (Fig. {\ref 7}(a)(b)) shows the trajectories and phase portrait of the system (2.2) for ? 1 = 0.7, ? 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.\par Case iii : when ? 1 = 0.9 in stable interval (0, ? 10 ), and take ? 2 > 0 as the parameter, [Re( d? d?2 ) ?2= ?2,?=i ? ] -1 = 0, the transversality condition is satisfied. So when shows the trajectories and phase portrait of the system (2.2) for ? 1 = 0.9, ? 2 = 1.06. It can be clearly seen that the system (2.2) will converge to the positive equilibrium point E * . And (Fig. {\ref 9}(a)(b)) shows the trajectories and phase portrait of the system (2.2) for ? 1 = 0.9, ? 2 = 1.09; we find the delayed system (2.2) has periodic solutions near the positive equilibrium point E * in this case.\par 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 ? 1 > ? 10 , ? 2 > ? 20 and ? 2 > ? 2 .\par The dynamic changes of the system ( 1 ) with different ?, q 1 and q 2 in the first 10 days, other parameter values and initial conditions are the same as Table \hyperref[tab_2]{2}. (a)(b) : In the case of q 1 = q 2 = 0, ? = 0 and ? = 10, the TPP concentration will fluctuate and the NTP concentration will barely change. (c)(d) : For ? = 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 ? = 10, TPP reached a higher flowering concentration, while NTP still maintained a lower concentration. (g)(h) : Based on (f), for ? = 1000, NTP and TPP concentrations are almost unchanged. (i)(j): for ? = 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. \section[{Notes}]{Notes}\par The long-term dynamics of the system (2.1), all other parameter values are the same as Table \hyperref[tab_2]{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 ? = 10, all populations survive and the system stabilizes to a limit cycle. (c)(d) : For ? = 0, 0 ? q 1 =q 2 ? 0.36, NTP and TPP can coexist. (e)(f): when we fix ? = 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. \section[{Notes}]{Notes}\par The behavior of the system(2.2) for ? 1 = 0,? 2 = 1 with other parameters chosen in (6.1).\par The behavior of the system(2.2) for ? 1 = 0,? 2 = 1.08 with other parameters chosen in (6.1).\par The behavior of the system(2.2) for ? 1 = 0.7,? 2 = 0 with other parameters chosen in (6.1). \section[{Notes}]{Notes}\par 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 \hyperref[tab_2]{2} and (6.1) to simulate several cases of double-delay systems in Matlab to verify all theoretical results.\begin{figure}[htbp] \noindent\textbf{1}\includegraphics[]{image-2.png} \caption{\label{fig_0}1}\end{figure} \begin{figure}[htbp] \noindent\textbf{}\includegraphics[]{image-3.png} \caption{\label{fig_1}}\end{figure} \begin{figure}[htbp] \noindent\textbf{2332}\includegraphics[]{image-4.png} \caption{\label{fig_3}with A 2 = 33 B 2 =}\end{figure} \begin{figure}[htbp] \noindent\textbf{}\includegraphics[]{image-5.png} \caption{\label{fig_4}}\end{figure} \begin{figure}[htbp] \noindent\textbf{20}\includegraphics[]{image-6.png} \caption{\label{fig_5}(4. 20 )}\end{figure} \begin{figure}[htbp] \noindent\textbf{}\includegraphics[]{image-7.png} \caption{\label{fig_6}}\end{figure} \begin{figure}[htbp] \noindent\textbf{2}\includegraphics[]{image-8.png} \caption{\label{fig_7}Fig. 2 :}\end{figure} \begin{figure}[htbp] \noindent\textbf{2}\includegraphics[]{image-9.png} \caption{\label{fig_8}? 2 < ? 2}\end{figure} \begin{figure}[htbp] \noindent\textbf{3}\includegraphics[]{image-10.png} \caption{\label{fig_9}Fig. 3 :}\end{figure} \begin{figure}[htbp] \noindent\textbf{}\includegraphics[]{image-11.png} \caption{\label{fig_10}}\end{figure} \begin{figure}[htbp] \noindent\textbf{4}\includegraphics[]{image-12.png} \caption{\label{fig_11}Fig. 4 :}\end{figure} \begin{figure}[htbp] \noindent\textbf{56789}\includegraphics[]{image-13.png} \caption{\label{fig_12}Fig. 5 :Fig. 6 :Fig. 7 :Fig. 8 :Fig. 9 :}\end{figure} \begin{figure}[htbp] \noindent\textbf{}\includegraphics[]{image-14.png} \caption{\label{figure14}}\end{figure} \begin{figure}[htbp] \noindent\textbf{}\includegraphics[]{image-15.png} \caption{\label{figure15}}\end{figure} \begin{figure}[htbp] \noindent\textbf{}\includegraphics[]{image-16.png} \caption{\label{figure16}}\end{figure} \begin{figure}[htbp] \noindent\textbf{}\includegraphics[]{image-17.png} \caption{\label{figure17}}\end{figure} \begin{figure}[htbp] \noindent\textbf{2} \par \begin{longtable}{} \end{longtable} \par \caption{\label{tab_2}Table 2 :}\end{figure} \backmatter \par Ref\par ), 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\par + /\{0\} be the optimal equilibrium. Now, we seek min S 0 (v) subject to the solution to Eq. {\ref (5.5)}.\par Using the maximum principle, we will get the adjoint variables ? 1 , ? 2 , ? 3 and ? 4 as\par (5.20)\par The adjoint variables ? 1 , ? Eq.(5.19) specifies a set of initial conditions for ? 1 , ? 2 , ? 3 and ? 4 , and Eq.(5.20) uses these initial conditions to determine the unique solution of ? 1 , ? 2 , ? 3 and ? 4 . Therefore, it is easy to obtain the optimal tax as follows:\par (5.22)\par 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:\par From the above analysis, we can easily observe the following points:\par (i) From Eqs.(5.7) and (5.11)-(5.13), we note that ? i e -?t , (i = 1, 2, 3, 4), where ? 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 ? ?.\par (ii) Considering the coexistence equilibrium point\par The fourth equation of Eq.(5.8) can be written as\par 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.\par (iii) From Eqs.(5.11) and (5.13), we can obtain \subsection[{Notes}]{Notes}\par The optimal solution of (5.5) forv = 0.867.\par 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.\par 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\par and the discounting factor ? = 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. {\ref 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.\par 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 {\ref 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 ?, 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. \subsection[{Notes}]{Notes}\par The behavior of the system(2.2) for ? 1 = 0.9,? 2 = 1.06 with other parameters chosen in (6.1).\par The behavior of the system(2.2) for ? 1 = 0.9,? 2 = 1.09 with other parameters chosen in (6.1).\par 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.\par 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. {\ref 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 ? 1 and the toxin onset delay ? 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(? 0 1 , ? 0 2 , ? * 10 and ? * 20 , respectively), but \begin{bibitemlist}{1} \bibitem[Dubey and Hussain ()]{b19}\label{b19} ‘A model for the allelopathic effect on two competing species’. B Dubey , J Hussain . \textit{Ecol. Model} 2000. 129 (23) p. . \bibitem[Samanta ()]{b20}\label{b20} ‘A stochastic two species competition model: Nonequilibrium fluctuation and stability’. G P Samanta . \textit{Int. J. Stoch. Anal} 2011. 2011 p. . \bibitem[Mitra and Flynn ()]{b2}\label{b2} ‘Accounting for variation in prey selectivity by zooplankton’. A Mitra , K J Flynn . \textit{Ecol. Model} 2006. 199 p. . \bibitem[Linhart et al. ()]{b3}\label{b3} ‘Avoidance of prey by captive coyotes punished with electric shock’. S B Linhart , J D Roberts , S A Shumake , R Johnson . \textit{Proceedings of the Vertebrate Pest Conference, escholarship}, (the Vertebrate Pest Conference, escholarship) 1976. p. . \bibitem[Sharma et al. ()]{b29}\label{b29} ‘Bifurcation behaviors analysis of a plankton model with multiple delays’. A K Sharma , A Sharma , K Agnihotri . \textit{Int. J. Biomath} 2016. 9 (06) p. 1650086. \bibitem[Xiao and Jennings ()]{b27}\label{b27} ‘Bifurcations of a ratio-dependent predatorprey system with constant rate harvesting’. D Xiao , L Jennings . \textit{Siam. Appl. Math} 2005. 65 p. . \bibitem[Clark ()]{b24}\label{b24} C W Clark . \textit{Bioeconomic Modelling and Fisheries Management}, (New York) 1985. USA) Wiley. \bibitem[Kar ()]{b33}\label{b33} ‘Conservation of a fishery through optimal taxation: a dynamic reaction model’. T Kar . \textit{Commun. Nonlinear Sci. Numer. Simul} 2005. 10 p. . \bibitem[Panday and Samanta ()]{b32}\label{b32} ‘Delay induced multiple stability switch and chaos in a predatorprey model with fear effect’. P Panday , S Samanta , N . \textit{Math. Comput. Simulation} 2020. 172 p. . \bibitem[Zhao et al. ()]{b9}\label{b9} ‘Dynamic behavior analysis of phytoplanktonzooplankton system with cell size and time delay’. Q Y Zhao , S T Liu , D D Tian . \textit{Chaos. Solitons. Fractals} 2018. 113 p. . \bibitem[Frost ()]{b8}\label{b8} ‘Effect of size and density of food particles on the feeding behaviour of the marine planktonic copepod Calanuspacificus’. B W Frost . \textit{Limnol. Oceanogr} 1972. 17 p. . \bibitem[Chattopadhyay ()]{b14}\label{b14} ‘Effect of toxic substances on a two-species competitive system’. J Chattopadhyay . \textit{Ecol. Model} 1996. 84 (13) p. . \bibitem[Pei et al. ()]{b22}\label{b22} ‘Evolutionary consequences of harvesting for a two-zooplankton one-phytoplankton system’. Y Pei , Y Lv , C Li . \textit{Appl. Math. Model} 2012. 36 (4) p. . \bibitem[Uye and Takamatsu ()]{b11}\label{b11} ‘Feeding interactions between planktonic copepods and red-tide flagellates from Japanese coastal waters’. S Uye , K Takamatsu . \textit{Mar. Ecol. Prog. Ser} 1990. 59 p. . \bibitem[Sotomayor ()]{b35}\label{b35} ‘Generic bifurcations of dynamical systems’. J Sotomayor . \textit{Dynamical systems}, 1973. Academic Press. \bibitem[Lv et al. ()]{b28}\label{b28} ‘Harvesting of a phytoplanktonzooplankton model’. Y Lv , Y Pei , S Gao , C Li . \textit{Nonlinear Anal: Real World Appl} 2010. 11 p. . \bibitem[Liu et al. ()]{b5}\label{b5} ‘Microzooplankton grazing and selective feeding during bloom periods in the Tolo Harbour area as revealed by HPLC pigment analysis’. X J Liu , C H Tang , C K Wong . \textit{J. Sea Res} 2014. 90 p. . \bibitem[Jang et al. ()]{b18}\label{b18} ‘Nutrient-phytoplankton-zooplankton models with a toxin’. S J Jang , J Baglama , J Rick . \textit{Math. Comput. Modelling} 2006. 43 (12) p. . \bibitem[Chakraborty et al. ()]{b23}\label{b23} ‘Optimal control of effort of a stage structured preypredator fishery model with harvesting’. K Chakraborty , S Das , T Kar . \textit{Nonlinear Anal. RWA} 2011. 12 (6) p. . \bibitem[Agnihotri and Kaur ()]{b13}\label{b13} ‘Optimal control of harvesting effort in a phytoplanktonzooplankton model with infected zooplankton under the influence of toxicity’. K Agnihotri , H Kaur . \textit{Math. Comput. Simul} 2021. 19 p. . \bibitem[Sarkar ()]{b34}\label{b34} ‘Optimal fishery harvesting rules under uncertainty’. S Sarkar . \textit{Resour. Energy Econ} 2009. 31 p. . \bibitem[Panja et al. ()]{b25}\label{b25} P Panja , S K Mondal , D K Jana . \textit{Effects of toxicants on Phytoplankton-Zooplankton-Fish dynamics and harvesting}, 2017. 104 p. . \bibitem[Ghorai et al. ()]{b4}\label{b4} ‘Preferential selection of zooplankton and emergence of spatiotemporal patterns in plankton population’. S Ghorai , B Chakraborty , N Bairagi . \textit{Chaos. Soliton. Fract} 2021. 153 p. 111471. \bibitem[Ref References Références Referencias]{b1}\label{b1} \textit{Ref References Références Referencias}, \bibitem[Mondal et al. ()]{b31}\label{b31} ‘Rich dynamics of non-toxic phytoplankton, toxic phytoplankton and zooplankton system with multiple gestation delays’. A Mondal , A K Pal , G P Samanta . \textit{Int. J. Dyn. Control} 2020. 8 p. . \bibitem[Khare et al. ()]{b16}\label{b16} ‘Role of toxin producing phytoplankton on a plankton ecosystem’. S Khare , O P Misra , J Dhar . \textit{Nonlinear Anal} 2010. 4 (3) p. . \bibitem[Porter ()]{b7}\label{b7} ‘Selective grazing and differential digestion of algae by zooplankton’. K G Porter . \textit{Nature} 1973. 244 p. . \bibitem[Zheng et al. ()]{b6}\label{b6} ‘Selective grazing of zooplankton on phytoplankton defines rapid algal succession and blooms in oceans’. Y L Zheng , X Gong , H W Gao . \textit{Ecol.l Modelling} 2022. 468 p. 109947. \bibitem[Zilverberg et al. ()]{b0}\label{b0} ‘Sensitivity of diet choices and environmental outcomes to a selective grazing algorithm’. C J Zilverberg , J Angerer , J Williams , L J Metz , K Harmoney . \textit{Ecol. Model} 2018. 390 p. . \bibitem[Zilverberg et al. ()]{b10}\label{b10} ‘Sensitivity of diet choices and environmental outcomes to a selective grazing algorithm’. C J Zilverberg , J Angerer , J Williams , L J Metz , K Harmoney . \textit{Ecol. Model} 2018. 390 p. . \bibitem[Zhang et al. ()]{b26}\label{b26} ‘Spatiotemporal pattern selection in a nontoxic-phytoplankton-toxicphytoplankton -zooplankton model with toxin avoidance effects’. F F Zhang , J M Sun , W Tian . \textit{App. Math. Comput} 2022. 423 p. 127007. \bibitem[Agnihotri and Kaur ()]{b30}\label{b30} ‘The dynamics of viral infection in toxin producing phytoplankton and zooplankton system with time delay’. K Agnihotri , H Kaur . \textit{Chaos Solitons. Fractals} 2019. 118 p. . \bibitem[Chakraborty et al. ()]{b15}\label{b15} ‘The role of avoidance by zooplankton for survival and dominance of toxic phytoplankton’. S Chakraborty , S Bhattacharya , U Feudel , J Chattopadhyay . \textit{Ecol. Complexity} 2012. 11 p. . \bibitem[Sole et al. ()]{b12}\label{b12} ‘The role of selective predation in harmful algal blooms’. J Sole , E Garcia-Ladona , M Estrada . \textit{J. Mar. Syst} 2006. 62 p. . \bibitem[Turriff et al. ()]{b17}\label{b17} ‘Toxin accumulation and feeding behaviour of the planktoniccopepod Calanus jinmarchicus exposed to the red -tide dinoflagellate Alexandrium excavatum’. N Turriff , J A Runge , A D Cembella . \textit{Mar. Biol} 1995. 123 p. . \bibitem[Roy and Chattopadhyay ()]{b21}\label{b21} ‘Toxin-allelopathy among phytoplankton species prevents competitive exclusion’. S Roy , J Chattopadhyay . \textit{J. Biol. Syst} 2007. 15 (01) p. . \end{bibitemlist} \end{document}