Dynamics Analysis of an Aquatic Ecological Model with Temperature Effect

Heng Yang¹, Qi Wang², Hengguo Yu^1,2*

¹School of Mathematics and Physics, Wenzhou University, Wenzhou, China.
²Key Laboratory for Subtropical Oceans & Lakes Environment and Biological Resources Utilization Technology of Zhejiang, Wenzhou University, Wenzhou, China.

DOI: 10.4236/jamp.2022.1010199                

Abstract

Under the control framework of algae bloom in eutrophic lakes and reservoirs based on biological manipulation, the temperature variable is introduced into ecological modeling to show that it is a necessary condition for the rapid occurrence of algal blooms, and an aquatic ecological model with temperature effect is proposed to describe dynamic relationship between algae and biological manipulation predator. The mathematical theory work mainly investigates the existence and stability of some equilibrium points and some critical conditions for the occurrence of transcritical bifurcation and Hopf bifurcation. The numerical simulation mainly shows the dynamic evolution process of bifurcation dynamics, which can not only verify the validity and feasibility of these theoretical works but also analyze the influence of some key parameters on dynamic behavior evolution. Furthermore, It is worth emphasizing that temperature plays an important role in the coexistence of algae and biological manipulation predators. Moreover, the coexistence mode of algae and biological manipulation predators is discovered by means of dynamic bifurcation evolution. Finally, it is hoped that these research results can provide some reference for the study of aquatic ecosystems.

Keywords

Algae, Biological Manipulation Predator, Temperature, Bifurcation, Coexistence Mode

1. Introduction

As we all know, the temperature is an important environmental factor that determines the growth of the algal population, is one of the important factors that affect the growth of algal cells, the composition and content of biological macromolecules in cells, and is also a key ecological factor that affects the growth, reproduction and population succession of aquatic plants [1] [2]. Furthermore, suitable temperature is a necessary condition for algal bloom outbreaks and also an important environmental factor for the replacement of the dominant algal population [3] [4]. Thus, it is important to study the effect of temperature on algal growth and provide the theoretical basis for the prevention and treatment of water eutrophication.

Algae is a ubiquitous photosynthetic organism; the algal growth rate can increase with increasing temperature up to a certain limit; this is because that temperature strongly influences the cellular chemical composition and uptake of nutrients and also plays a significant role in algal growth [5]. The paper [6] pointed out that the ability to model algal productivity under transient conditions of temperature was critical for assessing the profitability and sustainability of full-scale algae cultivation outdoors. The paper [7] investigated the effect of temperature and irradiance on the growth and reproduction of the green macroalga and gave that a suitable temperature range over 21˚C - 29˚C was more favorable for growth and reproduction. The paper [8] proposed that Microcystis aeruginosa had a wide range of adaptation to temperature, and the optimal growth temperature was 25˚C - 30˚C. The paper [9] showed that the temperature condition that was most conducive to algal growth was 25˚C, and the optimal condition for algal toxins release was 28˚C. The paper [10] inquired into the effects of temperature on the germination of micro-propagules via laboratory experiments and indicated that sea temperature played a significant role in the germination of green algae. The paper [11] sought to elucidate the effects of temperature on algal growth rates, biomass accumulation, fatty acid production and composition and pointed out that temperature significantly impacted the overall productivity of algal biofuel systems by influencing species growth rates and fatty acid production. The paper [12] studied the effects of different temperatures and illumination time on algal growth and obtained that the value of the algal growth rate constant was reduced to 0.812d⁻¹ by lowering the water temperature to 16˚C. In conclusion, the algal population has strong adaptability to the environment, and temperature is one of the important ecological factors affecting algal growth.

In an aquatic ecosystem, the ecological effects between algae and biological manipulation predators are mutual, mainly including the harm of algal blooms to fish and the feeding and regulation of algae by biological manipulation predators [13]. The use of biological manipulation predators to control the excessive algal growth in eutrophic lakes was proposed and gradually attracted attention after the eutrophication of lakes became more and more common, and certain results in the control of cyanobacteria blooms have been achieved in some lakes at home and abroad [14] [15]. The paper [16] pointed out that silver carp completely eliminated cyanobacteria Microcystis by size and biovolume reduction. The paper [17] deemed that silver and bighead carps were just suitable for controlling cyanobacteria bloom not total algae biomass, and the application of the fish for algal control in water supply reservoirs should be cautious. The paper [18] believed that it was feasible to use silver carp and bighead carp to control Microcystis in eutrophic water. Furthermore, using ecological models to explore the dynamic relationship between algae and biological manipulation predators has also developed rapidly. The paper [19] showed that Microcystis aeruginosa aggregation could effectively control the dynamic feeding behavior of filter-feeding fish and provide shelter from predators. The paper [20] obtained that the filter-feeding fish population could be a crucial factor in controlling the proliferation of the algae population based on an algae-fish model. Therefore, fully understanding the ecological function relationship between algae and biological manipulation predator is a premise for better implementing algae control strategy with biological manipulation predators and ensuring algae control effect.

In this paper, firstly, the temperature variable is introduced to build an aquatic ecological model which can describe the dynamic relationship between algae and biological manipulation predators. Secondly, mathematical theory works are implemented to obtain some critical threshold conditions by investigating the evolution process of some specific dynamic properties. Finally, the numerical simulation works are carried out to show the evolution process of dynamic properties and the influence mechanism of temperature variables. Generally speaking, the main purpose of this paper is to use an aquatic ecological model with temperature effect to explore the coexistence mode between algae and biological manipulation predators and reveal how temperature variable affect their dynamic relationships.

2. Aquatic Ecological Model

At present, the temperature has been regarded as a key factor for monitoring and predicting algal blooms. Furthermore, the reproduction rate of algae has an important relationship with temperature, and the appropriate temperature is a necessary condition for the rapid propagation and growth of the algae population. In order to better investigate the effect of temperature on the dynamic relationship between algae and biological manipulation predator, we will propose an aquatic ecological model with temperature effect, which can be described as follow:

$(\begin{array}{l}\frac{\text{d}N}{\text{d}T}=r_1\left(\frac{R-R_{\min }}{R_{\text{ref}}-R_{\min }}\right)N\left(1-\frac{N}{k_1}\right)-\left(u_1+s+m_1\right)N-\frac{aNP}{b+N},\\ \frac{\text{d}P}{\text{d}T}=r_{2}P+\frac{e\left(1-\alpha \right)aNP}{b+N}-m_{2}P,\end{array}$ (2.1)

where $N\left(T\right)$ and $P\left(T\right)$ are density of algae population and biological manipulation predator (silver carp, bighead carp and anodonta woodiana elliptica) respectively. $r_1$ is maximum growth rate of algae population at reference temperature $R_{\text{ref}}$ (it is considered to be 25˚C). $k_1$ is the maximum environmental capacity for algae population. R represents actual temperature during the test run, which can be either a constant or a function of time t. $R_{\min }$ is the lowest temperature value when the algal growth is 0, and it is generally considered to be 15˚C. $u_1$ is a respiratory rate, $m_1$ is a nonpredation-induced mortality of algae population, s is a sedimentation rate, $r_2$ is intrinsic growth rate of biological manipulation predator, $m_2$ is mortality rate of biological manipulation predator, e is energy conversion rate, a is capture rate of biological manipulation predator, b is a semi saturation coefficient, and $\alpha$ is assimilated food of catabolic loss during predation period.

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

$\begin{array}{l}x=\frac{N}{k_1},c=u_1+s+m_1,g=\frac{R-R_{\min }}{R_{\text{ref}}-R_{\min }},t=r_{1}gT,y=P,\\ m=\frac{c}{r_{1}g},n=\frac{a}{r_{1}gb},p=\frac{k_1}{b},u=\frac{r_2}{r_{1}g},\beta =\frac{m_2}{r_{1}g},v=\frac{e\left(1-\alpha \right)a{k}_1}{r_{1}gb}.\end{array}$

Then the model (2.1) can be rewritten as

$(\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(1-x\right)-mx-\frac{nxy}{1+px},\\ \frac{\text{d}y}{\text{d}t}=uy+\frac{vxy}{1+px}-\beta y.\end{array}$ (2.2)

For model (2.2), we first discuss the existence and stability of all possible equilibrium points and explore the existence of a limit cycle with some key constraints. Then we give some critical conditions to demonstrate the occurrence of transcritical bifurcation and Hopf bifurcation. Finally, some numerical simulation results are implemented to verify the validity and feasibility of theoretical results. At the same time, through some numerical simulation results, the influence of temperature on dynamic relationship between algae and biological manipulation predators will be explored, and then ecological evolution significance represented by bifurcation dynamic evolution behavior is also given.

3. Equilibrium Points and Their Stability

In this section, some preliminary results shall be presented, including the existence and stability of all possible equilibrium points of the model (2.2).

We will consider the following equation to explore all possible equilibrium points of the model (2.2),

$(\begin{array}{l}F\left(x,y\right)=x\left(1-x\right)-mx-\frac{nxy}{1+px}=0,\\ G\left(x,y\right)=uy+\frac{vxy}{1+px}-\beta y=0.\end{array}$ (3.1)

It is easy to acquire that the model (2.2) always has trivial equilibrium point $E_0\left(\mathrm{0,0}\right)$, and a biological manipulation predator extinction equilibrium point $E_1\left(1-m\mathrm{,0}\right)$ if $1>m$. In view of the biological significance and the characteristics of the model (2.2), the interior equilibrium points are conditional.

In order to explore the existence and stability of all possible equilibrium points of the model (2.2), we will gradually give the following Theorems 1 - 6.

Theorem 1. The model (2.2) has a unique interior equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ if and only if the following conditions hold,

$0<x_{\ast }<1-m,x_{\ast }=\frac{\beta -u}{v-p\beta +pu}.$

Proof: From the previous analysis, we know that $G\left(x,y\right)=0$ has two distinct real roots, that is $y_1=0$ or $u+\frac{vx}{1+px}-\beta =0$, thus there are $y_1=0$ and $x_{\ast }=\frac{\beta -u}{v-p\beta +pu}$. If we take $y_1=0$ into $F\left(x,y\right)=0$, there are $x_1=0$ or $x_2=1-m$. Similarly, we take $x_{\ast }=\frac{\beta -u}{v-p\beta +pu}$ into $F\left(x,y\right)=0$, we can obtain $1-x-m-\frac{ny}{1+px}=0$ and $y_{\ast }=\frac{\left(1-x_{\ast }-m\right)\left(1+p{x}_{\ast }\right)}{n}$. According to the biological significance of the internal equilibrium point, $\frac{\beta -u}{v-p\beta +pu}>0$ and $\frac{\left(1-x_{\ast }-m\right)\left(1+p{x}_{\ast }\right)}{n}>0$ must be established, hence the existence equivalence condition is $0<\frac{\beta -u}{v-p\beta +pu}<1-m$. This ends the proof.

Theorem 2. The model (2.2) has a boundary equilibrium point $E_0\left(\mathrm{0,0}\right)$ and it is always unstable.

1) If $u>\beta$, $E_0$ is an unstable node; 2) if $u<\beta$, then $E_0$ is a saddle.

Proof: The Jacobian matrix of the model (2.2) evaluated at $E_0$ is

$J_{E_0}=\left(\begin{array}{cc}1-m& 0\\ 0& u-\beta \end{array}\right).$

The eigenvalues of $J_{E_0}$ are ${\lambda }_1=1-m$, ${\lambda }_2=u-\beta$. Based on the biological significance of the model (2.2), the value m must be less than 1. so $\lambda >0$. Thereby, if $u>\beta$, then ${\lambda }_2>0$, and $E_0$ is an unstable node; if $u<\beta$, then ${\lambda }_2<0$, and $E_0$ is a saddle. This completes the proof.

Theorem 3. The model (2.2) always has a boundary equilibrium point $E_1\left(1-m\mathrm{,0}\right)$

1) if $p>\frac{v}{\beta -u}-\frac{1}{1-m}$, $E_1$ is a stable node; 2) if $p<\frac{v}{\beta -u}-\frac{1}{1-m}$, then $E_1$ is a saddle.

Proof: The Jacobian matrix of the model (2.2) evaluated at $E_1$ is

$J_{E_1}=\left(\begin{array}{cc}m-1& -\frac{n-nm}{1+p-pm}\\ 0& u-\beta +\frac{v-vm}{1+p-pm}\end{array}\right)\mathrm{.}$

The eigenvalues of $J_{E_1}$ are ${\lambda }_1=m-1$, ${\lambda }_2=u-\beta +\frac{v-vm}{1+p-pm}$. Since ${\lambda }_1<0$, so when $p<\frac{v}{\beta -u}-\frac{1}{1-m}$, there is ${\lambda }_2>0$, then $E_1$ is a saddle. When $p>\frac{v}{\beta -u}-\frac{1}{1-m}$, there is ${\lambda }_2<0$, then $E_1$ is a stable node. This completes the proof.

Theorem 4. Under the condition $F\left(x_{\ast }\right)<0$, the boundary equilibrium point $E_1$ is globally asymptotically stable, where $x_{\ast }=\frac{\beta -u}{up+v-\beta p}$.

Proof: This theorem can be derived by flow analysis. When $F\left(x_{\ast }\right)<0$, the model (2.2) has no internal equilibrium point, and has only two boundary equilibrium points $E_0$ and $E_1$. It is easy to know that $E_0$ is always unstable and is a saddle, its unstable manifold is x-axis. Thus, we can divide the positive quadrant into the following three regions,

$R_1=\left\{\left(x,y\right)\in R_{+}^2|0<x\le m,y>0\right\},$

$R_2=\left\{\left(x,y\right)\in R_{+}^2|m<x<1-m,0<y<G\left(x\right)\right\},$

$R_3=\left\{\left(x,y\right)\in R_{+}^2|x>m,y>G\left(x\right)\right\}.$

Here $y<G\left(x\right)$ is derived from the first equation of the Equations (3.1). Considering the biological significance of the model (2.2), we need only discuss the region $R_2$ and $R_3$. In fact, solutions that start in region $R_3$ will eventually enter region $R_2$ by crossing the algae population nullcline vertically (downwards) in finite time, for the reason that there is no internal equilibrium point in $R_3$ and $\frac{\text{d}x}{\text{d}t}<0$ and $\frac{\text{d}y}{\text{d}t}<0$. Once in region $R_2$, it is clear to see that solutions are trapped, because $\frac{\text{d}x}{\text{d}t}>0$ and $\frac{\text{d}y}{\text{d}t}<0$, so the solutions in $R_2$ cannot cross the algae population nullcline, hence we must have $\left(x\left(t\right)\mathrm{,}y\left(t\right)\right)\to \left(1-m\mathrm{,0}\right)$ as $t\to +\infty$. This ends the proof.

Theorem 5. Under the condition that the theorem 1 exists, the equilibrium point $E_{\ast }$ is locally asymptotically stable if $Tr\left(J_{E_{\ast }}\right)<0$, and is unstable if $Tr\left(J_{E_{\ast }}\right)>0$.

Proof: The Jacobian matrix of the model (2.2) evaluated at $E_{\ast }$ is given by

$J_{E_{\ast }}=\left(\begin{array}{cc}1-2x_{\ast }-m-\frac{n{y}_{\ast }}{{\left(1+p{x}_{\ast }\right)}^2}& -\frac{n{x}_{\ast }}{1+p{x}_{\ast }}\\ \frac{v{y}_{\ast }}{{\left(1+p{x}_{\ast }\right)}^2}& 0\end{array}\right).$

The determinant and the trace of matrix $J_{E_{\ast }}$ are given by

$Det\left(J_{E_{\ast }}\right)=\left(\frac{n{x}_{\ast }}{1+p{x}_{\ast }}\right)\left(\frac{v{y}_{\ast }}{{\left(1+p{x}_{\ast }\right)}^2}\right),$

$Tr\left(J_{E_{\ast }}\right)=1-2x_{\ast }-m-\frac{n{y}_{\ast }}{{\left(1+p{x}_{\ast }\right)}^2}.$

It is easy to check that $Det\left(J_{E_{\ast }}\right)$ is always positive. Hence, if $Tr\left(J_{E_{\ast }}\right)<0$, then the equilibrium point $E_{\ast }$ is locally asymptotically stable; if $Tr\left(J_{E_{\ast }}\right)>0$, then the equilibrium point $E_{\ast }$ is unstable. If we can substitute the expression of $E_{\ast }$ into the trace of the Jacobian matrix, we cannot directly derive the sign of $Tr\left(J_{E_{\ast }}\right)$ because of algebraic complexity of the expression; therefore, we will calculate the values of $Tr\left(J_{E_{\ast }}\right)$ by using numerical simulation in Section 5.

Theorem 6. If $\frac{pn{y}_{\ast }}{1+p{x}_{\ast }}<1$ is satisfied, then the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ is globally asymptotically stable.

Proof: To prove the global asymptotic stability of the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$, the following Lyapunov function was constructed:

$V\left(x,y\right)=\left(x-x_{\ast }-x_{\ast }\ln \frac{x}{x_{\ast }}\right)+B\left(y-y_{\ast }-y_{\ast }\ln \frac{y}{y_{\ast }}\right).$

Based on the basic properties of the function, it can be concluded that $V\left(x,y\right)$ is continuous for all $x>0$ and $y>0$, and is computationally available:

$\frac{\partial V}{\partial x}=1-\frac{x_{\ast }}{x},\frac{\partial V}{\partial y}=B\left(1-\frac{y_{\ast }}{y}\right).$

Therefore, the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ is only extreme value of the $V\left(x\mathrm{,}y\right)$ function in the positive quadrant,

$\underset{x\to 0}{\lim }V\left(x,y\right)=\underset{y\to 0}{\lim }V\left(x,y\right)=\underset{x\to +\infty }{\lim }V\left(x,y\right)=\underset{y\to +\infty }{\lim }V\left(x,y\right)=+\infty .$

For all the $x>0$ and $y>0$, we can get

$V\left(x,y\right)>V\left(x_{\ast },y_{\ast }\right)=0.$

Then, we can obtain

$\begin{array}{c}\frac{\text{d}V}{\text{d}t}=\frac{\text{d}x}{\text{d}t}-\frac{x_{\ast }}{x}\frac{\text{d}x}{\text{d}t}+B\left(\frac{\text{d}y}{\text{d}t}-\frac{y_{\ast }}{y}\frac{\text{d}y}{\text{d}t}\right)\\ =\left(x-x_{\ast }\right)\left(1-x-m-\frac{ny}{1+px}\right)+B\left(y-y_{\ast }\right)\left(u-\beta +\frac{vx}{1+px}\right).\end{array}$

Furthermore, we know that the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ satisfies the following equation:

$\begin{array}{l}1-x_{\ast }-m-\frac{n{y}_{\ast }}{1+p{x}_{\ast }}=0,\\ u-\beta +\frac{v{x}_{\ast }}{1+p{x}_{\ast }}=0.\end{array}$

Combined with the above two formulas for the simplification, we can obtain

$\begin{array}{c}\frac{\text{d}V}{\text{d}t}=\left(x-x_{\ast }\right)\left(x_{\ast }+\frac{n{y}_{\ast }}{1+p{x}_{\ast }}-x-\frac{ny}{1+px}\right)+B\left(y-y_{\ast }\right)\left(\frac{vx}{1+px}-\frac{v{x}_{\ast }}{1+p{x}_{\ast }}\right)\\ =\left(\frac{np{y}_{\ast }}{\left(1+p{x}_{\ast }\right)\left(1+px\right)}-1\right){\left(x-x_{\ast }\right)}^2-\frac{n\left(1+p{x}_{\ast }\right)\left(x-x_{\ast }\right)\left(y-y_{\ast }\right)}{\left(1+p{x}_{\ast }\right)\left(1+px\right)}+\frac{Bv\left(x-x_{\ast }\right)\left(y-y_{\ast }\right)}{\left(1+p{x}_{\ast }\right)\left(1+px\right)}\mathrm{.}\end{array}$

By choosing

$B=\frac{n\left(1+p{x}_{\ast }\right)}{v}\mathrm{,}$

if $\frac{pn{y}_{\ast }}{1+p{x}_{\ast }}<1$, we can obtain

$\frac{\text{d}V}{\text{d}t}=-{\left(x-x_{\ast }\right)}^2\left(1-\frac{pn{y}_{\ast }}{\left(1+px\right)\left(1+p{x}_{\ast }\right)}\right)<-{\left(x-x_{\ast }\right)}^2\left(1-\frac{pn{y}_{\ast }}{1+p{x}_{\ast }}\right)<0.$

According to the above formula, it is obvious that the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ makes $\frac{\text{d}V}{\text{d}t}=0$ and for all other $\left(x\mathrm{,}y\right)$ has $\frac{\text{d}V}{\text{d}t}<0$, so $V\left(x\mathrm{,}y\right)$ satisfies the global stability theorem of the Lyapunov function. Thus the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ is globally asymptotically stable when $E_{\ast }$ is feasible and the implicit condition $\frac{pn{y}_{\ast }}{1+p{x}_{\ast }}<1$ is satisfied. This ends the proof.

Theorem 7. When the following two conditions are true, there exists at least one limit cycle of the model (2.2).

1) $\frac{\left(1-x_{\ast }-m\right)\left(1+p{x}_{\ast }\right)}{n}>0$ ; 2) $\beta >u$.

Proof: From the conclusion of the Theorem 5, we know that the internal equilibrium point $E_{\ast }$ may be an unstable focus. Now we will prove Theorem 7 by constructing an invariant region $\Omega$, which consists of the following line $L_1\mathrm{,}{L}_2$, and $x,y$ axis,

$L_1:x=1-m,L_2:x+\frac{ny}{v}-M=0.$

Through the derivation calculation, we can get

$\frac{\text{d}{L}_1}{\text{d}t}={\left[x\left(1-x\right)-mx-\frac{nxy}{1+px}\right]|}_{x=1-m}=\frac{-n\left(1-m\right)y}{1+p\left(1-m\right)},$

$\begin{array}{c}\frac{\text{d}{L}_2}{\text{d}t}={\left[x\left(1-x\right)-mx+\frac{n\left(u-\beta \right)y}{v}\right]|}_{y=\frac{v\left(M-x\right)}{n}}\\ =-x^2+\left(1-m-u+\beta \right)x+\left(u-\beta \right)M.\end{array}$

It is easy to verify that $\frac{\text{d}{L}_1}{\text{d}t}<0$ for $0<x<1-m$ and $y>0$. Furthermore, if ${\left(1-m-u+\beta \right)}^2-4\left(\beta -u\right)M<0$ when M is a sufficiently large positive number, then $\frac{\text{d}{L}_2}{\text{d}t}<0$ for $0<x<1-m$, which implies that as long as $\beta >u$ holds, there must be a large enough positive number M to make ${\left(1-m-u+\beta \right)}^2-4\left(\beta -u\right)M<0$ hold. Thus, the model (2.2) has at least one limit cycle by Poincare-Bendixson Theorem [21] [22]. This ends the proof.

4. Bifurcation Analysis

It is well known that the evolution process of bifurcation dynamics has important biological significance in the process of population dynamics. Therefore, we will explore some bifurcation dynamics behaviors of the model (2.2) and give some threshold conditions for specific bifurcation dynamics of the model (2.2).

4.1. Transcritical Bifurcation

Here we will prove that the model (2.2) undergoes a transcritical bifurcation at $p=p_{TC}=\frac{v}{\beta -u}-\frac{1}{1-m}$. We recall from Theorem 3, the boundary equilibrium $E_1$ loses its stability at $p=p_{TC}$ and one of the eigenvalues of $J_{E_1}$ is zero; it implies that the equilibrium point $E_1$ becomes non-hyperbolic, so it is possible to undergo a transcritical bifurcation around $E_1$.

Theorem 8. The model (2.2) undergoes a transcritical bifurcation when $p=p_{TC}$, where $p_{TC}=\frac{v}{\beta -u}-\frac{1}{1-m}$.

Proof: We use Sotomayor’s theorem [23] [24] to prove that the model (2.2) undergoes a transcritical bifurcation with a bifurcation parameter p. If $p=p_{TC}$, then $Det\left(J_{E_1}\right)=0$, which means that the Jacobian matrix $J_{E_1}$ has a zero eigenvalue. Now, let V and W be the eigenvectors corresponding to the zero eigenvalue of $J_{E_1}$ and $J_{E_1}^{\text{T}}$ respectively. After a simple calculation they can be given by

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

Therefore,

$V=\left(\begin{array}{c}{v}_1\\ v_2\end{array}\right)=\left(\begin{array}{c}\frac{n-nm}{\left(1+p-pm\right)\left(m-1\right)}\\ 1\end{array}\right),$

$W=\left(\begin{array}{c}{w}_1\\ w_2\end{array}\right)=\left(\begin{array}{c}0\\ 1\end{array}\right).$

Due to

$F_p\left(E_1;p_{TC}\right)=\left(\begin{array}{c}{F}_{1p}\\ F_{2p}\end{array}\right)={\left(\begin{array}{c}\frac{ny{x}^2}{{\left(1+px\right)}^2}\\ -\frac{vy{x}^2}{{\left(1+px\right)}^2}\end{array}\right)}_{\left(E_1,p_{TC}\right)}=\left(\begin{array}{c}0\\ 0\end{array}\right),$

$D{F}_p\left(E_1;p_{TC}\right)V={\left(\begin{array}{l}{F}_{1px}{F}_{1py}\hfill \\ F_{2px}{F}_{2py}\hfill \end{array}\right)}_{\left(E_1;p_{TC}\right)}\left(\begin{array}{l}{v}_1\hfill \\ v_2\hfill \end{array}\right)=\left(\begin{array}{l}\frac{n{\left(1-m\right)}^2}{{\left(1+p\left(1-m\right)\right)}^2}\hfill \\ -\frac{v{\left(1-m\right)}^2}{{\left(1+p\left(1-m\right)\right)}^2}\hfill \end{array}\right),$

$\begin{array}{l}{D}^2{F}_p\left(E_1;p_{TC}\right)\left(V,V\right)\\ ={\left(\begin{array}{l}\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\hfill \\ \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\hfill \end{array}\right)}_{\left(E_1;p_{TC}\right)}=\left(\begin{array}{l}\frac{2p\left(1-m\right)n^2}{{\left(1+p-pm\right)}^3}\hfill \\ -\frac{2v}{{\left(1+p-pm\right)}^3}\hfill \end{array}\right).\end{array}$

Furthermore, we can obtain

$W^{\text{T}}{F}_p\left(E_1;p_{TC}\right)=\left(0,1\right)\left(\begin{array}{l}0\hfill \\ 0\hfill \end{array}\right)=0,$

$W^{\text{T}}\left[D{F}_P\left(E_1;p_{TC}\right)V\right]=\left(0,1\right)\left(\begin{array}{l}\frac{n{\left(1-m\right)}^2}{{\left(1+p\left(1-m\right)\right)}^2}\hfill \\ -\frac{v{\left(1-m\right)}^2}{{\left(1+p\left(1-m\right)\right)}^2}\hfill \end{array}\right)=-\frac{v{\left(1-m\right)}^2}{{\left(1+p\left(1-m\right)\right)}^2}\ne 0,$

$W^{\text{T}}\left[D^2{F}_p\left(E_1;p_{TC}\right)\left(V,V\right)\right]=\left(0,1\right)\left(\begin{array}{l}\frac{2p\left(1-m\right)n^2}{{\left(1+p-pm\right)}^3}\hfill \\ -\frac{2v}{{\left(1+p-pm\right)}^3}\hfill \end{array}\right)=-\frac{2v}{{\left(1+p-pm\right)}^3}\ne 0.$

Thus, based on Sotomayor’s theorem we can deduce that the model (2.2) undergoes a transcritical bifurcation as the parameter p passes through a critical threshold $p_{TC}$. This completes the proof.

4.2. Hopf Bifurcation

From the analysis of Theorem 5, we know that if $Tr\left(J_{E_{\ast }}\right)<0$, $E_{\ast }$ is locally asymptotically stable; if $Tr\left(J_{E_{\ast }}\right)>0$, then $E_{\ast }$ is unstable. Hence, we can easily conclude that the interior equilibrium point $E_{\ast }$ may lose its stability by Hopf bifurcation under sufficient conditions. Considering p as a Hopf bifurcation control parameter, the Hopf bifurcation threshold $p=p_H$ can be detected by using $Tr\left(J_{E_{\ast }}\right)=0$. When the value of the parameter p passes from one side of $p=p_H$ to the other side, the stability property of $E_{\ast }$ changes and the periodic orbits can be generated. Thus, we can yield the following theorem.

Theorem 9. Based on Theorem 1, the model (2.2) undergoes a Hopf bifurcation around $E_{\ast }$ when $p=p_H$.

Proof: In order to ensure the existence of Hopf bifurcation, we need to verify the transversality condition; then we have

$\frac{\text{d}}{\text{d}p}{\left[Tr\left(J_{E_{\ast }}\right)\right]}_{p=p_H}={\left[-2\frac{\text{d}{x}_{\ast }}{\text{d}p}+\frac{\left(1+p-pm\right)\frac{\text{d}{x}_{\ast }}{\text{d}p}+\left(1-x_{\ast }-m\right)x_{\ast }}{{\left(1+p{x}_{\ast }\right)}^2}\right]}_{p=p_H}\ne 0,$

and

$\frac{\text{d}{x}_{\ast }}{\text{d}p}=\frac{{\left(u-\beta \right)}^2}{{\left(v-p\beta +pu\right)}^2}\mathrm{.}$

This can guarantee the existence of Hopf bifurcation around $E_{\ast }$.

Next we will compute the first Lyapunov number to discuss the stability of limit cycle. Firstly, we will transform $E_{\ast }$ to the origin by the transformation $\left(x\mathrm{,}y\right)=\left(x-x_{\ast }\mathrm{,}y-y_{\ast }\right)$, then we can get

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

where

${\alpha }_{10}=1-2x_{\ast }-m-\frac{n{y}_{\ast }}{{\left(1+p{x}_{\ast }\right)}^2},{\alpha }_{01}=-\frac{n{x}_{\ast }}{1+p{x}_{\ast }},{\alpha }_{20}=-2+\frac{2pn{y}_{\ast }}{{\left(1+p{x}_{\ast }\right)}^3},$

${\alpha }_{11}=-\frac{n}{{\left(1+p{x}_{\ast }\right)}^2},{\alpha }_{30}=-\frac{6p^{2}n{y}_{\ast }}{{\left(1+p{x}_{\ast }\right)}^4},{\alpha }_{21}=\frac{2pn}{{\left(1+p{x}_{\ast }\right)}^3},$

${\alpha }_{02}={\alpha }_{12}={\alpha }_{03}=0,$

and

${\beta }_{10}=\frac{v{y}_{\ast }}{{\left(1+p{x}_{\ast }\right)}^2},{\beta }_{01}=u-\beta +\frac{v{x}_{\ast }}{1+p{x}_{\ast }},{\beta }_{11}=\frac{v}{{\left(1+p{x}_{\ast }\right)}^2},$

${\beta }_{20}=\frac{-2pv{y}_{\ast }}{{\left(1+p{x}_{\ast }\right)}^3},{\beta }_{21}=\frac{-2pv}{{\left(1+p{x}_{\ast }\right)}^3},{\beta }_{30}=\frac{6p^{2}v{y}_{\ast }}{{\left(1+p{x}_{\ast }\right)}^4},$

${\beta }_{02}={\beta }_{12}={\beta }_{03}=0,$

and $P\left(x_m\mathrm{,}{y}_m\right)\mathrm{,}Q\left(x_m\mathrm{,}{y}_m\right)$ are power series with terms $x^i{y}^j\left(i+j\ge 4\right)$.

Therefore, the first Lyapunov number $l$ is given by

$\begin{array}{c}l=\frac{-3\pi }{2{\alpha }_{01}{\Delta }^{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)\left[3\left({\beta }_{10}{\beta }_{03}-{\alpha }_{01}{\alpha }_{30}\right)+2{\alpha }_{10}\left({\alpha }_{21}+{\beta }_{12}\right)+\left({\beta }_{10}{\alpha }_{12}-{\alpha }_{01}{\beta }_{21}\right)\right]\},\end{array}$

where

$\Delta ={\alpha }_{10}{\beta }_{01}-{\alpha }_{01}{\beta }_{10}>0.$

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 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 are given, and the critical conditions for inducing specific bifurcation dynamics of the model (2.2) are analyzed, which can provide a theoretical basis for subsequent numerical simulation work. Furthermore, it should also be emphasized that the key parameter p can seriously affect dynamic evolution characteristics of the model (2.2).

5. Numerical Simulations and Results

Now, we will investigate dynamic properties of the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ and explore the effect of some parameters on dynamic relationship between $x_{\ast }$ and $y_{\ast }$. From Theorem 1, we can see that the relation expression of x and y is $y=\frac{\left(1-x-m\right)\left(1+px\right)}{n}$ and $x=\frac{\beta -u}{v-p\beta +pu}$, thus some numerical simulations are given in Figure 1 with $m=0.2$, $n=0.9$, $u=0.2$, $v=0.4$,

Figure 1. (a) Dynamic relationship between population density x, y and parameter p value; (b) Dynamic relationship between population density x and parameter p value.

$\beta =0.25$. It is easy to find from Figure 1(b) that the density of algae $x_{\ast }$ increases with the increase of parameter p-value, and the growth is slow in the early stage and extremely fast in the later stage, which implies that when the value of parameter p exceeds a certain critical threshold, the density of algae $x_{\ast }$ will become larger and larger. Furthermore, it is obvious to know from Figure 1(a) that when the value of parameter p exceeds a certain critical threshold and is determined, the density of biological manipulation predator $y_{\ast }$ is a concave function of the density of algae $x_{\ast }$, and $y_{\ast }$ can get a maximum value when the value of $x_{\ast }$ is $\frac{\left(1-m\right)p-1}{2p}$, which must satisfy $0<\frac{\left(1-m\right)p-1}{2p}<1-m$, this to say $p>\frac{1}{1-m}$. Therefore, it must be emphasized that the dynamic properties of the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ and the dynamic relationship between population $x_{\ast }$ and $y_{\ast }$ mainly depends on key parameters m and p.

In order to verify the validity and feasibility of Theorem 5 and 7, the stability of the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ and the existence of the limit cycle are numerically simulated with $p=5.12$ and $p=3.8$. It is easy to find from Figure 2(a) that the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ is stable, which means that algae and biological manipulation predators can form a constant steady-state coexistence mode. Furthermore, it is obvious to know from Figure 2(b) that model (2.2) has a limit cycle, which indicates that algae and biological manipulation

Figure 2. (a) The stability of the internal equilibrium point $E_{\ast }$ with $m=0.2$, $n=0.9$, $u=0.2$, $v=0.4$, $\beta =0.25$, $p=5.12$ ; (b) The existence of limit cycle with $m=0.2$, $n=0.9$, $u=0.2$, $v=0.4$, $\beta =0.25$, $p=3.8$.

predators can form a periodic oscillation coexistence mode. In a word, the value of parameter p has an important influence on coexistence mode of algae and biological manipulation predators.

In order to better understand how the value of parameter p affects the dynamic behavior evolution of the model (2.2), we give a bifurcation diagram of the model (2.2) in Figure 3. It can be seen clearly from Figure 3 and Figure 4 that if the value of parameter p is larger than $p_{TC}=6.75$, the model (2.2) has only two boundary equilibrium points $E_0\left(\mathrm{0,0}\right)$ and $E_1\left(1-m\mathrm{,0}\right)$, $E_1\left(1-m\mathrm{,0}\right)$ is stable and $E_0\left(\mathrm{0,0}\right)$ is unstable, which implies that biological manipulation predator will eventually approach extinction and algae will eventually approach the maximum biomass state, that is to say, biological manipulation predator and algae cannot form a final coexistence mode. If the value of parameter p is between $\left(p_{PH}\mathrm{,}{p}_{TC}\right)$, the model (2.2) goes through a transcritical bifurcation, which will induce the model (2.2) to have an internal equilibrium point, and this internal equilibrium point is asymptotically stable, detailed numerical simulation results are shown in Figure 4(b) and Figure 4(a). Therefore, it can be said that the transcritical bifurcation induces the formation of a constant steady-state coexistence mode between biological manipulation predator and algae. If the value of parameter p gradually decreases and is less than the critical threshold

Figure 3. Bifurcation diagram of the model (2.2) by taking p as bifurcation parameter with $m=0.2$, $n=0.9$, $u=0.2$, $v=0.4$, $\beta =0.25$, Bifurcation diagram with respect to parameter $p_H=4.5543$, $p_{TC}=6.75$, where the red circle is a hollow dot, that means the boundary equilibrium point $E_0\left(\mathrm{0,0}\right)$ is not exist, the pink and green lines implicit the critical values for the Hopf bifurcation and transcritical bifurcation respectively. In addition, the blue circle and blue asterisk separately denote the Hopf (PH) and transcritical (TC) bifurcation point.

value $p_{PH}=4.5543$, the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ loses stability and a stable limit cycle appears. In other words, the model (2.2) undergoes a Hopf bifurcation; the numerical dynamic evolution process is shown in Figure 5. Therefore, it is worth pointing out that the Hopf bifurcation can produce a periodic oscillation coexistence mode between biological manipulation predator and algae. Thus, the numerical simulation results not only prove the validity and feasibility of the theoretical derivation, but also directly show that the value of key parameter p seriously affects bifurcation dynamics evolution characteristics of the model (2.2).

Figure 4. The phase portraits of the model (2.2), $E_1$ has different dynamics with the value of parameter p varying. (a) the interior equilibrium point $E_{\ast }$ exists and it is a stable node, boundary equilibrium $E_1$ is a saddle, and $E_0$ is an unstable node when $p=6.65<p_{TC}$. (b) the interior equilibrium point $E_{\ast }$ coincides with boundary equilibrium $E_1$, which is a saddle-node and the parabolic sector is on the upper half plane, and $E_0$ is a saddle when $p=6.75=p_{TC}$. (c) $E_1$ is a stable node, $E_0$ is a saddle, and there is no interior equilibrium point when $p=6.85>p_{TC}$. The green curves represent the stable or unstable orbits, and the red points are some equilibrium points.

Figure 5. The phase portraits of the model (2.2), where the interior equilibrium point $E_{\ast }$ has different dynamics with the value of parameter p varying. (a) Stable periodic orbits bifurcate through Hopf bifurcation around $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ with $p=p_H=4.5543$. (b) Local amplification of (a) for $\left(x\mathrm{,}y\right)\in \left[\mathrm{0.26,0.32}\right]\times \left[\mathrm{1.28,1.35}\right]$. (c) $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ is a spiral source point with $p=4.5<p_H$. (d) $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$ is locally asymptotically stable with $p=4.6>p_H$.

To explore in detail how temperature affects the coexistence mode of biological manipulation predators and algae, and analyze the advantages of the model (2.1), some numerical simulation results are shown in Figures 6-8. As we all know, in the laboratory algal-predator culture test, in order to better maintain the growth of algal population, we usually conduct the culture test at a constant temperature degrees. However, in subtropical reservoirs, the temperature has a periodic change in behavior with time, which is not a constant. Therefore, we will conduct numerical simulation with the state of laboratory culture temperature $R=25$, $R=32$ and the state of field natural temperature $R=25+10\sin t$ with $r_1=0.6$, $r_2=0.2$, $k_1=5$, $u_1=0.05$, $s=0.05$, $m_1=0.1$, $m_2=0.25$, $\alpha =0.2$, $a=0.4$, $b=2.5$, $e=0.7$, $R_{\text{ref}}=25$ and $R_{\min }=15$. It is obvious to find from Figure 6 and Figure 7 that algae and biological manipulation predators can

Figure 6. Time series diagram of the model (2.1) with $R=25$.

Figure 7. Time series diagram of the model (2.1) with $R=25+10\sin t$.

Figure 8. Time series diagram of the model (2.1) with $R=32$.

form a constant steady-state coexistence mode under the constant temperature assumption, however, algae and biological manipulation predator can form a periodic oscillation coexistence mode under the periodic change temperature assumption, In addition, for algae and biological manipulation predator, the center of the periodic oscillation amplitude is approximately a constant steady state value, which shows that the temperature expression does not affect the growth average biomass of algae and biological manipulation predator, but it will seriously affect their growth dynamics. Moreover, it is easy to see from Figure 8 that algae and biological manipulation predators have a periodic oscillation coexistence mode as the value of temperature parameter R increases to 32, and the biomass of biological manipulation predators has been greatly increased. These results not only show that the model (2.1) has experienced a Hopf bifurcation as the value of parameter R increases, but also indicate that the increase in temperature is conducive to rapid growth of biological manipulation predators.

Based on the numerical simulation results, can clearly indicate that the results of theoretical derivation are effective and feasible. Furthermore, it should also be emphasized that temperature not only affects the dynamic evolution characteristics of the model (2.1), but also affects the biomass level of biological manipulation predator. Moreover, the model (2.2) has specific bifurcation dynamic behaviors (transcritical bifurcation and Hopf bifurcation) under the influence of the value of key parameter p; these two bifurcation dynamics behaviors lead to a constant steady-state coexistence mode and a periodic oscillation coexistence mode of algae and biological manipulation predator respectively.

6. Conclusions

Under the conceptual framework of biological control of cyanobacteria in eutrophic lakes and reservoirs, based on the fact that temperature is an extremely important factor in determining ecology, which has an important relationship with algae proliferation rate, an aquatic ecological model with temperature effect is proposed to explore the coexistence modes of algae and biological manipulation predator and investigate how temperature affects their dynamic evolution. Suppose temperature parameter R is a constant variable, which can approximately describe ecological culture system of algae and biological manipulation predator under laboratory conditions if temperature parameter R is a periodic function variable, which mainly represents the natural ecosystem of algae and biological manipulation predators in a naturally eutrophic lake.

Based on dynamic population theory, some threshold conditions are given to guarantee the existence and stability of all possible equilibrium points, and some critical conditions for the occurrence of transcritical bifurcation and Hopf bifurcation are also deduced. Furthermore, some key parameters affecting the dynamic evolution characteristics of the model (2.2) are found through theoretical derivation and numerical simulation. All in all, these results are the theoretical basis for subsequent numerical simulation work and abstractly display the influence of some parameters on the dynamic evolution of the model (2.2).

Through the numerical simulation test on dynamic behaviors of the model (2.1), the influence mechanism of temperature on the stable succession of aquatic ecosystem is discovered in Figure 6 and Figure 8, the coexistence mode of algae and biological manipulation predator can change from a constant steady-state mode to a periodic oscillation mode with the temperature increasing gradually, which also indirectly indicates that the appropriate temperature range is one of the key factors for algae and biological manipulation predator to form a stable coexistence mode, and the periodic oscillation coexistence mode is more favorable to control the growth rate of algae population by biological manipulation. Based on the bifurcation dynamics evolution analysis of the model (2.2), it is worth pointing out that transcritical bifurcation can induce the appearance of the internal equilibrium point $E_{\ast }\left(x_{\ast }\mathrm{,}{y}_{\ast }\right)$, which represents the coexistence of algae and biological manipulation predator in a periodic oscillation mode, and completely changes the dynamic coexistence nature of algae and biological manipulation predator. Furthermore, when the value of control parameter p decreases and falls below a critical threshold, the coexistence mode of algae and biological manipulation predator has changed fundamentally again, periodic oscillation coexistence mode will replace the constant steady-state coexistence mode through a Hopf bifurcation. These results directly show that the value of key parameter p plays an important role in the bifurcation dynamic behavior evolution of the model (2.2). In general, some theoretical and numerical results obtained in this paper can provide a certain theoretical basis for the formation of a healthy and stable aquatic ecosystem and also provide certain numerical support for the feasibility of biological manipulation technology.

In the follow-up research works, firstly, we will introduce Arrhenius exponential temperature function and partial normal distribution temperature function into ecological modeling and investigate the impact of different temperature function manifestations on the dynamic behavior of the model (2.1). Secondly, we will continue to deepen the environmental impact factors of such aquatic ecological models and then further explore the influence of various environmental factors on the dynamic relationship between algae and biological manipulation predators. Finally, we will further explore the dynamic pattern behavior of the model (2.1) with the help of these papers [25] [26] [27]. In a word, all these results are expected to be useful in studying the dynamic behavior of aquatic ecosystems.

Acknowledgements

This work was supported by the key International Cooperation Projects of China (Grant No: 2018YFE0103700), the National Natural Science Foundation of China (Grant No: 31570364 and 61871293), key R&D Program Projects in Zhejiang Province of China (Grant No: 2021C03166), by the Science and Technology Program of Wenzhou, China (Grant No.X20210097).

Conflicts of Interest

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

References

[1]	Jorgensen, S.E. and Bendoeicchio, G. (2001) Fundamentals of Ecological Modelling. 3rd Edition, Elsevier Science BV, The Netherlands.
[2]	Ma, Q., Deng, C.N. and Guo, F.F. (2018) Effects of Different Temperature on Growth of Chlorella and Microcystis meruginosa and Chlorophy II Fluorescence. Journal of Zhongzhou University, 35, 108-112.
[3]	Eppley, R.E. (1972) Temperature and Phytoplankton Growth in the Sea. Fishery Bulletin, 70, 1063-1085.
[4]	Robarts, R.D. and Zohary, T. (2010) Temperature Effects on Photosynthetic Capacity, Respiration, and Growth Rates of Bloom Forming Cyanobacteria. New Zealand Journal of Marine and Freshwater Research, 21, 391-399. https://doi.org/10.1080/00288330.1987.9516235
[5]	Singh, S.P. and Singh, P. (2015) Effect of Temperature and Light on the Growth of Algae Species: A Review. Renewable and Sustainable Energy Reviews, 50, 434-444. https://doi.org/10.1016/j.rser.2015.05.024
[6]	Bechet, Q., Shilton, A. and Guieysse, B. (2013) Modeling the Effects of Light and Temperature on Algae Growth: State of the Art and Critical Assessment for Productivity Predaiction Outdoor Cultivation. Biotechnology Advances, 31, 1648-1663. https://doi.org/10.1016/j.biotechadv.2013.08.014
[7]	Deng, Y.Y., Tang, X.R., Huang, B.X., et al. (2012) Effect of Temperature and Irradiance on the Growth and Reproduction of the Green Macroalga, Chaetomorpha valida (Cladophoraceae, Chlorophyta). Journal of Applied Phycology, 24, 927-933. https://doi.org/10.1007/s10811-011-9713-0
[8]	Chen, J.Z., Liu, Z.L., Li, X.M., et al. (2010) Effects of Temperature, pH, Nitrogen and Phosphorus on Growth of Microcytis aeruginosa. Oceanologia ET Limnologia Sinica, 41, 714-718.
[9]	Wu, R., Cui, L.F., Lyu, S., et al. (2010) The Influence of Temperature and Illumination on the Microcystis Production. Enviromental Science & Technology, 33, 33-36.
[10]	Song, W., Peng, K.Q., Xiao, J., et al. (2015) Efeects of Temperature on the Germination of Gree Algae Micropropagules in Coastal Waters of the Subei Shoal China. Estuarine, Coastal and Shelf Science, 163, 63-68. https://doi.org/10.1016/j.ecss.2014.08.007
[11]	Nalley, J.O., O’Donnell, D.R. and Litchman, E. (2018) Temperature Effects on Growth Rates and Fatty Acid Content in Freshwater Algae and Cyanobacteria. Algal Research, 35, 500-507. https://doi.org/10.1016/j.algal.2018.09.018
[12]	Zhang, L., Wang, C., Zhou, Y., et al. (2021) Research on Inhibition Effects of Low Temperature and Short Illumination Time on Algal Growth. Journal of Shanghai University of Engineering Science, 35, 48-52.
[13]	Gu, X.H., Li, H.M., Mao, Z.G., et al. (2021) Ecological Interaction between Cyanobacterial Blooms and Freshwater Fish (in Chinese). China Science Bulletin, 66, 2649-2662. https://doi.org/10.1360/TB-2020-0905
[14]	Wang, Y., Liu, L.S., Fang, Y.D., et al. (2009) Research Progress on Controlling Lake Eutrophication by Biological Manipulation. Progress in Natural Science, 19, 1296-1301. (In Chinese) https://doi.org/10.1016/j.pnsc.2009.03.009
[15]	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
[16]	Mueller, C.R., Eversole, A.G., Arnold, G., et al. (2004) Effect of Silver Carp Hypophthalmichthys molitrix and Freshwater Mussel Elliptio complanata Filtration on the Phytoplankton Community of Partitioned Aquaculture System Units. Journal of the World Aquaculture Society, 35, 372-382. https://doi.org/10.1111/j.1749-7345.2004.tb00101.x
[17]	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.
[18]	Guo, L., Wang, Q., Xie, P., et al. (2015) A Non-Classical Biomanipulation Experiment in Gonghu Bay of Lake Taihu: Control of Microcystis Blooms Using Silver and Bighead Carps. Aquaculture Research, 46, 2211-2224. https://doi.org/10.1111/are.12375
[19]	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
[20]	Zhuang Z.Y., Yan J.Y., Xie X.Y., et al. (2022) Population Interaction Dynamics Analysis of an Algae-Fish System. Applied Mathematics, 13, 544-565. https://doi.org/10.4236/am.2022.136035
[21]	Zhang, Z.F., Huang, W.Z. and Dong, Z.X. (1992) Qualitative Theory of Differential Equation. Science Press, Beijing.
[22]	Hu, D.P. and Cao, H.J. (2017) Stability and Bifurcation Analysis in a Predator-Prey System with Michaelis-Menten Type Predator Harvesting. Nonlinear Analysis: Real World Applications, 33, 58-82. https://doi.org/10.1016/j.nonrwa.2016.05.010
[23]	Jiang, J. and Song, Y.L. (2014) Delay-Induced Bogdanov-Takens Bifurcation in a Leslie-Gower Predator-Prey Model with Nonmonotonic Functional Response. Communications in Nonlinear Science and Numerical Simulation, 19, 2454-2465. https://doi.org/10.1016/j.cnsns.2013.11.020
[24]	Perko, L. (2001) Differential Equations and Dynamical Systems Differential Equations and Dynamical Systems. Springer, New York. https://doi.org/10.1007/978-1-4613-0003-8
[25]	Sambath, M., Balachandran, K. and Guin, L.N. (2018) Spatiotemporal Patterns in a Predator-Prey Model with Cross-Diffusion Effect. International Journal of Bifurcation and Chaos, 28, Article ID: 1830004. https://doi.org/10.1142/S0218127418300045
[26]	Guin, L.N., Mondal, B. and Chakravarty, S. (2020) Cross-Diffusion-Driven Pattern Formation and Selection in a Modified Leslie-Gower Predator-Prey Model with Fear Effect. Journal of Biological Systems, 28, 27-64. https://doi.org/10.1142/S0218339020500023
[27]	Guin, L.N., Das, E. and Sambath, M. (2020) Pattern Formation Scenario via Turing Instability in Interacting Reaction-Diffusion Systems with Both Refuge and Nonlinear Harvesting. Journal of Applied Nonlinear Dynamics, 9, 1-21. https://doi.org/10.5890/JAND.2020.03.001