# 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 mechanisms of selective grazing include prey morphology (size, color, shape, and colony formation), intestinal genetic strains, and poisonous chemicals released by prey [6][7][8][9][10][11][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][14][15][16][17][18][19][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][22][23][24][25][26][27]. Since time delay is widely studied in the phytoplankton-zooplankton model [28][29][30][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.

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):
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 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,(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 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,(2.2)
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.(2.3)
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.

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:
dX dt = F (X),(3.1)
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 ? ? ? ? ? ? ? .

# Notes

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 have
N (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},
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). (3.2)
Then, by differentiating (3.2) concerning t, we obtain
dW 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
, therefore, applying a theorem on differential inequalities [34], we obtain
0 ? W ? ? ? + W (N (0),T (0),Z(0)) e ?t
, let t ? +?, W (N, T, Z) ? ? ? . So, all solutions of system (2.3) enter the region
D = {(N, T, Z) ? R 3 + : 0 ? W (N, T, Z) ? ? ? }.(3.3)
This shows that every solution of the system is bounded. 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 = ( N , T , 0), where
N = ? 1 k 2 -k 1 ? 1 ? 2 -1 - q 1 k 1 E r 1 , T = ? 2 k 1 -k 2 ? 1 ? 2 -1 ;
(v)TPP-free equilibrium E 4 = ( N , 0, Z), where
N = (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 ;
(vi)the interior equilibrium, E * = (N * , T * , Z * ), where
T * = 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 ;
and N * can be obtained from
r 2 (p 2 + T * + ?N * )(k 2 -T * -? 2 N * ) -w 2 k 2 Z * = 0. (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 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 is
k 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 * )
> 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 2 E. Therefore, it is a saddle point and hence always unstable.

(ii) The eigenvalues of the TPP and zooplankton-free equilibrium
E 1 = (k 1 , 0, 0) are -r 1 -q 1 E, r 2 (1 -k1?2 k2 ) and c1w1k1 p1+k1 -d -q 2 E. When c 1 ? 1 -d -q 2 E
? 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 equilibrium
E 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 .
(iv) The eigenvalues of the zooplankton-free equilibrium
E 3 = ( N , T , 0) are c1w1 N p1+ N -c2w2 T p2+ T +? N -d -q 2 E
, ? 1 and ? 2 , where ? 1 and ? 2 are the roots of the equation
? 2 + b1 ? + c1 = 0,(1)
b) Equilibrium points and their stability

Notes
where 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 +?
N , b1 > 0 and c1 > 0. If the above conditions are satisfied, E 3 is LAS.

(v) The eigenvalues of the TPP-free equilibrium
E 4 = ( N , 0, Z) are r 2 (1 -?2 N k2 ) -w2 Z p2+?
N , ?1 and ?2 , where ?1 and ?2 are the roots of the equation
? 2 -(ã 2 + b2 )? + ã2 b2 + c2 = 0,(2)
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+?
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.

(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.(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,
where
D 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 ).
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).

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 ? 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 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 :
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 [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.


# Ref

(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.

(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 )
, 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).

(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
), 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 )
, 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).

(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.

Existence and stability conditions of the equilibrium points.


# Equilibrium Existence conditions

Stability conditions E 0 = (0, 0, 0) Always exist Always unstable
E 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:

# Notes

The existence and stability of these equilibrium points are summarized in Table 1 andFig 
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
), 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 > 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 .

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
} and ( r2(k2-?2 N )
k2 < w2 Z p2+?
N ) hold simultaneous, we can find the bistability of E 2 and E 4 (Fig. 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
} 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+?
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 (?) 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 ? 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+?
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 ?, 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.

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
dU (t) dt = F (U (t), U (t -? 1 ), U (t -? 2 )),(4.1)
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 )].

# Notes

Next, assuming ?
1 (t) = N (t) -N * , ? 2 (t) = T (t) -T * , ? 3 (t) = Z(t)
-Z * at the positive equilibrium point, and linearizing the system (2.2), we can obtain
d 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 ) ? ? ,(4.2)
where
L = ?F ?U (t) E * , M = ?F ?U (t -? 1 ) E * , S = ?F ?U (t -? 2 ) E * .
We linearize the system(2.2) about positive equilibrium E * = (N * , T * , Z * ), and get
dU (t) dt = LU (t) + M U (t -? 1 ) + SU (t -? 2 ),(4.3)
Fig. 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 (?) ? ? ? ? ,
where N 1 , T 1 , Z 1 are small perturbations around the equilibrium point E * = (N * , T * , Z * ). We have
l 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 ) .
The characteristic equation for the linearized system (2.2) is obtained as
D(?, ? 1 , ? 2 ) ? P (?) + Q(?)e -??1 + R(?)e -??2 = 0,(4.4)
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.
In this case, Section 3 covers the analysis of the system when ? 1 = ? 2 = 0.

Case (2):
? 1 = 0, ? 2 > 0.
In this case, the characteristic equation(4.4) becomes
D(?, ? 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,(4.5)
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)
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)
Simplifying Eq.(4.7) and substituting ? 2 = , the above equation can be written as
?( ) ? 3 + a 2 2 + a 1 + a 0 = 0,(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 ) -2C 0 C 2 -C 2 1 , a 0 = -C 2 0 .
(H1):
a 2 > 0, a 0 > 0, a 2 a 1 -a 0 > 0.
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 ? 2 > 0 when (H1) holds.

(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)
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)
For the transversality condition, differentiating Eq.(4.5) with respect to ? 2 , we get
d? 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 .
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 .
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 ,(4.11)
here
A = 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 .
From this, we can get
sgn[Re( d? d? 2 ) ?2=?20,?=i?0 ] -1 = sgn[AD + BC].
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.

(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 .


# Case (3):
? 1 > 0, ? 2 = 0.
In this case, the characteristic equation(4.4) becomes as follows
D(?, ? 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. (4.12)
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)
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)
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,(4.15)
where
b 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 .
Theorem 4.1.


# Notes
(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 ? 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 ?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?0
2 ) sin( ?0 ? 1 ) -B 1 ?0 cos( ?0 ? 1 ). (4.16)

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)
For the transversality condition, differentiating Eq.(4.13) with respect to ? 1 , we get
d? 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 .
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 .
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 ,
where
M R = -3 ?0 2 + A 1 + C 1 , M I = 2(A 2 + C 2 ) ?0 , N R = (B 0 ?0 -B 2?0
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 * ,(4.18)
Balancing Coexistence: Ecological Dynamics and Optimal Tax Policies in a Dual Phytoplankton-Zooplankton System Influenced by Toxin Avoidance and Harvesting


# Notes

here
A * = 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 .
From this, we can get
[Re( d? d? 1 ) ?1=?10,?=i? ] -1 = sgn[A * D * + B * C * ].
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.

(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. (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 ).
(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 ). (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 ).

# Notes

-? 2 cos(?? 2 ) + ? 1 sin(?? 2 ) = ? 6 -? 5 cos(?? 1 ) + ? 4 sin(?? 1 ), (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 ?.
Without losing generality, the Eq.( 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, .
.., 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 ,
where
M R = -3 ? 2 + A 1 , M I = 2A 2 ?, N R = (C 0 ? -C 1 ? 2 -C 2 ? 3 ) sin( ?? 2 )
,
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 ,(4.27)
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 .

# Notes

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 .

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.

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) = ?Q(t), 0 ? ? ? 1, dQ dt = I(t) -?Q(t), Q(0) = Q 0 .(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 ? 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],
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?[(u 1 -v)q 1 N + (u 2 -v)q 2 Z -C], 0 ? ? < 1.
(5.2)

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 get
dE dt = E{??[(u 1 -v)q 1 N + (u 2 -v)q 2 Z -C] -?}.(5.3)
From this we can get
sgn[Re( d? d? 2 ) ?2= ?2,?=i ? ] -1 = sgn[EH + F G].
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.

Theorem 4.3.


# V.

Optimal Tax Policy


# Notes

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[(u 1 -v)q 1 N + (u 2 -v)q 2 Z -C] + E[v(q 1 N ) + v(q 2 N )],
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
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 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] -?}.
(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 = +? 0 E(t)e -?t [u 1 q 1 N + u 2 q 2 Z -C]dt,(5.5)
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 by
H = (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] -?},(5.6)
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 [9] ?H ?v = -? 4 E??(q 1 N + q 2 Z) = 0 ? ? 4 = 0.

(5.7)

Now, the adjoint equations are
d? 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 )
,
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].(5.8)
Now start with Eqs.(5.8) and (5.7), using the equilibrium equation we have
d? 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 ]
,
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 ]
,
d? 3 dt = -u 2 q 2 Ee -?t + ? 1 ( w 1 N p 1 + N ) + ? 2 ( w 2 T p 2 + T + ?N
).

(5.9)

Using the second and third equations of Equation (5.9) from the fourth equation of Equation (5.8), we can obtain
d?1 dt = M 1 e -?t + M 2 ? 1 + M 3 ? 2 ,
where
M 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 .
The solution of this linear equation is
? 1 = N 0 e -M2t - M 1 e -?t M 2 + ? - M 3 ? 2 M 2 .
(5.10)

Using the same method as above, we can get
? 3 = I 0 e H2t - H 1 e -?t H 2 + ? ,(5.11)
where
H 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 .
Identically
d? 2 dt = R 1 e -?t + R 2 ? 2 ,(5.12)
where
R 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 .
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 .
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 + ? ).
(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 ].
(5.14)

Because of the computational complexity, our optimal equilibrium solution can be expressed as
T * = [(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 * ).(5.15)
N * available from the following equation
r 2 (k 2 -T * -? 2 N * )(p 2 + T * + ?N * ) -w 2 k 2 Z * = 0.
(5.16)

E * available from the following equation
r 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 .
(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
[ (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).
(5.18)

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  


# VI. Numerical Simulations


# Notes

In Fig. 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 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. 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. 3(c)(d)). Based on the values of q 1 and q 2 of (fig. 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. 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. 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. 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 ? from 0 to 10. First, we take ? = 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 ? 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 ? = 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 ? = 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)). 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 follows
r 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. (6 1)
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 ? 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. 5(a)(b)), the positive equilibrium E * is locally asymptotically stable, the positive equilibrium E * is unstable when ? 2 > ? 20 (Fig. 6(a)(b)), when ? 2 = ? 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 ? 1 = 0, ? 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 ? 1 = 0, ? 2 = 1.08. In this case, the delay system (2.2) has a periodic solution near the positive equilibrium point (E * ).

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. 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.

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. 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.

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 .

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 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. 


# Notes

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 ? = 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.  


# Notes

The behavior of the system(2.2) for ? 1 = 0,? 2 = 1 with other parameters chosen in (6.1).

The behavior of the system(2.2) for ? 1 = 0,? 2 = 1.08 with other parameters chosen in (6.1).

The behavior of the system(2.2) for ? 1 = 0.7,? 2 = 0 with other parameters chosen in (6.1). 


# Notes

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![Year 2023 III. Dynamical Behavior of Non-Delayed Model a) Positivity and boundedness of the solution Lemma 3.1.](image-2.png "1")
![Ecological Dynamics and Optimal Tax Policies in a Dual Phytoplankton-Zooplankton System Influenced by Toxin Avoidance and Harvesting](image-3.png "")
2332![-(l 33 + l 22 -l 11 ), A 1 = l 11 l 22 + l 11 l 33 + l 22 l 33 -l 12 l 21 , A 0 = -l 11 l 22 l 33 + l 12 l 21 l -m 33 , B 1 = -l 11 m 33 -l 22 m 33 -l 13 m 31 , B 0 = +l 13 l 22 m 31 + l 11 l 22 m 33 + l 12 l 21 m 33 -l 12 l 23 m 31 , C 2 = -s 33 , C 1 = -l 13 s 31 + l 11 s 33 -l 22 l 23 s 32 -l 22 s 33 , C 0 = l 11 s 33 + l 11 l 23 s 32 + l 12 l 21 s 33 + l 13 l 22 s 31 -l 12 l 23 s 31 -l 13 l 21 s 32 .](image-4.png "with A 2 = 33 B 2 =")
![Ecological Dynamics and Optimal Tax Policies in a Dual Phytoplankton-Zooplankton System Influenced by Toxin Avoidance and Harvesting Notes](image-5.png "")
20![Putting u = 0 in Eqs.(4.19) and (4.20), we obtain](image-6.png "(4. 20 )")
![Ecological Dynamics and Optimal Tax Policies in a Dual Phytoplankton-Zooplankton System Influenced by Toxin Avoidance and Harvesting](image-7.png "")
2![Fig. 2:](image-8.png "Fig. 2 :")
2![Balancing Coexistence: Ecological Dynamics and Optimal Tax Policies in a Dual Phytoplankton-Zooplankton System Influenced by Toxin Avoidance and Harvesting a) Time series analysis b) Double time delay analysis 8(a)(b)), the positive equilibrium E * is locally asymptotically stable, the positive equilibrium E * is unstable when ? 2 > ? 2 (Fig.9(a)(b)), when ? 2 = ? 2 , the system undergoes Hopf bifurcation around the positive equilibrium E * . (Fig.8(a)(b))](image-9.png "? 2 < ? 2")
3![Fig. 3:](image-10.png "Fig. 3 :")
![Ecological Dynamics and Optimal Tax Policies in a Dual Phytoplankton-Zooplankton System Influenced by Toxin Avoidance and Harvesting](image-11.png "")
4![Fig. 4:](image-12.png "Fig. 4 :")
56789![Fig. 5:](image-13.png "Fig. 5 :Fig. 6 :Fig. 7 :Fig. 8 :Fig. 9 :")
![](image-14.png "")
![](image-15.png "")
![](image-16.png "")
![](image-17.png "")
2
		
		
			

Ref

), 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

+ /{0} be the optimal equilibrium. Now, we seek min S 0 (v) subject to the solution to Eq. (5.5).

Using the maximum principle, we will get the adjoint variables ? 1 , ? 2 , ? 3 and ? 4 as

(5.20)

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:

(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:

From the above analysis, we can easily observe the following points:

(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 ? ?.

(ii) Considering the coexistence equilibrium point

The fourth equation of Eq.(5.8) can be written as

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


## Notes

The optimal solution of (5.5) forv = 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

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. 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.

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 ?, 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. 


## Notes

The behavior of the system(2.2) for ? 1 = 0.9,? 2 = 1.06 with other parameters chosen in (6.1).

The behavior of the system(2.2) for ? 1 = 0.9,? 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 ? 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 
			
			

				


* 
	
		Sensitivity of diet choices and environmental outcomes to a selective grazing algorithm
		
			CJZilverberg
		
		
			JAngerer
		
		
			JWilliams
		
		
			LJMetz
		
		
			KHarmoney
		
	
	
		Ecol. Model
		
			390
			
			2018
		
	




* 
	
		Ref References Références Referencias
		
	




* 
	
		Accounting for variation in prey selectivity by zooplankton
		
			AMitra
		
		
			KJFlynn
		
	
	
		Ecol. Model
		
			199
			
			2006
		
	




* 
	
		Avoidance of prey by captive coyotes punished with electric shock
		
			SBLinhart
		
		
			JDRoberts
		
		
			SAShumake
		
		
			RJohnson
		
	
	
		Proceedings of the Vertebrate Pest Conference, escholarship
		the Vertebrate Pest Conference, escholarship
		
			1976
			
		
	




* 
	
		Preferential selection of zooplankton and emergence of spatiotemporal patterns in plankton population
		
			SGhorai
		
		
			BChakraborty
		
		
			NBairagi
		
	
	
		Chaos. Soliton. Fract
		
			153
			111471
			2021
		
	




* 
	
		Microzooplankton grazing and selective feeding during bloom periods in the Tolo Harbour area as revealed by HPLC pigment analysis
		
			XJLiu
		
		
			CHTang
		
		
			CKWong
		
	
	
		J. Sea Res
		
			90
			
			2014
		
	




* 
	
		Selective grazing of zooplankton on phytoplankton defines rapid algal succession and blooms in oceans
		
			YLZheng
		
		
			XGong
		
		
			HWGao
		
	
	
		Ecol.l Modelling
		
			468
			109947
			2022
		
	




* 
	
		Selective grazing and differential digestion of algae by zooplankton
		
			KGPorter
		
	
	
		Nature
		
			244
			
			1973
		
	




* 
	
		Effect of size and density of food particles on the feeding behaviour of the marine planktonic copepod Calanuspacificus
		
			BWFrost
		
	
	
		Limnol. Oceanogr
		
			17
			
			1972
		
	




* 
	
		Dynamic behavior analysis of phytoplanktonzooplankton system with cell size and time delay
		
			QYZhao
		
		
			STLiu
		
		
			DDTian
		
	
	
		Chaos. Solitons. Fractals
		
			113
			
			2018
		
	




* 
	
		Sensitivity of diet choices and environmental outcomes to a selective grazing algorithm
		
			CJZilverberg
		
		
			JAngerer
		
		
			JWilliams
		
		
			LJMetz
		
		
			KHarmoney
		
	
	
		Ecol. Model
		
			390
			
			2018
		
	




* 
	
		Feeding interactions between planktonic copepods and red-tide flagellates from Japanese coastal waters
		
			SUye
		
		
			KTakamatsu
		
	
	
		Mar. Ecol. Prog. Ser
		
			59
			
			1990
		
	




* 
	
		The role of selective predation in harmful algal blooms
		
			JSole
		
		
			EGarcia-Ladona
		
		
			MEstrada
		
	
	
		J. Mar. Syst
		
			62
			
			2006
		
	




* 
	
		Optimal control of harvesting effort in a phytoplanktonzooplankton model with infected zooplankton under the influence of toxicity
		
			KAgnihotri
		
		
			HKaur
		
	
	
		Math. Comput. Simul
		
			19
			
			2021
		
	




* 
	
		Effect of toxic substances on a two-species competitive system
		
			JChattopadhyay
		
	
	
		Ecol. Model
		
			84
			13
			
			1996
		
	




* 
	
		The role of avoidance by zooplankton for survival and dominance of toxic phytoplankton
		
			SChakraborty
		
		
			SBhattacharya
		
		
			UFeudel
		
		
			JChattopadhyay
		
	
	
		Ecol. Complexity
		
			11
			
			2012
		
	




* 
	
		Role of toxin producing phytoplankton on a plankton ecosystem
		
			SKhare
		
		
			OPMisra
		
		
			JDhar
		
	
	
		Nonlinear Anal
		
			4
			3
			
			2010
		
	




* 
	
		Toxin accumulation and feeding behaviour of the planktoniccopepod Calanus jinmarchicus exposed to the red -tide dinoflagellate Alexandrium excavatum
		
			NTurriff
		
		
			JARunge
		
		
			ADCembella
		
	
	
		Mar. Biol
		
			123
			
			1995
		
	




* 
	
		Nutrient-phytoplankton-zooplankton models with a toxin
		
			SJJang
		
		
			JBaglama
		
		
			JRick
		
	
	
		Math. Comput. Modelling
		
			43
			12
			
			2006
		
	




* 
	
		A model for the allelopathic effect on two competing species
		
			BDubey
		
		
			JHussain
		
	
	
		Ecol. Model
		
			129
			23
			
			2000
		
	




* 
	
		A stochastic two species competition model: Nonequilibrium fluctuation and stability
		
			GPSamanta
		
	
	
		Int. J. Stoch. Anal
		
			2011
			
			2011
		
	




* 
	
		Toxin-allelopathy among phytoplankton species prevents competitive exclusion
		
			SRoy
		
		
			JChattopadhyay
		
	
	
		J. Biol. Syst
		
			15
			01
			
			2007
		
	




* 
	
		Evolutionary consequences of harvesting for a two-zooplankton one-phytoplankton system
		
			YPei
		
		
			YLv
		
		
			CLi
		
	
	
		Appl. Math. Model
		
			36
			4
			
			2012
		
	




* 
	
		Optimal control of effort of a stage structured preypredator fishery model with harvesting
		
			KChakraborty
		
		
			SDas
		
		
			TKar
		
	
	
		Nonlinear Anal. RWA
		
			12
			6
			
			2011
		
	




* 
	
		
			CWClark
		
		Bioeconomic Modelling and Fisheries Management
		New York
		
			USA) Wiley
			1985
		
	




* 
	
		
			PPanja
		
		
			SKMondal
		
		
			DKJana
		
	
	
		Effects of toxicants on Phytoplankton-Zooplankton-Fish dynamics and harvesting
		
			2017
			104
			
		
	




* 
	
		Spatiotemporal pattern selection in a nontoxic-phytoplankton-toxicphytoplankton -zooplankton model with toxin avoidance effects
		
			FFZhang
		
		
			JMSun
		
		
			WTian
		
	
	
		App. Math. Comput
		
			423
			127007
			2022
		
	




* 
	
		Bifurcations of a ratio-dependent predatorprey system with constant rate harvesting
		
			DXiao
		
		
			LJennings
		
	
	
		Siam. Appl. Math
		
			65
			
			2005
		
	




* 
	
		Harvesting of a phytoplanktonzooplankton model
		
			YLv
		
		
			YPei
		
		
			SGao
		
		
			CLi
		
	
	
		Nonlinear Anal: Real World Appl
		
			11
			
			2010
		
	




* 
	
		Bifurcation behaviors analysis of a plankton model with multiple delays
		
			AKSharma
		
		
			ASharma
		
		
			KAgnihotri
		
	
	
		Int. J. Biomath
		
			9
			06
			1650086
			2016
		
	




* 
	
		The dynamics of viral infection in toxin producing phytoplankton and zooplankton system with time delay
		
			KAgnihotri
		
		
			HKaur
		
	
	
		Chaos Solitons. Fractals
		
			118
			
			2019
		
	




* 
	
		Rich dynamics of non-toxic phytoplankton, toxic phytoplankton and zooplankton system with multiple gestation delays
		
			AMondal
		
		
			AKPal
		
		
			GPSamanta
		
	
	
		Int. J. Dyn. Control
		
			8
			
			2020
		
	




* 
	
		Delay induced multiple stability switch and chaos in a predatorprey model with fear effect
		
			PPanday
		
		
			SSamanta
		
		
			N
		
	
	
		Math. Comput. Simulation
		
			172
			
			2020
		
	




* 
	
		Conservation of a fishery through optimal taxation: a dynamic reaction model
		
			TKar
		
	
	
		Commun. Nonlinear Sci. Numer. Simul
		
			10
			
			2005
		
	




* 
	
		Optimal fishery harvesting rules under uncertainty
		
			SSarkar
		
	
	
		Resour. Energy Econ
		
			31
			
			2009
		
	




* 
	
		Generic bifurcations of dynamical systems
		
			JSotomayor
		
	
	
		Dynamical systems
		
			Academic Press
			1973