BIFURCATION CONTROL OF A MINIMAL MODEL OF MARINE PLANKTON INTERACTION WITH MULTIPLE DELAYS

Plankton blooms and its control is an intriguing problem in ecology. To investigate the oscillatory nature of blooms, a two-dimensional model for plankton species is considered where one species is toxic phytoplankton and other is zooplankton. The delays required for the maturation time of zooplankton, the time for phytoplankton digestion and the time for phytoplankton cells to mature and release toxic substances are incorporated and the delayed model is analyzed for stability and bifurcation phenomena. It proves that periodic plankton blooms can occur when the delay (the sum of the above three delays) changes and crosses some threshold values. The phenomena described by this mechanism can be controlled through the toxin release rates of phytoplankton. Then, a delay feedback controller with the coefficient depending on delay is introduced to system. It concludes that the onset of the bifurcation can be delayed as negative feedback gain (the decay rate) decreases (increases). Some numerical examples are given to verify the effectiveness of the delay feedback control method and the existence of crossing curve. These results show that the oscillatory nature of blooms can be controlled by human behaviors. Mathematics Subject Classification. 34K18, 34D20, 92D25, 93B52, 93D15. Received October 10, 2020. Accepted February 13, 2021.

alternative method of improving our knowledge of the physical and biological processes relating to plankton ecology. Mathematical modeling of plankton populations went back to the work of Fleming et al. [10,17,31,34]. Subsequent researchers used the models of varying levels of complexity. These studies take advantage of various factors related to marine ecosystems and the different characteristics of plankton.
The role of toxins in the occurrence and termination of phytoplankton blooms had been discussed [5,32,36]. Zooplankton grazing also plays an important role in the initial stages of outbreaks. Some studies had shown that many copepods had reduced grazing and fecundity in the presence of HABs species. In the last few decades, the decrease of grazing pressure of zooplankton due to the release of toxic substances by phytoplankton was one of the research hotspots [4,6,7,36]. Chattopadhyay et al. [6] provided a mathematical model for TPP and zooplankton interaction and explored the role of TPP behind HABs. In [7], the authors analyzed a model based on the assumption that the toxin release process followed a discrete change, τ being the discrete time required for phytoplankton cells to mature and release toxic substances. Authors studied the oscillation of phytoplankton and zooplankton populations and the stability conditions of oscillation behavior to explain the cyclic nature of bloom dynamics. The adverse effect of HABs is clear, but the control of such problems is under investigation.
In order to describe the dynamic interaction between phytoplankton and zooplankton, many zooplankton population models were established from the improved predator-prey model [18,21,35,40]. The dynamic relationship between predator and prey, due to its ubiquity and importance, has long been and will continue to be one of the dominant themes in ecology and mathematical ecology. Pioneering work on predator-prey interactions went back to the mid-1920s [27,38]. In the study of biological phenomena, many factors affect the dynamic properties of biological and mathematical models. Functional response is one of the common nonlinear factors. The interaction between predator and prey is also another factor in population dynamics. An interesting example is when the predator's growth function is different from the prey's function. Moreover, the growth term of predator is not only described by a function of predator density, but by the ratio of predator to prey [24][25][26]. Especially, Friedman et al. [11,16] proposed a predator-prey system in the following form ẋ(t) = x(t)g(x(t), δ) − y(t)p(x(t)), y(t) = y(t)q y(t)/x(t) , (1.1) where x(t) and y(t) are the population densities of the prey and predator at time t. g(x, δ) is a continuously differentiable function that describes the specific growth rate of prey in the absence of predator satisfying g(0, δ) = r > 0, g(δ, δ) = 0, g x (δ, δ) < 0, g x (x, δ) ≤ 0, and g δ (x, δ) > 0 for all x > 0. The logistic growth with g(x, δ) = r 1 (1 − x/δ) as the prototype and satisfying all hypotheses, δ is the carrying capacity of prey. r 1 is the intrinsic growth rate of prey. Choosing q(y/x) = r 2 (1 − y/γx), it became known as the Leslie-Gower model, which measured the decrease in the number of predator due to rare foods. r 2 denotes the intrinsic growth rate of predator. γ represents the quality of food provided by a prey that translates into a predator. The functional response of predator to prey p(x) describes the change of prey density attacked by each predator with the change of prey density in unit time. It is continuously differentiable, such that p(0) = 0. Hsu and Huang [16] used Dulac criterion to construct an appropriate Lyapunov function with all biological admissible parameters when p(x) = mx, and obtained the globally asymptotic stability of the uniquely positive equilibrium of system (1.1). Beretta and Kuang [3] studied the Leslie-Gower model by considering the maturation time and ratio dependence of predator. In addition, Ma [29] studied a Leslie-Gower system with two delays, and obtained the existence of global Hopf bifurcation. In addition, in [23,30,39], authors considered a delayed prey-predator system that combined an improved Leslie-Gower functional response with Holling-II functional response. However, we know that not all the species in the natural world are suitable for Holling-II functional response, and some species are suitable for Holling-III functional response, and that some others suitable for Ivlev type or other equivalent forms. In particular, Holling-III response which is suitable for the vertebral predators exists universally in population dynamics [9,19,37], and so on. Hopf bifurcation is widely used to obtain more information about the periodic solution properties of nonlinear systems near steady state. In particular, bifurcation control is a very essential and effective method to eliminate oscillation. The delayed feedback control (DFC) approach has received extensive attention since it was proposed by Pyragas in 1992 [33]. The delay feedback controller is composed of continuous linear feedback of the difference between the current state and the delay state of the system. By this method, a controller can be designed to suppress or reduce some of the existing bifurcation dynamics of a given nonlinear system, and some satisfactory dynamical behaviors can be obtained [8,20,28,41]. In recent years, DFC method has been widely used in mechanical, electronic communication and biological fields.
The interactions of phytoplankton and zooplankton basically follow the prey-predator relations. In this article, we improve system (1.1) to describe TPP and zooplankton interactions and combine the following functional forms: g(x, δ) = r 1 (1 − x/δ), p(x) = mf (x) and q(y/x) = r 2 (1 − y/γx), that is, Leslie-Gower type system with generalized functional response and two delays is investigated: where x(t) and y(t) are the population densities of TPP and zooplankton. TPP may be either Enteromorphs linza, the unicellular green alga, Chlorella vulgaris, or Noctiluca scintillans, etc. Chlorella vulgaris produces an autotoxin which has the ability to regulate the growth of its own population and also inhibits the growth of Asterionella formosa and Nitzschia frustulum. The zooplankton population may be Daphnia species or Paracalanus. Since the earliest days of phytoplankton ecology, nutrients have been invoked as one of the variables controlling phytoplankton community structure and biomass. In system (1.2), the nutrient uptake by phytoplankton is reflected by the parameter r 1 . τ 1 denotes the maturation time of zooplankton, while τ 2 represents the time required for phytoplankton digestion and the time required for phytoplankton cells to mature and release toxic substances. ρ is the ratio of phytoplankton producing toxic substances per unit biomass. c is the half saturation constant. This term ρx(t−τ2)y(t) c+x(t−τ2) represents the death of zooplankton caused by the toxins of phytoplankton with Holling II type. In most of cases, f (x) is chosen by x (Holling I type), x/(b + x) (Holling II type), x 2 /(b + x 2 ) (Holling III type), 1 − e −ax (Ivlev type) or some equivalent forms [15]. In this paper, we give the general assumptions of the above forms for f (x): (H) f (0) = 0, f (x) > 0 for x ≥ 0 and lim x→∞ f (x) = 1. Hence, our model represents a significant generalization of those previously analyzed.
In the real world, it is sometimes necessary to control a population at a reasonable level, otherwise the population may lead to the increase, decrease or even extinction of some populations. In terms of biological system control, the current research focus is to achieve feedback control by changing the structure of the biological community and increasing the feeding pressure of prey. For example, an effective way to eliminate algal blooms is to introduce appropriate fish (whitefish, etc.) which usually feed on plankton, so that algal blooms can be controlled. Based on this, a feedback controller with time-delay correlation coefficients is designed to control the bifurcation of system. Improving the ideal of Pyragas [33], we add a time-delayed force ke −dσ [y(t) − y(t − σ)] to predator equation of system (1.2) and obtain the following system: Parameters k and σ are used as human control parameters. The parameter k represents the control intensity per unit time, and the parameter σ represents the age of control object. The case k > 0 can be regarded as reproduction rate, which means that the number of predators of age σ will be artificially increased. On the contrary, the case k < 0 is regarded as the capture rate, which means it would artificially reduce the number of predators with an age of σ. The case k = 0 means no human control. The significance of the control is that by artificially increasing or decreasing the number of predators of the σ age, we can modify the bifurcation characteristics of nonlinear system to obtain some specific dynamical behaviors. Note that the strength of feedback control is in the form of ke −dσ , and the function decreases exponentially with delay σ. This means that the feedback effects of past states diminish over time. Mathematically, the larger the d, the smaller the factor ke −dσ . Therefore, d is called the decay rate and k the feedback gain. In fact, it can combine τ 1 and τ 2 by setting u(t) = x(t − τ 2 ) and v(t) = y(t), then system (1.3) becomes where τ = τ 1 + τ 2 . By biological meaning, the initial conditions are chosen as where (ϕ 1 (θ), ϕ 2 (θ)) ∈ C([− τ , 0], R 2 + ) and τ = max{τ, σ}. The fundamental existence uniqueness theorem (see, for example, [13]) ensures that solutions of (1.4) exist uniquely for all time.
Equilibria of system (1.4) are given by the solutions of equations The zooplankton-free equilibrium is readily identified given by E 1 (δ, 0). The possible positive equilibrium and Hence, F (u) and G(u) must intersect in a 0 < u * < δ. Hence system (1.4) exists a uniquely positive equilibrium E * (u * , v * ) under the conditions (H) and ρ ≤ r 2 . For biological system, the positive equilibrium is what we're interested in. Hence, in this paper, its main goal is to study the dynamic behaviors of the positive equilibrium. We always assume that the condition (H) and ρ ≤ r 2 hold, i.e., the toxin release rate of phytoplankton is no more than the intrinsic growth rate of zooplankton. Otherwise, two populations can not coexist. It now perturbs system (1.3) around E * and obtains the following linearized system of differential equations: where The characteristic equation corresponding to the Jacobian matrix of of system (1.6) around (0, 0) is In the following, the root distribution of (1.7) is examined to determine the stability of E * . The remaining sections of this article is organized as follows. In Section 2, the conditions for the Hopf bifurcation occurrence of uncontrolled systems are derived by taking delay τ as the bifurcation parameter. In Section 3, it analyzes the Hopf bifurcation of the controlled system and demonstrates the influence of feedback on bifurcation. The stability changes of equilibrium in (σ, τ ) plane can be obtained by using the crossing curve method. In Section 4, the theoretical conclusions are verified by numerical examples and simulations. Finally, some conclusions are given.

Bifurcation analysis in the uncontrolled system
In this section, it will investigated the dynamic behaviors of uncontrolled system only including one delay τ . (2.1) It can easily obtain the positivity and boundedness of solutions for system (2.1). The set Ω The system (2.1) has the same equilibrium with system (1.4). It is easy to see that for ρ > r 2 (c + δ)/δ, the system will exhibit stable behavior around zooplankton-free equilibrium E 1 (δ, 0). Ecologically, this result indicates that for large toxin release rate of phytoplankton, the zooplankton species will become extinct. In the following, it investigates the dynamic behavior of the positive equilibrium E * of system (2.1). In this case, the characteristic equation (1.7) becomes where A 4 (0) = −r 2 v * /γu * . By Routh-Hurwitz criterion, when τ = 0, the non-delayed system exhibits stable behavior for Furthermore, when A 1 + A 4 (0) = 0 and A 1 A 4 (0) + A 2 A 3 > 0, the ODE system can occur Hopf bifurcation. Substituting iω (ω > 0) into (2.2), then ω satisfies the following equations: Squaring the both of sides of (2.3) and adding, it has ) has a uniquely positive root denoted ω + and Let λ(τ ) = β(τ ) + iω(τ ) satisfying β(τ j ) = 0 and ω(τ j ) = ω + . Derivating the two sides of (2.2) with respect to τ , it has ) has no positive root and E * is locally asymptotically stable for all τ ≥ 0; , then E * is locally asymptotically stable for τ ∈ [0, τ 0 ) and unstable for τ > τ 0 , where τ 0 is the smallest value for which there is a solution to (2.2) with real part zero. Further as τ increases through τ 0 , E * bifurcates into small amplitude of periodic solutions.

Bifurcation for the controlled system
Now let's go back to system (1.4) and examine the effects of the control term on the system. Here we rewrite system (1.4) as Next, we talk about system (3.1) in two cases: two equal delays and opposite aspect.
At the same time, the characteristic equation (1.7) becomes In the following, the existence of purely imaginary roots λ = iw (w = w(τ ) > 0) for (3.3) is investigated. The equation (3.3) is an exponential polynomial of λ, whose coefficients depend on the delay of τ . Beretta and Kuang [2] established a geometric method to determine the existence of purely imaginary roots when the coefficients contain delays. It can easily verify the following relations: then it has finite roots; (v) If w > 0 exists satisfying F(w, τ ) = 0, then it is continuous and differentiable in τ . Let iw (w = w(τ ) > 0) be the root of (3.3), then w satisfies whose roots are given by where Then for all τ ∈ I τ , w satisfies (3.7) and w is not definited if τ / ∈ I τ . For τ ∈ I τ , let θ(τ ) ∈ [0, 2π) be defined by , and the maps τ n (τ ) : Furthermore, it can introduce the continuous and differentiable functions S n (τ ): S n (τ ) := τ − τ n (τ ), τ ∈ I τ , n ∈ N.

The case σ = τ
In this part, two kinds of methods are used to investigate the stability of positive equilibrium.
Furthermore, it can compute the property of Hopf bifurcation to determine the bifurcating direction and stability of bifurcating periodic solution using the methods in [14]. Here we leave it out.

Crossing curve methods
The results stated in Theorem 3.7 clearly show the stability of equilibrium E * when σ and τ change. However, if we set σ in the unstable region, then there may exist no τ 0 such that E * is unstable in τ ∈ [0, τ 0 ), it is stable in τ > τ 0 . That is to say, no information is given on couples (σ, τ ) that generate a stable or an unstable stationary state. In order to overcome this concern, an effective approach was the one proposed by Gu et al. [12] through the use of the stability crossing curve, which were defined as the curves that separated the stable and unstable regions in the (σ, τ ) plane. However, the system in [12] did not depend on the delay in the coefficients, hence, these results could not be used straightway to (2.2). Recently, An et al. [1] improved these results in [12] to the system with the coefficients depending on the delay. In the following, we will use the methods in [1] and [22] to give the couples (σ, τ ) that generate a stable or an unstable equilbrium.
In the following, we will identify the admissible values of (σ, τ ) ∈ R + × R + such that λ = iω (ω > 0) is a zero of (3.17). The curve formed by all such points is referred to the "crossing curve".
By (C5), the endpoints (ω e ,σ i ± (ω e )) of every component of curve C must lie on the boundary of Ω and can be divided into three types: Type A: (3.27) holds for (ω, σ ω ) = (ω e ,σ i ± (ω e )). Type B: (3.28) holds for (ω, σ ω ) = (ω e ,σ i ± (ω e )). Type C: (3.29) holds for (ω, σ ω ) = (ω e ,σ i ± (ω e )). For convenience, Type AA indicates that both end points of the component are Type A. Based on the type of the end point, the each component of curve C falls into the following six categories: Type AA, Type BB, Type CC, Type AB, Type AC and Type BC. (c) A spiral-like curve spreading out along τ -axis or a series of spiral-like curves oriented along τ -axis; and each of these curves approaching ∞ in the direction of σ-axis; (d) Truncated curves of one of the above three cases.

Numerical simulation
We have gained analytical understanding of possible dynamics of the phytoplankton-zooplankton system to some extent. We now perform some numerical simulations with some hypothetical set of parameters for better understanding of our analytical findings.  it can obtain τ 0 = 1.6806 and ω 0 = 0.8334. Choosing τ = 1 < τ 0 , τ = 1.69 . = τ 0 , τ = 2 > τ 0 , τ = 3 > τ 0 , respectively, the phase plots are shown in Figure 1(A)-(D). It can furthermore confirm the global existence of periodic solution (see Fig. 2), which leaves it for future research.
Furthermore, under the parameters (4.1), it can obtain the (ρ, τ 0 ) plots (see Fig. 3(A)-(B)). From the plots, it can see that as ρ increases from 0 to 2.7104, the first critical value of Hopf bifurcation keeps increasing, that is, the bifurcation occurrence is delayed. Only in a small interval (2.7104, 2.8123), ρ leads to the first Hopf bifurcation value earlier. When ρ > 2.8123, the first critical value of Hopf bifurcation is still delayed. These show that the rate of toxic substances produced by phytoplankton ρ can delay the occurrence of Hopf bifurcation and reduce the oscillation of the system.

Example 2: the controlled system (3.2)
Choosing f (u) = u, k = −2, d = 0.1 and other parameters are the same with (4.1). By computing, it gets τ 0 = 2.2895 > 1.6806. It can obtain (τ, k) plot (see Fig. 4(A)-(B)), in which the area I is no Hopf bifurcation region while II-IV areas are the regions in which Hopf bifurcation occurs. It can see that the stable interval changes from [0, 1.6806) for system (2.1) (the red line segment in τ -axis) to [0, 2.2895) for system (3.2) and the green line segment is the increased stability interval after the control term is introduced (see Fig. 4(B)). Choosing  = τ 0 , τ = 5 > τ 0 , respectively, the phases are sown in Figure 5(A)-(D). This indicates the validity of the DFC method to postpone the occurrence of Hopf bifurcation.
From Figure 4(A), it can see that if the value of the negative feedback gain k decreases, the stability interval of the positive equilibrium will become larger, that is, reducing k is conducive to the stability of system, which can also see as shown in the following figure. Fixed d = 0.1 and choosing k = 0, −1, −2, −4, respectively, it finds that the stable interval of the positive equilibrium increases as k decreases (see Fig. 6). This shows that proper capture can help to control blooms. The curve C and crossing curves are shown in Figure 7. It can obtain the feasible region Ω for (σ, ω), which is enclosed by the curve f 0 (ω, σ) = 0, f 2 (ω, σ) = 0 and ω-axis in Figure 7(A). It can see that the components of curve C are Type AC. The curve C on Ω lead to a series of spiral-like crossing curves on (σ, τ ) plane shown in Figure 7(B). It can also calculate the crossing direction in the light of Theorem 3.12 and present the final result in Figure 7(B), that is, when (σ, τ ) changes along the direction of the arrows, the characteristic roots passes through the imaginary axis from left to right. Choosing (τ, σ) = (1, 2) and (τ, σ) = (1, 5), it can obtain the phase plots in Figure 8(A)-(B). In addition, it also can choose the other parameters to see the changes of crossing curves as curve C. Choosing the following parameters: The feasible region and curve C are shown in Figure 9(A) (under parameter (4.3)) and Figure 9(B) (under parameter (4.4)). It can obtain the feasible region Ω for (σ, ω), which is possibly enclosed by the curve f 0 (ω, σ) = 0, f 2 (ω, σ) = 0, ω-axis and σ-axis (see Fig. 9(B)). It can see that the curve C on Ω still lead to a series of spirallike crossing curves on (σ, τ ) plane shown in Figures 10(A) (under parameter (4.3)) and Figure 10

Conclusion
This paper attempts to search a suitable mechanism by which it can explain the oscillatory succession of planktonic blooms and also tried to find out a possible mechanism for controlling blooms. To explain these, we study a minimal model for phytoplankton-zooplankton interaction relationship and Hopf bifurcation control problem. It is observed that bloom phenomena described by this mechanism can be controlled by two ways, (i) through toxin producing phytoplankton (TPP) and (ii) through introduction of human control. From our analytical study we have observed that enhancing the rate of toxin production it can increase the region of stability of the system and as a result can control the periodicity of planktonic blooms. The DFC method was proposed in order to retard the Hopf bifurcation occurrence. By choosing the delay as a bifurcation parameter  and using stability analysis, we determined the conditions for Hopf bifurcation occurrence for the uncontrolled and controlled system. By selecting appropriate control parameters, the proposed controller can effectively postpone the onset of the Hopf bifurcation. This shows that DFC method is effective for delay systems for retard the Hopf bifurcation occurrence. In addition, using the crossing curve methods, it can obtain the stable changes of equilibrium in (σ, τ ) plane. Simulation results have demonstrated the correctness of the theoretical analysis. By simulations, it can also know that increasing the release rate of toxins can delay the occurrence of Hopf bifurcation. Our analytical findings and supportive numerical simulations indicate that the delay differential equation system can capture the oscillatory behavior usually exhibited by various phytoplankton species, for certain choice of parametric values. The combination of delays as well as the age of control object the responsible  factors for the oscillation in phytoplankton population densities. Before ending the article we would like to mention that all the possible mechanism existing in the planktonic world may not be captured in a simple model. Moreover, the model proposed here can be extended in different ways, for example, introduction of stochastic perturbation around the coexisting steady state, turbulent diffusion, etc.