Dynamic analysis of a new aquatic ecological model based on physical and ecological integrated control

: Within the framework of physical and ecological integrated control of cyanobacteria bloom, because the outbreak of cyanobacteria bloom can form cyanobacteria clustering phenomenon, so a new aquatic ecological model with clustering behavior is proposed to describe the dynamic relationship between cyanobacteria and potential grazers. The biggest advantage of the model is that it depicts physical spraying treatment technology into the existence pattern of cyanobacteria, then integrates the physical and ecological integrated control with the aggregation of cyanobacteria. Mathematical theory works mainly investigate some key threshold conditions to induce Transcritical bifurcation and Hopf bifurcation of the model (2 . 1), which can force cyanobacteria and potential grazers to form steady-state coexistence mode and periodic oscillation coexistence mode respectively. Numerical simulation works not only explore the influence of clustering on the dynamic relationship between cyanobacteria and potential grazers, but also dynamically show the evolution process of Transcritical bifurcation and Hopf bifurcation, which can be clearly seen that the density of cyanobacteria decreases gradually with the evolution of bifurcation dynamics. Furthermore, it should be worth explaining that the most important role of physical spraying treatment technology can break up clumps of cyanobacteria in the process of controlling cyanobacteria bloom, but cannot change the dynamic essential characteristics of cyanobacteria and potential grazers represented by the model (2 . 1), this result implies that the physical spraying treatment technology cannot fundamentally eliminate cyanobacteria bloom. In a word, it is hoped that the results of this paper can provide some theoretical support for the physical and ecological integrated control of cyanobacteria bloom.


Introduction
In eutrophic lakes, cyanobacterial Microcystis colonies usually aggregate and float upwards to form nuisance mucilaginous Microcystis blooms, the development of mucilaginous cyanobacterial Microcystis blooms is a serious environmental and ecological problem, and information on the bloom-formation mechanism has been lacking until now [1]. Based on their long-term research work on cyanobacteria bloom in Taihu Lake, Qin et al. [2] believed that the outbreak of cyanobacteria bloom actually experienced four stages: 1) cyanobacteria cell proliferation stage; 2) cyanobacteria cell groups formation stage; 3) cyanobacteria cell groups floating stage; 4) cyanobacteria bloom outbreak stage. At the same time, they also deemed that in the process of cyanobacteria cell aggregation, the collision probability of cell groups could increase, and smaller cell masses were easier to bond, which could form larger groups and accelerate the emergence of cyanobacteria bloom. Furthermore, Chen et al. [3] pointed out that the aggregation and migration of Microcystis had many ecological advantages, which could promote the formation of Microcystis bloom. Kong et al. [4] proposed that the most direct reason for the formation of cyanobacteria bloom was the massive aggregation, growth, migration and floating to the water surface stage of Microcystis. Obviously, the formation of cyanobacteria cell aggregation is one of the keys to the outbreak of cyanobacteria bloom. Therefore, how to prevent and disperse cyanobacteria cell groups is one of the key control measures to prevent cyanobacteria bloom.
Cyanobacterial blooms are commonly referred to as harmful algal blooms (HABs) due to their vast ecological impacts and ability to generate metabolites and taste, a variety of in-lake/reservoir control measures are implemented to reduce the abundance of nuisance cyanobacteria biomass [5]. Kibuye [5,6] gave a critical review on operation and performance of source water control strategies for cyanobacterial bloom: chemical control methods, mechanical/physical control methods and biological control methods. And they pointed out that each control method had site-specific strengths, limitations, and ecological impacts, which should hint that none of the reviewed control strategies can provide a comprehensive solution to cyanobacterial blooms. In view of the cyanobacteria blooms in Taihu Lake, the competent authorities mainly adopt the comprehensive method of combining mechanical control method and biological control method, this is because that Taihu Lake is a drinking water source. Although most mechanical and biological control strategies can offer long-term control, their application can be cost-prohibitive and treatment efficacy is influenced by source water geometry and continual nutrient inputs from external sources [6]. Therefore, how to coordinate mechanical control methods with biological control methods is very important, and it is also one of the problems worthy of our in-depth study.
In order to curb the scale of cyanobacteria bloom in Taihu Lake, the management department of Taihu Lake should realize the transformation from focusing on source control and pollution interception to paying equal attention to source control, pollution interception and ecological regulation [7]. At the same time, there are many kinds of cyanobacteria in Taihu Lake, usually Microcystis [8], and allowing cyanobacteria cells to aggregate is one of the necessary conditions for cyanobacteria bloom [9]. Furthermore, once the cyanobacteria bloom breaks out, the regeneration from cell turnover and nutrient recycling by closely associated heterotrophic bacteria and microzooplankton grazers can help sustain bloom biomass, so the key biological factors to control cyanobacteria bloom must include grazing of zooplankton (possibly benthos and fish) [10]. However, some literatures have studied the dynamic relationship between cyanobacteria cell aggregation and zooplankton/fish grazing. Yang et al. [11] investigated the effect of grazing by different sorts of zooplankton on the induction of defensive morphology in the cyanobacterium Microcystis aeruginosa, and pointed out that cyanobacterium Microcystis aggregation could effectively block zooplankton predation and thus increase the survival of Microcystis aeruginosa. Burkert et al. [12] inquired into the interactions between the cyanobacteria and the potential grazer, pointed out that the predator came in direct contact with Microcystis, aggregates were formed. Chow-Fraser et al. [13] had evidence that grazers in oligotrophic lakes exert a greater impact on algae than those in eutrophic lakes. To summarise, cyanobacteria species can also benefit by their tendencies to congregate as large filamentous and colonial colonies, which reduces zooplankton predation.
At present, the cyanobacteria bloom research team of Wenzhou University is trying to use a novel comprehensive method to deal with the problem of cyanobacteria bloom in Taihu Lake, which is a combination of mechanical control method and biological control method. Firstly, in the initial stage of cyanobacteria bloom, a part of cyanobacteria agglomerates were broken by physical jet treatment device (a cyanobacteria bloom treatment ship is developed by Wenzhou University), its purpose is to destroy the aggregation stage of cyanobacteria cells and improve the chance of individual cyanobacteria cells being grazed. Secondly, some streamline ecosystems and food web structures will be reconstructed to create unfavorable conditions for cyanobacteria, some filter-feeding organisms (such as zooplankton, mussels, bivalves and fish) can be added to feed on cyanobacteria [14][15][16]. Finally, with the help of the internal operation characteristics of aquatic ecosystem, we expect that aquatic ecosystem can quickly change from algae bloom ecosystem to algae-zooplankton-fish cycle ecosystem. However, whether or not cyanobacterial bloom can be successfully controlled, the novel comprehensive method can influence changes in competition, predation, and behavioral response to predation and competition in the restored habitat. Therefore, after the cyanobacteria bloom is treated by physical jet technology, the restoration and increase of zooplankton and filter-feeding fish how to affect the self purification ability and balance of aquatic ecosystem are worthy of our further study, whether algal population can coexist with zooplankton or filter-feeding fish in the range of relatively low algal population is also worth exploring.
With the rapid development of numerical analysis and simulation technology, dynamic models have gradually become one of the powerful tools for studying natural science phenomena and problems, which can promote the relevant scientific problems to obtain some excellent research results, such as relevant literatures [17][18][19][20][21][22][23][24][25][26][27]. In references [17][18][19] applied reaction-diffusion equation to examine the calcium diffusion in the cells, investigated some nonlinear dynamics problems and got some interesting results. In reference [20] explored some dynamical behaviors of forced KdV equation to describe the free surface critical flow over a hole, and some meaningful results were given. In reference [21] inquired into chaotic dynamics of a fractional order HIV-1 model involving AIDS-related cancer cells, some meaningful results showed that order of the fractional derivative had a significant effect on the dynamics process. The paper [22] studied the dynamics of infectious diseases within a human host and in the population, and gave some meaningful theoretical and numerical simulation results. In reference [23] gave some good results to verify the effectiveness of linear feedback control for chaos. In reference [24] investigated a newly developed model for Hepatitis-B infection in sense of the Atangana-Baleanu Caputo fractional-order derivative, and obtained some good research conclusions. In references [25][26][27] looked into multiple bifurcations of discrete-time dynamic model, and some important theoretical and numerical results are obtained. In a word, the research on dynamic models has been developed rapidly, and a large number of excellent research results have emerged. Now, the disaster of aquatic ecosystem caused by cyanobacteria bloom is one of the hot issues in the world. After the cyanobacteria bloom is treated by chemical control methods, mechanical control methods and biological control methods, the essential characteristics of the internal self recovery and strengthening mechanism of aquatic ecosystem is one of the problems that need to be studied urgently. Furthermore, the algae bloom research team of Wenzhou University has built a physical spraying treatment ship, which mainly deals with the algae blooms in subtropical lakes and reservoirs, and reconstructs the biological food chain for ecological algae control after treatment. Thus, in order to understand the internal self recovery and strengthening mechanism of aquatic ecosystem, we must deeply explore the interactions between cyanobacteria (Microcystis aeruginosa) and potential grazer after cyanobacteria bloom treatment, and make clear the coexistence mode and change trend of cyanobacteria (Microcystis aeruginosa) and potential grazers based on water ecological model and dynamic simulation. Therefore, we believe that the originality of this paper is to illustrate the dynamic evolution relationship between cyanobacteria and potential grazers by means of the bifurcation dynamic evolution process of the model (2.1), and give the evolutionary characteristics of the coexistence mode of cyanobacteria and potential grazers.

Ecological mathematical modeling
Bertalanffy [28] proposed a method to study the biological system using mathematical model in 1932, this approach is a combination of coordination, order and purpose, which can form three basic ideas of studying the biological system, namely: trophic-level analysis, system perspective and dynamic view. Since then, the research on ecological model and its related dynamics has developed unprecedentedly, and a large number of literature and works have appeared, such as fundamentals of ecological modelling [29], a nitrogen-phosphorus-blue-green algae model [30], some nitrogen and phosphorus dose-response models [31], a mathematical model including three types of zooplankton and two types of fish as well as two types of algae and nutrients [32], an aquatic bacteria-algae amensalism model [33], an ecological model considering two types of phytoplankton, three types of zooplankton, planktivorous fish, detritus and dissolved matters [34], an aquatic algae-fish model [35], one dimensional hydrodynamic-ecological model [36], a two-component model with nutrients and cyanobacteria [37], a four-component model including nutrient, unicellular algae, colonial algae and herbivorous zooplankton [38], an algal bloom mathematical model involving zooplankton-phytoplankton population [39]. Obviously, the ecological model has been developed rapidly in water eutrophication and algal bloom, and some meaningful theoretical results have been obtained. Although there is a certain understanding of the mechanism and control effect of cyanobacteria bloom in recent years, there are still many unsolved mysteries in the outbreak and control of cyanobacteria bloom.
In order to make clear the coexistence mode and change trend of cyanobacteria (Microcystis aeruginosa) and potential grazers (mussels, bivalves, and silver carp) after cyanobacteria bloom is treated by comprehensive control strategy, which can help us explore the interactions between cyanobacteria and potential grazers, and explicit the effectiveness and feasibility of the comprehensive control strategy, we will construct a new aquatic ecological model based on the comprehensive control strategy and internal operation law of aquatic ecosystem. Therefore, some model assumptions and cyanobacteria bloom background will be given as follows: 1) The growth environment of cyanobacteria and potential grazers is assumed to be a semi closed aquatic ecosystem, species biomass is always evenly distributed in space and changed instantaneously with time t. x(t) and y(t) represent respectively the biomass of total cyanobacteria (unicellular cyanobacteria and colonial cyanobacteria) and potential grazers. After cyanobacteria bloom is treated by mechanical control method (physical jet processing device), total cyanobacteria will be divided into colonial cyanobacteria α 1 x(t) and unicellular cyanobacteria α 2 x(t) with α 1 + α 2 = 1. If α 1 = 0 and α 2 = 1, this shows that cyanobacteria population mainly exists in the form of unicellular cyanobacteria, and the cyanobacteria population has no aggregation behavior and can not form a group aggregation state. If α 1 = 1 and α 2 = 0, this indicates that cyanobacteria bloom has erupted in a large area, especially serious, and all cyanobacteria populations are clustered, which is an extreme state of existence and almost does not exist in reality.
2) Because when unicellular cyanobacteria gather into groups, the light utilization rate of algae cells can be enhanced [40], the salt stress resistance and high light stress of algal cells can be improved [41,42], the ability to avoid the invasion of viruses, bacteria and algae phages can be strengthened [43], nutrient and inorganic carbon absorption capacity can be greatly promoted [44,45], the risk of predation can be greatly reduced [46]. In light of the above, cyanobacteria aggregation can find more suitable niche, obtain greater ecological advantages and enhance growth rate. Therefore, we assume that the growth kinetics function of cyanobacteria α 2 ) and maximum environmental capacity k, where r 1 is a median growth rate with α 1 = 1 2 and α 2 = 1 2 , r min is a growth rate without aggregation, and, r max is a growth rate with overall aggregation state.
3) During cyanobacteria blooms, cyanobacteria can induce morphological changes to colonial forms in order to avoid zooplankton grazing, this mechanism can be an adaptive reaction for survival that has developed through evolutionary processes [47,48]. For example, Microcystis exists as colonial forms, the relevant research results are in the papers [35,38,43]. However, unicellular cyanobacteria do not have this protective mechanism. Hence, we use Holling type II functional responses a 2 b+α 2 x(t) to describe the predation mechanism of potential grazers on unicellular cyanobacteria, this is because the predation process of potential grazers not only take time to capture, but also depend on the abundance of unicellular cyanobacteria, where a 2 is a grazing coefficient and b is a half-saturation constant. At the same time, we use mathematical functional a 1 (α 1 x(t) − g( α 1 α 2 ))y(t) to describe the predation mechanism of potential grazers on colonial cyanobacteria, a 1 is a grazing coefficient and predation function g( α 1 α 2 ) avoided by colonial cyanobacteria, and where m is a median biomass of colonial cyanobacteria (which can never be hunted) with α 1 = 1 2 and α 2 = 1 2 , 0 indicates that unicellular cyanobacteria do not have this ability, m max is a maximum biomass of colonial cyanobacteria (which can never be hunted) with overall aggregation state.
4) It is widely known that potential grazers can not eat only cyanobacteria, and must forage other species as alternate foods. Thus, the growth kinetics function of potential grazers is b+α 2 x(t) with alternate coefficient r 2 and absorption conversion coefficients e 1 and e 2 . Furthermore, cyanobacteria and potential grazers have their own other losses (such as natural death), we might as well set them to m 1 x(t) and m 2 y(t) with loss coefficient m 1 , m 2 .
Based on the above assumptions and analysis, we can construct a new aquatic ecological control model and two special aquatic ecological models with non-negative initial conditions x(0) ≡ x 0 ≥ 0 and y(0) ≡ y 0 ≥ 0, which can be described as following: The model (2.1) represents the mixed state of unicellular cyanobacteria and colonial cyanobacteria, the biomass of unicellular cyanobacteria and colonial cyanobacteria is controlled by critical parameter α 1 or α 2 because of α 1 + α 2 = 1. If the value of α 1 belongs to (0, 0.2], the model (2.1) represents cyanobacteria cell proliferation stage. If the value of α 1 belongs to (0.2, 0.4], the model (2.1) represents cyanobacteria cell groups formation stage. If the value of α 1 belongs to (0.4, 0.6], the model (2.1) represents cyanobacteria cell groups floating stage. If the value of α 1 belongs to (0.6, 1), the model (2.1) represents cyanobacteria bloom outbreak stage. Overall, the model (2.1) can describe that the outbreak of cyanobacteria bloom in Taihu Lake actually experienced four stages [2]. The models (2.2) and (2.3) represent the stage of unicellular cyanobacteria and colonial cyanobacteria respectively, neither of these two states will exist in reality, therefore, we only study them as two limit states.
With the help of bifurcation dynamic analysis and dynamic simulation of the model (2.1), the main purpose of this paper is to explore how the aquatic ecosystem can recover and strengthen itself after the comprehensive treatment of cyanobacteria bloom. To solve this problem, we firstly make a qualitative analysis of the model (2.1), deduce some critical threshold conditions to ensure the existence and stability of some equilibrium points, and give some critical conditions under which the model can induce Transcritical and Hopf bifurcation. Secondly, we will comprehensively investigate how comprehensive management measures affect the dynamic characteristics of cyanobacteria bloom and evaluate the effectiveness of comprehensive management measures in controlling cyanobacteria bloom. Finally, we carry out relevant dynamic experiments to verify the feasibility of the theoretical derivation results and identify the possible coexistence mode of cyanobacteria and potential grazers, then we can analyze interaction mechanism between cyanobacteria and potential grazers in the process of cyanobacteria bloom outbreak and control.
3. Superiority analysis of the proposed model (2.1) In order to investigate the superiority of the model (2.1), the dynamic relationship between the biomass of cyanobacteria and potential grazers was qualitatively analyzed and numerically simulated with  From the model (2.2), we can find that the relation expression of cyanobacteria x and potential grazers y is given by It is easy to see from Figure 1(a) that the dynamic relationship of cyanobacteria x and potential grazers y is a concave function as the value of cyanobacteria x increases, and there is a unique maximum value, and the maximum value of cyanobacteria population is still low. Unfortunately, the value range of cyanobacteria x is relatively small, and the biomass of potential grazers y becomes negative with the increase of cyanobacteria x. Obviously, this dynamic result can not better reflect the running law of cyanobacteria bloom control test in the laboratory. From the model (2.3), we can find that the relation expression of cyanobacteria x and potential grazers y is given by It can clearly find from Figure 1(b) that the dynamic relationship of cyanobacteria x and potential grazers y is a convex function as the value of cyanobacteria x increases, and there is a unique minimum value 0, which is not in line with the wild growth law of cyanobacteria population because cyanobacteria population will never die out. Furthermore, it is even more regrettable that the value of potential grazers y approaches positive infinity with the increase of cyanobacteria x value, which completely violates the dynamic relationship between cyanobacteria x and potential grazers y in the process of cyanobacteria bloom. From the model (2.1), we can find that the relation expression of cyanobacteria x and potential grazers y is given by It is worth pointing out that potential grazers y is monotonically increasing and has supremum r 1 as the value of cyanobacteria x increases to positive infinity. Furthermore, the better result is that cyanobacteria x approaches a small fixed value when potential grazers y approaches positive infinity, this result is more in line with the natural law (see Figure 1(c)). Moreover, what is more gratifying is that there is a unique minimum value, through which we can regulate the recyclability of cyanobacteria bloom test, which is also our favorite routine manipulation. In a word, through the comparative analysis of numerical simulation results, it can be seen that the model (2.1) can more accurately describe the dynamic relationship between cyanobacteria and potential grazers under the mechanical and ecological integrated control, so we mainly study the relevant dynamic problems of the model (2.1).

Existence and stability of equilibrium point
The existence and local stability of all possible equilibrium points will be investigated in detail. The equilibrium point is a kind of special solution of the model (2.1), which can determine some coexistence modes of cyanobacteria and potential grazers, and it is also the basis of subsequent bifurcation dynamics research. Thus, the existence and local stability of all possible equilibrium points will be futher investigated in this section. Now, we consider the algae isocline vertically and potential grazers isocline horizontally, respectively: (2.4) By solving the above equations with zero dash, there are two classes of equilibrium points: boundary equilibrium points E 0 and E 1 , the coordinates of internal equilibria E * 1 and E * 2 . It is obvious that the equilibrium points are the intersections of these nullclines. Thus, we easily see that the model (2.1) possesses two boundary equilibrium points E 0 (0, 0) and E 1 (x 1 , 0). For the possible positive equilibria, we only need consider the positive solutions of the following equations: For positive equilibrium points, the algae population x must satisfy m α 2 < x < k, Let ∆(x) denote the discriminant of the first equation of (2.5) and express ∆(x) in terms of, i.e, Hence, we have With respect to the conditions of the equilibrium points, we obtain the following theorem: Theorem 1. The equilibrium points of the model (2.1) are as follows: the model(2.1) always has a boundary equilibrium point E 0 (0, 0). If α 1 α 2 > m 1 r 1 , then the model (2.1) has a boundary equilibrium point Theorem 2. For the positive equilibrium point, we have: 1) When m ≤ I 1 , the model (2.1) has no positive equilibrium point E * 1 (x * 1 , y * 1 ), 2) When r 2 − m 2 + a 2 e 2 > 0 and r 2 − m 2 < 0, if max{0, I 1 } < m < min{I 2 , I 3 , I 4 } or max{0, I 1 , I 2 } < m < min{I 3 , I 4 }, then the model (2,1) has a unique positive equilibrium point E * 1 (x * 1 , y * 1 ), 3) When r 2 − m 2 + a 2 e 2 < 0, if max{0, I 1 , I 3 } < m < min{I 2 , I 4 } or max{0, I 1 , I 2 , I 3 } < m < I 4 , then the model (2.1) has a unique positive equilibrium point E * Proof. Let the solution of the nullclines g(x) = Ax 2 + Bx + C = 0, we know that the model (2.1) has the positive roots only if r 2 − m 2 − a 1 e 1 α 1 m α 2 < 0, which means m ≥ I 1 . At the same time, we also show that if m ≥ I 1 , then ∆ > 0 and C < 0, which means that g(x) has a unique positive root E * 1 (x * 1 , y * 1 ). Hence, if g( m α 2 ) < 0 and g(k) > 0, there is a positive real root in m α 2 < x < k. Due to the fact that ∆ > 0, hence, if g( m α 2 ) < 0 and g(k) > 0, we have the conditions with I 3 and I 4 . If g( m α 2 ) < 0, it follows that r 2 − m 2 + a 2 e 2 > 0 for m < I 3 or r 2 − m 2 + a 2 e 2 < 0 for m > I 3 . If m < I 4 , we have f (k) > 0. Next, we just discuss the symbols of B.
Theorem 5. The conditions for the existence of the internal equilibrium point E * 1 (x * 1 , y * 1 ) are shown in Theorem 2. In this paper, it has been proved that there is and only one internal equilibrium point E * 1 (x * 1 , y * 1 ). If T r(J E * 1 (x * 1 , y * 1 )) < 0, then the internal equilibrium point E * 1 (x * 1 , y * 1 ) is locally asymptotically stable; otherwise, the internal equilibrium point E * 1 (x * 1 , y * 1 ) is unstable. Proof. The Jacobian matrix of the internal equilibrium point E * 1 is The determinant of the matrix J E * 1 is ).
The trace of the matrix J E * 1 is It is obvious that all terms of the product of determinants are positive, hence, we can judge that the determinants of the Jacobian matrix of E * 1 is always positive. Thus, when T r(J E * 1 ) < 0, the internal equilibrium point E * 1 has two negative eigenvalues, which means that the internal equilibrium point E * 1 is locally asymptotically stable. Otherwise, both characteristic roots are positive, and the internal equilibrium point E * 1 is unstable.

Local bifurcation analysis
In order to probe some possible coexistence evolution modes of cyanobacteria and potential grazers during the outbreak of cyanobacteria bloom, and clarify how the sheltering effect of colonial cyanobacteria affects the coexistence mode of cyanobacteria and potential grazers, we choose m as the control parameter to study the bifurcation dynamic behavior of the model (2.1), mainly including Transcritical bifurcation and Hopf bifurcation. If the model (2.1) can have Transcritical bifurcation, which implies that the cyanobacteria and potential grazers can coexist in a certain range at least, and there will be no extinction of at least one population with the advance of time. If the model (2.1) can have Hopf bifurcation, which shows that cyanobacteria and potential grazers can coexist in a periodic oscillation mode in a certain range, and the two populations can continue to survive with the advance of time. Therefore, bifurcation dynamics of the model (2.1) has important significance in population ecology and is worthy of our in-depth study.

Hopf bifurcation
In the view of the proof of Theorem 5, we can know that Hopf bifurcation can occur at the positive equilibrium point E * 1 . At the same time, the positive equilibrium point E * 1 can change its stability when the value of m passes through the critical magnitude. Thus, we summarize our findings in the following theorem. Proof. The characteristic equation of matrix J E * 1 is λ 2 − T r(J E * 1 )λ + Det(J E * 1 ) = 0, and three crosssectional conditions of Hopf bifurcation occurs, i.e., 1) T r(J E * 1 ) m=m HP = 0, 2) Det(J E * 1 ) m=m HP > 0, 3) dT r(J E * 1 ) dα | m=m HP = 0. Previously, we can show that Det(J E * 1 ) m=m HP > 0. And when m = m HP , we have T r(J E * 1 ) m=m HP = 0.
In conclusion, the model (2.1) can induce a Hopf bifurcation when m = m HP . Next, we will discuss the stability of the limit cycle by computing the first Lyapunov number. Translate the equilibrium point E * 1 to the origin by using the following transformation: where a 10 = a 1 α 1 y * 1 + And then Q(x m , y m ), P(x m , y m ) are power series in (x m , y m ) with terms x i m y j m , which can satisfy i + j ≥ 4. Therefore, the first Lyapunov number to determine the stability of limit cycle is given by the formula It is well known that the limit cycle generated by the Hopf bifurcation is unstable if l 1 > 0 and stable if l 1 < 0. Nevertheless, the expression of the first Lyapunov number is too complex to judge the positive and negative. Therefore, we will give some numerical solutions in Section 6.

Numerical simulation analysis
Under the research framework of theoretical derivation, we will not only verify the feasibility and effectiveness of the theoretical results, but also explore how clustering affects the dynamic relationship between cyanobacteria and potential grazers, thus the parameter m and β = α 1 α 2 are selected as control variables for numerical simulation with r 1 = 0.56, r 2 = 0.2, k = 5, a 1 = 0.6, a 2 = 0.4, b = 0.15, m 1 = 0.18, m 2 = 0.28, e 1 = 0.6, e 2 = 0.28, it should be noted that all parameter values are estimated according to biological significance and theoretical conditions. Furthermore, β = α 1 α 2 represents the proportion of colonial cyanobacteria and unicellular cyanobacteria, which can also indicate the outbreak degree of cyanobacteria bloom. If β = α 1 α 2 > 1 indicates that the manifestation of cyanobacteria is mainly colonial cyanobacteria, which is not conducive to the capture of potential grazers. If β = α 1 α 2 < 1 indicates that the manifestation of cyanobacteria is mainly unicellular cyanobacteria, which is conducive to the capture of potential grazers. Moreover, it should be pointed out that the predation behavior of potential grazers on cyanobacteria is essentially the implementation of ecological control technology.
When we take β = 0.5385 with α 1 = 0.35 and α 2 = 0.65, which represents that the cyanobacteria in the whole ecosystem are in the bloom period, showing the phenomenon of aggregation. In other words, the current situation is not conducive to potential grazers to prey on cyanobacteria. Based on numerical simulation work, the bifurcation diagram of the model (2.1) with control parameter m has been shown in Figure 2(a). It is easy to find that the model (2.1) will experience a Hopf bifurcation and Transcritical bifurcation with the value of m changing within [0, 1.8]. When the value of m gradually decreases from 1.8 to m TC = 1.42282, the model (2.1) will occur a Transcritical bifurcation (Figure 3), the boundary equilibrium point E 1 will change from a stable state to an unstable state, which suggests that potential grazers will not always be close to extinction.
And when the value of m is between (m HP , m TC ), the model (2.1) has an asymptotically stable internal equilibrium point E * 1 (seeing Figure 3(b)), which indicates that cyanobacteria and potential grazers can coexist in a stable equilibrium state. If the value of m gradually decreases to m HP = 0.2957018, the model (2.1) will occur Hopf bifurcation, the internal equilibrium point will change from a stable state to an unstable state, and a stable limit cycle will appear, the detailed numerical simulation results are shown in Figures 4 and 5. Furthermore, we can get that the Lyapunov exponent is −2.57010 × 10 7 , which can further verify that the limit cycle is stable. Moreover, it should be explained that the dynamic characteristics of the model (2.1) have changed substantially, which means    When we take β = 1.083 with α 1 = 0.52 and α 2 = 0.48, which represents that the cyanobacteria in the whole ecosystem are still in the bloom period, but the current situation is quite disadvantageous to potential grazers to prey on cyanobacteria because unicellular cyanobacteria account for 0.48. It is easy to discover from Figure 2(b) that the model (2.1) has almost the same dynamic evolution process.  When the value of m gradually decreases from 2 to m TC = 1.74941, a Transcritical bifurcation will occur ( Figure 6). When the value of m continues to decrease to m HP = 0.125, the instability of the equilibrium point can produce a stable limit cycle, the detailed dynamic results are shown in Figures 7  and 8. At the same time, the Lyapunov exponent is −5.15158 × 10 14 , which can show the stability of the limit cycle. From the essence of system dynamics evolution, whether the value of β is greater than 1 or less than 1, the model (2.1) has exactly the same dynamic evolution characteristics, that is, Transcritical bifurcation and Hopf bifurcation. That is to say, the ratio of colonial cyanobacteria to unicellular cyanobacteria does not affect the evolution characteristics of the dynamic state of the model (2.1). Moreover, it is also worth emphasizing that the value of parameter m seriously affects the dynamic relationship and coexistence mode between cyanobacteria and potential grazers.  Based on the comparative analysis of numerical simulation results, a Transcritical bifurcation will occur at m TC = 1.42282 when the value of β is less than 1, however a Transcritical bifurcation will occur at m TC = 1.74941 when the value of β is greater than 1. This result means that during the outbreak of cyanobacteria bloom, the earlier the clustered cyanobacteria population is crushed by physical spray treatment technology, the more opportunities there are to change the dynamic relationship between cyanobacteria and potential grazers, and provide an entry point for the coexistence of cyanobacteria and potential grazers. However, the critical value of inducing Hopf bifurcation is just the opposite.
When the value of β is less than 1, Hopf bifurcation can be induced at m = 0.2957018, which suggests that potential grazers actively prey on unicellular cyanobacteria to obtain sufficient food, which can lead to a large increase in the biomass of potential grazers. This increased predation ability can lead to essential changes in the dynamic relationship between them, and then play a role in controlling the large-scale outbreak of cyanobacteria bloom. However, when the value of β is greater than 1, Hopf bifurcation can be induced at m = 0.14124142, which means that the number of unicellular cyanobacteria can not basically maintain the growth needs of potential grazers, and a large number of colonial cyanobacteria need to be prey, which will cause colonial cyanobacteria to be gradually broken. This result also leads to a stable equilibrium between cyanobacteria and potential grazers in the range of m. Therefore, from the perspective of controlling cyanobacteria bloom by physical spray treatment technology, the colonial cyanobacteria turns into unicellular cyanobacteria after the cyanobacteria bloom is treated by physical spray, which is conducive to potential grazers to prey, and then lead to the coexistence mode of potential grazers and cyanobacteria in a stable state.

Conclusions
Within the framework of physical and ecological integrated control and management of cyanobacteria bloom, according to the growth characteristics and comprehensive management characteristics of wild cyanobacteria population, we mainly focus on the impact of physical spraying treatment technology on the structural characteristics of cyanobacteria population and the predation effect of potential grazers, and give some ecological modeling assumptions. On this basis, a aquatic ecological dynamic new model is proposed to describe the dynamic relationship between cyanobacteria and potential grazers. At the same time, it is worth emphasizing from Figure 1 that the biggest advantage of this new model is that it integrates three kinds of dynamic action modes, of which the sub model (2.1) is the most distinctive and has not been studied by predecessors. Furthermore, the biggest advantage of the model (2.1) is that it can perfectly describe the dynamic relationship between cyanobacteria and potential grazers under the physical and ecological integrated control, most commendably, it can control cyanobacteria bloom from the perspective of minimum value, and the consumption of potential predators is limited.
Underlying the analysis of mathematical theory, we first explore the existence and stability of the equilibrium point of the model (2.1), and give some critical conditions of some critical parameters. Then, we obtain some key threshold conditions to ensure that the model (2.1) has specific dynamic behaviors, such as Transcritical bifurcation and Hopf bifurcation. These theoretical results are the basis for our numerical simulation analysis, which provides theoretical support for in-depth exploration of the dynamic relationship between cyanobacteria and potential grazers.
Based on the results of numerical analysis, the dynamic behavior type and dynamic evolution process of the model (2.1) are clarified, focusing on the Transcritical bifurcation and Hopf bifurcation. It is easy to find from Figure 2 that when cyanobacteria bloom break out and gather into clusters, we implement a physical and ecological comprehensive control strategy for cyanobacterial bloom, the dynamic relationship between cyanobacteria and potential grazers suddenly changes by Transcritical bifurcation, which can form coexistence mode in a stable equilibrium state, the detailed dynamic simulation results are shown in Figures 3 and 6. After then, the biomass of cyanobacteria and potential grazers gradually decrease with the reduction of agglomeration, finally, a stable periodic oscillation coexistence mode is formed to maintain a good cycle of ecological algae control strategy, these dynamic simulation results are shown in Figures 4, 5, 7 and 8. Furthermore, it is more meaningful to obtain from Figures 3 and 6 that Transcritical bifurcation represents evolution process of steady-state coexistence mode, and it is the most ecologically significant result form Figures 4 and 7 that Hopf bifurcation represents evolution process of periodic oscillation coexistence mode. In addition, it should be emphasized that the most important role of physical spraying treatment technology is to break up clumps of cyanobacteria in the process of controlling cyanobacteria bloom, which will not change the dynamic essential characteristics of cyanobacteria and potential grazers represented by the model (2.1).
Based on theoretical and numerical research results, the main results of this paper are as follows: 1) The population density of cyanobacteria decreases gradually as the value of key parameter m decreases, that is to say, the implementation of physical spraying treatment technology can effectively control the population density of cyanobacteria. 2) Transcritical bifurcation can induce steate-state coexistence mode between cyanobacteria and potential grazers, and Hopf bifurcation can force cyanobacteria and potential grazers to form periodic oscillation coexistence mode. 3) The evolutionary characteristics of bifurcation dynamics suggest that physical and ecological integrated control technique can effectively inhibit the growth of cyanobacteria, and promote the formation of a stable coexistence mode between cyanobacteria and potential grazers within a controllable range.
In general, some results have been obtained by proposing a new aquatic ecological model, however, the work of this paper still has some limitations, such as: 1) the characteristics of toxicity released by cyanobacteria population were not considered in the modeling, 2) cyanobacteria migration behavior is not mentioned, 3) modeling assumptions make that the application of the model (2.1) has certain limitations. Thus there is still much work to be further studied, for example, how does the release of cyanobacteria toxins affects the predatory ability of potential grazers during physical spray treatment, effect of cyanobacteria migration behavior on algae control effect of physical spraying technology,etc. In a word, it is hoped that the research results of this paper will provide some theoretical support for the application of physical and ecological comprehensive management of cyanobacteria bloom.