Stability and bifurcation analysis of the Bazykin’s predator-prey ecosystem

Abstract: In the paper, stability and bifurcation behaviors of the Bazykin’s predator-prey ecosystem with Holling type II functional response are studied theoretically and numerically. Mathematical theory works mainly give some critical threshold conditions to guarantee the existence and stability of all possible equilibrium points, and the occurrence of Hopf bifurcation and Bogdanov-Takens bifurcation. Numerical simulation works mainly display that the Bazykin’s predator-prey ecosystem has complex dynamic behaviors, which also directly proves that the theoretical results are effective and feasible. Furthermore, it is easy to see from numerical simulation results that some key parameters can seriously affect the dynamic behavior evolution process of the Bazykin’s predator-prey ecosystem. Moreover, limit cycle is proposed in view of the supercritical Hopf bifurcation. Finally, it is expected that these results will contribute to the dynamical behaviors of predator-prey ecosystem.


Introduction
As we all know, the classic Lotka-Volterra predator-prey system has been used to simulate predatory phenomena in nature, and its impact on the researches of mathematical biology and ecology will be roughly equivalent to the atomic bomb effect. In 1965, Holling proposed several different functional responses to characterize the dynamic predator-prey relationship between populations, which could describe the predator population how to transform the captured prey population into its own growth ability. Generally speaking, functional responses in predator-prey systems mainly depend on many internal and external key factors, such as the densities of prey and predator. On the other hand, preydependent functional responses are an important role in mathematical ecology, especially the dynamics of predator-prey systems. In recent decades, many scholars have done a lot of research on the predatorprey function response, and have made some excellent research results. Now, it is briefly summarized to enrich the population dynamics modeling system.
(i) The Holling type I functional response [1] is where c is a positive constant.
(iii) The Holling type III or the p = 2 S-type functional response [1,4,5] is where α and β are positive constants. The generalized Holling type III or sigmoidal is Φ(x) = mx 2 where m and a are positive constants, b is a constant. When b = 0, the function Φ(x) can reduce to above Holling type III functional response. The S-type functional response with index p is Φ(x) = αx p β+x p , where α, β and p are positive constants. (iv) The generalized Holling IV or the Monod-Haldane functional response [1,[4][5][6][7] is where m and a are all positive constants, but b is a constant. When b = 0, it is called the Holling type IV functional response or the simplified Holling type IV functional response [8,9].
(v) The Beddington-DeAngelis type functional response [4,[10][11][12][13] is where α, a, b and c are positive constants, which is originally and independently introduced by Beddington and DeAngelis [10,11]. At the same time, it is similar to the Holling type II functional response incorporating an extra term cy in denominator, which can describe mutual interference among predators [12,13]. This functional response has some same qualitative features as the ratio-dependent form, but can avoid some singular behaviors of ratio-dependent models at low densities [12].
(vii) The Crowley-Martin type functional response [19] is Φ(x, y) = αx a + bx + cy + exy , where α, a, b, c and e are positive constants. It involves the interference among individual predator engaged in searching and handling the preys. While in [20][21][22], the authors particularly take the Crowley-Martin functional response in the new type of Φ(x, y) = αx (1+ax)(1+by) , which is proposed by Bazykin and is an immense breakthrough of the Holling type II and Beddington-DeAngelis functional responses.
In this paper, we will continually concentrate on a detailed discussion in the well-known Bazykin's predator-prey ecosystem with the Holling type II functional response and interspecific densityrestricted effect on the predator, which is also a variation of Volterra's classical predator-prey model and is expressed in the form of following ordinary differential equations (ODEs) [23,24]: Here the functions x = x(t) and y = y(t) represent densities of the prey population and predator population at time t, respectively. All positive parameters have practically biological meanings: r 1 denotes the intrinsic growth rate of the prey population, K 1 represents the carrying capacity of the environment, a is the half-saturation constant; α is the search efficiency of predator for prey, m 1 and m 2 are the mortality rate of the prey and predator population, e is the biomass conversion, d is the intra-specific competition coefficient. The specific growth term r 1 x(1 − x K 1 ) governs the increase of the prey in the lack of predator. The square non-linear term term dy 2 , denotes intrinsic decrease of the predator, and represents interspecific density-restricted effect on the predator. Excluding the dy 2 , the system (1.1) is based on the classical Gause type predator-prey system, which expresses the following form [25]:ẋ = xg(x, k) − yp(x), (1.2a) y = y(−d + cq(x)), (1.2b) where g(x, k) is a continues function for x > 0, p(x) is a functional response of predators, which is called Holling-type-II predator-prey model as well. Such ordinary differential system of predator-prey populations is familiar to the Lotka-Volterra system, in which populations have the addition of damping terms(or self inhabit). The positivity of solutions with respect to initial condition x(0) > 0 and y(0) > 0 is easy to prove and we omit its proof. This system with positive and bounded solutions is also well behaved as we intuits from the biological significance. Bazykin [24] wholly discussed the stability of equilibria, global existence of limit cycles, global attractivity of such equilibria, Hopf and codimension 2 bifurcations. In [26], analytical description and alteration of local stability were given. Here the authors investigated a familiar Lotka-Volterra system, in which the populations have self-inhibit(the addition of damping term) for global stability and existence of limit cycles [27]:ẋ = x(1 − k 1 x − k 2 x 2 ) − xy 1 + ax , (1.3a) y = rxy 1 + ax − y(δ 0 + δ 1 y), (1.3b) where specific growth rate governs the growth of the prey in the absence of predator and it has an increase (or decrease) intrinsic rate on the predator. They already proved the existence of two limit cycles with the help of idea from the Poincare-Bendixson theory. Obviously, at special case k 2 = 0, we note that the system (1.3) can reduce to above system (1.1), which was analyzed by [28] and [29] for stability of equilibria, Hopf bifurcations, global attraction and codimension two bifurcations as well. While in the respect of global behavior, the system (1.1) was investigated by [30]. In [31], for a particular form of the system (1.1) with a modified Holling type II functional response β(x−m) incorporating a constant prey refuge m, the authors therein gave sufficient conditions to guarantee the global stability of the positive equilibrium and uniqueness of a stable limit cycle. In [32], the authors revealed a rich bifurcation structure, including supercritical and subcritical Bogdanov-Takens bifurcation.
Based on classical biological manipulation theory, the Bazykin's predator-prey ecosystem can be used to explore the dynamic relationship between Microcystis aeruginosa and filter-feeding fish from the perspective of population dynamics, where x(t) and y(t) represent respectively the density of Microcystis aeruginosa and filter-feeding fish (bighead carp and silver carp), the growth kinetics function of Microcystis aeruginosa x is r 1 x(1 − x K 1 ) with intrinsic growth rate r 1 and maximum environmental capacity K 1 . The grazing function of filter-feeding fish y is αxy 1+ax with capture coefficient α and density restriction coefficient a. Furthermore, the parameter m 1 and m 2 are natural mortality of Microcystis aeruginosa and filter-feeding fish, the parameter d and e are internal competition coefficient and energy conversion rate of filter feeding fish. In order to deeply explore the dynamic relationship between Microcystis aeruginosa and filter-feeding fish, it is necessary to investigate some bifurcation dynamic behaviors of the Bazykin's predator-prey ecosystem, thus we mainly focus on the stability and bifurcation of the Bazykin's predator-prey ecosystem in this paper. Firstly, we investigate the existence and stability of hyperbolic equilibrium point and non-hyperbolic equilibrium point, and conveniently study the cusp of condimension 3. Secondly, we explore Hopf bifurcation and Bogdanov-Takens bifurcation in detail, and give some sufficient threshold conditions. Finally, for Hopf bifurcation dynamics, we especially analyze the limit cycle via a perturbation procedure and canonical transformation. Moreover, we think that these mathematical analysis results can provide a theoretical basis for numerical simulation, which can give some biological interpretation for Hopf bifurcation and Bogdanov-Takens bifurcation, hence, we mainly study stability and bifurcation dynamics of the Bazykin's predator-prey ecosystem from the perspective of mathematical analysis, other biological significance issues will be completed in the follow-up work.

Existence and stability analysis of equilibrium point
All solutions of the system (1.1) are non-negative and bounded with initial conditions x(0) 0, y(0) 0, thus it is namely dissipative in the first quadrant R +2 of 2 dimensional space R 2 and well-defined on the closed domain R 2 + = R +2 . Furthermore, the system is uniformly bounded with lim sup t→+∞ x(t) M 1 , lim sup t→+∞ y(t) M 2 , in which two positive constants M 1 and M 2 only depend on parameters r 1 , K 1 , α, a, m 1 , e, m 2 and d [33]. In other words, the system (1.1) is confined in the domain {(x, y)|0 z M + , for any > 0}, (2.1) with z = ex + y and a constant M > 0. Moreover, the system (1.1) is permanent if the value of all parameters can satisfy with λ ∈ (0, 1) [34].

Boundary equilibrium point
It is clear that the system (1.1) admits two biological boundary equilibria E 0 := (0, 0) and shows that E 0 is a hyperbolic asymptotically stable node(unstable node) when r 1 < m 1 (r 1 > m 1 ), while it is a stable node(non-hyperbolic attractor) with only one zero eigenvalue if r 1 = m 1 (see the Theorem 7.1 in [35]). The Jacobian matrix at point E 2 is then, the hyperbolic point E 2 is an asymptotically stable node(unstable node) when αex 2 a+x 2 < m 2 ( αex 2 a+x 2 > m 2 ). When αex 2 a+x 2 = m 2 , E 2 is also a stable node(non-hyperbolic attractor) with characteristic direction tan(θ) = (m 1 −r 1 )e m 2 under the polar coordinate transformation. Now, we use an example to verify the stability of the equilibrium point E 0 and E 2 with r 1 = 0.6, a = 1.5, α = 0.5, e = 0.6, K 1 = 20 and d = 0.1. It is easy to find from Figure 1(a) that E 0 is a stable node with characteristic direction θ = 0 when r 1 = m 1 and m 2 = 0.06. Furthermore, it is obvious to see from Figure 1(b) that E 2 is a stable node with characteristic direction tan(θ) = (m 1 −r 1 )e m 2 when αex 2 a+x 2 = m 2 and m 1 = 0.3. Figure 1. (a) E 0 is a stable node; (b) E 2 is a stable node.

Hyperbolic equilibrium point
At first, an interior equilibrium E * = (x * , y * ) of the system (1.1) always satisfies following algebraic equations Furthermore, from Sylvester's resultants in polynomial Equation (2.4), components x * and y * must be positive roots of third-order polynomial equations(cubic equations) p(x) = 3 i=0 a i x i = 0 and q(y) = 3 i=0 b i y i = 0, respectively. Here the coefficients are listed as follows: Then, we define these complicated expressions which is a discriminant of above cubic equations p(x) = 0 and q(y) = 0 for later use respectively [26]. The Eq (2.4) also implies that such interior equilibrium E * does not exist when one of conditions holds: (i) r 1 m 1 ; (ii) αe m 2 . The rest of our paper always assume the necessary existence condition r 1 > m 1 and αe > m 2 . If condition 0 < am 2 αe−m 2 < x 2 holds, it is cleat that E * always exists. Thus we define the Jacobian matrix of the system (1.1) at an interior equilibrium E * = (x * , y * ) as (2.5) The trace, determinant and discriminant of the matrix J(E * ) are denoted as and ∆ * = ∆ * (E * ) := A 2 1 − 4A 2 , respectively. Then, by using the Perron's theorems and the Routh-Hurwitz's criteria, we have following local stability of a hyperbolic equilibrium E * in general cases: (a) If A 1 < 0 and (a1) A 2 > 0, ∆ * 0, then E * is an asymptotically stable node; (a2) A 2 > 0, ∆ * < 0, then E * is an asymptotically stable focus; (a3) A 2 < 0, then E * is a saddle point; (b) If A 1 = 0 and (b1) A 2 > 0, then E * is a center or a focus; (b2) A 2 < 0, then E * is a saddle point; (c) If A 1 > 0, then E * is unstable and (c1) ∆ * = 0, then E * is a node; (c2) ∆ * < 0, then E * is a focus; (c3) ∆ * > 0 and A 2 > 0, then E * is a node; (c4) ∆ * > 0 and A 2 < 0, then E * is a saddle point.
A non-hyperbolic point E * is a stable(unstable) node if A 1 < 0(A 1 > 0) and A 2 = 0. The nilpotent E * is probable a cusp of codimension at least 2, which can ensure potential Bogdanov-Takens bifurcation when A 1 = A 2 = 0.
Obviously, when r 1 m 1 , the equilibrium point E 0 is globally asymptotically stable, which can be proved by using a Lyapunov function V = ex + y. Similarly, when r 1 > m 1 and αe < m 2 , the point E 0 is unstable and E * does not exist, thus E 2 is globally asymptotically stable. For equilibrium point E * , we define a positive definite Lyapunov function Now, along solutions of the system (1.1), differentiate V with regard to time t to obtain Substituting the value of dx dt and dy dt from the system (1.1), we can get Thus, it is obvious that if αy * a(a+x * ) < r 1 K 1 , then dV dt 0. This equality holds if and only if (x, y) = E * , i.e., the equilibrium point E * is globally asymptotically stable.
Under a generalized condition αy * a(a+x * ) r 1 K 1 , the hyperbolic point E * is a locally asymptotically stable focus or node since A 1 < 0 and A 2 > 0. Hence we assume y * = a(a+x * )r 1 ρ αK 1 with the introducing of a control variable ρ ∈ (0, 1], components of the corresponding equilibrium point E * = (x * , y * ) and restricted parameter d are (2.7) Thus, we can obtain Theorem 2.1, which can guarantee that the equilibrium point E * is globally asymptotically stable.
Theorem 2.1. If the system (1.1) has a unique interior equilibrium point E * with αy * a(a+x * ) < r 1 K 1 , then the equilibrium point E * is globally asymptotically stable (hyperbolic focus or node).
At the same time, by using the Bendixson-Dulac criteria and a proper Dulac function B(x, y) = 1 xy or B(x, y) = 1 (x+c)y with a positive constant c > 0, then theorem 2.2 can also guarantee that a unique equilibrium point E * is globally asymptotically stable on account of the non-existence of closed orbits and limit cycles.
Theorem 2.2. If the system (1.1) has a unique interior equilibrium E * , and 0 < am 2 αe−m 2 < x 2 < a, then the equilibrium E * is globally asymptotically stable.
In order to verify feasibility of Theorem 2.1 and 2.2, we will give some numerical simulations. If we take r 1 = 0.6, m 1 = 0.2, m 2 = 0.1, α = 0.5, a = 1.5, K 1 = 20, e = 0.6 and ρ = 0.9, the calculation shows that the values of these parameters can satisfy with the condition of Theorem 2.1, the equilibrium point E * is a globally asymptotically stable node, which can be seen in Figure 2(a). If we take r 1 = 0.6, m 1 = 0.2, m 2 = 0.1, α = 0.5, a = 1.5, e = 0.6, d = 0.05 and x 2 = aρ(ρ = 0.9), then Theorem 2.2 is true, thus the equilibrium point E * is a globally asymptotically stable point, which can be found from Figure 2(b). In a word, the equilibrium point E * is globally asymptotically stable under certain conditions. In this section, we mainly consider existence and stability of interior equilibrium in special cases, which can ensure the occurrence of Hopf bifurcation and Bogdanov-Takens bifurcation. In order to simplify A 1 as zero or get potential Hopf bifurcation, we take some parameters in a special case as where λ ∈ (0, 1) and µ = 1 + 8λμ > 1 are two control variables for later use to scale parameters. Then the parameter d is constrained by d = 4α(µ+1)ϕ d (λ,µ) a(µ+3) 2 ψ d (λ,µ) , where auxiliary functions are It is quite clear that we can derive a little complicated expressions of a required interior equilibrium point is also an auxiliary function.
Now, we will continually use transformations (III): u = p + a 02 pq, v = q − a 20 p 2 ; (IV): p = w, q = z − c 11 wz and (V): w = x 1 + 1 2 f 02 x 2 1 , z = y 1 + f 02 x 1 y 1 to derive a standard form with discriminants d 1 and d 2 :ẋ hence it also support above conclusions and we can obtain the Theorem 2.3. Theorem 2.3. As we take the value of parameters under the condition (2.8), the system (1.1) admits an interior equilibrium point E (2) * with zero trace. (i) When µ m < µ < µ 1 , the equilibrium point E (2) * is a stable multiple focus with multiplicity one. (ii) When µ = µ 1 , the equilibrium point E (2) * is a cusp of codimension 2 (BT bifurcation point). (iii) When µ > µ 1 , the equilibrium point E (2) * is a saddle point.
Here we can take the value λ = 2 3 to get some interesting result, the first positive(meaningful) root of Eq (2.9) is µ 1 = µ 1, Solving out the Eq (2.4), we have three possible interior equilibrium point and auxiliary functions are With the techniques in Calculus, we know: √ 157 3 ≈ 12.687, the function ϕ s (µ) is negative and monotonically decreasing.
(ii) when µ = µ + s , the function ϕ s (µ) owns a negative (local) minimum. (iii) when µ > µ + s , the function ϕ s (µ) is monotonically increasing and has a unique positive root, which is denoted as µ s , where Here the equilibrium point E (2) 4 is a stable multiple focus with multiplicity one when µ < µ 1 , while it becomes a saddle point when µ > µ 1 .
Case 2. When µ = µ 1 , there exist two interior equilibrium points, including an asymptotically stable node E (1) is a cusp of codimension 2. Furthermore, it is worth noting that the two equilibrium point E (2) 4 and E (3) 4 can coincide with each other.
Case 3. When µ = µ s or Φ s (µ) = 0, there exist two interior equilibria, that is to say, above two interior equilibrium point E (1) 4 and E (3) 4 can coincide with each other (but we still denote it as E (1)  4 ), where .
< 0 and Seeing [35] in detail, the interior equilibrium E (2) 4 is still a stable multiple focus with multiplicity one.
In order to verify the feasibility of theoretical derivation, we take r 1 = 0.6 α = 0.5, a = 1.5 and e = 0.6, then some numerical simulations are implemented. Figure 3 depicts the curves of functions ϕ A 2 (λ, µ), ϕ d (λ, µ) and ψ d (λ, µ), which can show the existence of key values. When µ = 20 < µ 1 , the unique equilibrium point E (2) 4 is a stable multiple focus with multiplicity one(see Figure 4(a)). When µ = 25 > µ 1 , there exist three interior equilibrium points including a stable node E (1) 4 , a saddle point E (2)   4 and an unstable node E (3) 4 (see Figure 4(b)). When µ = µ 1 , a cusp of codimension 2 E (2) 4 and a stable node E (1) 4 will occur (see Figure 4(c)). When µ = µ s , there exist two equilibrium points including a stable multiple focus E (2) 4 with multiplicity one and a stable node E (1) 4 (see Figure 5(a)). When µ = 21.5, there exist three equilibrium points including a stable multiple focus E (2) 4 with multiplicity one, a saddle point E (3) 4 and a stable node E (1) 4 (see Figure 5(b)). In a word, there are several kinds of internal equilibrium points with different characteristics in the system (1.1).
When µ = µ 1, 5 9 , the system (1.1) owns two separate interior equilibrium points E (1) ). The equilibrium point E (1) 5 is an unstable node since  1) is apparently equivalent to a new system in the standard form of Eq (2.14) with discriminants 0. (2.21) The equilibrium point E (2) 5 is just a cusp of codimension 2 due to the non-degeneracy condition d 1 d 2 0. Hence we can obtain Theorem 2.4.
In order to verify the feasibility of theoretical derivation, we will give some numerical simulations. Figure 6(b) is the curves of functions ϕ A 2 (red) and ϕ d (blue) with λ = 5 9 , which mainly displays the threshold value of control parameters. And then, we take r 1 = 1, a = 1.5, α = 0.5 and e = 0.6, the system (1.1) exists a multiple stable focus with multiplicity one (see Figure 7(a)), a unique multiple unstable focus with multiplicity one and a limit cycle (see Figure 7(b)), a cusp of codimension 2 and an unstable node (see Figure 7(c) ), a unique stable weak focus of order 2 (see Figure 7(d)). In a word, the system (1.1) has different internal equilibrium points with the value change of key parameters.  5 and an unstable node E (1) 5 when µ = µ 1, 5 9 ; (d) A unique stable weak focus of order 2 when µ = µ σ .

Cusp of codimension 3
In the following, we will investigate a cusp of codimension 3 in the system (1.1). First of all, translating the equilibrium point E * = (x * , y * ) to the origin O via a transformation (I): Next, following the technique in above subsection and making a linear transformation (II): the above system can be rewritten as the forṁ (2.24) From the Lemma 1 in [38], the system (2.24) is equivalent to the system in standard forṁ after some nonsingular transformations in the neighborhood of O, where d 1 = b 20 , d 2 = b 11 + 2a 20 are discriminates. Solving out an degenerate equilibrium E 3 = ( a(2αe+m 2 ) αe−m 2 , 12e 2 αa 2αe+m 2 ) from the condition d 1 = 0, which also satisfies A 1 = 0, A 2 = 0 and p x = ∆ x = p y = ∆ y = 0, thus we have parameters r 1 , K 1 and d (suppose they are positive) with restrictions (2.26) These restrictions (2.26) can deduce another discriminate d 2 = (8αe+7m 2 )(m 2 −αe) 18aαe 2 0, i.e., the degeneracy condition d 1 d 2 = 0, which also yields that the equilibrium point E 3 is a cusp of codimension at least 3. Indeed, the degenerate equilibrium point E 3 is a codimension 3 Bogdanov-Takens singularity(focus or center) after some nonsingular transformations. Finally, the existence of the equilibrium point E 3 can be seen from Figure 8 in details with parameters a = 1.5, α = 0.5, m 1 = 0.6, m 2 = 0.06 and e = 0.6. At the same time, we can obtain the Theorem 2.5.
Theorem 3.1. (Hopf bifurcation curve) For the equilibrium point E (2) 4 with µ m < µ < µ 1, 2 3 (A 1 = 0 and A 2 > 0), when parameter δ varies in a small neighbourhood of the origin in parameter plane, the Hopf bifurcation curve of the system (3.1) is defined by (3.2) (notice the range of parameter w) and the approximation is a straight line with slope k in a small neighbourhood of the origin. Furthermore, the curve divides a small neighbourhood of the origin in the parameter plane into two regions I and II, in which dynamical behaviors of the system (3.1) can be exhibited.
Theorem 3.2. (Hopf bifurcation curve) For the equilibrium point E (2) 5 with µ ∈ (µ m , µ 1, 5 9 )\{µ σ }, when parameter δ varies in a small neighbourhood of the origin in parameter plane, the Hopf bifurcation curve of the system (3.1) is defined by (3.4) (notice the range of parameter w) and the approximation is a straight line with slope k in a small neighbourhood of the origin. Furthermore, the curve divides a small neighbourhood of the origin in the parameter plane into two regions I and II, in which dynamical behaviors of the system (3.1) can be exhibited.
In order to verify the feasibility of the Theorems 3.1 and 3.2, we will give some numerical simulations. For the equilibrium point E (2) 4 with r 1 = 0.6 α = 0.5, a = 1.5 , e = 0.6 and µ = 10, the Hopf bifurcation curve with w in parameter plane is and the straight line of the approximate representation of bifurcation curve H p is δ 2 ≈ −0.704575δ 1 . Figure 9(a) depicts the Hopf bifurcation curve, Figure 10(a),(b) depict phase diagrams of asymptotically stable focus and unstable focus(with a limit cycle) corresponding to δ 1 = δ 2 = 0.001 (in region I) and δ 1 = δ 2 = −0.001(in region II) respectively. For the equilibrium point E (2) 5 with r 1 = 0.6 α = 0.5, a = 1.5, e = 0.6 and µ = 8 in the case (C2), the bifurcation curve Hp is analytically formulated by such solutions of the system (3.1), i.e., H p = {δ | δ satisfies (3.1)}, which can be seen from Figure 9(b), and a Hopf bifurcation curve Hp(red) is defined by then the curve Hp divides a parameter plane into separate two regions. Figure 11

Bogdanov-Takens bifurcation of codimension 2
We firstly recall the system (1.1) with a cusp E (2) 4 of codimension 2 when parameter conditions λ = 2 3 and µ = µ 1, 2 3 hold, in other words, we can start with an unfolding systeṁ by introducing bifurcation parameters e and m 2 . Naturally, a parameter vector δ = (δ 1 , δ 2 ) is in a sufficiently small neighbourhood of the origin O in the parameter plane. By using the transformation (I): x = X + x (2) 4 , y = Y + y (2) 4 and expanding such system in a power series around the origin, it can be rewritten asẊ where b 00 (0, 0) = 0. Secondly, we will use a transformation (II): u = X, v = F 1 (X, Y), the system (3.9) can be reduced to a new systemu > 0. Then we will construct a transformation (IV): w = p, z = (1 − f 02 (δ)p)q, dt = (1 − f 02 (δ)p)dτ, one can rewrite above system as(the symbol τ is denoted as t) Similarly, we omit the expressions of coefficients h i j (δ) although they are expressed iteratively by f i j (δ), and h 20 (0, 0) = f 20 (0, 0) > 0, h 11 (0, 0) = d 11 (0, 0) < 0. That is to say, there is a small neighbourhood of the origin such that h 20 (δ) is positive and h 11 (δ) is negative when δ falls in this neighbourhood. where the symbol τ is still denoted as t, and two discriminants are This system (3.8) is indeed a generic family unfolding at the codimension 2 cusp E (2) 4 according to the rank of a Jacobian matrix or the nonzero Jacobian determinant Therefore we obtain local approximated representations of saddle-node (SN), Hopf (H) and homoclinic (HL) bifurcation curves up to second-order with slope k BT = (1+ √ 17)α 6 ≈ 0.853851α > 0 around the origin for the system (3.8) [39]. These bifurcation curves can divide the parameter plane into several regions, which can exhibit separately dynamical behaviors.
(i) The saddle-node bifurcation curve is formulated by  bifurcation parameters e and m 2 , in a small neighbourhood of the equilibrium point E (2) 4 , the system undergoes an attracting Bogdanov-Takens bifurcation of codimension 2 when the value of parameter δ varies in such sufficiently small neighbourhood of the origin. Furthermore, this system is a generic family unfolding at the cusp E (2) 4 of codimension 2. Here we will take r 1 = 0.6, a = 1.5, α = 0.5 and e = 0.6, then Figure 12 depicts the saddle-node, Hopf and homoclinic bifurcation curves in different colors, which can show the existence of critical thresholds.
(ii) When δ 1 = 0.01, δ 2 = 0 (δ lies in positive δ 1 axis) or δ falls in region I (the region between saddle-node bifurcation curve S N 2 and homoclinic bifurcation curve), there exist three interior equilibrium points, where two interior equilibrium points are bifurcated from a stable node in (viii), which can be seen from Figure 13 (iii) When (δ 1 , δ 2 ) ≈ (0.01, 4.265486 × 10 −3 ) or δ lies in homoclinic bifurcation curve, there exist three interior equilibrium points and a homoclinic loop.
(iv) When (δ 1 , δ 2 ) ≈ (0.01, 4.267370 × 10 −3 ) or δ falls in region II(the region between homoclinic bifurcation curve and Hopf bifurcation curve), there exist three interior equilibrium points, including an unstable focus, a saddle point and a stable node, which can be seen from Figure 14(a),(b).
(vi) When (δ 1 , δ 2 ) ≈ (0.01, 4.271301 × 10 −3 ) or δ falls in region III(the region between Hopf bifurcation curve and saddle-node bifurcation curve S N 1 ), there exist three interior equilibrium points, including a stable focus which is unstable in case (iv), which can be seen from Figure 15(a). However, by combining the case (iv), it can ensure the potential Hopf bifurcation, but the homoclinic loop is broken.
(viii) When (δ 1 , δ 2 ) ≈ (0.01, 8.546695 × 10 −3 ) or δ falls in region IV(the region on the left hand side of saddle-node bifurcation curve), there exists a unique stable node, which can be seen from Figure  15   However, if we choose m 2 and d as Bogdanov-Takens bifurcation parameters and rewrite the original system in the unfolding form (3.1) with µ = µ 1, 2 3 : where δ 1 and δ 2 are sufficiently small parameters and the vector δ = (δ 1 , δ 2 ) is in a small neighbourhood of the origin O as well. Following the procedures above and the values of parameters, the unfolding system (3.1) is also a generic family unfolding at the codimension 2 Bogdanov-Takens cusp E (2) 4 ac- cording to the nonzero Jacobian: > 0.
and we have representations of saddle-node, Hopf and homoclinic bifurcation curves with slope k BT = 1− √ 17 6ae ≈ − 0.520518 ae < 0 around the origin. Furthermore, it should be noticed that the slope k can be viewed as the limiting case of the slope (3.3) when µ → µ 1, 2 3 . At the same time, when δ lies on the Hopf bifurcation curve H, for instance, (δ 1 , δ 2 ) ≈ (1 × 10 −5 , −5.783264 × 10 −6 ), the Hopf bifurcation can not undergo; when δ lies on the homoclinic bifurcation curve HL, for instance, (δ 1 , δ 2 ) ≈ (1 × 10 −5 , −5.783278 × 10 −6 ), the homoclinic loop does not exist. Moreover, it should be noticed more that these two cases both occur owing to d 2 > 0 or the minus of (3.2). While the saddle-node bifurcation curve up to second-order can be formulated by For the saddle-node bifurcation curve, when δ lies on the region I (the left hande side of the SN curve), there exist three interior equilibrium points. When δ lies on the region II (the right hand side of the SN curve), there exists a unique interior equilibrium point. When δ lies on the saddle-node bifurcation curve S N 1 , there exist three interior equilibrium points. When δ lies on the saddle-node bifurcation curve S N 2 , there exists a unique interior equilibrium point. All the detailed results can be seen in the Figure 16 for saddle-node bifurcation curve in this novel phenomenon with r 1 = 0.6, a = 1.5, α = 0.5 and e = 0.6.  Similarly, for the case (C2) with λ = 5 9 and Bogdanov-Takens bifurcation parameters m 2 and d, according to a Jacobian matrix with rank 2 and (i) The saddle-node bifurcation curve is formulated by (3.20) where the half curves S N 1 (S N 2 ) is the "right" ("left") part of curve S N in the forth(second) quadrant, respectively.
(iv) When δ lies in Hopf bifurcation curve, there exist three interior equilibrium points, including an unstable node, a saddle point, and a non-hyperbolic equilibrium (focus or center) with zero trace and positive determinant, which can ensure potential Hopf bifurcation.
(vi) When δ lies in homoclinic bifurcation curve, there exist three interior equilibrium points and a homoclinic loop, including an unstable node, a saddle point and a stable focus.

Conclusions
Within the framework of predator-prey ecological dynamics, this paper mainly discusses the dynamic properties of the Bazykin's predator-prey ecosystem with Holling type II functional response and interspecific density-restricted effects on the predators, including the existence and stability of all possible equilibrium points, Hopf bifurcation and Bogdanov-Takens bifurcation.
Aiming at all possible equilibrium points of the Bazykin's predator-prey ecosystem, mathematical theory works mainly investigate the existence and stability of boundary equilibrium point, hyperbolic equilibrium point, non-hyperbolic equilibrium point and cusp of condimension 3, and then give some corresponding threshold conditions of some key parameters, for example, by introducing control variables λ and µ, the stability analysis and type of the interior equilibrium point E (2) * ( * = 4, 5) are ascertainable and clear in detail, and this equilibrium point E (2) * is a multiple focus with multiplicity one(or cusp of codimension 2 or saddle point) when µ m < µ < µ 1 (or µ = µ 1 or µ > µ 1 ). At the same time, it is easy to find from numerical simulation works that all the equilibrium points derived from theoretical derivation always exist and have corresponding stability states, which indirectly prove that the Bazykin's predator-prey ecosystem has different intrinsic dynamic properties with the value change of critical parameters, which also represents that predator population and prey population have different coexistence modes. For Hopf bifurcation and Bogdanov-Takens bifurcation, mathematical theory works mainly investigate Hopf bifurcation and Bogdanov-Takens bifurcation about the equilibrium point E (2) 4 and E (2) 5 , and give some threshold conditions of some key parameters to ensure the occurrence of Hopf bifurcation and Bogdanov-Takens bifurcation. With the help of numerical simulation, the formulated Hopf bifurcation curve when µ m < µ < µ 1 and Bogdanov-Takens bifurcation of codimension 2 when µ = µ 1 are both presented, which can directly verify the validity and feasibility of theoretical derivation, and indirectly explain that the Bazykin's predator-prey ecosystem has complex bifurcation dynamic evolution process with the value change of critical parameters, such as saddle-node, Hopf and homoclinic bifurcation. Furthermore, it is obvious to find that the Bazykin's predator-prey ecosystem can exhibit different dynamical behaviours when parameter vector (δ 1 , δ 2 ) varies in different regions in a small neighbourhood of the origin O in the parameter plane. However, it is also worth noting that the Hopf and homoclinic bifurcations do not exist if we choose m 2 and d as bifurcation parameters in (C1). Moreover, we specifically investigate the limit problem by comparing the Runge-Kutta 45 method with perturbation procedure.
In the follow-up research works, based on the research results of this paper and biological manipulation theory, we will further explore the dynamic relationship between Microcystis aeruginosa and filter-feeding fish by using bifurcation dynamic analysis and explain the biological significance of the Bazykin's predator-prey ecosystem. At the same time, due to the difference of nutrient load, biological composition and hydrodynamic conditions in different water bodies, the dynamic relationship between filter feeding fish and Microcystis aeruginosa on the basic of the Bazykin's predator-prey ecosystem still needs to be further studied in the experiment.
In summary, all results in this paper further show that the system (1.1) proposed by the paper [24] has more abundant dynamic behavior, and vigorously develop the dynamic properties of the Bazykin's predator-prey ecosystem. Furthermore, the results of bifurcation and stability can be helpful to better understand the interaction mechanism between prey population and predator population in natural real ecosystem.