Global Dynamics of Modeling Flocculation of Microorganism

Abstract

From a biological perspective, a dynamic model describing the cultivation and flocculation of a microorganism that uses two different kinds of nutrients (carbon source and nitrogen source) is proposed. For the proposed model, there always exists a boundary equilibrium, i.e., $Rhodopseudomonaspalustris$ -free equilibrium. Furthermore, under additional conditions, the model also has five positive equilibria at most, i.e., the equilibria for which carbon source, nitrogen source, $Rhodopseudomonaspalustris$ and flocculants are coexistent. The phenomena of backward and forward bifurcations are extensively discussed by using center manifold theory. The global stability of the boundary equilibrium of the proposed model is deeply investigated. Moreover, the local stability of the positive equilibrium and the uniform persistence of the proposed model are discussed. Under additional conditions, the global stability of the positive equilibrium is studied. Some control strategies are given by the theoretical analysis. Finally, some numerical simulations are performed to confirm the correctness of the theoretical results.

1. Introduction

Photosynthetic bacteria, which are common microorganisms in the natural environment, have been applied in the field of environmental protection, such as in the treatment of sewage, domestic wastewater and the bioremediation of sediment mud polluted with organic matter (see, for example, [1,2,3,4]). On the other hand, photosynthetic bacteria can produce relatively large amounts of physiologically-active substances, such as vitamin $B_{12}$ , ubiquinone (coenzyme $Q10$ ), 5-aminolevulinic acid (ALA) and RNA (see, for example [5]). In particular, vitamin $B_{12}$ has been used in treating anemia and as an eye lotion. Recently, applications as health food supplements have received considerable attention. Coenzyme $Q10$ has been used in treating heart diseases for many years. Further, coenzyme $Q10$ has been used not only as a medicine, but also as some food supplements, because of its physiological activities. One of the developments of ALA applications is in the area of photodynamic diagnosis. RNA is an attractive source of 5${}^{\prime }$-ribonucleotides for use as a flavor enhancer in the food industry. In recent years, the production of RNA has been used as a dietary source of pyrimidine for human immune functions (see, for example [6,7]).

Some photosynthetic bacteria, such as $Rhodopseudomonaspalustris$, are extensively used in the production of lycopene, aquaculture, and so on [8]. It can use sunlight, inorganic and organic compounds for energy. Further, $Rhodopseudomonaspalustris$ can have practical value for removing microcystin from the water body during algal blooms [9]. It can also degrade 2,4,6-trinitrotoluene (TNT), which has negative effects on the human body and aquatic life, resulting in a major threat to drinking and irrigation water supplies, as well as the recreational use of surface waters worldwide. Moreover, $Rhodopseudomonaspalustris$ is regarded as the most promising microbial system for the biological production of hydrogen, which has been extensively developed because of its high-energy content and clean product after combustion [10]. However, the concentration of $Rhodopseudomonaspalustris$ is very low under anaerobic light culture conditions (see, for example, [11,12]). Therefore, cost-efficient harvesting of $Rhodopseudomonaspalustris$ is a new challenge. In order to harvest $Rhodopseudomonaspalustris$ from the liquid, it is necessary to flocculate the single cells into large cell aggregates. Flocculation is a chemically-based separation process that requires less energy than centrifugation and ultrafiltration and, thus, is regarded as the most promising means for degrading microorganisms. Since algal toxins of blooms have happened occasionally in recent years, the problems of degrading microorganisms have received wide attentions (see, for example, [13,14,15,16]).

Flocculants are a kind of important water treatment reagent, which can be divided into organic flocculants and inorganic flocculants according to the chemical compositions [17]. Although organic flocculants, such as polyacrylamide, are frequently used in wastewater treatment and industrial downstream processes because of their high efficiency, some of them are not easily degraded in nature [18,19], and some of the monomers derived from synthetic polymers are harmful to the human body (see, for example, [20,21]). To solve these environmental problems, inorganic flocculants are increasingly being seen as an alternative in the settlement of microorganisms, more specifically in wastewater treatment owing to their inexpensive and nontoxic characteristics. Thus, inorganic flocculants may be used as nontoxic, cost-effective and widely-available flocculants for harvesting $Rhodopseudomonaspalustris$ (see, for example, [22,23,24,25,26,27] and the references therein).

Mathematical models have played an important role in better understanding microbiology and population biology (see, for example, [28,29]). In recent years, the dynamics of the chemostat models has received considerable attentions (see, for example, [29]). The article of Smith and Waltman has played an important role in the development of the chemostat models [30]. From then on, much research on the chemostat models has been extensively studied by many authors. A model describing two populations of microorganisms competing for one single limiting nutrients was proposed in [31]. Later, the model was extended to an arbitrary number of populations in [32,33]. These models, which were studied in articles [31,32,33], have proven that they all include a competitive exclusion effect. In the articles [34,35,36,37,38,39], some further developments have been performed on the chemostat models to place the relevant models in a naturally more sensible manner.

Let $S\left(t\right)$ and $X\left(t\right)$ denote the concentration of the nutrients and $Rhodopseudomonaspalustris$, respectively, in the culture vessel at time t. $P\left(t\right)$ denotes the concentration of inorganic flocculants, which are used for harvesting $Rhodopseudomonaspalustris$ (see Figure 1). The constants $S^0>0$ and $P^0>0$ denote the input concentration of the nutrients and flocculants, respectively. For simplicity, it is assumed that the input of the nutrients and flocculants is continuous. The constant $D>0$ is the dilution rate of the chemostat.

The constant $\tau \ge 0$ denotes the time delay involved in the conversion of nutrients to $Rhodopseudomonas$ $palustris$. The flocculation rate of microorganisms is assumed to be a bilinear mass-action function response $m_{1}X\left(t\right)P\left(t\right)$, where $m_1\ge 0$ is the per capita contact rate. At the same time, flocculants produce loss or consumption [17], and the loss rate of flocculants is also assumed to be a bilinear mass-action function response $m_{2}X\left(t\right)P\left(t\right)$, where $m_2\ge 0$ is constant. Thus, in [40], the following dynamic model has been proposed:

$$\left\{\begin{array}{c}{\displaystyle \frac{dS\left(t\right)}{dt}=(S^0-S\left(t\right))D-r_{1}S\left(t\right)X\left(t\right),}\hfill \\ {\displaystyle \frac{dX\left(t\right)}{dt}=rS(t-\tau )X(t-\tau )-DX\left(t\right)-m_{1}X\left(t\right)P\left(t\right),}\hfill \\ {\displaystyle \frac{dP\left(t\right)}{dt}=(P^0-P\left(t\right))D-m_{2}X\left(t\right)P\left(t\right).}\hfill \end{array}\right.$$

(1)

In Model (1), $r\ge 0$ and $r_1\ge 0$ are constants, and the bilinear mass-action uptake function $S\left(t\right)X\left(t\right)$ has been used.

It should be mentioned here that, the analysis reveals that Model (1) proposed in [40] exhibits the phenomenon of backward bifurcation for the existence of positive equilibria. Moreover, the local stability properties of the equilibria have been dealt with in detail.

People found that the influence of different nutrients has played an important role in the culture of microorganisms. In order to take this into consideration, appropriate combinations of nutrients are considered in chemostat models. Models with two competitors and two perfectly-complemented growth-limiting nutrients are studied in [41,42]. Local asymptotic conditions for the equilibria are derived.

When there is a microorganism to compete for two or more resources, it may become necessary to consider how the resources, once consumed, interact to promote growth. In [41], the authors employ consumer needs to provide a criterion to classify resources and classify resources as perfectly complementary, perfectly substitutable or imperfectly substitutable. Perfectly-complementary resources are different essential substances that must be taken together. In this case, each substance fulfills different functions with respect to the growth of microorganisms. For example, carbon source and nitrogen source may be complementary for the growth of bacterium.

Motivated by the papers mentioned above, in this paper, we further consider a dynamic model describing the cultivation and flocculation of $Rhodopseudomonaspalustris$, and the nutrients presented in [40] will be divided into carbon source and nitrogen source, which are perfectly complementary in the culture of $Rhodopseudomonaspalustris$ (see, for example, [5,41,42,43,44]). We assume that the growth of $Rhodopseudomonas$ $palustris$ is always co-limited by carbon and nitrogen for all possible nutrient conditions. We do not consider the case where only one of these nutrients limits growth, for example, in environmental scenarios of high carbon, but very low nitrogen loads, the growth of bacteria may be purely nitrogen limited.

Let $C\left(t\right)$ and $N\left(t\right)$ denote the concentration of carbon source and nitrogen source, respectively, in the culture vessel at time t. The constants $C^0>0$ and $N^0>0$ denote the input concentration of carbon source and nitrogen source, respectively. We assume that the conversion of nutrients to microorganism biomass occurs instantly. That is, the time delay τ in Model (1) equals zero. Hence, we have the following dynamic model describing the cultivation and flocculation of a microorganism:

$$\left\{\begin{array}{c}{\displaystyle \frac{dC\left(t\right)}{dt}=(C^0-C\left(t\right))D-\frac{X\left(t\right)}{{\delta }_1}{r}_1{\mu }_1\left(C\left(t\right)\right){\mu }_2\left(N\left(t\right)\right),}\hfill \\ {\displaystyle \frac{dN\left(t\right)}{dt}=(N^0-N\left(t\right))D-\frac{X\left(t\right)}{{\delta }_2}{r}_2{\mu }_1\left(C\left(t\right)\right){\mu }_2\left(N\left(t\right)\right),}\hfill \\ {\displaystyle \frac{dX\left(t\right)}{dt}=r{\mu }_1\left(C\left(t\right)\right){\mu }_2\left(N\left(t\right)\right)X\left(t\right)-DX\left(t\right)-m_{1}X\left(t\right)P\left(t\right),}\hfill \\ {\displaystyle \frac{dP\left(t\right)}{dt}=(P^0-P\left(t\right))D-m_{2}X\left(t\right)P\left(t\right).}\hfill \end{array}\right.$$

(2)

In Model (2), the parameters D, r, $m_1$, $m_2$ and $P^0$ are the same as Model (1). The term $r{\mu }_1\left(C\left(t\right)\right){\mu }_2\left(N\left(t\right)\right)$ is the growth rate of $Rhodopseudomonas$ $palustris$, and the terms $r_1{\mu }_1\left(C\left(t\right)\right){\mu }_2\left(N\left(t\right)\right)$ and $r_2{\mu }_1\left(C\left(t\right)\right){\mu }_2\left(N\left(t\right)\right)$ represent the quantity of the decreasing of the carbon source and nitrogen source, respectively, where $r_1$ and $r_2$ are non-negative constants; the functions ${\mu }_1\left(C\left(t\right)\right)$ and ${\mu }_2\left(N\left(t\right)\right)$ are nonnegative and continuous for $C\left(t\right)\ge 0$, $N\left(t\right)\ge 0$. For the simplicity of the theoretical analysis, in this paper, the functions ${\mu }_1\left(C\left(t\right)\right)$ and ${\mu }_2\left(N\left(t\right)\right)$ are chosen as Monod-type functions, i.e.,

$${\mu }_1\left(C\left(t\right)\right)=\frac{C\left(t\right)}{K_1+C\left(t\right)},{\mu }_2\left(N\left(t\right)\right)=\frac{N\left(t\right)}{K_2+N\left(t\right)},$$

where $K_1>0$ and $K_2>0$ are the half-saturation constants with respect to the carbon source and nitrogen source, respectively. ${\delta }_i$ $(i=1,2)$ are yield coefficients (see, for example, [29]), which are defined as:

$${\delta }_i=\frac{massoforganismformed}{massofsubstrateconsumed},(i=1,2).$$

Therefore, the dynamic Model (2) can be rewritten in the following form:

$$\left\{\begin{array}{c}{\displaystyle \frac{dC\left(t\right)}{dt}=(C^0-C\left(t\right))D-\frac{r_{1}C\left(t\right)N\left(t\right)X\left(t\right)}{{\delta }_1(K_1+C\left(t\right))(K_2+N\left(t\right))},}\hfill \\ {\displaystyle \frac{dN\left(t\right)}{dt}=(N^0-N\left(t\right))D-\frac{r_{2}C\left(t\right)N\left(t\right)X\left(t\right)}{{\delta }_2(K_1+C\left(t\right))(K_2+N\left(t\right))},}\hfill \\ {\displaystyle \frac{dX\left(t\right)}{dt}=\frac{rC\left(t\right)N\left(t\right)X\left(t\right)}{(K_1+C\left(t\right))(K_2+N\left(t\right))}-DX\left(t\right)-m_{1}X\left(t\right)P\left(t\right),}\hfill \\ {\displaystyle \frac{dP\left(t\right)}{dt}=(P^0-P\left(t\right))D-m_{2}X\left(t\right)P\left(t\right).}\hfill \end{array}\right.$$

(3)

It is convenient to introduce dimensionless variables. In particular, we define:

$$\overline{C}=\frac{C}{C^0},\overline{N}=\frac{N}{N^0},\overline{X}=X,\overline{P}=\frac{P}{P^0},{\overline{K}}_1=\frac{K_1}{C^0},{\overline{K}}_2=\frac{K_2}{N^0},\overline{t}=tD,$$

$${\overline{r}}_1=\frac{r_1}{{\delta }_{1}D{C}^0},{\overline{r}}_2=\frac{r_2}{{\delta }_{2}D{N}^0},\overline{r}=\frac{r}{D},{\overline{m}}_1=\frac{m_1{P}^0}{D},{\overline{m}}_2=\frac{m_2}{D},$$

and still denote $\overline{C},\overline{N},\overline{X},\overline{P}$, ${\overline{K}}_1$, ${\overline{K}}_2$, $\overline{t}$, ${\overline{r}}_1$, ${\overline{r}}_2$, $\overline{r}$, ${\overline{m}}_1$ and ${\overline{m}}_2$ with $C,N,X,P$, $K_1$, $K_2$, t, $r_1$, $r_2$, r, $m_1$ and $m_2$, then Model (3) becomes:

$$\left\{\begin{array}{c}{\displaystyle \frac{dC\left(t\right)}{dt}=1-C\left(t\right)-\frac{r_{1}C\left(t\right)N\left(t\right)X\left(t\right)}{(K_1+C\left(t\right))(K_2+N\left(t\right))},}\hfill \\ {\displaystyle \frac{dN\left(t\right)}{dt}=1-N\left(t\right)-\frac{r_{2}C\left(t\right)N\left(t\right)X\left(t\right)}{(K_1+C\left(t\right))(K_2+N\left(t\right))},}\hfill \\ {\displaystyle \frac{dX\left(t\right)}{dt}=\frac{rC\left(t\right)N\left(t\right)X\left(t\right)}{(K_1+C\left(t\right))(K_2+N\left(t\right))}-X\left(t\right)-m_{1}X\left(t\right)P\left(t\right),}\hfill \\ {\displaystyle \frac{dP\left(t\right)}{dt}=1-P\left(t\right)-m_{2}X\left(t\right)P\left(t\right).}\hfill \end{array}\right.$$

(4)

According to the biological considerations, the initial condition of Model (4) is given as:

$$\begin{array}{c}\hfill C\left(0\right)=C_0\ge 0,N\left(0\right)=N_0\ge 0,X\left(0\right)=X_0\ge 0,P\left(0\right)=P_0\ge 0,\end{array}$$

(5)

where the constants $C_0,N_0,X_0$ and $P_0$ represent the initial concentration of the carbon source, nitrogen source, $Rhodopseudomonaspalustris$ and flocculants respectively.

The purpose of this paper is to tackle the existence of backward and forward bifurcations by using center manifold theory and to investigate the global stability properties of the two classes of equilibria by constructing the suitable Lyapunov functions.

The organization of the paper is as follows. The global existence, nonnegativity and boundedness of the solutions of Model (4) are investigated in Section 2. In Section 3, the existence of the equilibria and the phenomena of backward and forward bifurcations are extensively discussed. In Section 4, the global stability of the boundary equilibrium of Model (4) is discussed by the stability theory of ordinary differential equations. Furthermore, we consider the local stability of positive equilibrium, the uniform persistence of Model (4) and the global asymptotic stability of the positive equilibrium in Section 5. In Section 6, some control strategies are given by the theoretical analysis. Some discussions are given in Section 7.

2. The Global Existence, Nonnegativity and Boundedness of Solutions

From the biological considerations, it is necessary to show that all of the solutions of Model (4) are nonnegative and bounded for all $t\ge 0$. By using the basic theory of ordinary differential equations [45] and some simple calculations, it is not difficult to show the following result.

Theorem 1. 

The solution $\left(C\right(t),N(t),X(t),P(t\left)\right)$ of Model (4) with the initial Condition (5) is existent, unique and nonnegative for all $t\ge 0$ and satisfies:

$$\underset{t\to \infty }{lim sup}C\left(t\right)\le 1,\underset{t\to \infty }{lim sup}N\left(t\right)\le 1,\underset{t\to \infty }{lim sup}M\left(t\right)\le \alpha ,$$

$$\underset{t\to \infty }{lim sup}P\left(t\right)\le 1,\underset{t\to \infty }{lim inf}C\left(t\right)\ge \underline{C},\underset{t\to \infty }{lim inf}N\left(t\right)\ge \underline{N},\underset{t\to \infty }{lim inf}P\left(t\right)\ge \underline{P},$$

where $M\left(t\right)=\frac{r}{2r_1}C\left(t\right)+\frac{r}{2r_2}N\left(t\right)+X\left(t\right)$, $\alpha =\frac{r}{2r_1}+\frac{r}{2r_2}$, $\underline{C}$ = $\frac{K_1(K_2+1)}{K_1(K_2+1)+r_1\alpha },$ $\underline{N}$ = $\frac{K_2(K_1+1)}{K_2(K_1+1)+r_2\alpha },$ $\underline{P}=\frac{1}{1+m_2\alpha }.$

Proof. 

From the theory of the local existence of solutions for ordinary differential equations, it can be obtained that the solution $\left(C\right(t),N(t),X(t),P(t\left)\right)$ of Model (4) is existent and unique for $t\in [0,\delta )$. Here, δ is some positive constant [29,45,46]. Furthermore, we also have that the solution $\left(C\right(t),N(t),X(t),P(t\left)\right)$ is nonnegative for $t\in [0,\delta )$.

We can easily show that $C\left(t\right)$, $N\left(t\right)$ and $P\left(t\right)$ are bounded on $t\in [0,\delta )$. Let us further show that $X\left(t\right)$ is also bounded on $t\in [0,\delta )$. For $t\ge 0$, define:

$$\begin{array}{c}\hfill M\left(t\right)=\frac{r}{2r_1}C\left(t\right)+\frac{r}{2r_2}N\left(t\right)+X\left(t\right).\end{array}$$

From Model (4), we obtain that, for $t\ge 0$,

$$\begin{array}{c}\hfill \dot{M}\left(t\right)\le \frac{r}{2r_1}+\frac{r}{2r_2}-M\left(t\right).\end{array}$$

(6)

From (6) and the well-known comparison principle, we have that $M\left(t\right)$ is also bounded for $t\in [0,\delta )$. Hence, by employing the continuation theorems of the solutions [29,45,46], the solution $\left(C\right(t),N(t),X(t),P(t\left)\right)$ is existent and unique for any $t\ge 0$. Similarly, the solution is nonnegative for any $t\ge 0$.

Thus, from the comparison principle, we have that:

$$\begin{array}{c}\hfill \underset{t\to \infty }{lim sup}C\left(t\right)\le 1,\underset{t\to \infty }{lim sup}N\left(t\right)\le 1,\underset{t\to \infty }{lim sup}M\left(t\right)\le \frac{r}{2r_1}+\frac{r}{2r_2}=\alpha ,\underset{t\to \infty }{lim sup}P\left(t\right)\le 1.\end{array}$$

By the first equation of Model (4), we obtain that, for $t\ge 0$,

$$\begin{array}{cc}\hfill \dot{C}\left(t\right)& =1-C\left(t\right)-\frac{r_{1}C\left(t\right)N\left(t\right)X\left(t\right)}{(K_1+C\left(t\right))(K_2+N\left(t\right))}\ge 1-\left(1+\frac{r_1\alpha }{K_1(K_2+1)}\right)C\left(t\right).\hfill \end{array}$$

Again, we can conclude that $\underset{t\to \infty }{lim inf}C\left(t\right)\ge \underline{C}.$ By using the technique similar above, we can show that $\underset{t\to \infty }{lim inf}N\left(t\right)\ge \underline{N},$ $\underset{t\to \infty }{lim inf}P\left(t\right)\ge \underline{P}.$ We complete the proof of Theorem 1. □

Theorem 2. 

The compact set:

$$\begin{array}{c}\hfill \mathsf{\Omega }=\{(C,N,X,P)\in R_{+}^4|\underline{C}\le C\le 1,\underline{N}\le N\le 1,0\le M\le \alpha ,\underline{P}\le P\le 1\}\end{array}$$

attracts all of the solutions of Model (4) and is positively invariant with respect to Model (4).

Proof. 

According to Theorem 1, it only needs to be proven that Ω is positively invariant with respect to Model (4). That is, it needs to be shown that $\underline{C}\le C\left(t\right)\le 1,\underline{N}\le N\left(t\right)\le 1,M\left(t\right)\le \alpha ,\underline{P}\le P\left(t\right)\le 1$ for any $t\ge 0$ if $\left(C\right(0),N(0),X(0),P(0\left)\right)\in \mathsf{\Omega }$. Let us show $M\left(t\right)\le \alpha$ for any $t\ge 0$.

In fact, if there exists some $t_1>0$, such that $M\left(t_1\right)>\alpha$, then $t_1^{*}=sup\left\{t\right|M\left(t\right)=\alpha ,t\in [0,t_1]\}$ is existent and $t_1^{*}\ge 0$. Hence, we obtain that $M\left(t_1^{*}\right)=\alpha$, $M\left(t_1\right)>\alpha$ and $M\left(t\right)>\alpha$ for $t\in (t_1^{*},t_1)$. By the Lagrange mean-value theorem, there exists some $t_2\in (t_1^{*},t_1)$, such that $\dot{M}\left(t_2\right)=\frac{M\left(t_1\right)-M\left(t_1^{*}\right)}{t_1-t_1^{*}}>0.$

On the other hand, from (6), we have that $\dot{M}\left(t_2\right)\le \alpha -M\left(t_2\right)<\alpha -\alpha =0,$ which is a contradiction. Thus, $M\left(t\right)\le \alpha$ for any $t\ge 0$. Therefore, from Model (4), we have that for any $t\ge 0$, $\dot{C}\left(t\right)\le 1-C\left(t\right),\dot{N}\left(t\right)\le 1-N\left(t\right),\dot{P}\left(t\right)\le 1-P\left(t\right)$, from which we easily have that $C\left(t\right)\le 1$, $N\left(t\right)\le 1$, $P\left(t\right)\le 1$ for any $t\ge 0$.

Next, we prove that $C\ge \underline{C}$ for any $t\ge 0$ if $\left(C\right(0)$, $N\left(0\right)$, $X\left(0\right)$, $P\left(0\right))\in \mathsf{\Omega }$.

In fact, if there exists some $t_3>0$, such that $C\left(t_3\right)<\underline{C}$, then $t_2^{*}=sup\left\{t\right|C\left(t\right)=\underline{C},t\in [0,t_3]\}$ is existent and $t_2^{*}\ge 0$. Hence, we have that $C\left(t_2^{*}\right)=\underline{C}$, $C\left(t_3\right)<\underline{C}$ and $C\left(t\right)<\underline{C}$ for $t\in (t_2^{*},t_3)$. By the Lagrange mean-value theorem, there exists some $t_4\in (t_2^{*},t_3)$, such that $\dot{C}\left(t_4\right)=\frac{C\left(t_3\right)-C\left(t_2^{*}\right)}{t_3-t_2^{*}}<0.$

On the other hand, from Model (4), we obtain that:

$$\begin{array}{cc}\hfill \dot{C}\left(t_4\right)& \ge 1-C\left(t_4\right)-\frac{r_{1}C\left(t_4\right)\alpha }{K_1(K_2+1)}>1-\left(1+\frac{r_1\alpha }{K_1(K_2+1)}\right)\underline{C}=0,\hfill \end{array}$$

which is a contradiction. Thus, $C\left(t\right)\ge \underline{C}$ for any $t\ge 0$. Therefore, from Model (4), we obtain that, for any $t\ge 0$,

$$\begin{array}{cc}\hfill \dot{N}\left(t\right)& \ge 1-\left(1+\frac{r_2\alpha }{K_2(K_1+1)}\right)N\left(t\right),\dot{P}\left(t\right)\ge 1-(1+m_2\alpha )P\left(t\right),\hfill \end{array}$$

from which we easily obtain that $\underline{N}\le N\left(t\right)$, $\underline{P}\le P\left(t\right)$, for any $t\ge 0$. This completes the proof of Theorem 2. □

3. The Existence of the Equilibria and Its Classification

Let $(C,N,X,P)$ be any equilibrium of Model (4). Then, $(C,N,X,P)$ satisfies the following nonlinear algebraic equations,

$$\left\{\begin{array}{c}{\displaystyle 1-C-\frac{r_{1}CNX}{(K_1+C)(K_2+N)}=0,}\hfill \\ {\displaystyle 1-N-\frac{r_{2}CNX}{(K_1+C)(K_2+N)}=0,}\hfill \\ {\displaystyle \frac{rCNX}{(K_1+C)(K_2+N)}-X-m_{1}XP=0,}\hfill \\ {\displaystyle 1-P-m_{2}XP=0.}\hfill \end{array}\right.$$

(7)

Model (4) always has the boundary equilibrium $E_0(1,1,0,1)$. The existence of $E_0$ indicates that, if there is no $Rhodopseudomonaspalustris$ to be added into the culture vessel at the beginning of the culture, the concentrations of the carbon source, nitrogen source and flocculants always maintain the constant values 1, 1 and 1, respectively. The equilibrium $E_0(1,1,0,1)$ is also called $Rhodopseudomonaspalustris$-free equilibrium.

Define the basic bifurcation parameter as:

$$R_0=\frac{r}{(K_1+1)(K_2+1)(m_1+1)}.$$

Let $(C^{*},N^{*},X^{*},P^{*})$ be any positive equilibrium of Model (4). From (7), we have that:

$$\left\{\begin{array}{c}{\displaystyle (K_1+C^{*})(K_2+N^{*})(1+\frac{m_1}{1+m_2{X}^{*}})-r{C}^{*}{N}^{*}=0,}\hfill \\ {\displaystyle P^{*}=\frac{1}{1+m_2{X}^{*}},}\hfill \\ {\displaystyle C^{*}=\frac{-r_1{m}_2{\left(X^{*}\right)}^2+(m_{2}r-r_1(m_1+1))X^{*}+r}{r(1+m_2{X}^{*})},}\hfill \\ {\displaystyle N^{*}=\frac{-r_2{m}_2{\left(X^{*}\right)}^2+(m_{2}r-r_2(m_1+1))X^{*}+r}{r(1+m_2{X}^{*})}.}\hfill \end{array}\right.$$

(8)

Clearly, $X^{*}$ should satisfy the following conditions:

$$\begin{array}{c}\hfill -r_1{m}_2{\left(X^{*}\right)}^2+\left(m_{2}r-r_1(m_1+1)\right)X^{*}+r>0,\\ \hfill -r_2{m}_2{\left(X^{*}\right)}^2+\left(m_{2}r-r_2(m_1+1)\right)X^{*}+r>0.\end{array}$$

(9)

Substituting the second, third and forth equations of (8) into the first equation gives a fifth order algebraic equation,

$$\begin{array}{cc}\hfill f\left(X^{*}\right)=& a{\left(X^{*}\right)}^5+b{\left(X^{*}\right)}^4+c{\left(X^{*}\right)}^3+d{\left(X^{*}\right)}^2+e{X}^{*}+f=0,\hfill \end{array}$$

where:

$$\begin{array}{cc}\hfill a=& r_1{r}_2{m}_2^3(1-r),\hfill \\ \hfill b=& r_1{m}_2^{3}r\left(r-(K_2+1)-\frac{r_2(m_1+2)}{m_2}\right)+3r_1{r}_2{m}_2^2(m_1+1)\hfill \\ & +r_2{m}_2^{3}r\left(r-(K_1+1)-\frac{r_1(m_1+1)}{m_2}\right),\hfill \\ \hfill c=& m_2^3{r}^2\left((K_1+1)(K_2+1)+\frac{(r_1+r_2)(m_1+3)}{m_2}-r\right)+q,\hfill \\ \hfill d=& m_2^{2}r{R}_0{(K_1+1)}^2{(K_2+1)}^2{(m_1+1)}^2(1-R_0)+d_1,\hfill \\ \hfill e=& r_{1}r(m_1+1)\left(r-(K_2+1)(m_1+1)\right)+r_{2}r(m_1+1)\left(r-(K_1+1)(m_1+1)\right)\hfill \\ & +m_2{r}^2\left((K_1+1)(K_2+1)(2m_1+3)-3r\right),\hfill \\ \hfill f=& r{R}_0{(K_1+1)}^2{(K_2+1)}^2{(m_1+1)}^2(1-R_0),\hfill \\ \hfill q=& r_1{r}_2{m}_2(3-r){(1+m_1)}^2-2q_{1}r(1+m_1)-r{m}_2^2(K_1{r}_2+K_2{r}_1),\hfill \\ \hfill q_1=& r_1{m}_2^2+r_1{r}_2{m}_2+r_2{m}_2^2+K_1{r}_2{m}_2^2+K_2{r}_1{m}_2^2,\hfill \\ \hfill d_1=& m_2^2{r}^2\left(2(K_1+1)(K_2+1)+\frac{(r_1+r_2)(2m_1+3)}{m_2}-2r\right)+h,\hfill \\ \hfill h=& -m_2^{2}r\left(2(m_1+1)+1\right)(r_2{K}_1+r_1{K}_2)+2r_1{r}_2{m}_2{(1+m_1)}^2\hfill \\ & -r_1{m}_2^{2}r-r_2{m}_2^{2}r-h_{1}r(1+m_1),\hfill \\ \hfill h_1=& 2r_1{m}_2^2+3r_1{r}_2{m}_2+r_2{m}_2^2+r_2{m}_2^2+K_2{m}_2^2.\hfill \end{array}$$

Let us consider the necessary condition for the existence of the positive equilibria of Model (4).

From the first equation in (8), we obtain the following function:

$$\begin{array}{cc}\hfill F\left(X^{*}\right)& =\frac{m_1(K_1+C^{*})(K_2+N^{*})}{1+m_2{X}^{*}}+K_1{K}_2+K_1{N}^{*}+K_2{C}^{*}+(1-r)C^{*}{N}^{*},\hfill \end{array}$$

which implies that $r>1$ is a necessary condition for a positive equilibrium to exist.

Using the methods similar to [47], we can give the sufficient conditions of the existence of the positive equilibria of Model (4). The following results (Theorem 3) follow from the various possibilities enumerated in

Table A1. The number of the positive roots of $f\left(X^{*}\right)$ for $R_0<1$ and $R_0>1$.

Cases	a	b	c	d	e	f	R0	m	n
1	−	+	+	+	+	+	$R_0<1$	1	1
	−	+	+	+	+	−	$R_0>1$	2	$0,2$
2	−	+	+	+	−	+	$R_0<1$	3	$1,3$
	−	+	+	+	−	−	$R_0>1$	2	$0,2$
3	−	+	+	−	+	+	$R_0<1$	3	$1,3$
	−	+	+	−	+	−	$R_0>1$	4	$0,2,4$
4	−	+	+	−	−	+	$R_0<1$	3	$1,3$
	−	+	+	−	−	−	$R_0>1$	2	$0,2$
5	−	+	−	+	+	+	$R_0<1$	3	$1,3$
	−	+	−	+	+	−	$R_0>1$	4	$0,2,4$
6	−	+	−	+	−	+	$R_0<1$	5	$1,3,5$
	−	+	−	+	−	−	$R_0>1$	4	$0,2,4$
7	−	+	−	−	+	+	$R_0<1$	3	$1,3$
	−	+	−	−	+	−	$R_0>1$	4	$0,2,4$
8	−	+	−	−	−	+	$R_0<1$	3	$1,3$
	−	+	−	−	−	−	$R_0>1$	2	$0,2$
9	−	−	+	+	+	+	$R_0<1$	1	1
	−	−	+	+	+	−	$R_0>1$	2	$0,2$
10	−	−	+	+	−	+	$R_0<1$	3	$1,3$
	−	−	+	+	−	−	$R_0>1$	2	$0,2$
11	−	−	+	−	+	+	$R_0<1$	3	$1,3$
	−	−	+	−	+	−	$R_0>1$	4	$0,2,4$
12	−	−	+	−	−	+	$R_0<1$	3	$1,3$
	−	−	+	−	−	−	$R_0>1$	2	$0,2$
13	−	−	−	+	+	+	$R_0<1$	1	1
	−	−	−	+	+	−	$R_0>1$	2	$0,2$
14	−	−	−	+	−	+	$R_0<1$	3	$1,3$
	−	−	−	+	−	−	$R_0>1$	2	$0,2$
15	−	−	−	−	+	+	$R_0<1$	1	1
	−	−	−	−	+	−	$R_0>1$	2	$0,2$
16	−	−	−	−	−	+	$R_0<1$	1	1
	−	−	−	−	−	−	$R_0>1$	0	0

(see Appendix A):

Theorem 3. 

(i)

Model (4) has a unique positive equilibrium if $R_0<1$ and Conditions (9) hold and whenever Cases 1, 9, 13, 15 and 16 in Table A1 are satisfied;

(ii)

Model (4) could have more than one positive equilibrium if $R_0<1$ and Conditions (9) hold and whenever Cases 2–8, 10–12 and 14 in Table A1 are satisfied;

(iii)

Model (4) could have five positive equilibria at most if $R_0>1$ and Conditions (9) hold and whenever Cases 1–16 in Table A1 are satisfied.

Hence, under suitable conditions, there may be at most five different positive roots for the fifth order algebraic equation. Let $X=X^{*}$ be any such positive root, which also satisfies Conditions $\left(9\right)$. Thus, from (8), $C=C^{*}>0$, $N=N^{*}>0$ and $P=P^{*}>0$ can be obtained. Therefore, Model (4) at most has five positive equilibria of the type of $E^{*}(C^{*},N^{*},X^{*},P^{*})$. The equilibrium $E^{*}(C^{*},N^{*},X^{*},P^{*})$ indicates that the carbon source, nitrogen source, $Rhodopseudomonaspalustris$ and flocculants may be coexistent for any time $t\ge 0$.

Remark 1. 

The existence of multiple positive equilibria of Model (4) when $R_0<1$ (shown in Table A1; see Appendix A) indicates the possibility of the existence of backward bifurcation (see, for example, [40,47,48]), where the stable boundary equilibrium co-exists with a stable positive equilibrium. This can be explored below via numerical simulations (see Figure 2a,b). A rigorous result can be obtained using center manifold theory [49]. The detailed proof is given in Appendix B.

From the analysis given in Appendix B, we have established the following result.

Theorem 4. 

(i)

If:

$$\begin{array}{c}\hfill m_1{m}_2>\frac{r{r}_1{K}_1}{{(K_1+1)}^3{(K_2+1)}^2}+\frac{r{r}_2{K}_2}{{(K_1+1)}^2{(K_2+1)}^3},\end{array}$$

then Model (4) undergoes a backward bifurcation at $R_0=1$.

(ii)

If:

$$\begin{array}{c}\hfill m_1{m}_2<\frac{r{r}_1{K}_1}{{(K_1+1)}^3{(K_2+1)}^2}+\frac{r{r}_2{K}_2}{{(K_1+1)}^2{(K_2+1)}^3},\end{array}$$

then Model (4) undergoes a forward bifurcation at $R_0=1$.

The existence of backward bifurcation implies that stable boundary equilibrium and stable positive equilibrium may be coexistent. In biology, this means that the basic bifurcation parameter $R_0$ is not the threshold value, which is used to determine whether $Rhodopseudomonaspalustris$ can be harvested successfully or not. In this case, more complicated dynamic properties may occur. There may exist a new threshold value, which is less than $R_0$ and used to determine whether $Rhodopseudomonaspalustris$ can be harvested successfully or not.

4. The Global Stability of the Boundary Equilibrium

Global stability properties of the equilibria $E_0$ or $E^{*}$ imply that the asymptotic properties of the carbon source, nitrogen source, $Rhodopseudomonaspalustris$ and flocculants in the culture vessel are not dependent on the initial values $C_0,N_0,X_0$ and $P_0$. For the global stability property of the boundary equilibrium $E_0$ of Model (4), we have the following result.

Theorem 5. 

If $R_0<1$, then the boundary equilibrium $E_0(1,1,0,1)$ of Model (4) is locally asymptotically stable. Further, if:

$$\begin{array}{c}\hfill R_0\le \frac{m_1\underline{P}+1}{m_1+1}<1,\end{array}$$

(10)

then the boundary equilibrium $E_0$ of Model (4) is globally asymptotically stable.

Proof. 

It is easy to show that the boundary equilibrium $E_0(1,1,0,1)$ of Model (4) is locally asymptotically stable by the characteristic equation of the linearization of Model (4). Next, we prove the global asymptotic stability of the boundary equilibrium $E_0$ of Model (4).

Since Ω is attractive and positively forward invariant for Model (4), hence it just considers Model (4) in Ω. Define:

$$\begin{array}{cc}\hfill V_1=& X.\hfill \end{array}$$

Apparently, $V_1(1,1,0,1)=0$ and $V_1$ is continuous on Ω. If $R_0\le \frac{m_1\underline{P}+1}{m_1+1}<1$, then the derivative of $V_1$ along the solutions of Model (4) is:

$$\begin{array}{cc}\hfill {\dot{V}}_1=& \frac{rC\left(t\right)N\left(t\right)X\left(t\right)}{(K_1+C\left(t\right))(K_2+N\left(t\right))}-X\left(t\right)-m_{1}X\left(t\right)P\left(t\right)\hfill \\ \hfill \le & (\frac{rC\left(t\right)N\left(t\right)}{(K_1+C\left(t\right))(K_2+N\left(t\right))}-1-m_1\underline{P})X\left(t\right)\hfill \\ \hfill \le & (m_1+1)(\frac{rC\left(t\right)N\left(t\right)}{(K_1+C\left(t\right))(K_2+N\left(t\right))}\frac{1}{m_1+1}-R_0)X\left(t\right)\hfill \\ \hfill \le & (m_1+1)(\frac{r}{(K_1+1)(K_2+1)}\frac{1}{m_1+1}-R_0)X\left(t\right)=0\hfill \end{array}$$

for any $t\ge 0$. Hence, $V_1$ is a Lyapunov function of Model (4) on Ω.

Define $E=\left\{(C,N,X,P)\right|(C,N,X,P)\in \mathsf{\Omega },{\dot{V}}_1=0\}.$ We have that:

$$E\subset \left\{\right(C,N,X,P\left)\right|(C,N,X,P)\in \mathsf{\Omega },X=0,\text{or}C=1\text{and}N=1\}.$$

Let $\widehat{M}$ be the largest set in E, which is invariant with respect to Model (4). Clearly, $\widehat{M}$ is not empty, since $(1,1,0,1)\in \widehat{M}.$ For any $(C_0,N_0,X_0,P_0)\in \widehat{M}$, let $\left(C\right(t),N(t),X(t),P(t\left)\right)$ be the solution of Model (4) with the initial Condition (5). From the invariance of $\widehat{M}$, we get $(C\left(t\right),N\left(t\right),X\left(t\right),P\left(t\right))\in \widehat{M}\subseteq E$ for any $t\in R$. Thus, we get, for each t, $X\left(t\right)=0$, or $C\left(t\right)=1$, and $N\left(t\right)=1$.

If for some $\tilde{t}$, $C\left(\tilde{t}\right)=1$ and $N\left(\tilde{t}\right)=1$, then we have that $\dot{C}\left(\tilde{t}\right)=\dot{N}\left(\tilde{t}\right)=0$. Hence, from the first or second equation of Model (4), we obtain that $X\left(\tilde{t}\right)=0$. Thus, for any $t\in R$, we have that $X\left(t\right)\equiv 0$.

Subsequently, from the first, second and forth equations of Model (4), we have that, for any $t\in R$,

$$\dot{C}\left(t\right)=1-C\left(t\right),\dot{N}\left(t\right)=1-N\left(t\right),\dot{P}\left(t\right)=1-P\left(t\right).$$

Furthermore, for any $t\in R$, we have that:

$$C\left(t\right)=1-(1-C\left(0\right))e^{-t},N\left(t\right)=1-(1-N\left(0\right))e^{-t},P\left(t\right)=1-(1-P\left(0\right))e^{-t}.$$

If $C\left(0\right)<1$ or $N\left(0\right)<1$ or $P\left(0\right)<1$, then $C\left(t\right)$, $N\left(t\right)$ and $P\left(t\right)$ become negative values $(t\to -\infty )$. Then, we obtain that $C\left(0\right)=1$, $N\left(0\right)=1$ and $P\left(0\right)=1$. Hence, we obtain that, for any $t\in R$, $C\left(t\right)=N\left(t\right)=P\left(t\right)\equiv 1.$ Therefore, we obtain $\widehat{M}=\left\{(1,1,0,1)\right\}.$ The classical Lyapunov–LaSalle invariance principle shows that $E_0$ is globally attractive. Since it has been shown that, if $R_0\le \frac{m_1\underline{P}+1}{m_1+1}<1$, the boundary equilibrium $E_0$ of Model (4) is locally asymptotically stable. Hence, the boundary equilibrium $E_0$ of Model (4) is globally asymptotically stable. □

Remark 2. 

Let us use some numerical simulations to check the correctness of the theoretical analyses. We set

$$r=1.5,K_1=0.36,K_2=0.3,m_1=1.5,m_2=4,r_1=0.96,r_2=1.0001.$$

By the analysis of Theorem 5, we obtain that the boundary equilibrium $E_0$ of Model (4) is globally asymptotically stable (see Figure 3a).

Furthermore, if the parameters are chosen as follows, $r=3.5,$ $K_1=0.36,$ $K_2$ $=0.3,$ $m_1=1.5,m_2=4,$ $r_1=0.96,r_2=1.0001,$ it is easy to check that the condition $R_0<1$ holds, but there are two positive equilibria, $E_1^{*}(C^{*},N^{*},$ $X^{*},P^{*})$ $\approx (0.4783,0.4565,1.5782,$ $0.1376)$ and $E_2^{*}(C^{*},N^{*},$ $X^{*},P^{*})$ $\approx (0.9100,0.9063,0.1741,$ $0.5895).$ Then, the backward bifurcation phenomenon is illustrated. Therefore, sufficient Condition (10) is reasonable.

5. Stability and Uniform Persistence

In this section, we give the local asymptotic stability of the positive equilibrium and uniform persistence of Model (4). Further, under additional conditions, the global asymptotic stability of the positive equilibrium is obtained by using some techniques of constructing Lyapunov functions.

5.1. The Local Asymptotic Stability of the Positive Equilibrium

For the local asymptotic stability of the positive equilibrium of Model (4), we have the following result.

Theorem 6. 

If the positive equilibrium $E^{*}(C^{*},N^{*},X^{*},P^{*})$ of Model (4) exists and the condition:

$$\begin{array}{c}\hfill m_1<\frac{K_1(1-C^{*})}{P^{*}{C}^{*}(K_1+C^{*})}+\frac{K_2(1-N^{*})}{P^{*}{N}^{*}(K_2+N^{*})}\end{array}$$

(11)

holds, then the positive equilibrium $E^{*}(C^{*},N^{*},X^{*},P^{*})$ of Model (4) is locally asymptotically stable.

Proof. 

For convenience, we assume that:

$$\begin{array}{cc}\hfill \Delta =& \frac{K_1{r}_1{N}^{*}{X}^{*}}{{(K_1+C^{*})}^2(K_2+N^{*})}+\frac{K_2{r}_2{C}^{*}{X}^{*}}{(K_1+C^{*}){(K_2+N^{*})}^2}.\hfill \end{array}$$

For the linearized system, the corresponding characteristic equation of Model (4) can be expressed as follows:

$$\begin{array}{c}\hfill {\lambda }^4+A_1{\lambda }^3+A_2{\lambda }^2+A_3\lambda +A_4=0,\end{array}$$

where:

$$\begin{array}{cc}\hfill A_1=& 3+m_2{X}^{*}+\Delta ,\hfill \\ \hfill A_2=& \Delta (m_1{P}^{*}+2+m_2{X}^{*})+1+\Delta +2(1+m_2{X}^{*})-m_1{m}_2{P}^{*}{X}^{*},\hfill \\ \hfill A_3=& \Delta (2m_1{P}^{*}+3+2m_2{X}^{*})-2m_1{m}_2{P}^{*}{X}^{*}+1+m_2{X}^{*},\hfill \\ \hfill A_4=& m_2{X}^{*}(\Delta -m_1{P}^{*})+(1+m_1{P}^{*})\Delta .\hfill \end{array}$$

According to the Routh-Hurwitz criterion, we need to show that:

$$\begin{array}{c}\hfill A_1{A}_2-A_3>0,A_3(A_1{A}_2-A_3)-A_1^2{A}_4>0,A_4>0.\end{array}$$

For the sake of simplicity, we define:

$$\begin{array}{cc}\hfill {\Delta }_1=& A_1,{\Delta }_2=A_1{A}_2-A_3,{\Delta }_3=A_3(A_1{A}_2-A_3)-A_1^2{A}_4,{\Delta }_4=A_4{\Delta }_3.\hfill \end{array}$$

By computation, it can be obtained that:

$$\begin{array}{cc}\hfill {\Delta }_1& =A_1>0,\hfill \\ \hfill {\Delta }_2& =m_2{X}^{*}(1+\Delta +m_2{X}^{*})(\Delta -m_1{P}^{*})+(2+2m_2{X}^{*}+\Delta )(3+m_2{X}^{*})\hfill \\ & +\Delta (6+m_1{P}^{*})+\Delta (m_2{X}^{*}+\Delta )(3+m_1{P}^{*}),\hfill \\ \hfill {\Delta }_3& =P_1(\Delta -m_1{P}^{*})+Q,\hfill \\ \hfill P_1& =2m_2{X}^{*}\Delta (1+\Delta )(3+m_1{P}^{*}+m_2{X}^{*})+4m_2{X}^{*}\left(1+\Delta (1+m_2{X}^{*})\right)\hfill \\ & +4m_2{X}^{*}(1+m_2{X}^{*})(m_2{X}^{*}+3)+2{\left(m_2{X}^{*}\right)}^2\Delta (1+m_2{X}^{*})\hfill \\ & +m_2{X}^{*}\Delta (1+m_2{X}^{*})(2m_1{P}^{*}+3+2m_2{X}^{*}),\hfill \\ \hfill Q& =a_1(1+m_1{P}^{*})+a_2(3+m_1{P}^{*}+m_2{X}^{*})+a_3(1+m_2{X}^{*})+a_4{m}_1{P}^{*}{m}_2{X}^{*}\Delta \hfill \\ & +2a_5{m}_1{P}^{*}{m}_2{X}^{*}+a_6{m}_2{X}^{*}\Delta +{\Delta }^2(2+3m_2{X}^{*})(3+2m_2{X}^{*})+14\Delta +4{\Delta }^2,\hfill \\ \hfill a_1& =2{\Delta }^3+8{\Delta }^2+7\Delta +m_2{X}^{*}\Delta (m_2{X}^{*}+7)+m_1{P}^{*}{\Delta }^2(\Delta +2),\hfill \\ \hfill a_2& =\Delta (1+\Delta )(1+m_2{X}^{*})+{\Delta }^3(2+m_1{P}^{*})+{\Delta }^2,\hfill \\ \hfill a_3& =2\left(1+(1+m_2{X}^{*})(3+m_2{X}^{*})\right),\hfill \\ \hfill a_4& =2{\left(m_2{X}^{*}\right)}^2+{(1+\Delta )}^2,\hfill \\ \hfill a_5& =4+{\left(m_2{X}^{*}\right)}^2+m_1{P}^{*}{m}_2{X}^{*}+2m_2{X}^{*},\hfill \\ \hfill a_6& =4m_2{X}^{*}+4\Delta +13,\hfill \\ \hfill {\Delta }_4& =A_4{\Delta }_3.\hfill \end{array}$$

If:

$$\Delta =\frac{K_1{r}_1{N}^{*}{X}^{*}}{{(K_1+C^{*})}^2(K_2+N^{*})}+\frac{K_2{r}_2{C}^{*}{X}^{*}}{(K_1+C^{*}){(K_2+N^{*})}^2}>m_1{P}^{*},$$

that is,

$$\begin{array}{c}\hfill m_1<\frac{K_1(1-C^{*})}{P^{*}{C}^{*}(K_1+C^{*})}+\frac{K_2(1-N^{*})}{P^{*}{N}^{*}(K_2+N^{*})},\end{array}$$

we obtain that ${\Delta }_1>0$, ${\Delta }_2>0$, ${\Delta }_3>0,$ and ${\Delta }_4>0.$ Hence, from the Routh-Hurwitz criterion, we obtain that the positive equilibrium $E^{*}(C^{*},N^{*},X^{*},P^{*})$ of Model (4) is locally asymptotically stable. This completes the proof of Theorem 6. □

5.2. Uniform Persistence

As pointed out in [29,45,50], uniform persistence is an important concept in the cultivation of microorganisms. From a biological perspective, one basic question about biological models involves the long-time survival of the species. To this end, we may want to know whether the constructed model is uniformly persistent with respect to one or more species. That is, the number of the species will keep positive and bounded away from zero for any positive time. In recent years, this kind of research has been particularly common in microorganisms (see, for example, [45]), especially in the simulation of harmful algae growth, where we usually want to predict if the harmful algae will go towards extinction or if blooms of the harmful algae will remain, which directly impacts human health and food webs in aquatic ecosystems. From a mathematical point of view, uniform persistence can give sufficient criteria for the existence of a positive equilibrium for the dissipative system [51].

Model (4) is said to be uniformly persistent if there are positive constants $n_1,n_2,n_3$, $n_4,N_1,N_2,N_3,N_4$, such that each positive solution $\left(C\right(t),N(t),X(t),P(t\left)\right)$ of Model (4) satisfies:

$$\begin{array}{c}\hfill n_1\le \underset{t\to \infty }{lim inf}C\left(t\right)\le \underset{t\to \infty }{lim sup}C\left(t\right)\le N_1,\\ \hfill n_2\le \underset{t\to \infty }{lim inf}N\left(t\right)\le \underset{t\to \infty }{lim sup}N\left(t\right)\le N_2,\\ \hfill n_3\le \underset{t\to \infty }{lim inf}X\left(t\right)\le \underset{t\to \infty }{lim sup}X\left(t\right)\le N_3,\\ \hfill n_4\le \underset{t\to \infty }{lim inf}P\left(t\right)\le \underset{t\to \infty }{lim sup}P\left(t\right)\le N_4.\end{array}$$

Now, we give a result on the uniform persistence of Model (4). To proceed, we introduce the following notation and terminology. Denote by $R\left(t\right)$ $(t\ge 0)$ the family of solution operators corresponding to Model (4). The ω-limit set $\omega \left(x\right)$ of x consists of y , such that there exists a sequence $t_n\to \infty$ as $n\to \infty$ with $R\left(t_n\right)x\to y$ as $n\to \infty$. Define:

$$G=\{(C,N,X,P)\in R_{+}^4\mid \underline{C}\le C\le 1,\underline{N}\le N\le 1,X\ge 0,\underline{P}\le P\le 1\},$$

$$G_0=\{(C,N,X,P)\in G\mid \underline{C}\le C\le 1,\underline{N}\le N\le 1,X>0,\underline{P}\le P\le 1\},$$

$$\partial G_0=G\backslash G_0,$$

$M_{\partial }=\{(C,N,X,P)\in \partial G_0\mid R\left(t\right)(C,N,X,P)$ satisfies Model (4) and $R\left(t\right)(C,N,X,P)\in \partial G_0,\forall t\ge 0\},$

$$\mathsf{\Omega }\left(M_{\partial }\right)={\cup }_{x\in M_{\partial }}\omega \left(x\right).$$

From Theorem 2, we obtain that G is positively invariant corresponding to $R\left(t\right)$. It is easy to see that $G_0$ is also positively invariant.

Clearly, $\partial G_0$ is relatively compact in G. Now, we show that $\mathsf{\Omega }\left(M_{\partial }\right)=\left\{(1,1,0,1)\right\}.$

In fact, $\left\{(1,1,0,1)\right\}\subseteq \mathsf{\Omega }\left(M_{\partial }\right).$ For any $\left(C\right(0),$ $N\left(0\right),$ $X\left(0\right),$ $P\left(0\right))$ ∈ $M_{\partial }$, it has that, for all $t\ge 0$, $X\left(t\right)\equiv 0$ and ${lim}_{t\to +\infty }C\left(t\right)=$ ${lim}_{t\to +\infty }N\left(t\right)=$ ${lim}_{t\to +\infty }P\left(t\right)=1.$ Thus, $\mathsf{\Omega }\left(M_{\partial }\right)=\left\{(1,1,0,1)\right\}$.

Theorem 7. 

If $R_0>1$, then Model (4) is uniformly persistent.

Proof. 

By Theorem 1, it can be obtained that nonnegative solutions of Model (4) are point dissipative. According to the definitions above, it suffices to show that $\partial G_0$ repels uniformly nonnegative solutions of Model (4). It is obvious that there is only one equilibrium $E_0$ in $M_{\partial }$.

We now show that $W^s\left(E_0\right)\cap G_0=\varnothing$. Assume $W^s\left(E_0\right)\cap G_0\ne \varnothing$, then there exists a positive solution of Model (4), such that $\underset{t\to \infty }{lim}(C\left(t\right),N\left(t\right),X\left(t\right),P\left(t\right))=(1,1,0,1)$. Assume:

$$\begin{array}{c}\hfill U(C\left(t\right),N\left(t\right),P\left(t\right))=\frac{rC\left(t\right)N\left(t\right)}{(K_1+C\left(t\right))(K_2+N\left(t\right))}-1-m_{1}P\left(t\right),\end{array}$$

then:

$$\begin{array}{c}\hfill U(1,1,1)=\frac{r}{(K_1+1)(K_2+1)}-1-m_1.\end{array}$$

Then, $r>(K_1+1)(K_2+1)(1+m_1)$ is equivalent to $U(1,1,1)>0$.

Since $\underset{t\to \infty }{lim}U(C\left(t\right),N\left(t\right),P\left(t\right))=U(1,1,1)>0$, there exists $T>0$, such that, for any $t\ge T$,

$$\begin{array}{c}\hfill U(C\left(t\right),N\left(t\right),P\left(t\right))>\frac{1}{2}U(1,1,1)>0.\end{array}$$

By the third equation of Model (4), we obtain that, for any $t\ge T$,

$$\begin{array}{cc}\hfill \dot{X}\left(t\right)& =\frac{rC\left(t\right)N\left(t\right)X\left(t\right)}{(K_1+C\left(t\right))(K_2+N\left(t\right))}-X\left(t\right)-m_{1}P\left(t\right)X\left(t\right)\hfill \\ & =U\left(C\right(t),N(t),P(t\left)\right)X\left(t\right)\hfill \\ & >\frac{1}{2}U(1,1,1)X\left(t\right).\hfill \end{array}$$

This implies that $X\left(t\right)\to +\infty$ as $t\to +\infty$, which leads to a contradiction. Thus, $W^s\left(E_0\right)\cap G_0=\varnothing$.

In the following, let us show that $E_0$ is factually isolated. That is, there exists some neighborhood U of $E_0$, such that $E_0$ is the largest invariant set in U. In fact, for sufficiently small positive constant ε, let us choose:

$$U=U\left(E_0\right)=\{(C,N,X,P)\in G|1-C<\epsilon ,1-N<\epsilon ,X<\epsilon ,1-P<\epsilon \}.$$

We show that $E_0$ is the largest invariant set of U for some $\epsilon .$

If not, for any sufficiently small ε, there exists some invariant set $W(W\subset U)$, such that $W\backslash E_0$ is not empty. Let $(C_0,N_0,X_0,P_0)\in W\backslash E_0$ and $\left(C\right(t),N(t),X(t),P(t\left)\right)$ be the solution of Model (4) with the initial function (5). Then, we have that $\left(C\right(t),N(t),X(t),P(t\left)\right)\in U$, $t\in (-\infty ,+\infty )$. From the third equation of Model (4), we obtain that, for $t\in (-\infty ,+\infty )$,

$$\begin{array}{cc}\hfill \dot{X}\left(t\right)& \ge \frac{r(1-\epsilon )(1-\epsilon )X\left(t\right)}{(K_1+1-\epsilon )(K_2+1-\epsilon )}-X\left(t\right)-m_{1}X\left(t\right).\hfill \end{array}$$

Since $r>(K_1+1)(K_2+1)(m_1+1)$, we can choose sufficiently small ε, such that:

$$r(1-\epsilon )(1-\epsilon )>(K_1+1-\epsilon )(K_2+1-\epsilon )(1+m_1).$$

If $X_0>0$, then we have $X\left(t\right)\to +\infty$ $(t\to +\infty )$, which contracts the boundedness of $X\left(t\right)$. Hence, we get $X_0=0$. From the third equation of Model (4), we obtain that $X\left(t\right)\equiv 0$. From the first, second and forth equations of Model (4), for $t\in (-\infty ,+\infty )$, we obtain that $\dot{C}\left(t\right)=1-C\left(t\right)$, $\dot{N}\left(t\right)=1-N\left(t\right)$ and $\dot{P}\left(t\right)=1-P\left(t\right)$. Hence, we must have that $C_0=N_0=P_0=1$. Therefore, $(C_0,N_0,X_0,P_0)=(1,1,0,1)$, which is a contradiction. Thus, we obtain that $E_0$ is factually isolated.

Clearly, $E_0$ is acyclic in $M_{\partial }$. By paper [52], $\partial G_0$ repels uniformly nonnegative solutions of Model (4). It then follows that Model (4) is persistent.

Define $p:X\to R_{+}$ by $p(C,N,X,P)=X,$ $(C,N,X,P)\in G$. Obviously, we obtain that $G_0=p^{-1}$ $(0,+\infty )$ and $\partial G_0=p^{-1}\left(0\right)$. Thus, by ([53], Theorem 3), we have $\underset{t\to \infty }{lim inf}(C\left(t\right),N\left(t\right),X\left(t\right),P\left(t\right))\ge (\eta ,\eta ,\eta ,\eta )$. It then follows that Model (4) is uniformly persistent. The proof of Theorem 7 is completed. □

According to [51,53], we easily obtain the following result.

Corollary 8. 

If $R_0>1$, then Model (4) has one global attractor $A_0$ and also has at least one positive equilibrium $E^{*}(C^{*},N^{*},X^{*},P^{*})\in A_0$.

5.3. The Global Asymptotic Stability of the Positive Equilibrium

In this subsection, the global asymptotic stability of the positive equilibrium $E^{*}(C^{*},N^{*},X^{*},P^{*})$ of Model (4) is studied.

Theorem 9. 

Assume that $R_0>1$ and Model (4) has a unique positive equilibrium $E^{*}(C^{*},N^{*},X^{*},P^{*})$. Let:

$$\begin{array}{cc}\hfill D_1& =r_1(2+m_1{P}^{*})(K_1+C^{*})(K_2+N^{*})-r^2{N}^{*},\hfill \\ \hfill D_2& =r_2(2+m_1{P}^{*})(K_1+C^{*})(K_2+N^{*})-r^2,\hfill \\ \hfill D_3& =b_{13}^2+b_{23}^2+4b_{33},\hfill \\ \hfill D_4& =4b_{33}{b}_{44}+(b_{23}^2+b_{13}^2)b_{44}+(b_{14}^2+b_{24}^2)b_{33}-b_{34}^2+b_{23}{b}_{24}{b}_{34}+b_{13}{b}_{14}{b}_{34}.\hfill \end{array}$$

If $D_1\ge 0$, $D_2\ge 0$, $D_3<0$ and $D_4>0$ hold, then the positive equilibrium $E^{*}(C^{*},N^{*},X^{*},P^{*})$ of Model (4) is globally asymptotically stable. Here:

$$\begin{array}{cc}\hfill b_{13}& =\frac{2r_1+r_1{m}_1{P}^{*}}{r}+\frac{r}{(K_1+1)(K_2+1)(K_1+C^{*})}-\frac{r{N}^{*}}{(K_1+C^{*})(K_2+N^{*})},\hfill \\ \hfill b_{14}& =\frac{r_1{m}_1\alpha }{r},b_{23}=\frac{2r_2+r_2{m}_1{P}^{*}}{r},b_{24}=\frac{r_2{m}_1\alpha }{r},b_{33}=-\frac{(r_1^2+r_2^2)(1+m_1{P}^{*})}{r^2},\hfill \\ \hfill b_{34}& =\frac{(r_1^2+r_2^2)m_1\alpha }{r^2}+m_1+m_2,b_{44}=-(1+m_2{X}^{*}).\hfill \end{array}$$

Proof. 

Let us consider the following Lyapunov function on $A_0$,

$$\begin{array}{cc}\hfill V_2=& \frac{1}{2}{\left(C-C^{*}+\frac{r_1}{r}(X-X^{*})\right)}^2+\frac{1}{2}{\left(N-N^{*}+\frac{r_2}{r}(X-X^{*})\right)}^2\hfill \\ & +X-X^{*}-X^{*}ln\frac{X}{X^{*}}+P-P^{*}-P^{*}ln\frac{P}{P^{*}}.\hfill \end{array}$$

Clearly, $V_2(C^{*},N^{*},X^{*},P^{*})=0$, and $V_2$ is positive definite with respect to $E^{*}(C^{*},N^{*}$ $,X^{*},P^{*})$. The derivative of $V_2$ along the solutions of Model (4) is:

$$\begin{array}{cc}\hfill {\dot{V}}_2=& \left(C\left(t\right)-C^{*}+\frac{r_1}{r}(X\left(t\right)-X^{*})\right)\left(\dot{C}\left(t\right)+\frac{r_1}{r}\dot{X}\left(t\right)\right)\hfill \\ & +\left(N\left(t\right)-N^{*}+\frac{r_2}{r}(X\left(t\right)-X^{*})\right)\left(\dot{N}\left(t\right)+\frac{r_2}{r}\dot{X}\left(t\right)\right)\hfill \\ & +\frac{X\left(t\right)-X^{*}}{X\left(t\right)}\dot{X}\left(t\right)+\frac{P\left(t\right)-P^{*}}{P\left(t\right)}\dot{P}\left(t\right)\hfill \\ \hfill =& \left(C\left(t\right)-C^{*}+\frac{r_1}{r}(X\left(t\right)-X^{*})\right)\left(1-C\left(t\right)-\frac{r_1}{r}X\left(t\right)-\frac{r_1{m}_{1}X\left(t\right)P\left(t\right)}{r}\right)\hfill \\ & +\left(N\left(t\right)-N^{*}+\frac{r_2}{r}(X\left(t\right)-X^{*})\right)\left(1-N\left(t\right)-\frac{r_2}{r}X\left(t\right)-\frac{r_2{m}_{1}X\left(t\right)P\left(t\right)}{r}\right)\hfill \\ & +(X\left(t\right)-X^{*})\left(\frac{rC\left(t\right)N\left(t\right)}{(K_1+C\left(t\right))(K_2+N\left(t\right))}-1-m_{1}P\left(t\right)\right)\hfill \\ & +\frac{P\left(t\right)-P^{*}}{P\left(t\right)}(1-P\left(t\right)-m_{2}X\left(t\right)P\left(t\right))\hfill \end{array}$$

for any $t\ge 0$.

Noting that:

$$\left\{\begin{array}{c}{\displaystyle 1-C^{*}-\frac{r_1{C}^{*}{N}^{*}{X}^{*}}{(K_1+C^{*})(K_2+N^{*})}=0,}\hfill \\ {\displaystyle 1-N^{*}-\frac{r_2{C}^{*}{N}^{*}{X}^{*}}{(K_1+C^{*})(K_2+N^{*})}=0,}\hfill \\ {\displaystyle \frac{r{C}^{*}{N}^{*}{X}^{*}}{(K_1+C^{*})(K_2+N^{*})}-X^{*}-m_1{X}^{*}{P}^{*}=0,}\hfill \\ {\displaystyle 1-P^{*}-m_2{X}^{*}{P}^{*}=0.}\hfill \end{array}\right.$$

Thus, we have that:

$$\begin{array}{cc}\hfill {\dot{V}}_2=& \left(C\left(t\right)-C^{*}+\frac{r_1}{r}(X\left(t\right)-X^{*})\right)\left(C^{*}-C\left(t\right)+\frac{r_1}{r}(X^{*}-X\left(t\right))+\frac{r_1{m}_1}{r}(X^{*}{P}^{*}-X\left(t\right)P\left(t\right))\right)\hfill \\ & +\left(N\left(t\right)-N^{*}+\frac{r_2}{r}(X\left(t\right)-X^{*})\right)\left(N^{*}-N\left(t\right)+\frac{r_2}{r}(X^{*}-X\left(t\right))+\frac{r_2{m}_1}{r}(X^{*}{P}^{*}-X\left(t\right)P\left(t\right))\right)\hfill \\ & +(X\left(t\right)-X^{*})\left(\frac{rC\left(t\right)N\left(t\right)}{(K_1+C\left(t\right))(K_2+N\left(t\right))}-\frac{r{C}^{*}{N}^{*}}{(K_1+C^{*})(K_2+N^{*})}\right)\hfill \\ & +m_1(X\left(t\right)-X^{*})(P^{*}-P\left(t\right))+\frac{P\left(t\right)-P^{*}}{P\left(t\right)}\left(P^{*}-P\left(t\right)+m_2(X^{*}{P}^{*}-X\left(t\right)P\left(t\right))\right)\hfill \end{array}$$

for any $t\ge 0$.

For any $t\ge 0$, we have that:

$$\begin{array}{c}\frac{rC\left(t\right)N\left(t\right)}{(K_1+C\left(t\right))(K_2+N\left(t\right))}-\frac{r{C}^{*}{N}^{*}}{(K_1+C^{*})(K_2+N^{*})}\hfill \\ =-\frac{r{C}^{*}{N}^{*}}{(K_1+C^{*})(K_2+N^{*})}+\frac{rC\left(t\right)N^{*}}{(K_1+C^{*})(K_2+N^{*})}-\frac{rC\left(t\right)N^{*}}{(K_1+C^{*})(K_2+N^{*})}\hfill \\ +\frac{rC\left(t\right)N\left(t\right)}{(K_1+C^{*})(K_2+N^{*})}-\frac{rC\left(t\right)N\left(t\right)}{(K_1+C^{*})(K_2+N^{*})}+\frac{rC\left(t\right)N\left(t\right)}{(K_1+C^{*})(K_2+N\left(t\right))}\hfill \\ -\frac{rC\left(t\right)N\left(t\right)}{(K_1+C^{*})(K_2+N\left(t\right))}+\frac{rC\left(t\right)N\left(t\right)}{(K_1+C\left(t\right))(K_2+N\left(t\right))}.\hfill \end{array}$$

Hence, it can be obtained that, for any $t\ge 0$,

$$\begin{array}{cc}\hfill {\dot{V}}_2=& \left(C\left(t\right)-C^{*}+\frac{r_1}{r}(X\left(t\right)-X^{*})\right)(C^{*}-C\left(t\right)+\frac{r_1}{r}(X^{*}-X\left(t\right))\hfill \\ & +\frac{r_1{m}_1{P}^{*}}{r}(X^{*}-X\left(t\right))+\frac{r_1{m}_{1}X\left(t\right)}{r}(P^{*}-P\left(t\right)))\hfill \\ & +\left(N\left(t\right)-N^{*}+\frac{r_2}{r}(X\left(t\right)-X^{*})\right)(N^{*}-N\left(t\right)+\frac{r_2}{r}(X^{*}-X\left(t\right))\hfill \\ & +\frac{r_2{m}_1{P}^{*}}{r}(X^{*}-X\left(t\right))+\frac{r_2{m}_{1}X\left(t\right)}{r}(P^{*}-P\left(t\right)))\hfill \\ & +(X\left(t\right)-X^{*})(-\frac{r{C}^{*}{N}^{*}}{(K_1+C^{*})(K_2+N^{*})}+\frac{rC\left(t\right)N^{*}}{(K_1+C^{*})(K_2+N^{*})}\hfill \\ & -\frac{rC\left(t\right)N^{*}}{(K_1+C^{*})(K_2+N^{*})}+\frac{rC\left(t\right)N\left(t\right)}{(K_1+C^{*})(K_2+N^{*})}-\frac{rC\left(t\right)N\left(t\right)}{(K_1+C^{*})(K_2+N^{*})}\hfill \\ & +\frac{rC\left(t\right)N\left(t\right)}{(K_1+C^{*})(K_2+N\left(t\right))}-\frac{rC\left(t\right)N\left(t\right)}{(K_1+C^{*})(K_2+N\left(t\right))}+\frac{rC\left(t\right)N\left(t\right)}{(K_1+C\left(t\right))(K_2+N\left(t\right))})\hfill \\ & -\frac{{(P\left(t\right)-P^{*})}^2}{P\left(t\right)}+\frac{m_2(P\left(t\right)-P^{*})}{P\left(t\right)}\left(X^{*}(P^{*}-P\left(t\right))+P\left(t\right)(X^{*}-X\left(t\right))\right)\hfill \\ & +m_1(X\left(t\right)-X^{*})(P^{*}-P\left(t\right)).\hfill \end{array}$$

Motivated by the papers mentioned in [50,54], the derivative of $V_2$ along the solutions of Model (4) is represented in quadratic form. Then, we obtain that:

$$\begin{array}{cc}\hfill {\dot{V}}_2=& -{(C\left(t\right)-C^{*})}^2-{(N\left(t\right)-N^{*})}^2+(\frac{2r_1+r_1{m}_1{P}^{*}}{r}+\frac{rC\left(t\right)N\left(t\right)}{(K_1+C\left(t\right))(K_2+N\left(t\right))(K_1+C^{*})}\hfill \\ & -\frac{r{N}^{*}}{(K_1+C^{*})(K_2+N^{*})})(C\left(t\right)-C^{*})(X^{*}-X\left(t\right))+(\frac{2r_2+r_2{m}_1{P}^{*}}{r}\hfill \\ & +\frac{rC\left(t\right)N\left(t\right)}{(K_2+N\left(t\right))(K_2+N^{*})(K_1+C^{*})}-\frac{rC\left(t\right)}{(K_1+C^{*})(K_2+N^{*})})(N\left(t\right)-N^{*})(X^{*}-X\left(t\right))\hfill \\ & -\frac{1+m_2{X}^{*}}{P\left(t\right)}{(P\left(t\right)-P^{*})}^2+\frac{r_1{m}_{1}X}{r}(C\left(t\right)-C^{*})(P^{*}-P\left(t\right))\hfill \\ & -\frac{(r_1^2+r_2^2)(1+m_1{P}^{*})}{r^2}{(X\left(t\right)-X^{*})}^2+\left(\frac{(r_1^2+r_2^2)m_{1}X\left(t\right)}{r^2}+m_1+m_2\right)(X^{*}-X\left(t\right))(P^{*}-P\left(t\right))\hfill \\ & +\frac{r_2{m}_{1}X\left(t\right)}{r}(N\left(t\right)-N^{*})(P^{*}-P\left(t\right))\hfill \end{array}$$

for any $t\ge 0$.

If the condition $D_1\ge 0$ holds, we obtain that, for any $t\ge 0$,

$$\begin{array}{c}\hfill \frac{2r_1+r_1{m}_1{P}^{*}}{r}+\frac{rC\left(t\right)N\left(t\right)}{(K_1+C\left(t\right))(K_2+N\left(t\right))(K_1+C^{*})}-\frac{r{N}^{*}}{(K_1+C^{*})(K_2+N^{*})}\ge 0.\end{array}$$

Furthermore, from the condition $D_2\ge 0$ and $C\left(t\right)\le 1$ for any $t\ge 0$, we obtain that:

$$\begin{array}{cc}\hfill \frac{2r_2+r_2{m}_1{P}^{*}}{r}+\frac{rC\left(t\right)N\left(t\right)}{(K_2+N\left(t\right))(K_2+N^{*})(K_1+C^{*})}-\frac{rC\left(t\right)}{(K_1+C^{*})(K_2+N^{*})}& \ge 0.\hfill \end{array}$$

These conditions ensure that the coefficients of $(C\left(t\right)-C^{*})(X\left(t\right)-X^{*})$ and $(N\left(t\right)-N^{*})(X\left(t\right)-X^{*})$ are nonnegative. Thus, we obtain that, for any $t\ge 0$,

$$\begin{array}{cc}\hfill {\dot{V}}_2\le & -{(C\left(t\right)-C^{*})}^2-{(N\left(t\right)-N^{*})}^2+(\frac{2r_1+r_1{m}_1{P}^{*}}{r}+\frac{r}{(K_1+1)(K_2+1)(K_1+C^{*})}\hfill \\ & -\frac{r{N}^{*}}{(K_1+C^{*})(K_2+N^{*})})\left|C\left(t\right)-C^{*}\right|\left|X^{*}-X\left(t\right)\right|\hfill \\ & +\left(\frac{2r_2+r_2{m}_1{P}^{*}}{r}\right)\left|N\left(t\right)-N^{*}\right|\left|X^{*}-X\left(t\right)\right|-(1+m_2{X}^{*}){(P\left(t\right)-P^{*})}^2\hfill \\ & +\frac{r_1{m}_1\alpha }{r}\left|C\left(t\right)-C^{*}\right|\left|P^{*}-P\left(t\right)\right|-\frac{(r_1^2+r_2^2)(1+m_1{P}^{*})}{r^2}{(X\left(t\right)-X^{*})}^2\hfill \\ & +\left(\frac{(r_1^2+r_2^2)m_1\alpha }{r^2}+m_1+m_2\right)\left|X^{*}-X\left(t\right)\right|\left|P^{*}-P\left(t\right)\right|+\frac{r_2{m}_1\alpha }{r}\left|N\left(t\right)-N^{*}\right|\left|P^{*}-P\left(t\right)\right|\hfill \\ \hfill =& Y^T\left(t\right)WY\left(t\right),\hfill \end{array}$$

where $Y\left(t\right)=(\left|C\left(t\right)-C^{*}\right|,\left|N\left(t\right)-N^{*}\right|,\left|X\left(t\right)-X^{*}\right|,\left|P\left(t\right)-P^{*}\right|)$. W is a symmetric $4\times 4$ matrix and $W={\left\{b_{ij}\right\}}_{1\le i,j\le 4}$ with:

$$W=\left(\begin{array}{cccc}{b}_{11}& {\displaystyle \frac{1}{2}{b}_{12}}& {\displaystyle \frac{1}{2}{b}_{13}}& {\displaystyle \frac{1}{2}{b}_{14}}\\ {\displaystyle \frac{1}{2}{b}_{12}}& b_{22}& {\displaystyle \frac{1}{2}{b}_{23}}& {\displaystyle \frac{1}{2}{b}_{24}}\\ {\displaystyle \frac{1}{2}{b}_{13}}& {\displaystyle \frac{1}{2}{b}_{23}}& b_{33}& {\displaystyle \frac{1}{2}{b}_{34}}\\ {\displaystyle \frac{1}{2}{b}_{14}}& {\displaystyle \frac{1}{2}{b}_{24}}& {\displaystyle \frac{1}{2}{b}_{34}}& b_{44}\end{array}\right).$$

Here, $b_{11}=-1,b_{12}=0,b_{22}=-1.$ The other parameters have the same definitions as in Theorem 9. Then, W is negative definite. By Lemma 6.2 provided in [50], we can obtain that the real quadratic form $Y^T\left(t\right)WY\left(t\right)$ is negative definite. Then, by using some classical analysis techniques of differential equations, the positive equilibrium $E^{*}(C^{*},N^{*},X^{*},P^{*})$ of Model (4) is globally asymptotically stable. □

Remark 3. 

Obviously, it is easy to obtain that $b_{11}=-1<0,$ and:

$$\begin{array}{c}\hfill \left|\begin{array}{cc}{b}_{11}& {\displaystyle \frac{1}{2}{b}_{12}}\\ {\displaystyle \frac{1}{2}{b}_{12}}& b_{22}\end{array}\right|=1>0.\end{array}$$

If the conditions $D_3<0$ and $D_4>0$ hold, we obtain that:

$$\begin{array}{c}\hfill \left|\begin{array}{ccc}{\displaystyle b_{11}}& {\displaystyle \frac{1}{2}{b}_{12}}& {\displaystyle \frac{1}{2}{b}_{13}}\\ {\displaystyle \frac{1}{2}{b}_{12}}& b_{22}& {\displaystyle \frac{1}{2}{b}_{23}}\\ {\displaystyle \frac{1}{2}{b}_{13}}& {\displaystyle \frac{1}{2}{b}_{23}}& b_{33}\end{array}\right|<0,\end{array}$$

and:

$$\begin{array}{c}\hfill \left|\begin{array}{cccc}{\displaystyle b_{11}}& {\displaystyle \frac{1}{2}{b}_{12}}& {\displaystyle \frac{1}{2}{b}_{13}}& {\displaystyle \frac{1}{2}{b}_{14}}\\ {\displaystyle \frac{1}{2}{b}_{12}}& b_{22}& {\displaystyle \frac{1}{2}{b}_{23}}& {\displaystyle \frac{1}{2}{b}_{24}}\\ {\displaystyle \frac{1}{2}{b}_{13}}& {\displaystyle \frac{1}{2}{b}_{23}}& b_{33}& {\displaystyle \frac{1}{2}{b}_{34}}\\ {\displaystyle \frac{1}{2}{b}_{14}}& {\displaystyle \frac{1}{2}{b}_{24}}& {\displaystyle \frac{1}{2}{b}_{34}}& b_{44}\end{array}\right|>0.\end{array}$$

Next, we illustrate that the conditions given in Theorem 9 are reasonable. We set $r=5.8,K_1=1,K_2=0.6,m_1=0.0001,m_2=1,r_1=6,r_2=10.$ It is easy to see that there exists a unique positive equilibrium $E^{*}(C^{*},N^{*},X^{*},$ $P^{*})\approx (0.6709,0.4515,0.3181,0.7587)$, and the conditions $D_1$–$D_4$ hold. By Theorem 9, the positive equilibrium of Model (4) is globally asymptotically stable (see Figure 3b).

In general, it is very difficult to obtain the global stability properties of the positive equilibrium of Model (4). In this paper, we have obtained the sufficient conditions to ensure the global stability properties of the positive equilibrium by constructing a suitable Lyapunov function. Because of the complexity of the conditions, it is difficult to grasp biological intuition. However, from the numerical simulations, we have found some interesting biological phenomena. The conditions $D_1$–$D_4$ can be satisfied if the flocculation rate of microorganisms ($m_1$) is small enough (see, Figure 3b). From a biological point of view, these conditions are reasonable. For more detailed biological considerations, we will leave it for further investigation.

6. Control Strategies

In this section, some control strategies are provided by suitable theoretical analysis. If $R_0<1$ holds, then, we have from Theorem 5 that the boundary equilibrium $E_0$ of Model (4) is locally asymptotically stable. We note that the condition $R_0<1$ is equivalent to the following inequality,

$$\begin{array}{c}\hfill R_0=\frac{r}{(K_1+1)(K_2+1)(m_1+1)}<1.\end{array}$$

(12)

All of the parameters r, $K_1$, $K_2$ and $m_1$ in (12) are defined as in Model (4). (12) can be further written as the following form,

$$\begin{array}{c}\hfill \frac{r}{(\frac{K_1}{C^0}+1)(\frac{K_2}{N^0}+1)(m_1{P}^0+D)}<1.\end{array}$$

(13)

Here, all of the parameters r, $K_1$, $K_2$, $m_1$, $C^0$, $N^0$, $P^0$ and D in (13) are defined as in Model (3).

In view of the biological meanings of the parameters in Model (3) and Condition $\left(13\right)$, Theorem 5 indicates that the concentration of $Rhodopseudomonaspalustris$ in the chemostat tends to zero, and the concentration of the carbon source, nitrogen source and flocculants may tend to the constant values $C^0$, $N^0$ and $P^0$, respectively, as time t increases, if one of the following two cases occurs: (a) reducing the absorption of $Rhodopseudomonaspalustris$ or the carbon input concentration, or the nitrogen input concentration; (b) improving the velocity or the flocculation effect or flocculant input concentration. These cases are reasonable, since they imply the insufficient sources for $Rhodopseudomonaspalustris$ to grow. Hence, in the environmental science field, it can be used to remove algae and heavy metals.

From Theorem 7, we have that Model (3) is uniformly persistent if $R_0>1$. This means that the concentrations of the carbon source, nitrogen source, $Rhodopseudomonas$ $palustris$ and flocculants in the chemostat may be ultimately maintained at some positive constant values, as time t increases, if one of the following two cases occurs: (a) improving the absorption of $Rhodopseudomonas$ $palustris$, or the carbon input concentration, or the nitrogen input concentration; (b) reducing the velocity or the flocculation effect or flocculant input concentration.

These control strategies can be performed by numerical simulations.

In the following, for convenience, we simulate the extinction or persistence of microorganism ($X\left(t\right)$) numerically by using (12) and Model (4).

If the parameters are chosen as in

Table A2. Parameter values used in the simulations of control strategies shown in Figure 3, Figure 4, Figure 5 and Figure 6.

Description	Parameter	Value
the growth rate of $Rhodopseudomonaspalustris$	r	5.8
the quantity of decreasing of carbon source	$r_1$	6
the quantity of decreasing of nitrogen source	$r_2$	10
the flocculation effect	$m_1$	1
the flocculation ratio	$m_2$	1
Michaelis–Menten constant of carbon	$K_1$	1
Michaelis–Menten constant of nitrogen	$K_2$	0.6

(see Appendix A), $Rhodopseudomonas$ $palustris$ in the chemostat will tend to extinction (see Figure 4a).

If the absorption of $Rhodopseudomonas$ $palustris$ (r) is improved from $r=5.8$ to $r=8$ and the other parameters are the same as Table A2, $Rhodopseudomonas$ $palustris$ in the chemostat will tend to be constant (see Figure 4b).

If the flocculation effect ($m_1$) is reduced from $m_1=1$ to $m_1=0.1$, $Rhodopseudomonas$ $palustris$ in the chemostat will tend to be constant (see Figure 5a).

If Michaelis–Menten constant of carbon ($K_1$) is reduced from $K_1=1$ to $K_1=0.5$, $Rhodopseudomonas$ $palustris$ in the chemostat will tend to be constant (see Figure 5b).

If Michaelis–Menten constant of nitrogen ($K_2$) is reduced from $K_2=1$ to $K_2=0.25$, $Rhodopseudomonas$ $palustris$ in the chemostat will tend to be constant (see Figure 6).

7. Discussion and Conclusions

In the paper, based on some biological considerations and chemostat models, a dynamic model governed by ordinary differential equations with four variables (carbon source, nitrogen source, $Rhodopseudomonas$ $palustris$ and flocculants) is presented. There is a boundary equilibrium and at most five positive equilibria for the proposed model. To give a theoretical analysis for the existence of all of the positive equilibria of Model (4), the method of the Descartes rule of signs is applied to the classifications of the positive roots of a fifth order algebraic equation.

The local and global stability properties of the boundary equilibrium of Model (4) have been studied in detail. An interesting phenomenon of backward and forward bifurcations is observed. That is, there may exist two positive equilibria even if the condition $R_0<1$ holds. Hence, sufficient Condition (10) to ensure the global stability of the boundary equilibrium is reasonable in mathematics.

The local stability of the positive equilibrium of Model (4) is also carried out. From Condition (11), we have that the positive equilibrium is locally asymptotically stable when the flocculation coefficient $m_1$ is small enough. Hence, Condition (11) is also reasonable in biology.

Uniform persistence of Model (4) has also been completely studied under the condition $R_0>1$. Uniform persistence has very important significance both in mathematics and biology, and it characterizes the long-term survival of some microorganisms [45].

Finally, some control strategies are provided by simple theoretical analysis. From Theorem 5, we have that $Rhodopseudomonas$ $palustris$ in the chemostat will tend to extinction if $R_0<1$. In this case, these control strategies can be applied to remove $Cyanobacteria$, which are well known to produce a variety of toxins and have serious harm on human health. From Theorem 7, we have that $Rhodopseudomonas$ $palustris$ in the chemostat will tend to be positive constant if $R_0>1$. In this case, these control strategies can be widely used for the collection of useful microorganisms.

It is well-known that the existence of time delays is inevitable in biology. For example, in the cultivation of microorganisms, there are always time delays in the process of transferring nutrients and the uptake of nutrients. Hence, chemostat models with time delays that account for the time lapsing between the uptake of nutrients by cells and the incorporation of these nutrients as biomass have been given much attention [40,55,56,57]. Based on Model (3), it may have the following more general form with time delays,

$$\left\{\begin{array}{c}{\displaystyle \frac{dC\left(t\right)}{dt}=(C^0-C\left(t\right))D-\frac{r_{1}C\left(t\right)N\left(t\right)X\left(t\right)}{{\delta }_1(K_1+C\left(t\right))(K_2+N\left(t\right))}+{\rho }_{1}X(t-\sigma ),}\hfill \\ {\displaystyle \frac{dN\left(t\right)}{dt}=(N^0-N\left(t\right))D-\frac{r_{2}C\left(t\right)N\left(t\right)X\left(t\right)}{{\delta }_2(K_1+C\left(t\right))(K_2+N\left(t\right))}+{\rho }_{2}X(t-\sigma ),}\hfill \\ {\displaystyle \frac{dX\left(t\right)}{dt}=\frac{r{e}^{-d_1\tau }C(t-\tau )N(t-\tau )X(t-\tau )}{(K_1+C(t-\tau ))(K_2+N(t-\tau ))}-DX\left(t\right)-d_{1}X\left(t\right)-m_{1}X\left(t\right)P\left(t\right),}\hfill \\ {\displaystyle \frac{dP\left(t\right)}{dt}=(P^0-P\left(t\right))D-m_{2}X\left(t\right)P\left(t\right).}\hfill \end{array}\right.$$

(14)

In Model (14), the constants ${\rho }_1\ge 0$ and ${\rho }_2\ge 0$ are the rate constants at which the carbon source and nitrogen source are recycled because of the death of $Rhodopseudomonas$ $palustris$. The constant $\sigma \ge 0$ is a fixed time during which the carbon source and nitrogen source are released completely from dead $Rhodopseudomonas$ $palustris$. The constant $\tau \ge 0$ denotes the time delay involved in the conversion of nutrients to $Rhodopseudomonas$ $palustris$. The factor $e^{-d_1\tau }$ is the probability constant at which $Rhodopseudomonas$ $palustris$ remains in the culture vessel during the conversion process. The theoretical analysis of Model (14) will be studied separately.