Finite element approximation of a population spatial adaptation model

  • In [18], Sighesada, Kawasaki and Teramoto presented a system of partial differential equations for modeling spatial segregation of interacting species. Apart from competitive Lotka-Volterra (reaction) and population pressure (cross-diffusion) terms, a convective term modeling the populations attraction to more favorable environmental regions was included. In this article, we study numerically a modification of their convective term to take account for the notion of spatial adaptation of populations. After describing the model, in which a time non-local drift term is considered, we propose a numerical discretization in terms of a mass-preserving time semi-implicit finite element method. Finally, we provied the results of some biologically inspired numerical experiments showing qualitative differences between the original model of [18] and the model proposed in this article.

    Citation: Gonzalo Galiano, Julián Velasco. Finite element approximation of a population spatial adaptation model[J]. Mathematical Biosciences and Engineering, 2013, 10(3): 637-647. doi: 10.3934/mbe.2013.10.637

    Related Papers:

    [1] Huanyi Liu, Hengguo Yu, Chuanjun Dai, Zengling Ma, Qi Wang, Min Zhao . Dynamical analysis of an aquatic amensalism model with non-selective harvesting and Allee effect. Mathematical Biosciences and Engineering, 2021, 18(6): 8857-8882. doi: 10.3934/mbe.2021437
    [2] Shengyu Huang, Hengguo Yu, Chuanjun Dai, Zengling Ma, Qi Wang, Min Zhao . Dynamics of a harvested cyanobacteria-fish model with modified Holling type Ⅳ functional response. Mathematical Biosciences and Engineering, 2023, 20(7): 12599-12624. doi: 10.3934/mbe.2023561
    [3] Moitri Sen, Malay Banerjee, Yasuhiro Takeuchi . Influence of Allee effect in prey populations on the dynamics of two-prey-one-predator model. Mathematical Biosciences and Engineering, 2018, 15(4): 883-904. doi: 10.3934/mbe.2018040
    [4] Fang Liu, Yanfei Du . Spatiotemporal dynamics of a diffusive predator-prey model with delay and Allee effect in predator. Mathematical Biosciences and Engineering, 2023, 20(11): 19372-19400. doi: 10.3934/mbe.2023857
    [5] Zhenliang Zhu, Yuming Chen, Zhong Li, Fengde Chen . Dynamic behaviors of a Leslie-Gower model with strong Allee effect and fear effect in prey. Mathematical Biosciences and Engineering, 2023, 20(6): 10977-10999. doi: 10.3934/mbe.2023486
    [6] Xiaofen Lin, Hua Liu, Xiaotao Han, Yumei Wei . Stability and Hopf bifurcation of an SIR epidemic model with density-dependent transmission and Allee effect. Mathematical Biosciences and Engineering, 2023, 20(2): 2750-2775. doi: 10.3934/mbe.2023129
    [7] Kunlun Huang, Xintian Jia, Cuiping Li . Analysis of modified Holling-Tanner model with strong Allee effect. Mathematical Biosciences and Engineering, 2023, 20(8): 15524-15543. doi: 10.3934/mbe.2023693
    [8] Minjuan Gao, Lijuan Chen, Fengde Chen . Dynamical analysis of a discrete two-patch model with the Allee effect and nonlinear dispersal. Mathematical Biosciences and Engineering, 2024, 21(4): 5499-5520. doi: 10.3934/mbe.2024242
    [9] Yue Xia, Lijuan Chen, Vaibhava Srivastava, Rana D. Parshad . Stability and bifurcation analysis of a two-patch model with the Allee effect and dispersal. Mathematical Biosciences and Engineering, 2023, 20(11): 19781-19807. doi: 10.3934/mbe.2023876
    [10] Manoj K. Singh, Brajesh K. Singh, Poonam, Carlo Cattani . Under nonlinear prey-harvesting, effect of strong Allee effect on the dynamics of a modified Leslie-Gower predator-prey model. Mathematical Biosciences and Engineering, 2023, 20(6): 9625-9644. doi: 10.3934/mbe.2023422
  • In [18], Sighesada, Kawasaki and Teramoto presented a system of partial differential equations for modeling spatial segregation of interacting species. Apart from competitive Lotka-Volterra (reaction) and population pressure (cross-diffusion) terms, a convective term modeling the populations attraction to more favorable environmental regions was included. In this article, we study numerically a modification of their convective term to take account for the notion of spatial adaptation of populations. After describing the model, in which a time non-local drift term is considered, we propose a numerical discretization in terms of a mass-preserving time semi-implicit finite element method. Finally, we provied the results of some biologically inspired numerical experiments showing qualitative differences between the original model of [18] and the model proposed in this article.


    Nowadays, great quantities of rivers and lakes around the world are more or less affected by algal bloom, which phenomenon has received more and more attention from aquatic ecologists and environmentalists[1,2]. Algae are the simplest vegetative organisms capable of photosynthesis, they have various kinds and are widely distributed. The outbreak of algal bloom will bring a lot of troubles to our life and environment. For examples, some algaes can secrete toxic substances, which can not be easily removed by conventional water treatment processes, then directly affect the quality and safety of water supply. A large number of algae populations will also consume too much oxygen, causing fish and other aquatic organisms to die due to the lackness of oxygen, which also further deteriorates water quality [3,4]. However, until now, the outbreak and control of algal bloom is still one of the world's difficult problems. Therefore, it is urgent for scientists to choose new ideas and methods to explore this world problem.

    Since the 1950s, as system models and mathematical methods have penetrated into the field of ecological research, a research method combining mathematical model and numerical simulation is formed, which can not only dynamically reflect the internal essential characteristics of real ecological problems, but also analyze their causality and find out the laws reflecting the internal mechanism according to their deep understanding[5,6,7]. Consequently, more and more scholars use mathematical models to describe ecological phenomena and use numerical simulation to predict the future trend, and then some excellent research results were obtained in these papers[8,9]. The papers[9,10,11,12,13,14] introduced refuge effect into predator-prey model to explore its impact on dynamic behaviors, and gave some meaningful results, such as: refuge effect has a crucial role in the evolution of dynamic behaviors[9,10,11], prey refuge can balance the relationship of predator-prey under the presence of habitat complexity[12], prey refuge can affect spatiotemporal patterns[13], Both nonlinear harvesting and refuge will influence the Turing instability [14]. the papers[15,16,17,18,19,20,21] added Allee effect to predator-prey model to explore its impact mechanism, and some good results were obtained, such as the system without Allee effect will always take a shorter time to get a steady-state solution than that with Allee effect[15], Allee effect can create or destroy the interior attractor [16], Allee effect in prey growth can reduce the complex dynamics[17], Allee effect will cause strong impact on population dynamics[18,19,20,21]. The papers[22,23,24] considered how harvest factors affect the dynamic relationship of predator-prey model, and obtained some important results. The paper[25] pointed out that different level of fear will induce differently, such as striped inhomogeneous distribution is induced under high level of fear, spotted inhomogeneous distribution is induced under low level of fear, the mixture of spotted and striped inhomogeneous distribution is induced under intermediate level of fear. The paper[8] declared that cross-diffusion was able to create stationary patterns, which could enrich the findings of pattern formation in an ecosystem. All in all, it is successful to explore the action mechanism of eco-environmental factors in predator-prey model by using mathematical model and numerical simulation. Therefore, we believe that this method is particularly useful to study the outbreak and control of algal bloom.

    From the perspective of fish controlling cyanobacteria bloom, the most important aim of the paper is to propose a modified algae and fish model based on an aquatic ecological model, which can make the dynamic relationship between algae and fish represented by the modified model closer to the dynamic relationship in the natural ecological environment. On this basis, the dynamic behavior of the modified aquatic ecological model is investigated in detail especially bifurcations, and some practical results are obtained.

    The paper[26] proposed an ecological model of aquatic organisms to express the algal aggregation based on $ 5 $ modeling assumptions, and revealed the dynamic relationship between algae and fish, which can be described by the following differential equations:

    $ {˙X=r1X(1Xk)Yα(Xm)a+Xm,˙Y=ηYα(Xm)a+Xm+r2Y(1Xk)dY,
    $
    (2.1)

    where $ X(t) $ and $ Y(t) $ represent the population densities of algae and fish respectively, the ecological significance of all parameters are detailed in the paper [26].

    However, the aquatic ecological model (2.1) has two shortcomings. One is that the Allee effect is not considered in the growth function of algae. As far as we know, since the population size of microbial species is often "huge" and most microorganisms reproduce asexually, there are few researches on microorganisms with Allee effect. However, some scholars have found the existing evidences of Allee effect in microbial population experiments. For example, Smith et al. [27] engineered Allee effects in Escherichia coli, and Qi et al. [28] detected Allee effect in experimental population of Vibrio fischeri in a specific experimental environment. The causes of Allee effect are various and different for different species. Compared with algae monomer, algal aggregation has some obvious survival advantages [29,30,31]. There are gaps between algal aggregation, which will increase buoyancy then promote photosynthesis better [29]. Algal aggregation is usually composed of thousands of cells, which is large enough that can effectively prevent zooplankton from preying, and better avoid the invasion of viruses, bacteria and algae phagocytes [30]. Algal aggregation can secrete more microcystins than monomer, which reduces the risk of predation [31]. Thus algal aggregation can make algae to find a more suitable niche, then promote the growth and reproduction. The individual reproduction of algae monomer with low aggregation degree will be relatively reduced, which will result in Allee effect. The other is that the function $ r_{2}Y(1-\frac{X}{k}) $ of algae monomer affecting the abundance of fish is not reasonable, this is because that the function includes the influence of algal aggregation. Thus, relatively speaking, the function $ r_{2}Y(1-\frac{X-M}{K}) $ is more reasonable.

    Based on the above analysis, we will construct a modified aquatic ecological model, which can be described by the following differential equations:

    $ {dXdT=R1X(1XK)(XN1)Yα1(XM)A+XM,dYdT=β1Yα1(XM)A+XM+R2Y(1XMK)DY,
    $
    (2.2)

    where $ N $ $ \left(0 < N < K\right) $ is Allee effect threshold, the prey population is doomed to extinction when the prey population density or size is below the threshold. It is worth mentioning that the term of feeding on the alternative prey $ {R_2}Y\left({1 - \frac{{X - M}}{K}} \right) $ in model (2.2) is more reasonable than it in model (2.1), this is because the aggregated part of algae $ M $ is always composed by thousands of units, which is too big for fish, hence fish can only turn to the monomer part of algae $ X-M $ for grazing.

    There are too many parameters of model (2.2), which will bring inconvenience to our follow-up research, in order to reduce the number of parameters, replacing it with the following variables

    $ x = \frac{X}{A}, \; y = \frac{{\alpha _1} Y}{A{R_1}}, \; t = {R_1}T, $

    then model (2.2) can be rewritten as the following dimensionless form:

    $ {dxdt=x(1xk)(xn1)xm1+xmy,dydt=βxm1+xmy+ry(1xmk)dy,
    $
    (2.3)

    where $ m = \frac{M}{A}, \; k = \frac{K}{A}, \; n = \frac{N}{A}, \; \beta = \frac{{\alpha _1}{\beta _1}}{R_1}, \; r = \frac{{{R_2}}}{{{R_1}}} $ and $ d = \frac{D}{{{R_1}}} $ are positive constants and apparently $ 0\leqslant m \leqslant k $. It is worth noting that the density of prey $ X $ must greater than the aggregation parameter $ M $, therefore $ x $ must satisfy $ m \leqslant x $ and $ m $ must be euqal to zero when $ x = 0 $. To illustrate how Allee effect affects the per-capita growth rate of algae population, we numerically show the change trend of per-capita growth rate under different values of Allee effect control parameter $ n $. It can be seen from Figure 1 that if the Allee effect is not considered, the per-capita growth rate of algae population decreases monotonically with the increase of $ x $ value, if the Allee effect is considered, the per-capita growth rate of algae population increases first and then decreases monotonically with the increase of $ x $ value, and if the algae population density is lower than an Allee threshold, the per-capita growth rate is negative, which implies that the algae population density can not support the large-scale reproduction of algae population, and will lead to the extinction of algae population. At the same time, it is also worth pointing out that the per-capita growth rate curve is concave, which can better describe the proliferation stage of algae population. Thus, considering Allee effect in the construction of algae ecological dynamics model is more in line with the actual ecological and environmental significance.

    Figure 1.  Per-capita growth rate of algae under the influence of different Allee effect parameters $ n $ without predator. Other parameters are as follows: $ k = 5 $, $ m = 0.5 $, $ \beta = 0.17 $, $ r = 0.15 $, $ d = 0.21 $.

    In this paper, firstly, we analyze the relevant dynamic properties of model (2.3), and obtain several conditions of critical threshold to ensure that model (2.3) can undergo transcritical, saddle-node, Hopf and B-T bifurcation. Secondly, we carry out numerical experiments on the dynamic behaviors of model (2.3), which can confirm the effectiveness of theoretical derivation and visually display the dynamic behavior evolution process. Finally, we reveal the quality of dynamic relationship between algae and fish from the view of population dynamics evolution according to results of the numerical simulation.

    Equilibrium points are the special solutions, which will exhibit rich properties of model (2.3), therefore, the existence and stability of all possible equilibrium points of model (2.3) will be discussed in this section, and we will also use the $ Poincare-Bendixson $ $ Theorem $ to confirm the existece of a limit cycle.

    It is obvious that model (2.3) has at most five possible equilibrium points: trivial extinction equilibrium point $ {E_0}(0, 0) $, two predator-free equilibrium points $ {E_1}\left({{k}, 0} \right) $ and $ {E_2}\left({n, 0} \right) $, two coexistence equilibrium points $ {E_1^*}\left({{x_1}, {y_1}} \right) $ and $ {E_2^*}\left({{x_2}, {y_2}} \right) $. The equilibrium points $ {E_1} $ and $ {E_2} $ exist unconditionally, and $ {E_0} $ exists when $ m \ne 1 $ from the mathematical perspective, while from the biological perspective it exists without condition. Where $ {x_1} $ and $ x_2 $ are the roots of the equation:

    $ - \frac{r}{k}{x^2} + Bx + C = 0, $

    with

    $ B = \beta - d + r + \frac{{2mr}}{k} - \frac{r}{k}, \; C = dm - d - m\beta + r - mr + \frac{{mr}}{k} - \frac{{{m^2}r}}{k}, $

    then we can obtain the concrete expressions of $ {x_1} $ and $ x_2 $ as

    $ {x_1} = \frac{{kB - k\sqrt \Delta }}{{2r}}, {x_2} = \frac{{kB + k\sqrt \Delta }}{{2r}}, $

    where

    $ \Delta = {B^2} + \frac{{4r}}{k}C = \frac{{{{\left( {d - r - \beta } \right)}^2}{k^2} - 2r\left( {d - r + \beta } \right)k + {r^2}}}{{{k^2}}}, $

    corresponding, the concrete expressions of $ {y_1} $ and $ y_2 $ are

    $ {y_1} = \frac{{{x_1}\left( {1 - \frac{{{x_1}}}{k}} \right)\left( {\frac{{{x_1}}}{n} - 1} \right)\left( {1 + {x_1} - m} \right)}}{{{x_1} - m}}, {y_2} = \frac{{{x_2}\left( {1 - \frac{{{x_2}}}{k}} \right)\left( {\frac{{{x_2}}}{n} - 1} \right)\left( {1 + {x_2} - m} \right)}}{{{x_2} - m}}. $

    The existence of coexistence equilibrium points is conditional, and the condition $ m \leqslant x_i \leqslant k, \; i = 1, 2 $ must be satisfied, this is out of the consideration of biological significance. It is obvious that $ y_i > 0 $ when $ x_i > max\left\{ {m, n} \right\} $. The specific existence conditions of coexistence equilibrium points can be viewed in the appendix. Coexistence equilibrium points have a decisive influence on the dynamic behavior, for an example, the fish will finally die out and model (2.3) will not be persistent if there exists no coexistence equilibrium point.

    The stability of equilibrium points can be obtained through the signs of the eigenvalues of Jacobian matrix, firstly we obtain the expression of Jacobian matrix about model (2.3) as

    $ {J_{E\left( {x, y} \right)}} = \left( {3nkx2+(2n+2k)x1y(1+xm)2xm1+xmβy(1+xm)2rykβxm1+xm+r(1xmk)d
    } \right), $

    then we have the following theorems about the stability of equilibrium points.

    Theorem 3.1. $ {E_0}\left(0, 0\right) $ exists when $ m \ne 1 $ and it is stable when $ r < d $, while $ {E_0} $ is a saddle when $ r > d $.

    Proof. When $ x = 0 $, $ m $ must be equal to zero, therefore the Jacobian matrix of $ {E_0} $ can be written as

    $ {J_{{E_0}(0, 0)}} = \left( {100rd
    } \right), $

    it is obvious that both the two eigenvalues of matrix $ {J_{{E_0}(0, 0)}} $, $ {\lambda _1} = -1 $ and $ {\lambda _2} = r-d $, are negative when $ r < d $, then $ E_0 $ is stable. However the two eigenvalues of $ {J_{{E_0}(0, 0)}} $ have opposite signs when $ r > d $, then $ E_0 $ is an unstable saddle.

    Theorem 3.2. $ E_1\left(k, 0\right) $ exists unconditionally, and $ {E_1} $ is stable when $ d > \beta \frac{{k - m}}{{1 + k - m}} + \frac{m}{k}r $, but $ {E_1} $ is an unstable saddle when $ d < \beta \frac{{k - m}}{{1 + k - m}} + \frac{m}{k}r $.

    Proof. The Jacobian matrix of $ {E_1} $ can be written as

    $ {J_{{E_1}\left( {k, 0} \right)}} = \left( {1knkm1+km0βkm1+km+mkrd
    } \right), $

    the two characteristic roots of matrix $ {J_{{E_1}(k, 0)}} $ are $ {\lambda _1} = {1 - \frac{k}{n}} $ and $ {\lambda _2} = {\beta \frac{{k - m}}{{1 + k - m}} + \frac{m}{k}r - d} $, the former is negative since the Allee effect threshold $ N $ is smaller than the maximum environmental capacity $ K $, then $ n < k $ is satisfied. Therefore, when $ d > \beta \frac{{k - m}}{{1 + k - m}} + \frac{m}{k}r $, the two eigenvalues of $ {J_{{E_0}(0, 0)}} $ are both negative, then $ {E_0} $ is stable, but $ {E_0} $ is an unstable saddle as $ d < \beta \frac{{k - m}}{{1 + k - m}} + \frac{m}{k}r $.

    Theorem 3.3. $ {E_2} $ allways exists as an unstable equilibrium point. When $ d < \beta \frac{{n - m}}{{1 + n - m}} + r\left({1 - \frac{{n - m}}{k}} \right) $, $ {E_2} $ is an unstable node or focus, otherwise $ {E_2} $ is an unstable saddle if $ d > \beta \frac{{n - m}}{{1 + n - m}} + r\left({1 - \frac{{n - m}}{k}} \right) $.

    Proof. The expression of the Jacobian matrix around the equilibrium point $ {E_2} $ is given by:

    $ {J_{{E_2}\left( {n, 0} \right)}} = \left( {1nknm1+nm0βnm1+nm+r(1nmk)d
    } \right). $

    Obviously, $ {J_{{E_2}(n, 0)}} $ has two characteristic roots $ {\lambda _1} = {1 - \frac{n}{k}} $ and $ {\lambda _2} = \beta \frac{{n - m}}{{1 + n - m}} + r\left({1 - \frac{{n - m}}{k}} \right) -d $. From the previous content that $ n < k $ then the root $ \lambda_1 $ is positive, which means $ E_2 $ is always unstable. Moreover, $ E_2 $ is an unstable node or focus when $ d < \beta \frac{{n - m}}{{1 + n - m}} + r\left({1 - \frac{{n - m}}{k}} \right) $, and a saddle when $ d > \beta \frac{{n - m}}{{1 + n - m}} + r\left({1 - \frac{{n - m}}{k}} \right) $.

    Next we will focus on the stability of the internal equilibrium points $ E_1^* $ and $ E_2^* $, the Jacobian matrixs of them can be written as

    $ {J_{E_i^*\left( {{x_i}, {y_i}} \right)}} = \left( {3nkxi2+(2n+2k)xi1yi(1+xim)2xim1+ximβyi(1+xim)2ryik0
    } \right), $

    the expression of characteristic equations of $ {J_{E_i^*\left({{x_i}, {y_i}} \right)}} $ is

    $ {\lambda ^2} - \left[ { - \frac{3}{{nk}}{x_i}^2 + \left( {\frac{2}{n} + \frac{2}{k}} \right){x_i} - 1 - \frac{{{y_i}}}{{{{\left( {1 + {x_i} - m} \right)}^2}}}} \right]\lambda + {y_i}\left[ {\frac{\beta }{{{{\left( {1 + {x_i} - m} \right)}^2}}} - \frac{r}{k}} \right]\frac{{{x_i} - m}}{{1 + {x_i} - m}} = 0, $

    where $ i = 1, 2 $, then we have the following two theorems, which are under the conditions of their existence.

    Theorem 3.4. Suppose $ {E_1^*}\left(x_1, y_1\right) $ exists, then $ {E_1^*} $ is an unstable saddle as long as $ r < \frac{k}{{k + 1}}\left({d - \beta + \sqrt \Delta } \right) $. As for $ r > \frac{k}{{k + 1}}\left({d - \beta + \sqrt \Delta } \right) $, $ {E_1^*} $ is asymptotically stable in a small domain when $ Tr\left(J_{E_1^*}\right) < 0 $, while $ {E_1^*} $ is an unstable node or focus when $ Tr\left(J_{E_1^*}\right) > 0 $.

    Proof. It is not hard to get the following equation through calculation

    $ {S_1} = {\frac{\beta }{{{{\left( {1 + {x_1} - m} \right)}^2}}} - \frac{r}{k}}{\text{ = }}\frac{{ - 2r\sqrt \Delta }}{{\sqrt \Delta k + \left( {d - r - \beta } \right)k - r}}, $

    then for the stability of equilibrium point $ E_1^* $, we have the following two situations

    (1) If $ 0 < r < \frac{k}{{k + 1}}\left({d - \beta + \sqrt \Delta } \right) $, then $ {S_1} < 0 $, which means $ Det\left(E_1^*\right) < 0 $ i.e., the two eigenvalues of $ J_{E_1^*} $ have opposite signs, therefore $ {E_1^*} $ is an unstable saddle.

    (2) If $ r > \frac{k}{{k + 1}}\left({d - \beta + \sqrt \Delta } \right) $, then $ {S_1} > 0 $, which means $ Det\left(E_1^*\right) > 0 $. Under this condition, $ {E_1^*} $ is asymptotically stable in a small domain when $ Tr\left(J_{E_1^*}\right) < 0 $ i.e., $ J_{E_1^*} $ has two negative eigenvalues or the real part of eigenvalues are negative, but $ {E_1^*} $ is an unstable focus or node when $ Tr\left(J_{E_1^*}\right) > 0 $ i.e., both the two eigenvalues of $ J_{E_1^*} $ are positive or have positive real part.

    Theorem 3.5. The equilibrium point $ {E_2^*} $ is an unstable saddle as long as it exists.

    Proof. Similar to the above theorem, we can obtain the following equation

    $ {S_2} = {\frac{\beta }{{{{\left( {1 + {x_2} - m} \right)}^2}}} - \frac{r}{k}} = \frac{{ - 2r\sqrt \Delta }}{{\sqrt \Delta k + \left( { - d{\text{ + }}r{\text{ + }}\beta } \right)k{\text{ + }}r}}. $

    It's obvious that $ {S_2} < 0 $ due to the existing condition of $ E_2^* $ including $ r > d $, which means $ Det\left(E_2^*\right) < 0 $, then $ J_{E_2^*} $ has a positive eigenvalue and a negative eigenvalue i.e., $ {E_2^*} $ is an unstable saddle under the existing conditions.

    In this subsection, we assume that $ {E_1^*} $ is the unique internal equilibrium point, which is unstable, then as for model (2.3) we have the following main theorem about the existence of limit cycle.

    Theorem 3.6. For model (2.3), there exists one limit cycle when $ m > n $ at least.

    Proof. Consider the lines $ {L_1}:x - m = 0 $, $ {L_2}:x - k = 0 $, we have

    $ \frac{{d{L_1}}}{{dt}} = {\left. {x\left( {1 - \frac{x}{k}} \right)\left( {\frac{x}{n} - 1} \right) - \frac{{ {x - m} }}{{1 + x - m}}y} \right|_{x = m}} > 0, $
    $ \frac{{d{L_2}}}{{dt}} = {\left. {x\left( {1 - \frac{x}{k}} \right)\left( {\frac{x}{n} - 1} \right) - \frac{{ {x - m} }}{{1 + x - m}}y} \right|_{x = k}} < 0. $

    Therefore the orbit of model (2.3) will across the lines $ {L_1} $ and $ {L_2} $ from left and right, respectively. Moreover, we define the second line $ {L_3}:x + \frac{y}{\beta } - B = 0 $, then we have

    $ \frac{{d{L_3}}}{{dt}} = x\left( {1 - \frac{x}{k}} \right)\left( {\frac{x}{n} - 1} \right) - rx\left( {1 - \frac{{x - m}}{k}} \right) + dx + \left[ {r - \frac{r}{k}\left( {x - m} \right) - d} \right]B, $

    here $ B $ is a positive constant and large enough. According to the existence of $ E_1^* $ from appendix we have $ r < d $, then $ \left[{r- \frac{r}{k}\left({x - m} \right) - d} \right]B $ is a quite small negative number for any specific $ x \in [m, k] $, meanwhile $ x\left({1 - \frac{x}{k}} \right)\left({\frac{x}{n} - 1} \right) - rx\left({1 - \frac{{x - m}}{k}} \right) + dx $ is bounded. Thus $ \frac{{d{L_3}}}{{dt}} < 0 $ for all $ x \in [m, k] $, which means that model (2.3) has at least one limit cycle according to $ Poincare-Bendixson $ $ Theorem $.

    In order to make the results more visualized, we fix the parameters $ k = 5, n = 0.2, r = 0.15, d = 0.21, m = 0.4 $ and $ \beta = 0.25 $, then $ E_1^* $ is the unique equilibrium point of model (2.3), which is unstable. At the same time, the condition $ \frac{{d{L_3}}}{{dt}} < 0 $ is satisfied when $ B = 120 $, so it is obvious to find from Figure 2 that there exists one limit cycle at least. The condition $ \frac{{d{L_1}}}{{dt}} > 0 $ means that the density values of the two populations will tend to the right side of the line $ L_1 $ when the initial density values lie on the line, $ \frac{{d{L_2}}}{{dt}} < 0 $ and $ \frac{{d{L_3}}}{{dt}} < 0 $ denote similar meanings. Furthermore, the limit cycle describes such a phenomenon in biology that neither of algae and fish will be extinct, but reach a state of periodic oscillation and dynamic coexistence.

    Figure 2.  The limit cycle of this model with k = 5, n = 0.2, r = 0.15, d = 0.21, m = 0.4, $ \beta $ = 0.25, B = 120, where $ E_1^* $ is the unique internal equilibrium point in the region and the lines L1–L3 are defined in the main text.

    In this section, the local bifurcation of model (2.3) will be discussed in detail. We not only consider the codimension one bifurcations, scuh as the transcritical, saddle-node and Hopf bifurcation but also explore codimension two bifurcation as B-T bifurcation. Here the transversality conditions for transcritical bifurcation and saddle-node bifurcation will be verified by the Sotomayor's theorem[32].

    Usually, the boundary equilibrium points are the main research object of transcritical bifurcation. In this subsection, the existence of transcritical bifurcation at the equilibrium point $ {E_1}(k, 0) $ is studied. From Theorem 2, we know that $ {E_1} $ is an unstable saddle when $ g\left(m \right) < 0 $, while $ {E_1} $ is stable when $ g\left(m \right) > 0 $, where

    $ g\left( m \right) = r{m^2} + \left( {\beta k - dk - rk - r} \right)m - \beta {k^2} + d{k^2} + dk, $

    letting $ {m_1} \leqslant {m_2} $ be the two possible solutions of equition $ g\left(m \right) = 0 $, we have

    $ {m_1} = \frac{{r + \left( {d + r - \beta } \right)k - k\sqrt \Delta }}{{2r}}, {m_2} = \frac{{r + \left( {d + r - \beta } \right)k + k\sqrt \Delta }}{{2r}}. $

    Thus, when $ {m_1} $ or $ {m_2} $ is in the set $ [0, k] $, the equilibrium point $ {E_1} $ will translate its stability as the value of $ m $ passes through $ {m_1} $ or $ {m_2} $. Furthermore, it should be attention that when $ m = {m_1} $ or $ m = {m_2} $, the equilibrium point $ {E_1} $ coincides with $ {E_1^*} $ if $ \beta k-dk-rk+2mr-r > 0 $, and coincides with $ {E_2^*} $ if $ \beta k-dk-rk+2mr-r < 0 $.

    Theorem 4.1. Model (2.3) will undergo two transcritical bifurcations.

    (a) Transcritical bifurcation takes place in model (2.3) at $ {E_1} $ when $ m = {m_1} $ with either of the following two cases is satisfied.

    case 1. $ \Delta > 0 $, $ \left({\beta - d} \right)k < d $, $ \left| {r + dk - \beta k} \right| < rk $.

    case 2. $ \Delta > 0 $, $ \left({\beta - d} \right)k < d $, $ r + dk - \beta k > rk $, $ d < r $.

    (b) Transcritical bifurcation takes place in model (2.3) at $ {E_1} $ when $ m = {m_2} $ with either of the following two cases is satisfied.

    case 1. $ \Delta > 0 $, $ r < d $, $ \left| {r + dk - \beta k} \right| < rk $.

    case 2. $ \Delta > 0 $, $ r < d < \left({\beta - d} \right)k $, $ r + dk - \beta k < - rk $.

    Proof. (a) $ {m_1} \in (0, k) $, when case 1 or case 2 is satisfied. And $ {E_1} $ shifts it's stability from stable to unstable with $ m $ passing through $ {m_1} $ from left to right. When $ m = {m_1} $, the Jacobion matrix at $ {E_1} $ can be expressed as

    $ {J_{\left( {{E_1};{m}} \right)}} = \left( {1knkm1+km00
    } \right), $

    we assume that $ V $ and $ W $ respectively are eigenvectors of $ {J_{\left({{E_1};m} \right)}} $ and $ J_{\left({{E_1};m} \right)}^T $ with respect to eigenvalue zero, which means.

    $ J_{_{\left( {{E_1};m} \right)}}V = 0, \; \; J_{\left( {{E_1};m} \right)}^TW = 0, $

    we can set

    $ V = \left( {v1v2
    } \right) = \left( {km1+km1kn
    } \right), \; \; W = \left( {w1w2
    } \right) = \left( {01
    } \right), $

    then

    $ {W^T}{F_m}({E_1};m) = \left( {01
    } \right){\left( {y(1+xm)2βy(1+xm)2+ryk
    } \right)_{\left( {{E_1}, m} \right)}} = \left( {01
    } \right)\left( {00
    } \right) = 0, $
    $ WT[DFm(E1;m)V]=(01)(2y(1+xm)31(1+xm)22βy(1+xm)3β(1+xm)2+rk)(E1,m)(km11+km11kn)=2r(nk)[k2Δ((drβ)kr)kΔ]nk[(drβ)krkΔ]2,
    $
    $ WT[D2F(E1;m)(V,V)]=(01)(f1xxf1xyf1yxf1yyf2xxf2xyf2yxf2yy)(E1;m)(v1v1v1v2v2v1v2v2)=4r(kn)[(drβ)k+rkΔ][k2Δ((drβ)kr)kΔ]nk[(drβ)krkΔ]3,
    $

    where

    $ {f^1} = x\left( {1 - \frac{x}{k}} \right)\left( {\frac{x}{n} - 1} \right) - \frac{{ {x - m}}}{{1 + x - m}}y, \; \; {f^2} = \beta \frac{{ {x - m}}}{{1 + x - m}}y + ry\left( {1 - \frac{{x - m}}{k}} \right) - dy, $

    due to $ \beta \ne 0 $ and $ n \ne k $, it's obvious that

    $ WT[DFm(E1;m)V]0,WT[D2F(E1;m)(V,V)]0.
    $

    Therefore, by the $ Sotomayor's $ theorem, transcritical bifurcation takes place in model (2.3) at $ {E_1} $ when $ m = {m_1} $.

    (b) We omitted the process since it is similiar to the proof of (a).

    The existing conditions of the two internal equilibrium points $ {E_1^*} $ and $ {E_2^*} $ are given in the appendix. When $ \beta $ is chosen as the bifurcation parameter, the collision of $ {E_1^*} $ and $ {E_2^*} $ may overlap as an equilibrium point $ {E_{SN}}\left({x_{SN}}, {y_{SN}} \right) $ when $ \Delta = 0 $. With the change of the value of $ \beta $, the value of $ \Delta $ will change, and when $ \Delta < 0 $, there is no internal equilibrium point $ E_{SN} $. The change of equilibrium point number is caused by the taking place of saddle-node bifurcation in model (2.3) when $ \beta = {\beta_{SN}} $, where

    $ {\beta _{SN}} = \frac{{dk - rk + r + 2\sqrt {dkr - k{r^2}} }}{k}, $

    corresponding, we have

    $ {x_{SN}} = \frac{{mr + \sqrt {rk\left( {d - r} \right)} }}{r}, $
    $ {y_{SN}} = \frac{{{x_{SN}}\left( {1 - \frac{{{x_{SN}}}}{k}} \right)\left( {\frac{{{x_{SN}}}}{n} - 1} \right)\left( {1 + {x_{SN}} - m} \right)}}{{{x_{SN}} - m}}. $

    Theorem 4.2. Saddle-node bifurcation takes place in model (2.3) when $ \beta = {\beta _{SN}} $ under the conditions of

    (1) $ d > r $,

    (2) $ \left({n - m} \right)r < \sqrt {rk\left({d - r} \right)} < \left({k - m} \right)r $.

    Proof. The equilibrium point $ {E_{SN}} $ exists under the above two conditions according to the Appendix. The Jacobian matrix at $ {E_{SN}} $ when $ \beta = {\beta _{SN}} $ can be written as

    $ {J_{{E_{SN}}}} = \left( {3nkx2SN+(2n+2k)xSN1ySN(1+xSNm)2xSNm1+xSNm00
    } \right). $

    Letting the eigenvectors of the zero eigenvalues of $ {J_{{E_{SN}}}} $ and $ {J_{{E_{SN}}}^T} $ are $ V $ and $ W $ respectively, where

    $ V = \left( {v1v2
    } \right) = \left( {xSNm1+xSNm3nkx2SN+(2n+2k)xSN1ySN(1+xSNm)2
    } \right), W = \left( {w1w2
    } \right) = \left( {01
    } \right), $

    then, we can get

    $ {W^T}{F_\beta }\left( {{E_{SN}};{\beta _{SN}}} \right) = \left( {01
    } \right)\left( {0xSNm1+xSNmySN
    } \right) = \frac{{{x_{SN}} - m}}{{1 + {x_{SN}} - m}}{y_{SN}} \ne 0, $
    $ WT[D2Fm(ESN;βSN)(V,V)]=(01)(f1xxf1xyf1yxf1yyf2xxf2xyf2yxf2yy)(ESN;βSN)(v1v1v1v2v2v1v2v2)=2βSNySN(mxSN)2(1+xSNm)50.
    $

    Clearly, the transversality condition for the taking place of saddle-node bifurcation at $ {E_{SN}} $ is satisfied when $ \beta = {\beta _{SN}} $. Therefore, it is obvious that the number of internal equilibrium point of model (2.3) changes from zero to two when the value of parameter $ \beta $ passes through $ \beta = {\beta _{SN}} $.

    From the analysis in the previous content, it's easy to conclude that the equilibrium point $ {E_1^*} $ has different stability under different restrictions of parameters, which may caused by Hopf bifurcation. In order to figure out how algal aggregation and Allee effect influence the dynamic behavior of model (2.3), $ m $ and $ n $ are chosen as the control parameter of Hopf bifurcation respectively, then we have the following two Theorems.

    Theorem 4.3. Hopf bifurcation takes place in model (2.3) around $ {E_1^*} $ at $ m = {m_{Hp}} $ when $ r > \frac{k}{{k + 1}}\left({d - \beta + \sqrt \Delta } \right) $ based on Theorem 3.4.

    Proof. As for matrix $ {J_{E_1^*}} $, the characteristic equation of it can be written as $ {\lambda ^2} - Tr\left({{J_{E_1^*}}} \right)\lambda + Det\left({{J_{E_1^*}}} \right) = 0 $, then a Hopf bifurcation takes place when $ m = {m_{Hp}} $ such that

    (1) $ Tr\left({{J_{E_1^*}}} \right) = 0 $,

    (2) $ Det\left({{J_{E_1^*}}} \right) > 0 $,

    (3) $ \frac{d}{{dm}}{\left. {Tr\left({{J_{E_1^*}}} \right)} \right|_{m = {m_{Hp}}}} \ne 0 $.

    When $ m = {m_{Hp}} $, $ Tr\left({{J_{E_1^*}}} \right) = 0 $ is set up, and $ Det\left({{J_{E_1^*}}} \right) > 0 $ is satisfied when $ r > \frac{k}{{k + 1}}\left({d - \beta + \sqrt \Delta } \right) $ according to Theorem 3.4. Therefore we only need to certify the transersality condition (3) to guarantee the changes of stability of $ E_1^* $ through Hopf bifurcation.

    $ ddmTr(JE1)|m=mHp={6x12n2knk+x1(1x1k)(x1n1)+[x31nk(1n+1k)x21+x1](2+4x14m)(x1m)2(1+x1m)2}|m=mHp,
    $

    The condition (3) is satisfied through our numerical simulation, then Hopf bifurcation takes place in model (2.3) at $ m = {m_{Hp}} $. To find out the stability of the limit cycle brought by Hopf bifurcation, the first Lyapunov number $ {l_1} $ at the equilibrium point $ {E_1^*} $ is going to be computed following.

    Firstly, translating the equilibrium point $ {E_1^*} $ to the origin by using the transformation $ x = {x_m}+{x_1} $, $ y = {y_m}+{y_1} $, then model (2.3) can be rewritten as

    $ {˙xm=a10xm+a01ym+a20x2m+a11xmym+a02y2m+a30x3m+a21x2mym+a12xmy2m+a03y3m+P1(xm,ym),˙ym=b10xm+b01ym+b20x2m+b11xmym+b02y2m+b30x3m+b21x2mym+b12xmy2m+b03y3m+P2(xm,ym),
    $

    where $ a_{10} $, $ a_{01} $, $ b_{10} $, $ b_{01} $ are the components of the Jacobian matrix at $ E_1^* $, we have $ a_{10}+b_{01} = 0 $ and $ Det = a_{10}b_{01}-a_{01}b_{10} > 0 $. The coefficients $ a_{ij} $ and $ b_{ij} $ are determined by

    $ a02=a12=a03=b01=b02=b03=b12=0,a10=3nkx21+2(1n+1k)x11y1(1+x1m)2,a01=x1m1+x1m,a20=3nkx1+1n+1k+y1(1+x1m)3,a11=1(1+x1m)2,a30=1nky1(1+x1m)4,a21=1(1+x1m)3,b10=βy1(1+x1m)2rky1,b20=βy1(1+x1m)3,b11=β(1+x1m)2rk,b30=βy1(1+x1m)4,b21=β(1+x1m)3.
    $

    $ P_1(x_{m}, y_{m}) $ and $ P_2(x_{m}, y_{m}) $ are the power series in $ (x_{m}, y_{m}) $ with terms $ x_{m}^{i}y_{m}^{j} $ satisfying $ i+j\geq4 $.

    The expression of the first Lyapunov number can be expressed by the formula:

    $ l1=3π2a01(Det)3/2{[a10b10(a211+a11b02+a02b11)+a10a01(b211+a20b11+a11b02)2a10b10(b202a20a02)2a10a01(a220b20b02)a201(2a20b20+b11b20)+(a01b102a210)(b11b02a11a20)+b210(a11a02+2a02b02)](a210+a01b10)[3(b10b03a01a30)+2a10(a21+b12)+(b10a12a01b21)]},=3πβy1[6x1nk+2k+2n+2y1(1+x1m)3+βy1(1+x1m)2rk]2(x1m)1/2[βy1(1+x1m)2y1rk]3/2(1+x1m)5/23π[3nk3y1(1+x1m)4β(1+x1m)3]2(x1m)1/2[βy1(1+x1m)3y1rk(1+x1m)]1/23π[3x1nk+1k+1n+y1(1+x1m)3]2(x1m)3/2[βy1(1+x1m)y1rk(1+x1m)]1/2.
    $

    The limit cycle around equilibrium point $ E_1^* $ caused by Hopf bifurcation is unstable if $ l_1 > 0 $, otherwise the limit cycle is stable if $ l_1 < 0 $. For the expression of the first Lyapunov number $ l_1 $ is too cumbersome to tell the sign, the accuracy of this Theorem will be verified in the Section 5.

    Theorem 4.4. Hopf bifurcation takes place in model (2.3) around $ {E_1^*} $ at $ n = {n_{Hp}} $ when $ r > \frac{k}{{k + 1}}\left({d - \beta + \sqrt \Delta } \right) $.

    Proof. This Theorem is similar to the above one, we only need to prove the third condition

    (3') $ \frac{d}{{dn}}{\left. {Tr\left({{J_{E_1^*}}} \right)} \right|_{n = {n_{Hp}}}} \ne 0 $.

    Through calculation we can obtain

    $ \frac{d}{{dn}}{\left. {Tr\left( {{J_{E_1^*}}} \right)} \right|_{n = {n_{Hp}}}} = {\left. {\left( {\frac{3}{{{n^2}k}}x_1^2 - \frac{2}{{{n^2}}}{x_1} + \frac{{x_1^2\left( {1 - \frac{{{x_1}}}{k}} \right)}}{{\left( {1 + {x_1} - m} \right)\left( {{x_1} - m} \right){n^2}}}} \right)} \right|_{n = {n_{Hp}}}} \ne 0, $

    then we complete the proof.

    Usually several main parameters will affect model (2.3) collectively, therefore it is significant to focus on the dynamic behavior caused by the combined parameters. We study a bifurcation of codimension two with parameters $ m $ and $ \beta $ in this section, we pay main attention to the B-T bifurcation, which is caused by the simultaneous occurrence of saddle-node bifurcation and Hopf bifurcation.

    Theorem 4.5. When choose two bifurcation parameters $ m $ and $ \beta $, $ m_{BT} $ and $ \beta_{BT} $ are the bifurcation threshold valves, which satisfy

    $ {\left. {Tr\left( {{J_{E_1^*}}} \right)} \right|_{\left( {{m_{BT}}, {\beta _{BT}}} \right)}} = 0, \; \; {\left. {Det\left( {{J_{E_1^*}}} \right)} \right|_{\left( {{m_{BT}}, {\beta _{BT}}} \right)}} = 0, $

    then a B-T bifurcation takes place near the equilibrium point $ E_1^* $ with changing parameters $ \left(m, \beta\right) $ near $ \left(m_{BT}, \beta_{BT}\right) $.

    Here we derive model (2.3) into a normal form of the B-T bifurcation to obtain the specific expressions of the saddle-node, Hopf and homoclinic bifurcation curve in a small domain near the B-T point.

    Substituting $ m = {m_{BT}} + {\xi _1} $, $ \beta = {\beta _{BT}} + {\xi _2} $ into model (2.3), where $ \xi_1 $ and $ \xi_2 $ stand for two small perturbations. Then we can obtain the following model

    $ {dxdt=x(1xk)(xn1)xmBTξ11+xmBTξ1y,dydt=(βBT+ξ2)xmBTξ11+xmBTξ1y+ry(1xmBTξ1k)dy,
    $
    (4.1)

    Through introducing new variables $ u = x-x_1 $ and $ v = y-y_1 $, the equilibrium point $ E_1^* $ is translated to the origin, we have

    $ {dudt=a00(ξ)+a10(ξ)u+a01(ξ)v+a20(ξ)u2+a11(ξ)uv+a02(ξ)v2,dvdt=b00(ξ)+b10(ξ)u+b01(ξ)v+b20(ξ)u2+b11(ξ)uv+b02(ξ)v2+P3(u,v,ξ),
    $
    (4.2)

    where

    $ {a_{00}}\left( \xi \right) = {x_1}\left( {1 - \frac{{{x_1}}}{k}} \right)\left( {\frac{{{x_1}}}{n} - 1} \right) - \frac{{x_1 - {m_{BT}} - {\xi _1}}}{{1 + x_1 - {m_{BT}} - {\xi _1}}}{y_1}, \; \; {a_{01}}\left( \xi \right) = - \frac{{{x_1} - {m_{BT}} - {\xi _1}}}{{1 + {x_1} - {m_{BT}} - {\xi _1}}}, $
    $ {a_{10}}\left( \xi \right) = - \frac{3}{{nk}}x_1^2 + 2\left( {\frac{1}{n} + \frac{1}{k}} \right){x_1} - 1 - \frac{{{y_1}}}{{{{\left( {1 + {x_1} - {m_{BT}} - {\xi _1}} \right)}^2}}}, \; \; {a_{11}}\left( \xi \right) = - \frac{1}{{{{\left( {1 + {x_1} - {m_{BT}} - {\xi _1}} \right)}^2}}}, $
    $ {a_{20}}\left( \xi \right) = - \frac{3}{{nk}}{x_1} + \frac{1}{n} + \frac{1}{k} + \frac{{{y_1}}}{{{{\left( {1 + {x_1} - {m_{BT}} - {\xi _1}} \right)}^3}}}, \; \; {b_{10}}\left( \xi \right) = \frac{{\left( {{\beta _{BT}} + {\xi _2}} \right){y_1}}}{{{{\left( {1 + {x_1} - {m_{BT}} - {\xi _1}} \right)}^2}}} - \frac{r}{k}{y_1}, $
    $ {a_{02}}\left( \xi \right) = {b_{02}}\left( \xi \right) = 0, \; \; {b_{20}}\left( \xi \right) = - \frac{{\left( {{\beta _{BT}} + {\xi _2}} \right){y_1}}}{{{{\left( {1 + {x_1} - {m_{BT}} - {\xi _1}} \right)}^3}}}, \; \; {b_{11}}\left( \xi \right) = \frac{{{\beta _{BT}} + {\xi _2}}}{{{{\left( {1 + {x_1} - {m_{BT}} - {\xi _1}} \right)}^2}}} - \frac{r}{k}. $
    $ {b_{00}}\left( \xi \right) = \left( {{\beta _{BT}} + {\xi _2}} \right)\frac{{{x_1} - {m_{BT}} - {\xi _1}}}{{1 + {x_1} - {m_{BT}} - {\xi _1}}}{y_1} + r{y_1}\left( {1 - \frac{{{x_1} - {m_{BT}} - {\xi _1}}}{k}} \right) - d{y_1}, $
    $ {b_{01}}\left( \xi \right) = \left( {{\beta _{BT}} + {\xi _2}} \right)\frac{{{x_1} - {m_{BT}} - {\xi _1}}}{{1 + {x_1} - {m_{BT}} - {\xi _1}}} + r\left( {1 - \frac{{{x_1} - {m_{BT}} - {\xi _1}}}{k}} \right) - d, $

    and $ P_3\left(u, v, \xi \right) $ is power series in $ \left(u, v\right) $ with terms $ u^iv^j $ satisfying $ i+j \geqslant 4 $, whose coefficients are depend on $ \xi_1 $ and $ \xi_2 $ smoothly.

    Then, in a small domain of the origin $ \left(0, 0\right) $, we take the following $ {C^\infty } $ change of coordinates:

    $ {n_1} = u, \; \; \; {n_2} = {a_{00}}\left( \xi \right) + {a_{10}}\left( \xi \right)u + {a_{01}}\left( \xi \right)v + {a_{20}}\left( \xi \right){u^2} + {a_{11}}\left( \xi \right)uv, $

    model (4.2) can be written as

    $ {dn1dt=n2,dn2dt=c00(ξ)+c10(ξ)n1+c01(ξ)n2+c20(ξ)n21+c11(ξ)n1n2+c02(ξ)n22+P4(n1,n2,ξ),
    $
    (4.3)

    where

    $ c10(ξ)=a01(ξ)b10(ξ)a10(ξ)b01(ξ)a00(ξ)b11(ξ)+a11(ξ)b00(ξ)a00(ξ)a11(ξ)b01(ξ)a01(ξ),c00(ξ)=a01(ξ)b00(ξ)a00(ξ)b01(ξ),c01(ξ)=a10(ξ)+b01(ξ)a00(ξ)a11(ξ)a01(ξ),c20(ξ)=a01(ξ)b20(ξ)a20(ξ)b01(ξ)a10(ξ)b11(ξ)+a11(ξ)b10(ξ)+2a00(ξ)a112(ξ)b01(ξ)a012(ξ),c11(ξ)=2a20(ξ)+b11(ξ)a10(ξ)a11(ξ)a01(ξ)+a00(ξ)a112(ξ)a012(ξ),c02(ξ)=a11(ξ)a01(ξ),
    $

    and $ P_4\left(n_1, n_2, \xi \right) $ is power series in $ \left(n_1, n_2\right) $ with terms $ n_1^in_2^j $ satisfying $ i+j \geqslant 4 $, whose coefficients are depend on $ \xi_1 $ and $ \xi_2 $ smoothly.

    In order to remove the term $ {c_{02}}\left(\xi \right)n_2^2 $ from model (4.3), we take a new time variable $ \tau $, which satisfies $ \left({1 - {c_{02}}\left(\xi \right){n_1}} \right)d\tau = dt $. Then letting $ p_1 = n_1 $, $ p_2 = n_2\left(1-c_{02}n_1\right) $ and rewriting $ \tau $ as $ t $ for the seek of briefness, we obtain

    $ {dp1dt=p2,dp2dt=η00(ξ)+η10(ξ)p1+η01(ξ)p2+η20(ξ)p21+η11(ξ)p1p2+P5(p1,p2,ξ),
    $
    (4.4)

    where

    $ η00(ξ)=c00(ξ),η10(ξ)=c10(ξ)2c00(ξ)c02(ξ),η01(ξ)=c01(ξ),η20(ξ)=c20(ξ)+c00(ξ)c202(ξ)2c02(ξ)c10(ξ),η11(ξ)=c11(ξ)2c01(ξ)c02(ξ),
    $

    and $ P_5\left(p_1, p_2, \xi \right) $ is power series in $ \left(p_1, p_2\right) $ with terms $ p_1^ip_2^j $ satisfying $ i+j \geqslant 4 $, whose coefficients are depend on $ \xi_1 $ and $ \xi_2 $ smoothly.

    Let

    $ {q_1} = {p_1} + \frac{{{\eta _{10}}\left( \xi \right)}}{{2{\eta _{20}}\left( \xi \right)}}, \; \; \; {q_2} = {p_2}, $

    under which, model (4.4) can be translated into the following one

    $ {dq1dt=q2,dq2dt=ϖ00(ξ)+ϖ01(ξ)q2+ϖ20(ξ)q21+ϖ11(ξ)q1q2+P6(q1,q2,ξ),
    $
    (4.5)

    where

    $ ϖ00(ξ)=η00(ξ)η210(ξ)4η20(ξ),ϖ01(ξ)=η01(ξ)η10(ξ)η11(ξ)2η20(ξ),ϖ20(ξ)=η20(ξ),ϖ11(ξ)=η11(ξ),
    $

    and $ {P_6}\left({{q_1}, {q_2}, \xi } \right) $ is power series in $ \left(q_1, q_2\right) $ with terms $ q_1^iq_2^j $ satisfying $ i+j \geqslant 4 $, whose coefficients are depend on $ \xi_1 $ and $ \xi_2 $ smoothly.

    Notice that the expression of $ {\varpi _{20}}\left(\xi \right) $ is too complex to determine the sign, two cases are considered follows.

    case 1:$ {\varpi _{20}}\left(\xi \right) > 0 $, when $ \xi_i $ $ \left(i = 1, 2\right) $ is small enough, then take the following new variables

    $ {\nu _1} = {q_1}, \; \; {\nu _2} = \frac{{{q_2}}}{{\sqrt {{\varpi _{20}}\left( \xi \right)} }}, \; \; \tau = t\sqrt {{\varpi _{20}}\left( \xi \right)}, $

    retaining $ t $ to denote $ \tau $, then model (4.5) can be rewritten as

    $ {dν1dt=ν2,dν2dt=θ00(ξ)+θ01(ξ)ν2+ν21+θ11(ξ)ν1ν2+P7(ν1,ν2,ξ),
    $
    (4.6)

    where

    $ {\theta _{00}}\left( \xi \right) = \frac{{{\varpi _{00}}\left( \xi \right)}}{{{\varpi _{20}}\left( \xi \right)}}, \; \; {\theta _{01}}\left( \xi \right) = \frac{{{\varpi _{01}}\left( \xi \right)}}{{\sqrt {{\varpi _{20}}\left( \xi \right)} }}, \; \; {\theta _{11}}\left( \xi \right) = \frac{{{\varpi _{11}}\left( \xi \right)}}{{\sqrt {{\varpi _{20}}\left( \xi \right)} }}, $

    and $ {P_7}\left({{\nu _1}, {\nu _2}, \xi } \right) $ is power series in $ \left(\nu_1, \nu_2\right) $ with terms $ \nu_1^i\nu_2^j $ satisfying $ i+j \geqslant 4 $, whose coefficients are depend on $ \xi_1 $ and $ \xi_2 $ smoothly.

    Assume $ {\theta _{11}}\left(\xi \right) \neq 0 $, then through introducing new variables

    $ x = \theta _{11}^2\left( \xi \right){\nu _1}, \; \; y = \theta _{11}^3\left( \xi \right){\nu _2}, \; \; \tau = \frac{t}{{{\theta _{11}}\left( \xi \right)}}, $

    and retaining $ t $ to represent $ \tau $, then model (4.6) can be rewritten as

    $ {dxdt=y,dydt=σ00(ξ)+σ01(ξ)y+x2+xy+P8(x,y,ξ),
    $
    (4.7)

    where

    $ {\sigma _{00}}\left( \xi \right) = {\theta _{00}}\left( \xi \right)\theta _{11}^4\left( \xi \right), \; \; \; {\sigma _{01}}\left( \xi \right) = {\theta _{01}}\left( \xi \right){\theta _{11}}\left( \xi \right), $

    and $ {P_8}\left({x, y, \xi } \right) $ is power series in $ \left(x, y\right) $ with terms $ x^iy^j $ satisfying $ i+j \geqslant 4 $, whose coefficients depend on $ \xi_1 $ and $ \xi_2 $ smoothly.

    case 2:$ {\varpi _{20}}\left(\xi \right) < 0 $, when $ \xi_i $ $ \left(i = 1, 2\right) $ is small enough, then take the following new variables

    $ {\nu _1^\prime} = {q_1}, \; \; {\nu _2^\prime} = \frac{{{q_2}}}{{\sqrt {-{\varpi _{20}}\left( \xi \right)} }}, \; \; \tau^\prime = t\sqrt {-{\varpi _{20}}\left( \xi \right)}, $

    retaining $ t $ to denote $ \tau^\prime $, then model (4.5) can be rewritten as

    $ {dν1dt=ν2,dν2dt=θ00(ξ)+θ01(ξ)ν2ν21+θ11(ξ)ν1ν2+P7(ν1,ν2,ξ),
    $
    (4.6')

    where

    $ {\theta _{00}^\prime}\left( \xi \right) = - \frac{{{\varpi _{00}}\left( \xi \right)}}{{{\varpi _{20}}\left( \xi \right)}}, \; \; {\theta _{01}^\prime}\left( \xi \right) = \frac{{{\varpi _{02}}\left( \xi \right)}}{{\sqrt {-{\varpi _{20}}\left( \xi \right)} }}, \; \; {\theta _{11}^\prime}\left( \xi \right) = \frac{{{\varpi _{11}}\left( \xi \right)}}{{\sqrt {-{\varpi _{20}}\left( \xi \right)} }}, $

    and $ {P_7^\prime}\left({{\nu _1^\prime}, {\nu _2^\prime}, \xi } \right) $ is power series in $ \left(\nu_1^\prime, \nu_2^\prime\right) $ with terms $ \nu_1^{\prime i}\nu_2^{\prime j} $ satisfying $ i+j \geqslant 4 $, whose coefficients depend on $ \xi_1 $ and $ \xi_2 $ smoothly.

    Supposing $ {\theta _{11}^\prime}\left(\xi \right) \neq 0 $, then we make the following transformation:

    $ x = -\theta _{11}^{\prime 2}\left( \xi \right){\nu _1^\prime}, \; \; y = \theta _{11}^{\prime 3}\left( \xi \right){\nu _2^\prime}, \; \; \tau^\prime = -\frac{t}{{{\theta _{11}^\prime}\left( \xi \right)}}, $

    retaining $ t $ to denote $ \tau^\prime $, then model $ (4.6^\prime) $ can be rewritten as

    $ {dxdt=y,dydt=σ00(ξ)+σ01(ξ)y+x2+xy+P8(x,y,ξ),
    $
    (4.7')

    where

    $ {\sigma _{00}^\prime}\left( \xi \right) = -{\theta _{00}^\prime}\left( \xi \right)\theta _{11}^{\prime 4}\left( \xi \right), \; \; \; {\sigma _{01}^\prime}\left( \xi \right) = -{\theta _{01}^\prime}\left( \xi \right){\theta _{11}^\prime}\left( \xi \right), $

    and $ {P_8^\prime}\left({x, y, \xi } \right) $ is power series in $ \left(x, y\right) $ with terms $ x^iy^j $ satisfying $ i+j \geqslant 4 $, whose coefficients depend on $ \xi_1 $ and $ \xi_2 $ smoothly.

    We retain $ \sigma_{00}\left(\xi \right) $ and $ \sigma_{01}\left(\xi \right) $ to denote $ \sigma_{00}^\prime\left(\xi \right) $ and $ \sigma_{01}^\prime\left(\xi \right) $ in $ (4.7^\prime) $ to reduce the number of situations to be discussed. In a small domain of $ \left(0, 0\right) $, the two transformations of models (4.7) and $ (4.7^\prime) $ are homeomorphisms, and $ \sigma_{00}\left(\xi \right) $, $ \sigma_{01}\left(\xi \right) $ are independent parameters when $ {\left| {\frac{{\partial \left({{\sigma _{00}}, {\sigma _{01}}} \right)}}{{\partial \left({{\xi _1}, {\xi _2}} \right)}}} \right|_{{\xi _1} = {\xi _2} = 0}} \ne 0 $. According to the results of [33,34], model (2.3) undergoes a B-T bifurcation when $ \xi = \left(\xi_1, \xi_2\right) $ is in a small doamin of the origin. The local expressions around the origin of the bifurcation curves can be expressed as follows ("+" denotes $ {\varpi _{20}}\left(\xi \right) > 0 $ and "$ - $" denotes $ {\varpi _{20}}\left(\xi \right) < 0 $) :

    (1) The expression of saddle-node bifurcation curve:

    $ SN = \left\{ {\left({{\xi _1}, {\xi _2}} \right):{\sigma _{00}}\left({{\xi _1}, {\xi _2}} \right) = 0, {\sigma _{01}}\left({{\xi _1}, {\xi _2}} \right) \ne 0} \right\}; $

    (2) The expression of Hopf bifurcation curve:

    $ Hp = \left\{ {\left({{\xi _1}, {\xi _2}} \right):{\sigma _{01}}\left({{\xi _1}, {\xi _2}} \right) = \pm \sqrt { - {\sigma _{00}}\left({{\xi _1}, {\xi _2}} \right)}, {\sigma _{00}}\left({{\xi _1}, {\xi _2}} \right) < 0} \right\}; $

    (3) The expression of homoclinic bifurcation curve:

    $ HL = \left\{ {\left({{\xi _1}, {\xi _2}} \right):{\sigma _{01}}\left({{\xi _1}, {\xi _2}} \right) = \pm \frac{5}{7}\sqrt { - {\sigma _{00}}\left({{\xi _1}, {\xi _2}} \right)}, {\sigma _{00}}\left({{\xi _1}, {\xi _2}} \right) < 0} \right\}. $

    In real life, ecological control of cyanobacteria bloom is one of the most effective methods, especially, controlling algae bloom by fish. Although fish can effectively graze algae, algae can prevent grazing through aggregation. Thus, to better understand the dynamic variation of fish and algae, some numerical simulation work needs to be implemented. Since model (2.3) contains too many parameters, we fix some of the parameters as follows for the seek of convenient: $ k = 5, \; n = 0.2, \; r = 0.15, \; d = 0.21. $

    On the premise of the fixed parameters previously, we fix the parameter $ m = 1.5 $ and let parameter $ \beta $ vary within a certain range as Figure 3(a). It is can be find from Figure 3(a) obviously that model (2.3) has abundant dynamic properties. When $ \beta = \beta_{SN} = 0.1748528137 $ and $ \beta = \beta_{TC} = 0.2121428571 $, there takes place a saddle-node bifurcation and transcritical bifurcation in model (2.3) respectively. Firstly the internal equilibrium points $ E_1^*\left(x_1, y_1\right) $ and $ E_2^*\left(x_2, y_2\right) $ are not exist when $ \beta < \beta_{SN} $. In the meantime, $ E_0 $ and $ E_2 $ are unstable, but $ E_1 $ is a stable node. With the value of $ \beta $ increasing and passing through $ \beta_{SN} $, there appear two additional internal equilibrium points $ E_1^*\left(x_1, y_1\right) $ and $ E_2^*\left(x_2, y_2\right) $, the former is a stable node while the latter is an unstable saddle, this process can be seen more clearly from Figure 4. Secondly, the internal equilibrium point $ E_2^*\left(x_2, y_2\right) $ will coincide with the boundary equilibrium point $ E_1\left(5, 0\right) $ with the value of $ \beta $ increasing and passing through $ \beta_{TC} $. And this collision of the two equilibrium points changes $ E_1\left(5, 0\right) $ from a stable node to an unstable saddle.

    Figure 3.  Bifurcation diagrams of model (2.3) with the previously fixed parameters. The red and bule curves represent the internal equilibrium points $ E_1^* $ and $ E_2^* $, respectively. And three horizontal lines stand for three other equilibrium points: $ E_0 $(green), $ E_1 $(cyan) and $ E_2 $(pink). Equilibrium points presented as solid curves are stable, dotted curves are unstable. Moreover, the vertical lines with labels 'SN', 'TC' and 'Hp' indicate that model (2.3) undergoes saddle-node, Transcritical and Hopf bifurcation here respectively. (a) fix m = 1.5 then vary $ \beta $. (b) fix $ \beta = 0.17488281 $ then vary m.
    Figure 4.  The process of saddle-node bifurcation of model (2.3) with parameters fixed as mentioned in the main text. (a) Internal equilibrium points $ E_1^* $ and $ E_2^* $ are not exist when $ 0.1738 = \beta < \beta_{SN} $. (b) When $ \beta = \beta_{SN} $, there exists a special internal equilibrium point, which is a saddle-node. (c) Exist two internal equilibrium points $ E_1^* $ and $ E_2^* $, the former is a stable node and the latter is an unstable saddle, when $ 0.1758 = \beta > \beta_{SN} $. (d) Partical enlarged view of the whole process of saddle-node bifurcation. The black point represents a saddle-node, red for $ E_1^* $ and blue for $ E_2^* $.

    Similar to the above approach, now we fix an additional parameter $ \beta = 0.17488281 $, then let parameter $ m $ vary within a certain range as Figure 3(b). At the beginning, that is to say when $ m < m_{Hp} = 1.385035335 $, both two internal equilibrium points $ E_1^* $ and $ E_2^* $ are unstable, the former is an unstable focus and the latter is an unstable saddle. Furthermore, it is worth mentioning that there is a limit cycle in a small neighborhood containing $ E_1^* $, which is surrendered by one of the unstable trajectory of saddle $ E_2^* $. At the same time, we can calculate that the first Lyapunov number is $ l_1 = -770.954244\pi < 0 $, which means that the limit cycle is stable. Thus, the unstable focus will become a center when the value of $ m $ gets greater and reaches $ m_{Hp} $. that is to say, $ E_1^* $ will turn into a stable focus, as the value of $ m $ continues to increase and exceeds $ m_{Hp} $. In a word, the detailed dynamic evolution process of Hopf bifurcation can be seen in Figure 5. Furthermore, it is easy to see from Figure 5(d) that if the value of $ m $ is less than $ m_{Hp} $, model (2.3) has limit cycles around $ E_1^* $, which are represented by the colorful circles in the diagram, if the value of $ m $ is larger than $ m_{Hp} $, model (2.3) has a stable equilibrium point. Moreover, it is also worth emphasizing that the amplitude of limit cycle is increasing as the value of $ m $ is decreasing, which means that if the aggregation area of algae is smaller, the area where algae and fish periodically oscillate and coexist is more wider. And with the passage of time, fish and algae will eventually reach a state of coexistence within a small range of $ m > m_{Hp} $.

    Figure 5.  The process of Hopf bifurcation according to the bifurcation parameter $ n $, with $ k = 5, \; r = 0.15, \; d = 0.21, \; m = 1.385035335, \; \beta = 0.17488281 $. (a) $ E_1^* $ is locally asymptotically stable when $ 0.1998 = n < n_{Hp} $. (b) $ E_1^* $ is a center when $ n = n_{Hp} = 0.2 $. (c) Through Hopf bifurcation with $ 0.2005 = n = n_{Hp} $, there exists a stable periodic orbits around the unstable focus $ E_1^* $. (d) Hopf bifurcation diagram repressenting stable $ E_1^* $ and stable limit cycles with various values of $ n $.

    In order to explore how Allee effect affect the dynamic behavior of model (2.3), we choose Allee threshold $ n $ as the control parameter to simulate the Hopf bifurcation with $ k = 5, \; r = 0.15, \; d = 0.21, \; m = 1.385035335, \; \beta = 0.17488281 $, the detailed results are shown in Figure 6. It is can see clearly that the equilibrium point $ E_1^* $ is a stable focus and $ E_2^* $ is an unstable saddle when $ 0.1998 = n < n_{Hp} = 0.2 $. That is to say, when the population densities are within a certain range, the algae and fish will coexist at equilibrium point $ E_1^* $. When the value of $ n $ is larger than $ n_{Hp} = 0.2 $, the equilibrium point $ E_1^* $ will lose stability, and a stable limit cycle will appear, which implies that a supercritical Hopf bifurcation takes place. Moreover, when the value of $ n $ increases, the amplitude of limit cycle will be larger and larger. Obviously, this situation represents that the periodic oscillation coexistence mode between algae and fish will gradually take shape when $ n $ exceeds the Hopf bifurcation threshold value. And in a certain small range, the larger the value of $ n $, the more conducive it is for the coexistence of periodic oscillation of algae and fish. At the same time, it is also worth mentioning that the value of Allee threshold $ n $ seriously affects the dynamic behavior of model (2.3).

    Figure 6.  The process of Hopf bifurcation according to the bifurcation parameter $ m $. (a) Through Hopf bifurcation with $ 1.3840 = m < m_{Hp} $, there exists a stable periodic orbits around the unstable focus $ E_1^* $. (b) $ E_1^* $ is a center when $ m = m_{Hp} = 1.385035335 $. (c) $ E_1^* $ is locally asymptotically stable when $ 1.3851 = m > m_{Hp} $. (d) Hopf bifurcation diagram repressenting stable $ E_1^* $ and stable limit cycles with various values of $ m $.

    Now, it can be found from Figure 3(b) that model (2.3) can go through two transcritical bifurcations as the value of $ m $ increase. When $ m_{Hp} < m < m_{TC1} = 3.547677046 $, the equilibrium points $ E_1 $ and $ E_1^* $ are stable node, the equilibrium point $ E_2^* $ is an unstable saddle. As the value of $ m $ passes through $ m_{TC1} $, a transcritical bifurcation takes place, which can cause the collision of saddle $ E_2^* $ with stable node $ E_1 $. This collision transforms the stability of boundary equilibrium point $ E_1 $, and makes it an unstable saddle. Then if the value of $ m $ increases greater than $ m_{TC2} = 3.622895828 $, model (2.3) occurs a transcritical bifurcation at $ E_1 $ again, but this time is the consequence of the collision of $ E_1 $ with $ E_1^* $, which can prompt saddle $ E_1 $ regain its stability and back to a node. These switching of stability with respect to two Transcritical bifurcations can be better explained through the phase portraits in Figure 7.

    Figure 7.  The bifurcation curves diagram of model (2.3). The small diagram inside is a partial enlargement of the big diagram.

    In order to study how the parameters $ m $ and $ \beta $ synergistically affect the dynamic behavior of model (2.3), the numerically simulation of B-T bifurcation with $ k = 5, \; n = 0.2, \; r = 0.15, \; d = 0.21 $ will be carried out. By calculation, we obtain $ m_{BT} = 1.374006280 $ and $ \beta_{BT} = 0.1748528137 $, and we have

    $ {\left| {\frac{{\partial \left( {{\sigma _{00}}, {\sigma _{01}}} \right)}}{{\partial \left( {{\xi _1}, {\xi _2}} \right)}}} \right|_{{\xi _1} = {\xi _2} = 0}}{\text{ = }}\left| {0.00009876 - 43590.5557695.31449003 - 7.25883992
    } \right|{\text{ = }}4.154811591\cdot 10^6 \ne 0, $
    $ \varpi_{20} = 0.19835682+0.30459173\xi_1+1.52158465\xi_2+0.41550653\xi_1^2+1.88040682\xi_1\xi_2, $
    $ \theta_{11} = -5.51455564+24.53907955\xi_1+21.15093188\xi_2+12.40764023\xi_1^2-100.4590524\xi_1\xi_2-121.6857571\xi_2^2. $

    Therefore, we know that the transformation of parameters is nonsingular, and $ \varpi_{20} > 0 $ and $ \theta_{00} \ne 0 $ for small $ \xi_1 $ and $ \xi_2 $. Moreover, the local expressions of bifurcation curves $ SN $, $ Hp $ and $ HL $ around the origin are revealed up to second-order approximately as:

    (1) The saddle-node bifuecation curve satisfies $ {\sigma _{01}}\left({{\xi _1}, {\xi _2}} \right) \ne 0 $, and

    $ SN = \left\{\left( {{\xi _1}, {\xi _2}} \right)|\; 0.00009876\xi_1-43590.55576\xi_2+0.06615614\xi_1^2+855595.4398\xi_1\xi_2+1003142.522\xi_2^2 \right\}, $

    (2) The Hopf bifurcation curve satisfies $ {\sigma _{00}}\left({{\xi _1}, {\xi _2}} \right) < 0 $ and

    $ Hp = \left\{ \left( {{\xi _1}, {\xi _2}} \right)|\; 0.00009876\xi_1-43590.55576\xi_2+9084.918166\xi_1^2+854211.6945\xi_1\xi_2+1003195.213\xi_2^2 \right\}, $

    (3) The homoclinic bifurcation curve satisfies $ {\sigma _{00}}\left({{\xi _1}, {\xi _2}} \right) < 0 $ and

    $ HL = \left\{ \left( {{\xi _1}, {\xi _2}} \right)|0.00009876\xi_1-43590.55576\xi_2+17806.37610\xi_1^2+852883.2991\xi_1\xi_2+1003245.796\xi_2^2 \right\}. $

    Meanwhile, the corresponding bifurcation curves is depicted in Figure 8, and the small image in Figure 8 is a partial enlargement of the saddle-node bifurcation curve, which should be exactly coincided with horizontal line $ \xi_2 = 0 $. Here we need to explain that the error in the image is inevitable through simulation, but this error does not affect the readability of the numerical simulation results. Obviously, these three curves will divide the visible area into four regions as I, II, III and IV. It is easy to see form Figure 9(a) that when $ \left(\xi_1, \xi_2\right) = \left(0, 0\right) $, model (2.3) has an unique internal equilibrium point, which is a cusp of codimension 2. If we fix $ \xi_1 = 0.01 $ for the seek of convenience, then we have the following results.

    Figure 8.  The process of Transcritical bifurcation of model (2.3) with parameters fixed as mentioned in the main text. (a) Internal equilibrium point $ E_1^* $ and boundary equilibrium point $ E_1 $ are stable node but $ E_2^* $ is an unstable saddle when $ 3.2 = m < m_{TC1} = 3.547677046 $. (b) When $ m_{TC1} < m = 3.6 < m_{TC2} = 3.622895828 $, $ E_1 $ lose its stability and becomes a saddle through colliding with $ E_2^* $. (c) $ E_1 $ regains its stability and becomes a node through colliding with $ E_1^* $ when $ 4 = m > m_{TC2} $. (d) Partical enlarged view of the whole process of Transcritical bifurcation. The black point represents equilibrium point $ E_1 $, red for $ E_1^* $ and blue for $ E_2^* $.
    Figure 9.  Phase portraits of model (4.1) as varying the values of $ \xi_1 $ and $ \xi_2 $ around $ \left(0, 0\right) $. (a) When $ \left(\xi_1, \xi_2\right) = \left(0, 0\right) $, there exists a cusp with codimension two. (b) Model (4.1) has no internal equilibrium point as $ \left(\xi_1, \xi_2\right) = \left(0.01, -0.0001\right) $ locating in region I. (c) Two internal equilibrium points appear through a saddle-node bifurcation, $ E_1^* $ is a stable focus and $ E_2^* $ is an unstable saddle, when $ \left(\xi_1, \xi_2\right) = \left(0.01, 0.00001\right) $ lies in region II. (d) A stable limit cycle arises surrounding $ E_1^* $ with the happening of Hopf bifurcation as $ \left(\xi_1\, xi_2\right) $ into region III. (e) When $ \left(\xi_1, \xi_2\right) = \left(0.01, 0.00005307\right) $ lies on the curve $ HL $, the limit cycle will get larger and approaches to the saddle $ E_2^* $ then becomes a homoclinic orbit. (f) The homoclinic break as $ \left(\xi_1, \xi_2\right) $ into region IV, $ E_1^* $ and $ E_2^* $ are unstable, the former is a focus while the latter is a saddle.

    (1) There is no internal equilibrium point when $ \xi_2 = -0.0001 $ in region I, it can see from Figure 9(b) that the fish will finally extinct as long as the algae population exceeds a certain threshold.

    (2) Model (4.1) has gone through a saddle-node bifurcation when $ \xi_2 = 0.00001 $, $ \left(\xi_1, \xi_2\right) $ can go through the curve $ SN $ into region II, there exists a stable focus $ E_1^* $ and an unstable saddle $ E_2^* $, which is shown in Figure 9(c).

    (3) With the value of $ \xi_2 $ varying from 0.00001 to 0.00003, $ \left(\xi_1, \xi_2\right) $ moves from region I to II and model (4.1) occurs a supercritical Hopf bifurcation (see Figure 9(d)), which deprives the stability of $ E_1^* $ and produces a stable limit cycle surrounding $ E_1^* $.

    (4) As the value of $ \xi_2 $ increases to 0.00005307, $ \left(\xi_1, \xi_2\right) $ is exactly locating on the curve $ HL $, the stable limit cycle grows and goes through the saddle $ E_2^* $ and converts to an unstable homoclinic orbit, which can be seen in Figure 9(e).

    (5) The homoclinic orbit disappears when $ \xi_2 = 0.00007 $ and $ \left(\xi_1, \xi_2\right) $ is in region IV, and there exists an unstable focus $ E_1^* $ and a saddle $ E_2^* $ (see Figure 9(f)).

    From the above numerical example, it can be seen that the aggregation effect and conversion rate have a great influence on model (2.3), and the two populations show a series of rich properties such as extinction, coexistence and periodic oscillation near the B-T bifurcation point.

    Based on the above numerical simulation analysis, we can know that the values of three key parameters $ m $, $ n $ and $ \beta $ have an important influence on the dynamic behaviors of model (2.3), which can not only essentially change the dynamic characteristics, but also affect the survival and extinction of algae and fish. It is clearly visible from Figures 3(a) and 4 that the value of key parameter $ \beta $ can lead to the taking place of a saddle-node bifurcationthe in model (2.3), which can promote algae and fish to form a new stable coexistence model. It must be stated from Figures 3(b), 5 and 7 that the value of the key parameter $ m $ is very important for the occurrence of subcritical Hopf bifurcation and transcritical bifurcation of model (2.3), which is not only related to whether algae and fish can survive for a long time, but also can urge algae and fish to form a stable periodic oscillation coexistence mode. It is more worthy of our clear understanding from Figures 8 and 9 that the synergistic action mechanism of key parameters $ m $ and $ \beta $ plays an important role in changing the dynamic characteristics and internal essential laws of model (2.3), which can not only impel algae and fish to form two new types of periodic oscillation coexistence modes (Limit cycle coexistence mode and homoclinic orbit coexistence mode), but also lead to the extinction of algae and fish, or the extinction of fish and the outbreak of algae bloom. All in all, it is necessary to consider Allee effect and aggregation effect in building aquatic ecological model to further ascertain the dynamic relationship between algae and fish.

    Although the dynamic behaviors of the predator-prey model with Allee effect in the single or both species have been studied by many mathematicians and biologists, there was relatively few literature on the combination of it with algal aggregation. Within the research framework of mutual restriction between algae and fish, we proposed a modified algae and fish model to probe bifurcation dynamic behaviors. Firstly, we studied all possible equilibrium points of model (2.3) and their stability, which could not only directly give the ideal coexistence mode of algae and fish, but also provide a theoretical basis for the subsequent discussion of bifurcation dynamics. Secondly, when there is only one internal equilibrium point in model (2.3), we gave a threshold condition about aggregation effect parameter $ m $ and Allee effect parameter $ n $ to ensure that model (2.3) has a limit cycle, this result implied that the relationship between aggregation effect parameter $ m $ and Allee effect parameter $ n $ was very important for whether algae and fish could exist periodic oscillation survival mode. Finally, by selecting the absorption-conversion rate $ \beta $, aggregation effect parameter $ m $ and Allee effect threshold $ n $ as control parameters respectively, we theoretically derived some key threshold conditions to compel that model (2.3) could undergo transcritical, saddle-node, Hopf and B-T bifurcation, these bifurcation dynamic behaviors could force the internal essential changes in the dynamic relationship and coexistence mode between algae and fish. Moreover, it had to be said that these theoretical results not only could summarize and develop the previous theoretical research results, but also further promote the rapid development of bifurcation dynamics in aquatic ecosystem.

    In order to verify the feasibility of the theoretical analysis results and visually explore dynamic relationship between algae and fish, a large number of bifurcation numerical simulation results were implemented. Through analysis and comparison, it could be seen that in model (2.3) aggregation effect parameter $ m $ and Allee effect parameter $ n $ played an important role in the occurrence and evolution of bifurcation dynamics, which also indirectly showed that the coexistence mode of algae and fish depended heavily on aggregation effect and Allee effect. Furthermore, through numerical analysis of B-T bifurcation behavior, it was worth us to clarify that value relationship of the absorption-conversion rate $ \beta $ and aggregation effect parameter $ m $ could promote that model (2.3) experienced saddle-node, Hopf and Homoclinic bifurcation, these bifurcation behaviors could not only force algae and fish to form three new coexistence modes, but also was the internal driving force to dynamically adjust their coexistence mode. In a word, it should be emphasized that when we used the mathematical ecological model to study the dynamic relationship between algae and fish, aggregation effect and Allee effect were one of the ecological and environmental factors, which needed to be considered urgently.

    The theoretical and numerical results of this paper can get the following four results: (1) If the dynamic behavior of transcritical bifurcation can occur in model (2.3), which means that algae and fish can change from a single population survival mode to a dual population sustainable survival mode. (2) If the dynamic behavior of saddle-node bifurcation can occur in model (2.3), which hints that algae and fish can change from a dual population unsteady coexistence mode to a dual population steady coexistence mode. (3) If the dynamic behavior of Hopf bifurcation can occur in model (2.3), which shows that algae and fish can change from a dual population steady-state coexistence mode to a dual population periodic oscillation mode. (4) If the dynamic behavior of B-T bifurcation can occur in model (2.3), which suggests that the coexistence mode of algae and fish is easily affected by eco-environmental factors, and the coexistence mode can change back and forth between an equilibrium point steady-state mode, a stable periodic oscillation mode and an unstable survival mode. Therefore, the four kinds of bifurcation dynamic behaviors occurred in model (2.3) can represent the dynamic change characteristics of the coexistence mode of algae and fish under the dynamic change of eco-environmental factors. In other words, the interaction mechanism between algae and fish in nature can be described and explained by the bifurcation dynamic behavior, dynamic evolution process and internal essential characteristics of model (2.3).

    In the follow-up research works, we will firstly deepen theoretical research of bifurcation dynamics in some modified ecological models, especially B-T bifurcation. Since neither algae nor fish will remain in a fixed space for many factors (mate choice, food supplies, population density, etc.), it is meaningful to consider their changes of spatial diffusibility[35,36]. Finally explore how to control algae bloom by fish under human external control. In a word, all the results are expected to be helpful in the process of study bifurcation dynamic behavior in some aquatic ecosystems.

    This work was supported by the National Natural Science Foundation of China (Grants No. 61871293 and No. 61901303), the National Key Research and Development Program of China (Grant No. 2018YFE0103700).

    All authors declare that there is no conflict of interest.

    Table .  Existence of internal equilibrium points.
    Equilibrium points Conditions
    $E_{1}^*$ ${\bf case \ 1}$: $\left({m - n} \right)\left({1 - k - m + n} \right) - k > 0$
    $\max \left\{ {0, \frac{{\left({\beta - d} \right){k^2} + \left({dm - \beta m - d} \right)k}}{{m\left({m - k - 1} \right)}}} \right\} < r < \min \left\{ {d, \frac{{k\left[{d- \left({d - \beta } \right)\left({m - n} \right)} \right]}}{{\left({m - n} \right)\left({1 - k - m + n} \right) - k}}} \right\}$
    ${\bf case \ 2}$: $\left({m - n} \right)\left({1 - k - m + n} \right) - k < 0$
    $\max \left\{ {0, \frac{{\left({\beta - d} \right){k^2} + \left({dm - \beta m - d} \right)k}}{{m\left({m - k - 1} \right)}}, \frac{{k\left[{d- \left({d - \beta } \right)\left({m - n} \right)} \right]}}{{\left({m - n} \right)\left({1 - k - m + n} \right) - k}}} \right\} < r < d$
    $E_{2}^*$ ${\bf case \ 1}$: $m>n$
    $d < r < \frac{{\left({\beta - d} \right){k^2} + \left({dm - \beta m - d} \right)k}}{{m\left({m - k - 1} \right)}}$
    ${\bf case \ 2}$: $m < n$
    (1).$\left({m - n} \right)\left({1 - k - m + n} \right) - k > 0$
    $\max \left\{ {0, \frac{{k\left[{d- \left({d - \beta } \right)\left({m - n} \right)} \right]}}{{\left({m - n} \right)\left({1 - k - m + n} \right) - k}}} \right\} < r < \frac{{\left({\beta - d} \right){k^2} + \left({dm - \beta m - d} \right)k}}{{m\left({m - k - 1} \right)}}$
    (2).$\left({m - n} \right)\left({1 - k - m + n} \right) - k < 0$
    $0 < r < \min \left\{ {\frac{{\left({\beta - d} \right){k^2} + \left({dm - \beta m - d} \right)k}}{{m\left({m - k - 1} \right)}}, \frac{{k\left[{d- \left({d - \beta } \right)\left({m - n} \right)} \right]}}{{\left({m - n} \right)\left({1 - k - m + n} \right) - k}}} \right\}$
    $E_{1}^*$, $E_{2}^*$ $\left({2n - 2m + 1} \right)r < k\left({\beta - d} \right)$
    ${\left({k + 1} \right)^2}{r^2}+2\left[{\left({\beta - d} \right){k^2} - \left({\beta + d} \right)k} \right]r + {\left({\beta - d} \right)^2}{k^2} > 0$
    ${\bf case \ 1}$: $\left({m - n} \right)\left({1 - k - m + n} \right) - k > 0$
    $0 < r < \min \left\{ {d, \frac{{k\left[{d- \left({d - \beta } \right)\left({m - n} \right)} \right]}}{{\left({m - n} \right)\left({1 - k - m + n} \right) - k}}, \frac{{\left({\beta - d} \right){k^2} + \left({dm - \beta m - d} \right)k}}{{m\left({m - k - 1} \right)}}, \frac{k}{{1 + k}}\left({\beta - d} \right)} \right\}$
    ${\bf case \ 2}$: $\left({m - n} \right)\left({1 - k - m + n} \right) - k < 0$
    $\max \left\{ {0, \frac{{k\left[{d- \left({d - \beta } \right)\left({m - n} \right)} \right]}}{{\left({m - n} \right)\left({1 - k - m + n} \right) - k}}} \right\} < r < \min \left\{ {d, \frac{{\left({\beta - d} \right){k^2} + \left({dm - \beta m - d} \right)k}}{{m\left({m - k - 1} \right)}}, \frac{k}{{1 + k}}\left({\beta - d} \right)} \right\}$
    $E_{1}^* = E_{2}^*$ $r = \frac{{k\left({dk - \beta k + d + \beta \pm \sqrt { - {\beta ^2}k + \beta dk + \beta d} } \right)}}{{{{\left({k + 1} \right)}^2}}}$
    $0 < r < \frac{k}{{1 + k}}\left({\beta - d} \right)$
    $\max \left\{ {0, \frac{{2n - 2m + 1}}{{\beta - d}}r} \right\} < k \leqslant \frac{d}{{\beta - d}}$

     | Show Table
    DownLoad: CSV
    [1] Numer. Math., 98 (2004), 195-221.
    [2] SIAM J. Math. Anal., 36 (2004), 301-322.
    [3] Commun. Part. Diff. Eqs., 32 (2007), 127-148.
    [4] Math. Z., 194 (1987), 375-396.
    [5] Math. Nachr., 195 (1998), 77-114.
    [6] Appl. Math. Comput., 218 (2011), 4587-4594.
    [7] Comput. Math. Appl., 64 (2012), 1927-1936.
    [8] RACSAM Rev. R. Acad. Cienc. Exactas Fs. Nat. Ser. A. Mat., 95 (2001), 281-295.
    [9] Numer. Math., 93 (2003), 655-673.
    [10] Banach Center Publ., 63 (2004), 209-216.
    [11] SIAM J. Math. Anal., 35 (2003), 561-578.
    [12] Nonlinear Anal., 12 (2011), 2826-2838.
    [13] Appl. Numer. Math., 59 (2009), 1059-1074.
    [14] Nonlinear Analysis TMA, 8 (1984), 1121-1144.
    [15] J. Diff. Eqs., 131 (1996), 79-131.
    [16] Adv. Math., Beijing, 25 (1996), 283-284.
    [17] J. Math. Biol., 9 (1980), 49-64.
    [18] J. Theor. Biol., 79 (1979), 83-99.
    [19] Nonlinear Analysis TMA, 21 (1993), 603-630.
  • This article has been cited by:

    1. Li Yin, Ying Xu, Desheng Kong, Juan Wang, Kaipian Shi, Yong Zhang, Huan He, Shaogui Yang, Lixiao Ni, Shiyin Li, Role of extracellular polymeric substances in resistance to allelochemical stress on Microcystis aeruginsosa and its mechanism, 2023, 41, 2096-5508, 2219, 10.1007/s00343-023-2318-z
    2. Shengyu Huang, Hengguo Yu, Chuanjun Dai, Zengling Ma, Qi Wang, Min Zhao, Dynamics of a harvested cyanobacteria-fish model with modified Holling type Ⅳ functional response, 2023, 20, 1551-0018, 12599, 10.3934/mbe.2023561
    3. I. R. Belshiasheeela, Mini Ghosh, The role of zooplankton in the growth of algal bloom: a mathematical study, 2024, 42, 0736-2994, 591, 10.1080/07362994.2024.2323535
    4. Chenyu Liang, Hangjun Zhang, Yancong Xu, Libin Rong, Bifurcation Analysis of a New Aquatic Ecological Model with Aggregation Effect and Harvesting, 2023, 33, 0218-1274, 10.1142/S0218127423501808
  • Reader Comments
  • © 2013 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(2909) PDF downloads(500) Cited by(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog