Dynamics Analysis of an Aquatic Ecological Model with Allee Effect

Lingzhen Zheng, Jueyuan Yan, Hengguo Yu^*

School of Mathematics and Physics, Wenzhou University, Wenzhou, China.

DOI: 10.4236/am.2022.1310052                

Abstract

In this paper, based on the dynamic relationship between algae and protozoa, an aquatic ecological model with Allee effect was established to investigate how some ecological environment factors affect coexistence mode of algae and protozoa. Mathematical derivation works mainly gave some key conditions to ensure the existence and stability of all possible equilibrium points, and to induce the occurrence of transcritical bifurcation and Hopf bifurcation. The numerical simulation works mainly revealed ecological relationship change characteristics of algae and protozoa with the help of bifurcation dynamics evolution process. Furthermore, it was also worth emphasizing that Allee effect had a strong influence on the dynamic relationship between algae and protozoa. In a word, it was hoped that the research results could provide some theoretical support for algal bloom control, and also be conducive to the rapid development of aquatic ecological models.

Keywords

Algae, Protozoa, Allee Effect, Bifurcation, Relationship

1. Introduction

As everyone knows, lake eutrophication is a natural process and a stage of lake evolution. However, a common phenomenon associated with lake eutrophication is that many phytoplankton, especially those with buoyancy or mobility, usually multiply in large numbers to form algal blooms, which can cause a series of serious water environment problems [1] [2]. Therefore, the cause, harm and control measures of algal blooms have become one of the important environmental issues concerned by the academic community.

There are many organisms in nature that can inhibit algal growth, mainly including: cyanobacteria virus (algaphage), alginolytic bacteria, protozoa, fungi and actinomycetes [3] [4] [5] [6]. At the same time, protozoa are an important link in aquatic food chain, many protozoa can eat algae, and some even take algae as their only food [7]. A large number of studies have shown that the decline of algal biomass is often accompanied by a sharp increase in the number of protozoa, considering the environmental adaptability, reproductive ability, algal control efficacy, host range, adaptability to host changes of various algal control biological factors, scholars believe that protozoa are a control factor with great application prospects [7] [8] [9]. The paper [10] mainly focused on the dispute and consensus of non classical biological manipulation technology dominated by silver carp and bighead carp. The paper [11] pointed out that silver and bighead carps were just suitable for controlling cyanobacteria bloom by comparison with the increasing of blue-green algae’s proportion and the forming of microcystis bloom within the enclosures without fish. The paper [12] constructed a new aquatic ecological model to understand the dynamic relationship between Microcystisaeruginosa and filter-feeding fish, which could indirectly show the algae control effect of filter-feeding fish. The paper [13] proposed an aquatic amensalism model to explore the inhibition mechanism of algicidal bacteria on algae. In general, the use of protozoa to control algal blooms is a brand-new control idea, which is worth our in-depth exploration.

The Allee effect is an ecological concept with roots that go back at least to the 1920s, and fifty years have elapsed since the last edition of a book by W.C. Allee, the father of this process in the paper [14]. The paper [15] pointed out that Allee effect is divided into weak Allee effect and strong Allee effect, weak Allee effect refers to the unit individual growth rate at low densities, increasing with population density and always positive, with the population showing a positive growth trend, strong Allee effect refers to the unit individual growth rate at low density, increasing with population density but negative below a critical value, when population density below that becomes negative and tends to extinction. The paper [16] investigated sufficient conditions for the existence of coexisting solutions from the strong and weak Allee effects. The paper [17] established a predation-prey system with Allee effects to study the stability analysis of nonspatial systems and obtain the existence of Hopf branching at coexisting equilibrium points and the stability of branching periodic solutions. The paper [18] shown that the model with strong Allee effect has at most two positive equilibrium point in the first quadrant, while the model with weak Allee effect has at most three positive equilibrium point in the first quadrant. The paper [19] discussed the impacts of Allee effect on co-existence, stability, bistability and bifurcations, and pointed out that the introduction of Allee effect could induce more rich dynamics and compel the model to be more sensitive to initial population densities. The paper [20] pointed out that the model with strong Allee effect could exhibit multiple stability in the first quadrant, and the model with weak Allee effect could undergo saddle knot bifurcation, Hopf bifurcation and Bogdanov-Takes bifurcation in the first quadrant. In short, with more and more examples of Allee phenomenon in natural ecology, more and more researchers pay attention to Allee effect, then more and more excellent achievements will appear in the near future.

Other arrangements of this paper are as follows: In the second section, an aquatic ecological model with Allee effect is built to describe the ecological relationship between algae and protozoa. In the third section, the existence and stability of all possible equilibrium points are studied. In the fourth section, the possible bifurcation dynamic behavior of the model (2.2) is mainly explored. In the fifth section, relevant dynamic simulation tests are carried out to verify the feasibility of theoretical results and demonstrate the evolution trend of population coexistence mode. In the sixth section, we mainly give the main conclusions and make some explanations.

2. Ecological Mathematical Modeling

At present, it is of special significance to apply ecological models to study the problem of biological algae control, this is because that ecological model can form three basic of studying biological system, namely: trophic level analysis, system perspective and dynamic view [21] [22], which can improve the understanding of the ecological interactions between populations and their dependence on internal and external conditions [23] [24]. The paper [12] proposed a new aquatic ecological model to explore the aggregation behavior of algae population. According to the modeling framework of this new aquatic ecological model, we will propose an aquatic ecological model to characterize the dynamic relationship between algae and protozoa (protozoa can feed on algae), which can be described as follow:

$\{\begin{array}{l}\frac{\text{d}N}{\text{d}T}=r_{1}N\left(1-\frac{N}{K_1}\right)\frac{N}{N+a}-\frac{{\alpha }_1\left(N-g\right)P}{c+N-g}-m_{1}N,\hfill \\ \frac{\text{d}P}{\text{d}T}=r_{2}P+\frac{{\beta }_1{\alpha }_1\left(N-g\right)P}{c+N-g}-m_{2}P,\hfill \end{array}$ (2.1)

where $N\left(T\right)$ and $P\left(T\right)$ are density of algae population and protozoa respectively, r₁ is maximum growth rate of algae population, K₁ is maximum environmental capacity for algae population, m₁ is mortality rate of algae population, a is Allee effect coefficient, r₂ is intrinsic growth rate of protozoa, c is saturation coefficient, m₂ is mortality rate of protozoa, ${\alpha }_1$ is capture rate of protozoa preying on algae, ${\beta }_1$ is energy conversion rate, and g is algal aggregation parameter.

For simplicity, we will replace the model (2.1) with the following variable:

$\begin{array}{l}N=cx,P=\frac{r_{1}cy}{{\alpha }_1},T=\frac{t}{r_1},p=\frac{K_1}{c},d=\frac{g}{c},\\ m=\frac{m_1}{r_1},q=\frac{a}{c},b=\frac{r_2}{r_1},e=\frac{{\beta }_1{\alpha }_1}{r_1},n=\frac{m_2}{r_1},\end{array}$

then the model (2.2) is obtained:

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=\frac{x^2}{x+q}\left(1-\frac{x}{p}\right)-\frac{\left(x-d\right)y}{1+x-d}-mx,\hfill \\ \frac{\text{d}y}{\text{d}t}=by+\frac{ey\left(x-d\right)}{1+x-d}-ny.\hfill \end{array}$ (2.2)

For the model (2.2), the existence and stability of all possible equilibrium points will firstly be discussed. Then some critical conditions are given to demonstrate the occurrence of transcritical bifurcation and Hopf bifurcation. Finally, some numerical simulations were implemented to not only verify the feasibility of the theoretical results, but also dynamically evolve ecological dynamic relationship between algae and protozoa, which can abstract out ecological evolution significance represented by bifurcation dynamic evolution behavior.

3. Existence and Stability of All Possible Equilibrium Points

In this section, we will explore the existence and stability of all possible equilibrium points of the model (2.2), which represents the special dynamic relationship between populations.

To obtain all possible equilibrium points of the model (2.2), we list the following equations from the model (2.2):

$\{\begin{array}{l}\frac{x^2}{x+q}\left(1-\frac{x}{p}\right)-\frac{\left(x-d\right)y}{1+x-d}-mx=0,\hfill \\ by+\frac{ey\left(x-d\right)}{1+x-d}-ny=0.\hfill \end{array}$ (2.3)

It is easy to find that the model (2.2) has five possible equilibrium points: $E_0\left(\mathrm{0,0}\right)$, $E_1\left(x_1\mathrm{,0}\right)$, $E_2\left(x_2\mathrm{,0}\right)$, $E_3\left(x_3\mathrm{,0}\right)$, $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$, where

$x_1=\frac{\left(1-m\right)p+\sqrt{\Delta }}{2}>0,x_2=\frac{\left(1-m\right)p-\sqrt{\Delta }}{2}>0,$

$x_3=\frac{p\left(1-m\right)}{2}>0,\Delta ={\left(m-1\right)}^2{p}^2-4mpq>0.$

According to the equation(2.3), the model(2.2) has one internal equilibrium point if $n>b$, $d>\frac{b-n}{b+e-n}$ and $x_2<x_{\ast }<x_2$, where

$x_{\ast }=d-\frac{b-n}{b+e-n},y_{\ast }=\frac{\left[x_{\ast }^2\left(p-x_{\ast }\right)-mp{x}_{\ast }\left(x_{\ast }+q\right)\right]\left(1+x_{\ast }-d\right)}{\left(x_{\ast }+q\right)\left(x_{\ast }-d\right)p}.$

Thus, we can give Theorem 1, which is mainly the critical condition for the existence of all possible equilibrium points.

Theorem 1 1) The boundary equilibrium point $E_0\left(\mathrm{0,0}\right)$ always exists.

2) The boundary equilibrium point $E_1\left(x_1\mathrm{,0}\right)$ and $E_2\left(x_2\mathrm{,0}\right)$ exist if and only if $q<\frac{{\left(1-m\right)}^{2}p}{4m}$ and $0<m<1$.

3) The boundary equilibrium point $E_3\left(x_3\mathrm{,0}\right)$ exists if and only if $q=\frac{{\left(1-m\right)}^{2}p}{4m}$ and $0<m<1$.

4) The internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ exists if and only if $n>b$, $d>\frac{b-n}{b-n+e}$ and $x_2<x_{\ast }<x_1$.

Because the stability of the equilibrium point is determined by the properties of the eigenvalues of its Jacobian matrix, the stability of each equilibrium point is discussed, thus we can get that the Jacobi matrix of the model (2.2) is

$J=\left[\begin{array}{cc}\frac{-2x^3+\left(p-3q\right)x^2+2pqx}{p{\left(x+q\right)}^2}-\frac{y}{{\left(1+x-d\right)}^2}-m& -\frac{x-d}{1+x-d}\\ \frac{ey}{{\left(1+x-d\right)}^2}& b+\frac{e\left(x-d\right)}{1+x-d}-n\end{array}\right].$

On the based of the Jacobi matrix, we can obtain Theorem 2-6, which mainly explore the types and stability of all equilibrium points.

Theorem 2 Under the premise of $n>b$, we have

1) If $b-n+e>0$ and $d>1$ or $e+b-n<0$ and $1<d<\frac{b-n}{e+b-n}$ hold, $E_0\left(\mathrm{0,0}\right)$ is a saddle.

2) If $d<1$ or $e+b-n<0$ and $d>\frac{b-n}{e+b-n}$ hold $E_0\left(\mathrm{0,0}\right)$ is a stable node.

Proof. The Jacobi matrix of the $E_0\left(\mathrm{0,0}\right)$ is:

$J_{E_0}=\left[\begin{array}{cc}-m& \frac{d}{1-d}\\ 0& b-\frac{ed}{1-d}-n\end{array}\right].$

Apparently, the Jacobi matrix of $J_{E_0}$ has two characteristic roots, ${\lambda }_1=-m<0$, ${\lambda }_2=b-\frac{ed}{1-d}-n$. Under the premise of $n>b$, it is easy to know that if $b-n+e>0$ and $d>1$ or $e+b-n<0$ and $1<d<\frac{b-n}{e+b-n}$ hold, the boundary equilibrium point $E_0\left(\mathrm{0,0}\right)$ is a saddle; if $d<1$ or $e+b-n<0$ and $d>\frac{b-n}{e+b-n}$ hold, the boundary equilibrium point $E_0\left(\mathrm{0,0}\right)$ is a stable node.

Theorem 3 Under the premise of $n>b$, we have

1) If $q<\frac{p{\left(m-1\right)}^2}{4m}$ and $\frac{e\left(x_1-d\right)}{1+x_1-d}>n-b$ hold, the boundary equilibrium point $E_1\left(x_1\mathrm{,0}\right)$ is a saddle.

2) If $q<\frac{p{\left(m-1\right)}^2}{4m}$ and $\frac{e\left(x_1-d\right)}{1+x_1-d}<n-b$ hold, the boundary equilibrium point $E_1\left(x_1\mathrm{,0}\right)$ is a stable node.

Proof. The Jacobi matrix of the boundary equilibrium point $E_1\left(x_1\mathrm{,0}\right)$ is:

$J_{E_1\left(x_1,0\right)}=\left[\begin{array}{cc}\frac{-2x_1^3+\left(p-3q\right)x_1^2+2pq{x}_1}{p{\left(x_1+q\right)}^2}-m& -\frac{x_1-d}{1+x_1-d}\\ 0& b+\frac{e\left(x_1-d\right)}{1+x_1-d}-n\end{array}\right].$

Then, the Jacobi matrix of $J_{E_1\left(x_1,0\right)}$ has two characteristic roots,

${\lambda }_1=\frac{-2x_1^3+\left(p-3q\right)x_1^2+2pq{x}_1}{p{\left(x_1+q\right)}^2}-m,{\lambda }_2=b+\frac{e\left(x_1-d\right)}{1+x_1-d}-n.$

Firstly, we will analyze the positive and negativity of ${\lambda }_1$. It is easy to know that $p{\left(x_1+q\right)}^2>0$ and $x_1$ satisfies $x_1^2+p\left(m-1\right)x_1+mpq=0$, we put it in ${\lambda }_1$ and sort it out.

${\lambda }_1=\frac{\left[-{\left(m-1\right)}^2{p}^2+\left(3m-1\right)pq\right]x_1+2mp{q}^2+m{p}^2\left(1-m\right)q}{p{\left(x_1+q\right)}^2}.$

According to $q<\frac{p{\left(m-1\right)}^2}{4m}$, and

$x_1=\frac{\left(1-m\right)p+\sqrt{\Delta }}{2}>0,$

therefore

$\left[-{\left(m-1\right)}^2{p}^2+\left(3m-1\right)pq\right]x_1+2mp{q}^2+m{p}^2\left(1-m\right)q<\frac{p^2{\left(1-m\right)}^2\left(m+1\right)}{8m}\left(-\sqrt{\Delta }\right)<0,$

so ${\lambda }_1<0$.

Secondly, we analyze the positive and negativity of ${\lambda }_2$. Through calculation, we can get that if $q<\frac{p{\left(m-1\right)}^2}{4m}$ and $\frac{e\left(x_1-d\right)}{1+x_1-d}>n-b$ hold, we have ${\lambda }_2>0$, then the boundary equilibrium point $E_1\left(x_1\mathrm{,0}\right)$ is a saddle; If $q<\frac{p{\left(m-1\right)}^2}{4m}$ and $\frac{e\left(x_1-d\right)}{1+x_1-d}<n-b$ hold, we have ${\lambda }_2<0$, then the boundary equilibrium point $E_1\left(x_1\mathrm{,0}\right)$ is a stable node.

Theorem 4 Under the premise of $n>b$, we have

1) If $q<\frac{p{\left(m-1\right)}^2}{4m}$, $m\ne \frac{1}{3}$ and $\frac{e\left(x_2-d\right)}{1+x_2-d}>n-b$ hold, the boundary equilibrium point $E_2\left(x_2\mathrm{,0}\right)$ is a saddle.

2) If $q<\frac{p{\left(m-1\right)}^2}{4m}$, $m\ne \frac{1}{3}$ and $\frac{e\left(x_2-d\right)}{1+x_2-d}<n-b$ hold, the boundary equilibrium point $E_2\left(x_2\mathrm{,0}\right)$ is a stable node.

Proof. The Jacobi matrix of the boundary equilibrium point $E_2\left(x_2\mathrm{,0}\right)$ is

$J_{E_2\left(x_2\mathrm{,0}\right)}=\left[\begin{array}{cc}\frac{-2x_2^3+\left(p-3q\right)x_2^2+2pq{x}_2}{p{\left(x_2+q\right)}^2}-m& -\frac{x_2-d}{1+x_2-d}\\ 0& b+\frac{e\left(x_2-d\right)}{1+x_2-d}-n\end{array}\right]\mathrm{.}$

Then, the Jacobi matrix of $J_{E_2\left(x_2\mathrm{,0}\right)}$ has two characteristic roots,

${\lambda }_1=\frac{-2x_2^3+\left(p-3q\right)x_2^2+2pq{x}_2}{p{\left(x_2+q\right)}^2}-m,{\lambda }_2=b+\frac{e\left(x_2-d\right)}{1+x_2-d}-n.$

Firstly, we will analyze the positive and negativity of ${\lambda }_1$, because $p{\left(x_2+q\right)}^2>0$, owing to $x_2$ satisfies $x_2^2+p\left(m-1\right)x_2+mpq=0$, we put it in ${\lambda }_1$ and sort it out, we can get

${\lambda }_1=\frac{\left[-{\left(m-1\right)}^2{p}^2+\left(3m-1\right)pq\right]x_2+2mp{q}^2+m{p}^2\left(1-m\right)q}{p{\left(x_2+q\right)}^2}.$

1) If ${\lambda }_1>0$, then $\left[-{\left(m-1\right)}^2{p}^2+\left(3m-1\right)pq\right]x_2>-\left(2mp{q}^2+m{p}^2\left(1-m\right)q\right)$.

a) If $\left[-{\left(m-1\right)}^2{p}^2+\left(3m-1\right)pq\right]>0$, then $\left(3m-1\right)q>{\left(m-1\right)}^{2}p$, $x_2>\frac{-\left[2mp{q}^2+m{p}^2\left(1-m\right)q\right]}{-{\left(m-1\right)}^2{p}^2+\left(3m-1\right)pq}$.

i) When $0<m<\frac{1}{3}$, $q<\frac{p{\left(m-1\right)}^2}{3m-1}<0$ contradicts $q>0$.

ii) When $\frac{1}{3}<m<1$, $q>\frac{p{\left(m-1\right)}^2}{3m-1}<0$ has no intersection with $0<q<\frac{p{\left(m-1\right)}^2}{4m}$, so there is no solution.

b) If $\left[-{\left(m-1\right)}^2{p}^2+\left(3m-1\right)pq\right]<0$, we have $\left(3m-1\right)q<{\left(m-1\right)}^{2}p$, $x_2<\frac{-\left[2mp{q}^2+m{p}^2\left(1-m\right)q\right]}{-{\left(m-1\right)}^2{p}^2+\left(3m-1\right)pq}$, also because $x_2=\frac{\left(1-m\right)p-\sqrt{\Delta }}{2}$, therefore, it must be satisfied that

$\frac{\left(1-m\right)p-\sqrt{\Delta }}{2}+\frac{2mp{q}^2+m{p}^2\left(1-m\right)q}{-{\left(m-1\right)}^2{p}^2+\left(3m-1\right)pq}<0.$

It can be deformed

$\frac{\left[\left(1-m\right)p-\sqrt{\Delta }\right]\cdot \left[-{\left(m-1\right)}^{2}p+\left(3m-1\right)q\right]+4m{q}^2+2mpq\left(1-m\right)}{2\cdot \left[-{\left(m-1\right)}^{2}p+\left(3m-1\right)q\right]}<0$

for $\left[-{\left(m-1\right)}^{2}p+\left(3m-1\right)q\right]<0$. So it needs to be proved that

$\left[\left(1-m\right)p-\sqrt{\Delta }\right]\cdot \left[-{\left(m-1\right)}^{2}p+\left(3m-1\right)q\right]+4m{q}^2+2mpq\left(1-m\right)>0,$

simplified the upper type

$\left(4mq-{\left(m-1\right)}^{2}p\right)\cdot \left(p+q\right)>0,$

because $\left(p+q\right)>0$, and $4mq-{\left(1-m\right)}^{2}p>0$, $q>\frac{p{\left(m-1\right)}^2}{4m}$ contradicts $0<q<\frac{p{\left(m-1\right)}^2}{4m}$, so ${\lambda }_1<0$.

2) If ${\lambda }_1<0$, then $\left[-{\left(m-1\right)}^2{p}^2+\left(3m-1\right)pq\right]x_2<-\left(2mp{q}^2+m{p}^2\left(1-m\right)q\right)$,

a) If $\left[-{\left(m-1\right)}^2{p}^2+\left(3m-1\right)pq\right]>0$, then $x_2<\frac{-\left[2mp{q}^2+m{p}^2\left(1-m\right)q\right]}{-{\left(m-1\right)}^2{p}^2+\left(3m-1\right)pq<0}$, it contradicts $x_2>0$.

b) If $\left[-{\left(m-1\right)}^2{p}^2+\left(3m-1\right)pq\right]<0$, we have $\left(3m-1\right)q<{\left(m-1\right)}^{2}p$, $x_2>\frac{-\left[2mp{q}^2+m{p}^2\left(1-m\right)q\right]}{-{\left(m-1\right)}^2{p}^2+\left(3m-1\right)pq}$, also because $x_2=\frac{\left(1-m\right)p-\sqrt{\Delta }}{2}$, therefore, it must be satisfied that

$\frac{\left(1-m\right)p-\sqrt{\Delta }}{2}+\frac{2mp{q}^2+m{p}^2\left(1-m\right)q}{-{\left(m-1\right)}^2{p}^2+\left(3m-1\right)pq}>0.$

It can be deformed

$\frac{\left[\left(1-m\right)p-\sqrt{\Delta }\right]\cdot \left[-{\left(m-1\right)}^{2}p+\left(3m-1\right)q\right]+4m{q}^2+2mpq\left(1-m\right)}{2\cdot \left[-{\left(m-1\right)}^{2}p+\left(3m-1\right)q\right]}>0$

for $\left[-{\left(m-1\right)}^{2}p+\left(3m-1\right)q\right]<0$. So it needs to be proved that

$\left[\left(1-m\right)p-\sqrt{\Delta }\right]\cdot \left[-{\left(m-1\right)}^{2}p+\left(3m-1\right)q\right]+4m{q}^2+2mpq\left(1-m\right)<0,$

simplified the upper type

$\left(4mq-{\left(m-1\right)}^{2}p\right)\cdot \left(p+q\right)<0,$

because $\left(p+q\right)>0$, $4mq-{\left(1-m\right)}^{2}p<0$, so $q<\frac{p{\left(m-1\right)}^2}{4m}$, therefore, this situation is true. On account of $\left(3m-1\right)q<{\left(m-1\right)}^{2}p$,

i) When $0<m<\frac{1}{3}$, so $q>\frac{p{\left(m-1\right)}^2}{3m-1}<0$ and $0<q<\frac{p{\left(m-1\right)}^2}{4m}$, they take the intersection to get $0<q<\frac{p{\left(m-1\right)}^2}{4m}$.

ii) When $\frac{1}{3}<m<1$, $q<\frac{p{\left(m-1\right)}^2}{3m-1}<0$ and $0<q<\frac{p{\left(m-1\right)}^2}{4m}$, they take the intersection to get $0<q<\frac{p{\left(m-1\right)}^2}{4m}$.

In a word, when $0<q<\frac{p{\left(m-1\right)}^2}{4m}$ and $m\ne \frac{1}{3}$ hold, we have ${\lambda }_1<0$.

Secondly, we analyze the positive and negativity of ${\lambda }_2$. When $q<\frac{p{\left(m-1\right)}^2}{4m}$ and $\frac{e\left(x_2-d\right)}{1+x_2-d}>n-b$ hold, the boundary equilibrium point $E_2\left(x_2\mathrm{,0}\right)$ is a saddle. when $q<\frac{p{\left(m-1\right)}^2}{4m}$ and $\frac{e\left(x_2-d\right)}{1+x_2-d}<n-b$ hold, the boundary equilibrium point $E_2\left(x_2\mathrm{,0}\right)$ is a stable node.

Theorem 5 Under the premise of $n>b$, we have

1) If $q<\frac{p{\left(m-1\right)}^2}{4m}$ and $\frac{e\left(x_3-d\right)}{1+x_3-d}>n-b$ hold, the boundary equilibrium point $E_3\left(x_3\mathrm{,0}\right)$ is a saddle.

2) If $q<\frac{p{\left(m-1\right)}^2}{4m}$ and $\frac{e\left(x_3-d\right)}{1+x_3-d}<n-b$ hold, the boundary equilibrium point $E_3\left(x_3\mathrm{,0}\right)$ is a stable node.

Proof. The Jacobi matrix of the boundary equilibrium point $E_3\left(x_3\mathrm{,0}\right)$ is

$J_{E_3\left(x_3,0\right)}=\left[\begin{array}{cc}\frac{-2x_3^3+\left(p-3q\right)x_3^2+2pq{x}_3}{p{\left(x_3+q\right)}^2}-m& -\frac{x_3-d}{1+x_3-d}\\ 0& b+\frac{e\left(x_3-d\right)}{1+x_3-d}-n\end{array}\right].$

It is easy to obtain that the Jacobi matrix of $J_{E_3\left(x_3,0\right)}$ has two characteristic roots,

${\lambda }_1=\frac{-2x_3^3+\left(p-3q\right)x_3^2+2pq{x}_3}{p{\left(x_2+q\right)}^2}-m,{\lambda }_2=b+\frac{e\left(x_3-d\right)}{1+x_3-d}-n.$

Firstly, we will analyze the positive and negativity of ${\lambda }_1$, because $x_3$ satisfies $x_3^2+p\left(m-1\right)x_3+mpq=0$ and $x_3^2=-p\left(m-1\right)x_3-mpq$, we put it in ${\lambda }_1$ and sort it out, we can get

${\lambda }_1=\frac{\left[-{\left(m-1\right)}^2{p}^2+\left(3m-1\right)pq\right]x_3+2mp{q}^2+m{p}^2\left(1-m\right)q}{p{\left(x_3+q\right)}^2}.$

Owing to $x_3=\frac{\left(1-m\right)p}{2}$, simplified $\left[-{\left(m-1\right)}^2{p}^2+\left(3m-1\right)pq\right]x_3+2mp{q}^2+m{p}^2\left(1-m\right)q$, we can get ${\lambda }_1=\frac{-{\left(1-m\right)}^4{p}^3}{p{\left(x_3+q\right)}^2}<0$, so ${\lambda }_1<0$.

Secondly, we analyze the positive and negativity of ${\lambda }_2$. When $q=\frac{p{\left(m-1\right)}^2}{4m}$ and $\frac{e\left(x_3-d\right)}{1+x_3-d}>n-b$ hold, the boundary equilibrium point $E_3\left(x_3\mathrm{,0}\right)$ is a saddle. When $q=\frac{p{\left(m-1\right)}^2}{4m}$, and $\frac{e\left(x_3-d\right)}{1+x_3-d}<n-b$ hold, the boundary equilibrium point $E_3\left(x_3\mathrm{,0}\right)$ is a stable node.

Theorem 6 Under the condition of the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$,

1) If $Det\left(J_{E_{\ast }}\right)>0$, $Tr\left(J_{E_{\ast }}\right)<0$ hold, the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ is a stable node(focus).

2) If $Det\left(J_{E_{\ast }}\right)>0$, $Tr\left(J_{E_{\ast }}\right)>0$ hold, the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ is an unstable node(focus).

Proof. The Jacobi matrix of the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ is,

$J_{E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)}=\left[\begin{array}{cc}\frac{-2x_{\ast }^3+\left(p-3q\right)x_{\ast }^2+2pq{x}_{\ast }}{p{\left(x_{\ast }+q\right)}^2}-\frac{y_{\ast }}{{\left(1+x_{\ast }-d\right)}^2}-m& -\frac{x_{\ast }-d}{1+x_{\ast }-d}\\ 0& b+\frac{e\left(x_{\ast }-d\right)}{1+x_{\ast }-d}-n\end{array}\right]$

where $x_{\ast }=d-\frac{b-n}{b+e-n}$, $y_{\ast }=\frac{\left[x_{\ast }^2\left(p-x_{\ast }\right)-mp{x}_{\ast }\left(x_{\ast }+q\right)\right]\left(1+x_{\ast }-d\right)}{\left(x_{\ast }+q\right)\left(x_{\ast }-d\right)p}$. The following characteristic equation is obtained as follow

$\begin{array}{l}{\lambda }^2-\left[\frac{-2x_{\ast }^3+\left(p-3q\right)x_{\ast }^2+2pq{x}_{\ast }}{p{\left(x_{\ast }+q\right)}^2}-\frac{y_{\ast }}{{\left(1+x_{\ast }-d\right)}^2}-m\right]\cdot \lambda \\ +\frac{x_{\ast }-d}{1+x_{\ast }-d}\cdot \frac{e{y}_{\ast }}{{\left(1+x_{\ast }-d\right)}^2}=0,\end{array}$

and

$Tr\left(J_{E_{\ast }}\right)={\lambda }_1+{\lambda }_2,Det\left(J_{E_{\ast }}\right)={\lambda }_1{\lambda }_2,$

here

$Det\left(J_{E_{\ast }}\right)=\frac{x_{\ast }-d}{1+x_{\ast }-d}\cdot \frac{e{y}_{\ast }}{{\left(1+x_{\ast }-d\right)}^2},$

for $y_{\ast }>0$. Then we can get

$q<\frac{-x_{\ast }^2-mp{x}_{\ast }+p{x}_{\ast }}{mp},$

for

$T=\frac{-2x_{\ast }^3+\left(p-3q\right)x_{\ast }^2+2pq{x}_{\ast }}{p{\left(x_{\ast }+q\right)}^2}-\frac{y_{\ast }}{{\left(1+x_{\ast }-d\right)}^2}-m.$

Thus if $Det\left(J_{E_{\ast }}\right)>0$, $Tr\left(J_{E_{\ast }}\right)<0$ hold, the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ is a stable node(focus); if $Det\left(J_{E_{\ast }}\right)>0$, $Tr\left(J_{E_{\ast }}\right)>0$ hold, the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ is an unstable node(focus).

4. Local Bifurcation Analysis

In this section, we will choose parameter d as a bifurcation control parameter to investigate the bifurcation dynamics evolution characteristics of the model (2.2), and give the threshold conditions for transcritical bifurcation and Hopf bifurcation of the model (2.2).

4.1. Transcritical Bifurcation

Theorem 7 1) The model (2.2) undergoes a transcritical bifurcation at the equilibrium point $E_1\left(x_1\mathrm{,0}\right)$ when $d=d_{TC1}=x_1+\frac{b-n}{b+e-n}$.

2) The model (2.2) undergoes a transcritical bifurcation at the equilibrium point $E_2\left(x_2\mathrm{,0}\right)$ when $d=d_{TC2}=x_2+\frac{b-n}{b+e-n}$.

3) The model (2.2) undergoes a transcritical bifurcation at the equilibrium point $E_3\left(x_3\mathrm{,0}\right)$ when $d=d_{TC3}=x_3+\frac{b-n}{b+e-n}$.

Proof:

1) On the basis of the Theorem 3, when $d=d_{TC1}=x_1+\frac{b-n}{b+e-n}$, the Jacobi matrix of the equilibrium point $E_1$ is

$J_{E_{TC1}}=\left[\begin{array}{cc}\frac{-2x_1^3+\left(p-3q\right)x_1^2+2pq{x}_1}{p{\left(x_1+q\right)}^2}-m& -\frac{x_1-d}{1+x_1-d}\\ 0& 0\end{array}\right]\mathrm{,}$

suppose V and W are eigenvectors of $J_{E_{TC1}}$ and $J_{E_{TC1}}^{\text{T}}$, then

$J_{E_{TC1}}V=0\cdot V,J_{E_{TC1}}^{\text{T}}W=0\cdot W.$

Then we can get

$\begin{array}{l}V=\left[\begin{array}{c}{v}_1\\ v_2\end{array}\right]=\left[\begin{array}{c}\frac{x_1-d}{1+x_1-d}\\ \frac{-2x_1^3+\left(p-3q\right)x_1^2+2pq{x}_1}{p{\left(x_1+q\right)}^2}-m\end{array}\right],\\ W=\left[\begin{array}{c}{w}_1\\ w_2\end{array}\right]=\left[\begin{array}{c}0\\ 1\end{array}\right].\end{array}$

Due to

$F_d\left(E_1;d_{TC1}\right)=\left[\begin{array}{c}{F}_{1d}\\ F_{2d}\end{array}\right]={\left[\begin{array}{c}\frac{y}{{\left(1+x-d\right)}^2}\\ \frac{-ey}{{\left(1+x-d\right)}^2}\end{array}\right]}_{\left(E_1;d_{TC1}\right)}=\left[\begin{array}{c}0\\ 0\end{array}\right],$

so

$\begin{array}{c}D{F}_d\left(E_1;d_{TC1}\right)V={\left[\begin{array}{cc}{F}_{1d_x}& F_{1d_y}\\ F_{2d_x}& F_{2d_y}\end{array}\right]}_{\left(E_1;d_{TC1}\right)}\left[\begin{array}{c}{v}_1\\ v_2\end{array}\right]\\ =\left[\begin{array}{c}\frac{1}{{\left(1+x_1-d\right)}^2}\left(\frac{-2x_1^3+\left(p-3q\right)x_1^2+2pq{x}_1}{p{\left(x_1+q\right)}^2}-m\right)\\ \frac{-e}{{\left(1+x_1-d\right)}^2}\left(\frac{-2x_1^3+\left(p-3q\right)x_1^2+2pq{x}_1}{p{\left(x_1+q\right)}^2}-m\right)\end{array}\right],\end{array}$

$\begin{array}{l}{D}^2{F}_d\left(E_1;d_{TC1}\right)\left(V,V\right)={\left[\begin{array}{c}\frac{{\partial }^2{F}_1}{\partial x^2}{v}_1{v}_1+2\frac{{\partial }^2{F}_1}{\partial x\partial y}{v}_1{v}_2+\frac{{\partial }^2{F}_1}{\partial y^2}{v}_2{v}_2\\ \frac{{\partial }^2{F}_2}{\partial x^2}{v}_1{v}_1+2\frac{{\partial }^2{F}_2}{\partial x\partial y}{v}_1{v}_2+\frac{{\partial }^2{F}_2}{\partial y^2}{v}_2{v}_2\end{array}\right]}_{\left(E_1;d_{TC1}\right)}\\ =\left[\begin{array}{c}\frac{-2x_1^3-6q{x}_1^2-6q^2{x}_1+2p{q}^2}{p{\left(x_1+q\right)}^3}{v}_1{v}_1-\frac{2}{{\left(1+x_1-d\right)}^2}{v}_1{v}_2\\ \frac{2e}{{\left(1+x_1-d\right)}^2}{v}_1{v}_2\end{array}\right].\end{array}$

Thus, we can reach the following conclusions:

$W^{\text{T}}{F}_d\left(E_1;d_{TC1}\right)=\left[0,1\right]\left[\begin{array}{c}0\\ 0\end{array}\right]=0,$

$\begin{array}{l}{W}^{\text{T}}\left[D{F}_d\left(E_1;d_{TC1}\right)V\right]\\ =\frac{-e}{{\left(1+x_1-d\right)}^2}\left(\frac{-2x_1^3+\left(p-3q\right)x_1^2+2pq{x}_1}{p{\left(x_1+q\right)}^2}-m\right)\ne 0,\end{array}$

$\begin{array}{l}{W}^{\text{T}}\left[D^2{F}_d\left(E_1;d_{TC1}\right)\left(V,V\right)\right]\\ =\frac{2e\left(x_1-d\right)}{{\left(1+x_1-d\right)}^3}\left(\frac{-2x_1^3+\left(p-3q\right)x_1^2+2pq{x}_1}{p{\left(x_1+q\right)}^2}-m\right)\ne 0.\end{array}$

According to Sotomayors theorem, when $d=d_{TC1}=x_1+\frac{b-n}{b+e-n}$, then the model (2.2) undergoes a transcritical bifurcation at the equilibrium point $E_1\left(x_1\mathrm{,0}\right)$.

2) On the basis of the Theorem 4, we next prove that the model (2.2) will undergo a transcritical bifurcation at the equilibrium point $E_2\left(x_2\mathrm{,0}\right)$. When $d=d_{TC2}=x_2+\frac{b-n}{b+e-n}$, the Jacobi matrix of the equilibrium point $E_2$ is

$J_{E_{TC2}}=\left[\begin{array}{cc}\frac{-2x_2^3+\left(p-3q\right)x_2^2+2pq{x}_2}{p{\left(x_2+q\right)}^2}-m& -\frac{x_2-d}{1+x_2-d}\\ 0& 0\end{array}\right],$

suppose V and W are eigenvectors of $J_{E_{TC2}}$ and $J_{E_{TC2}}^{\text{T}}$, then

$J_{E_{TC2}}V=0\cdot V,J_{E_{TC2}}^{\text{T}}W=0\cdot W,$

then

$V=\left[\begin{array}{c}{v}_1\\ v_2\end{array}\right]=\left[\begin{array}{c}\frac{x_2-d}{1+x_2-d}\\ \frac{-2x_2^3+\left(p-3q\right)x_2^2+2pq{x}_2}{p{\left(x_2+q\right)}^2}-m\end{array}\right],W=\left[\begin{array}{c}{w}_1\\ w_2\end{array}\right]=\left[\begin{array}{c}0\\ 1\end{array}\right].$

Owing to

$F_d\left(E_2;d_{TC2}\right)=\left[\begin{array}{c}{F}_{1d}\\ F_{2d}\end{array}\right]={\left[\begin{array}{c}\frac{y}{{\left(1+x-d\right)}^2}\\ \frac{-ey}{{\left(1+x-d\right)}^2}\end{array}\right]}_{\left(E_2;d_{TC2}\right)}=\left[\begin{array}{c}0\\ 0\end{array}\right],$

and

$\begin{array}{c}D{F}_d\left(E_2;d_{TC2}\right)V={\left[\begin{array}{cc}{F}_{1d_x}& F_{1d_y}\\ F_{2d_x}& F_{2d_y}\end{array}\right]}_{\left(E_2;d_{TC2}\right)}\left[\begin{array}{c}{v}_1\\ v_2\end{array}\right]\\ =\left[\begin{array}{c}\frac{1}{{\left(1+x_2-d\right)}^2}\left(\frac{-2x_2^3+\left(p-3q\right)x_2^2+2pq{x}_2}{p{\left(x_2+q\right)}^2}-m\right)\\ \frac{-e}{{\left(1+x_2-d\right)}^2}\left(\frac{-2x_2^3+\left(p-3q\right)x_2^2+2pq{x}_2}{p{\left(x_2+q\right)}^2}-m\right)\end{array}\right],\end{array}$

$\begin{array}{l}{D}^2{F}_d\left(E_2;d_{TC2}\right)\left(V,V\right)={\left[\begin{array}{c}\frac{{\partial }^2{F}_1}{\partial x^2}{v}_1{v}_1+2\frac{{\partial }^2{F}_1}{\partial x\partial y}{v}_1{v}_2+\frac{{\partial }^2{F}_1}{\partial y^2}{v}_2{v}_2\\ \frac{{\partial }^2{F}_2}{\partial x^2}{v}_1{v}_1+2\frac{{\partial }^2{F}_2}{\partial x\partial y}{v}_1{v}_2+\frac{{\partial }^2{F}_2}{\partial y^2}{v}_2{v}_2\end{array}\right]}_{\left(E_2;d_{TC2}\right)}\\ =\left[\begin{array}{c}\frac{-2x_2^3-6q{x}_2^2-6q^2{x}_2+2p{q}^2}{p{\left(x_2+q\right)}^3}{v}_1{v}_1-\frac{2}{{\left(1+x_2-d\right)}^2}{v}_1{v}_2\\ \frac{2e}{{\left(1+x_2-d\right)}^2}{v}_1{v}_2\end{array}\right],\end{array}$

we can obtain

$W^{\text{T}}{F}_d\left(E_2;d_{TC2}\right)=\left[0,1\right]\left[\begin{array}{c}0\\ 0\end{array}\right]=0,$

$\begin{array}{l}{W}^{\text{T}}\left[D{F}_d\left(E_2;d_{TC2}\right)V\right]\\ =\frac{-e}{{\left(1+x_2-d\right)}^2}\left(\frac{-2x_2^3+\left(p-3q\right)x_2^2+2pq{x}_2}{p{\left(x_2+q\right)}^2}-m\right)\ne 0,\end{array}$

$\begin{array}{l}{W}^{\text{T}}\left[D^2{F}_d\left(E_2;d_{TC2}\right)\left(V,V\right)\right]\\ =\frac{2e\left(x_2-d\right)}{{\left(1+x_2-d\right)}^3}\left(\frac{-2x_2^3+\left(p-3q\right)x_2^2+2pq{x}_2}{p{\left(x_2+q\right)}^2}-m\right)\ne 0.\end{array}$

According to Sotomayors theorem, when $d=d_{TC2}=x_2+\frac{b-n}{b+e-n}$, then the model (2.2) undergoes a transcritical bifurcation at the equilibrium point $E_2\left(x_2\mathrm{,0}\right)$.

3) On the basis of the Theorem 5, we will prove that the model (2.2) will undergo a transcritical bifurcation at the equilibrium point $E_3\left(x_3\mathrm{,0}\right)$. When $d=d_{TC3}=x_3+\frac{b-n}{b+e-n}$ holds, the equilibrium point $E_1\left(x_1\mathrm{,0}\right)$ and $E_2\left(x_2\mathrm{,0}\right)$ will coincide as an equilibrium point, here the Jacobian matrix at $E_3\left(x_3\mathrm{,0}\right)$ is

$J_{E_{TC3}}=\left[\begin{array}{cc}\frac{-2x_3^3+\left(p-3q\right)x_3^2+2pq{x}_3}{p{\left(x_3+q\right)}^2}-m& -\frac{x_3-d}{1+x_3-d}\\ 0& 0\end{array}\right],$

suppose V and W are eigenvectors of $J_{E_{TC3}}$ and $J_{E_{TC3}}^{\text{T}}$, then

$J_{E_{TC3}}V=0\cdot V,J_{E_{TC3}}^{\text{T}}W=0\cdot W,$

so

$V=\left[\begin{array}{c}{v}_1\\ v_2\end{array}\right]=\left[\begin{array}{c}\frac{x_3-d}{1+x_3-d}\\ \frac{-2x_3^3+\left(p-3q\right)x_3^2+2pq{x}_3}{p{\left(x_3+q\right)}^2}-m\end{array}\right],W=\left[\begin{array}{c}{w}_1\\ w_2\end{array}\right]=\left[\begin{array}{c}0\\ 1\end{array}\right].$

Because of

$F_d\left(E_3;d_{TC3}\right)=\left[\begin{array}{c}{F}_{1d}\\ F_{2d}\end{array}\right]={\left[\begin{array}{c}\frac{y}{{\left(1+x-d\right)}^2}\\ \frac{-ey}{{\left(1+x-d\right)}^2}\end{array}\right]}_{\left(E_3;d_{TC3}\right)}=\left[\begin{array}{c}0\\ 0\end{array}\right],$

then, we can obtain

$\begin{array}{c}D{F}_d\left(E_3;d_{TC3}\right)V={\left[\begin{array}{cc}{F}_{1d_x}& F_{1d_y}\\ F_{2d_x}& F_{2d_y}\end{array}\right]}_{\left(E_3;d_{TC3}\right)}\left[\begin{array}{c}{v}_1\\ v_2\end{array}\right]\\ =\left[\begin{array}{c}\frac{1}{{\left(1+x_3-d\right)}^2}\left(\frac{-2x_3^3+\left(p-3q\right)x_3^2+2pq{x}_3}{p{\left(x_3+q\right)}^2}-m\right)\\ \frac{-e}{{\left(1+x_3-d\right)}^2}\left(\frac{-2x_3^3+\left(p-3q\right)x_3^2+2pq{x}_3}{p{\left(x_3+q\right)}^2}-m\right)\end{array}\right],\end{array}$

$\begin{array}{l}{D}^2{F}_d\left(E_3;d_{TC3}\right)\left(V,V\right)={\left[\begin{array}{c}\frac{{\partial }^2{F}_1}{\partial x^2}{v}_1{v}_1+2\frac{{\partial }^2{F}_1}{\partial x\partial y}{v}_1{v}_2+\frac{{\partial }^2{F}_1}{\partial y^2}{v}_2{v}_2\\ \frac{{\partial }^2{F}_2}{\partial x^2}{v}_1{v}_1+2\frac{{\partial }^2{F}_2}{\partial x\partial y}{v}_1{v}_2+\frac{{\partial }^2{F}_2}{\partial y^2}{v}_2{v}_2\end{array}\right]}_{\left(E_3;d_{TC3}\right)}\\ =\left[\begin{array}{c}\frac{-2x_3^3-6q{x}_3^2-6q^2{x}_3+2p{q}^2}{p{\left(x_3+q\right)}^3}{v}_1{v}_1-\frac{2}{{\left(1+x_3-d\right)}^2}{v}_1{v}_2\\ \frac{2e}{{\left(1+x_3-d\right)}^2}{v}_1{v}_2\end{array}\right].\end{array}$

Thus, we can get the following conclusions:

$W^{\text{T}}{F}_d\left(E_3;d_{TC3}\right)=\left[0,1\right]\left[\begin{array}{c}0\\ 0\end{array}\right]=0,$

$\begin{array}{l}{W}^{\text{T}}\left[D{F}_d\left(E_3;d_{TC3}\right)V\right]\\ =\frac{-e}{{\left(1+x_3-d\right)}^2}\left(\frac{-2x_3^3+\left(p-3q\right)x_3^2+2pq{x}_3}{p{\left(x_3+q\right)}^2}-m\right)\ne 0,\end{array}$

$\begin{array}{l}{W}^{\text{T}}\left[D^2{F}_d\left(E_3;d_{TC3}\right)\left(V,V\right)\right]\\ =\frac{2e\left(x_3-d\right)}{{\left(1+x_3-d\right)}^3}\left(\frac{-2x_3^3+\left(p-3q\right)x_3^2+2pq{x}_3}{p{\left(x_3+q\right)}^2}-m\right)\ne 0.\end{array}$

According to Sotomayors theorem, when $d=d_{TC3}=x_3+\frac{b-n}{b+e-n}$, then the model (2.2) undergoes a transcritical bifurcation at the equilibrium point $E_3\left(x_3\mathrm{,0}\right)$.

4.2. Hopf Bifurcation

According to the Theorem 6, the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ can lose its stability, hence the model (2.2) may occur a Hopf bifurcation under certain conditions.

Theorem 8 Under the conditions of the Theorem 6, the internal equilibrium point $E_{\ast }$ can change its stability when the controlling parameter d passes through a critical value $d=d_{Hp}$, then the model (2.2) will undergo a Hopf bifurcation, where ${Tr\left(J_{E_{\ast }}\right)|}_{d=d_{Hp}}=0$.

Proof: To determine the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ can change its stability by a Hopf bifurcation, we need to prove the cross-sectional condition of Hopf bifurcation:

$Tr\left(J_{E_{\ast }}\right)=\frac{-2x_{\ast }^3+\left(p-3q\right)x_{\ast }^2+2pq{x}_{\ast }}{p{\left(x_{\ast }+q\right)}^2}-\frac{y_{\ast }}{{\left(1+x_{\ast }-d\right)}^2}-m=0,$

${\frac{\text{d}}{\text{d}q}\left[Tr\left(J_{E_{\ast }}\right)\right]|}_{d=d_{Hp}}={\left(\frac{x_{\ast }^3-2pq{x}_{\ast }+3q{x}_{\ast }^2}{p{\left(x_{\ast }+q\right)}^2}+\frac{x_{*}^3-p{x}_{\ast }^2}{p\left(1+x_{\ast }-d\right){\left(x_{\ast }+q\right)}^2\left(x_{\ast }-d\right)}\right)|}_{d=d_{Hp}}\ne 0,$

hence the model (2.2) can occur a Hopf bifurcation at $q=q_{Hp}$.

Next, we discuss the stability of the limit cycle by computing the first Lyapunov coefficient of the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$. Translating the origin of coordinates at this equilibrium point through the following transformation $x=x_d-x_{\ast }$, $y=y_d-y_{\ast }$,

we can obtain

$\begin{array}{l}\{\begin{array}{l}{\stackrel{\dot{}}{x}}_d={\alpha }_{10}{x}_d+{\alpha }_{01}{y}_d+{\alpha }_{20}{x}_d^2+{\alpha }_{11}{x}_d{y}_d+{\alpha }_{02}{y}_d^2+{\alpha }_{30}{x}_d^3+{\alpha }_{21}{x}_d^2{y}_d+{\alpha }_{12}{x}_d{y}_d^2+{\alpha }_{03}{y}_d^3+P\left(x_d,y_d\right)\\ {\stackrel{\dot{}}{y}}_d={\beta }_{10}{x}_d+{\beta }_{01}{y}_d+{\beta }_{20}{x}_d^2+{\beta }_{11}{x}_d{y}_d+{\beta }_{02}{y}_d^2+{\beta }_{30}{x}_d^3+{\beta }_{21}{x}_d^2{y}_d+{\beta }_{12}{x}_d{y}_d^2+{\beta }_{03}{y}_d^3+Q\left(x_d,y_d\right)\end{array}\hfill \end{array}$

where

$\begin{array}{l}{\alpha }_{10}=\frac{-pq{x}_{*}-2q{x}_{*}^2+x_{*}^3}{p\left(q-x_{*}\right)}+\frac{y_{*}}{{\left(1-x_{*}-d\right)}^2}-m,{\alpha }_{01}=\frac{x_{*}+d}{1-x_{*}-d},\\ {\alpha }_{20}=\frac{q\left(p+2x_{*}\right)}{p{\left(q-x_{*}\right)}^2}+\frac{x_{*}^2}{{\left(q-x_{*}\right)}^3}-\frac{y_{*}}{{\left(1-x_{*}-d\right)}^3},{\alpha }_{11}=\frac{-1}{{\left(1-x_{*}-d\right)}^3},\\ {\alpha }_{30}=\frac{q^2-pq-2x_{*}^3}{p{\left(q-x_{*}\right)}^3}+\frac{x_{*}^2\left(q-p-2x_{*}\right)}{p{\left(q-x_{*}\right)}^4}+\frac{y_{*}}{{\left(1-x_{*}-d\right)}^4},{\alpha }_{21}=\frac{1}{{\left(1-x_{*}-d\right)}^3},\\ {\alpha }_{02}={\alpha }_{12}={\alpha }_{03}=0,\end{array}$

and

$\begin{array}{l}{\beta }_{10}=\frac{-e{y}_{*}}{{\left(1-x_{*}-d\right)}^2},{\beta }_{01}=b-n-\frac{e\left(x_{*}+d\right)}{1-x_{*}-d},{\beta }_{20}=\frac{e{y}_{*}}{{\left(1-x_{*}-d\right)}^3},\\ {\beta }_{11}=\frac{e}{{\left(1-x_{*}-d\right)}^2},{\beta }_{30}=\frac{-e{y}_{*}}{{\left(1-x_{*}-d\right)}^4},{\beta }_{21}=\frac{-e}{{\left(1-x_{*}-d\right)}^3},\\ {\beta }_{02}={\beta }_{12}={\beta }_{03}=0,\end{array}$

then $P\left(x_d\mathrm{,}{y}_d\right)\mathrm{,}Q\left(x_d\mathrm{,}{y}_d\right)$ are power series in $\left(x_d\mathrm{,}{y}_d\right)$ with terms $x_d^i{y}_d^j$ satisfying $i+j\ge 4$.

Thus, the first Lyapunov coefficient is

$\begin{array}{c}l=\frac{-3\pi }{2{\alpha }_{01}{\Delta }^{\frac{3}{2}}}\{[{\alpha }_{10}{\beta }_{10}\left({\alpha }_{11}^2+{\alpha }_{11}{\beta }_{02}+{\alpha }_{02}{\beta }_{11}\right)+{\alpha }_{10}{\alpha }_{01}\left({\beta }_{11}^2+{\alpha }_{20}{\beta }_{11}+{\alpha }_{11}{\beta }_{02}\right)\\ +{\beta }_{10}^2\left({\alpha }_{11}{\alpha }_{02}+2{\alpha }_{02}{\beta }_{02}\right)-2{\alpha }_{10}{\beta }_{10}\left({\beta }_{02}^2-{\alpha }_{20}{\alpha }_{02}\right)-2{\alpha }_{10}{\alpha }_{01}\left({\alpha }_{20}^2-{\beta }_{20}{\beta }_{02}\right)\\ -{\alpha }_{01}^2\left(2{\alpha }_{20}{\beta }_{20}+{\beta }_{11}{\beta }_{20}\right)+\left({\alpha }_{01}{\beta }_{10}-2{\alpha }_{10}^2\right)\left({\beta }_{11}{\beta }_{02}-{\alpha }_{11}{\alpha }_{20}\right)]\\ -\left({\alpha }_{10}^2+{\alpha }_{01}{\beta }_{10}\right)[3\left({\beta }_{10}{\beta }_{03}-{\alpha }_{01}{\alpha }_{30}\right)+2{\alpha }_{10}\left({\alpha }_{21}+{\beta }_{12}\right)+\left({\alpha }_{12}{\beta }_{10}-{\alpha }_{01}{\beta }_{21}\right)]\}\end{array}$

$\begin{array}{c}=\frac{-3\pi }{2{\alpha }_{01}{\Delta }^{\frac{3}{2}}}[{\alpha }_{10}{\beta }_{10}{\alpha }_{11}^2+{\alpha }_{10}{\alpha }_{01}\left({\beta }_{11}^2+{\alpha }_{20}{\beta }_{11}\right)-2{\alpha }_{10}{\alpha }_{01}{\alpha }_{20}^2\\ -{\alpha }_{01}^2\left(2{\alpha }_{20}{\beta }_{20}+{\beta }_{11}{\beta }_{20}\right)-{\alpha }_{11}{\alpha }_{20}\left({\alpha }_{01}{\beta }_{10}-2{\alpha }_{10}^2\right)\\ -\left({\alpha }_{10}^2+{\alpha }_{01}{\beta }_{10}\right)\left(-3{\alpha }_{01}{\alpha }_{30}+2{\alpha }_{10}{\alpha }_{21}-{\alpha }_{01}{\beta }_{21}\right)].\end{array}$

If $l<0$, the limit cycle is stable; if $l>0$, the limit cycle is unstable. However, the expression for Lyapunov number $l$ is rather cumbersome, we cannot directly judge the sign of it, so we will give some numerical simulation results in section 5.

Based on the mathematical theory, the existence and stability threshold conditions of all possible equilibrium points of the model (2.2) are deduced, and some critical conditions for inducing transcritical bifurcation and Hopf bifurcation of the model (2.2) are explored, which can provide a theoretical basis for some numerical simulation works. Furthermore, it is also worth pointing out that the key parameter d has a serious effect on bifurcation dynamics of the model (2.2).

5. Simulation Analysis and Results

In order to verify the validity of theoretical results, find some key control parameters that can induce bifurcation dynamics of the model (2.2), and explore ecological interaction between algae and protozoa, some numerical simulations are given with parameter values $n=0.4$, $b=0.2$, $e=0.6$, $m=0.1$, $q=0.6$ and $p=2$. From the equations 2.3, we can obtain that the dynamic relationship between algae x and protozoa y is $y=\frac{\left[x^2\left(p-x\right)-mpx\left(x+q\right)\right]\left(1+x-d\right)}{\left(x+q\right)\left(x-d\right)p}$. It is easy to find from Figure 1(b) that only when the algae density is greater than the value of d, the protozoa density can be positive, so the initial value of algae density in numerical simulation is larger than the value of d. At the same time, the density of protozoa y can reach a limit value and a maximum value within the range of algae x density, which implies that there may be an oscillatory coexistence mode between algae x and protozoa y. Furthermore, it is obvious to know from Figure 1(a) that the dynamic relationship between algae x and protozoa y is affected by the value of parameter d, thus, we can select parameter d as a control parameter of dynamic evolution process of the model (2.2).

The bifurcation dynamic evolution processes of the model (2.2) are shown in Figures 2-5. It is clearly visible from Figure 2 that if the value of d is greater than a critical value $d_{TC}=1.23$, the boundary equilibrium point $E_1$ is locally asymptotically stable, the model (2.2) has no internal equilibrium point. However, if the value of d is less than a critical value $d_{TC}=1.23$, the boundary equilibrium point $E_1$ loses stability and a new internal equilibrium point appears, this process implies that a transcritical bifurcation occurs, the detailed dynamic results are shown in Fi.3. Therefore, if the value of d is within the interval $\left(d_{Hp}\mathrm{,}{d}_{TC}\right)$, the model (2.2) has a stable internal equilibrium point, that is, algae and protozoa have a steady-state coexistence mode. As the value of d gradually decreases and is lower than a key value $d_{Hp}=0.1882$, the internal equilibrium point loses stability and a limit cycle appears, which implies that the model (2.2) has a Hopf bifurcation dynamic behavior, the dynamic evolution process of Hopf bifurcation is shown in Figure 4 and Figure 5, which shows that algae and protozoa coexist in a periodic oscillation mode. At the same time, because the first Lyapunov coefficient is $-0.046126845270109456077\Pi$, this limit cycle is stable. Furthermore, it is also worth emphasizing that when the value of parameter d gradually decreases from 1.5 to 0, the model (2.2) will undergo transcritical bifurcation and Hopf bifurcation successively, which means that the coexistence mode of algae and protozoa has changed fundamentally, from a protozoan extinction mode to a steady-state coexistence mode, and finally to a stable periodic oscillatory coexistence mode. Therefore, it is worth pointing out that the value of parameter d seriously affects the coexistence of algae and protozoa.

Figure 1. (a) Dynamic relationship between algae x, protozoa y and parameter d value; (b) Dynamic relationship between algae x and protozoa y with $d=0.1$.

In order to investigate the influence mechanism of Allee effect on dynamic behavior of the model (2.2), we will select parameter q as a control parameter for relevant dynamic simulation experiments. It is relatively clear from Figure 6 and Figure 7 that the model (2.2) has a constant steady state and a stable periodic

Figure 2. Bifurcation diagram of the model (2.2), here the red line indicates that the internal equilibrium point $x_{\ast }$ changes with the parameter d, the blue line and yellow line stand for the boundary equilibrium point $x_0$ and $x_1$ with the parameter d value changing, respectively. Here $x_0$ and $x_1$ represent the boundary equilibrium point $E_0\left(\mathrm{0,0}\right)$ and $E_1\left(x_1\mathrm{,0}\right)$. The solid curve shows that the equilibrium point is stable, the cyan dotted curve shows that the equilibrium point is unstable, and the vertical dot dotted line indicates a critical value of the equilibrium point which can induce bifurcation. More detailed, H_P and TC are some critical values for the Hopf bifurcation and transcritical bifurcation. Besides, the red dots are solid points, where the boundary equilibrium point $\left(\mathrm{0,0}\right)$ does not exist.

Figure 3. (a) If $d=d_{TC1}=1.23066238$, then a transcritical bifurcation occurs, where the boundary equilibrium point $E_1\left(x_1\mathrm{,0}\right)$ and the equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ coincide; (b) $E_1\left(x_1\mathrm{,0}\right)$ is a saddle when $d=1.2<d_{TC1}$, which can separate an internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ ; (c) $E_1\left(x_1\mathrm{,0}\right)$ is a stable node when $d_{TC1}<d=1.25$, and the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ does not exist.

Figure 4. (a) If $d=d_{Hp}=0.18821122$, then a Hopf bifurcation occurs at the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ and can generate a periodic solution; (b) Local magnification plot of (a); (c) $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ is an unstable node (or spiral source)if $d=0.16<d_{Hp}$ ; (d) $E_{\mathrm{<mo>\; *\; </mo>}}\left(x_{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}_{\mathrm{<mo>\; *\; </mo>}}\right)$ is a stable node (or spiral source) if $d_{Hp}<d=0.2$.

Figure 5. (a) Dynamic evolution diagram of Hopf bifurcation based on the change of parameter d value; (b) Local magnification plot of (a) when $d>d_{Hp}=0.1882$.

Figure 6. (a) Time series of algae x with $q=1.5$ ; (b) Time series of protozoa y with $q=1.5$ ; (c) Phase diagram of algae x and protozoa y with $q=1.5$.

Figure 7. (a) Time series of algae x with $q=1.85$ ; (b) Time series of protozoa y with $q=1.85$ ; (c) Phase diagram of algae x and protozoa y with $q=1.85$.

oscillation state, when the value of parameterq are 1.5 and 1.85 respectively, that is to say, the model (2.2) experiences a Hopf bifurcation dynamic behavior with the increase of parameter q value. Furthermore, this simulation result also indirectly shows that the size of Allee effect seriously affects the coexistence mode of algae and protozoa.

Based on the numerical simulation analysis, we first know that the model (2.2) has complex bifurcation dynamic behavior, mainly including transcritical bifurcation and Hopf bifurcation. Secondly, the algal population density gradually decreases with the transition from transcritical bifurcations to Hopf bifurcations. Finally, the value of key parameters of Allee effect seriously affects the coexistence mode of algae and protozoa.

6. Conclusions and Remarks

In this paper, based on the dynamic relationship between algae and protozoa, an aquatic ecological model with Allee effect was established to explore bifurcation behavior and investigate how Allee effect affects the coexistence mode of algae and protozoa. Some key conditions were given to ensure the existence and stability of all possible equilibrium points, and induce the model (2.2) to have transcritical bifurcation and Hopf bifurcation, which were theoretical basis for subsequent numerical simulation and the necessary conditions for parameter estimation value.

Through numerical simulation, we can see that the model (2.2) has complex bifurcation dynamics. It can be seen from Figure 3 and Figure 4 that transcritical bifurcation could make algae and protozoa to transform from a protozoa gradual extinction coexistence mode to a steady-state coexistence mode, and Hopf bifurcation could force algae and protozoa to transform from a constant steady-state coexistence mode to a stable periodic oscillation coexistence mode, which implied that the ecological relationship between algae and protozoa had changed substantially. Furthermore, it was easy to know from Figure 6 and Figure 7 that the ecological relationship between algae and protozoa could change from a constant steady state to a periodic oscillatory steady state with the increase of key parameters of Allee effect, which showed that Allee effect seriously affected coexistence mode of algae and protozoa.

Based on the theoretical analysis and numerical simulation results, it is worth pointing out that the algal aggregation behavior can change the coexistence mode of algae and protozoa, and the greater the algae population aggregation intensity, the more adverse to the permanent survival of protozoa, this research result is consistent with the fact that algae bloom is not conducive to the survival of protozoa. Furthermore, it should also be emphasized that algae population has Allee effect mechanism with small key value, which is conducive to form periodic oscillation coexistence mode of algae and protozoa.

Although some theoretical and numerical simulation results have been obtained in this study, there are still some deficiencies that need our follow-up research, such as: 1) the natural growth mode of protozoa is too simple, and the logistic growth function needs to be studied subsequently; 2) the influence of hydrodynamics on algae and protozoa should be considered in modeling dynamic process. However, it is hoped that the research results of this paper can play a theoretical supporting role in the study of aquatic ecological model.

Acknowledgements

This work was supported by scientific and technological innovation activity plan for college students in Zhejiang Province (new talent plan) (Grant No: 2021R429015) and the National Natural Science Foundation of China (Grant No: 31570364).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Yang, L.Y. and Xiao, L. (2011) Outbreak, Harm and Control of Lake Cyanobacteria Bloom. Science Press, Beijing.
[2]	Ryding, S.O. and Rast, W. (1989) Control of Eutrophication of Lakes and Reservoirs. Parthenon Publishing Group, Carnforth.
[3]	Zhao, Y.J. and Shi, Z.L. (1999) Blue-Green Algal Viruses(Cyanophages). Viro-Logica Sinica, 14, 100-104.
[4]	Zhao, Y.J. and Liu, Y.D. (1996) Possible Microbial Control on the Adverse Impacts of Algae-Current Information about the Relationship between Algae and Microbes. Acta Hydrobiologica Sinica, 20, 173-177.
[5]	Dryden, R.C. and Wright, S.J. (1987) Predation of Cyanobacteria by Protozoa. Canadian Journal of Microbiology, 33, 471-482. https://doi.org/10.1139/m87-080
[6]	Sigee, D.C. (1999) Biological Control of Cyanabacteria: Principles and Possibilities. Hydrobiologia, 395, 161-172. https://doi.org/10.1023/A:1017097502124
[7]	Liu, X.Y., Shi, M., Liao, Y.H., et al. (2005) Protozoa Capable Grazing on Cyanobacteria and Its Biological Control of the Algae Blooming. Acta Hydrobiologica Sinica, 29, 456-461.
[8]	Yamamoto, Y. (1981) Observation on the Occurrence of Microbial Agents Which Cause Lysis of Blue-Gree Algae in Lake Kasumigaura. Japanese Journal of Limnology, 42, 20-27. https://doi.org/10.3739/rikusui.42.20
[9]	Canter, H.M. Heaney, S.I. and Lund, J.W. (1990) The Ecological Significance of Grazing on Planktonic Populations of Cyanobacteria by the Ciliate Nassula. New Phytologist, 144, 247-263. https://doi.org/10.1111/j.1469-8137.1990.tb00397.x
[10]	Liu, G.Q. and Zhang, Z. (2016) Controlling the Nuisance Algae by Silver and Bighead Carps in Eutrophic Lakes: Disputes and Consensus. Journal of Lake Sciences, 28, 463-475. (In Chinese) https://doi.org/10.18307/2016.0301
[11]	Wang, S., Wang, Q.S., Zhang, L.B., et al. (2009) Large Enclosures Experimental Study on Algal Control by Silver Carp and Bighead. China Environmental Science, 29, 1190-1195.
[12]	Li, X.X., Yu, H.G., Dai, C.J., et al. (2021) Bifurcation Analysis of a New Aquatic Ecological Model with Aggregation Effect. Mathematics and Computers in Simulation, 190, 75-96. https://doi.org/10.1016/j.matcom.2021.05.015
[13]	Liu, H.Y., Yu, H.G., Dai, C.J., et al. (2021) Dynamical Analysis of an Aquatic Amensalism Model with Non-Selective Harvesting and Allee Effect. Mathematical Bioscience and Engineering, 18, 8857-8882. https://doi.org/10.3934/mbe.2021437
[14]	Franck, C., Ludek, B. and Joanna, G. (2008) Allee Effects in Ecology and Conservation. Oxford University Press, Oxford.
[15]	Taylor, C.N. and Hastings. A. (2005) Allee Effects in Biological Invasions. Ecology Letters, 8, 895-908. https://doi.org/10.1111/j.1461-0248.2005.00787.x
[16]	Cao, Q. and Li, Y.L. (2021) Differdivergent Solution of a Class of Prey-Prey Models of Baits with Allee Effect. Journal of Engineering Mathematics, 38, 377-388.
[17]	Surendar, M.S. and Sambath. M. (2021) Qualitative Analysis for a Phytoplankton-Zooplankton Model with Allee Effect and Holling Type II Response. Discontinuity, Nonlinearity and Complexity, 10, 1-18. https://doi.org/10.5890/DNC.2021.03.001
[18]	Arancibia-Ibarra, C. and Flores, J. (2021) Dynamics of a Leslie-Gower Predator-Prey Model with Holling Type II Functional Response, Allee Effect and a Generalist Predator. Mathematics and Computers in Simulation, 188, 1-22. https://doi.org/10.1016/j.matcom.2021.03.035
[19]	Tripathi, J.P., Mandal, P.S., Poonia, A. and Bajiya, V.P. (2021) A Widespread Interaction between Generalist and Specialist Enemies: The Role of Intraguild Predation and Allee Effect. Applied Mathematical Modelling, 89, 105-135. https://doi.org/10.1016/j.apm.2020.06.074
[20]	Yang, Q.Q. (2021) Analysis of a Class of Predation Models with Strong and Weak Allee Effects. Northwestern University, Evanston.
[21]	Yu, H.G., Zhao, M. and Agarwal, R.P. (2014) Stability and Dynamics Analysis of Time Delayed Eutrophication Ecological Model Based upon the Zeya Reservoir. Mathematics and Computers in Simulation, 97, 53-67. https://doi.org/10.1016/j.matcom.2013.06.008
[22]	Bertalanffy, L.V. (1932) Theoretische Biologise. Springer, Berlin.
[23]	Yu, H.G., Zhao, M., Wang, Q. and Agarwal, R.P. (2014) A Focus on Long-Run Sustainability of an Impulsive Switched Eutrophication Controlling System Based upon the Zeya Reservoir. Journal of the Franklin Institute, 351, 487-499. https://doi.org/10.1016/j.jfranklin.2013.08.025
[24]	Malchow, H., Petrovskii, S. and Medvinsky, A. (2011) Pattern Formation in Models of Plankton Dynamics: A Synthesis. Oceanologica Acta, 24, 479-487. https://doi.org/10.1016/S0399-1784(01)01161-6