# Dynamical Analysis of a Continuous Stirred-Tank Reactor with the Formation of Biofilms for Wastewater Treatment

Authors: Karen López Buriticá; Simeón Casanova Trujillo; Carlos D. Acosta; Héctor A. Granada Diaz
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

This paper analyzes the dynamics of a system that models the formation of biofilms in a continuous stirred-tank reactor (CSTR) when it is utilized for wastewater treatment. The growth rate of the microorganisms is modeled using two different kinetics, Monod and Haldane kinetics, with the goal of studying the influence of each in the system. The equilibrium points are identified through a stability analysis, and the bifurcations found are characterized.


---

## Body

## 1. Introduction

Water is the most important natural resource without which life could not exist. It is the most abundant liquid on Earth, covering 71% of the Earth’s surface, of which only 3% is fresh water. Water use results in a decrease in water quality, and serious environmental deterioration results from directly returning water to the environment; therefore, in recent years, the construction of increasingly more efficient treatment plants has gained vital importance. Consequently, the role of biological processing in wastewater treatment has increased considerably where biofilm reactors are one of the most interesting [1]. Biofilms are thin layers of microorganisms that adhere to a solid surface; they can develop on almost any type of surface exposed to an aqueous medium, and they are utilized in wastewater treatment to eliminate and oxidize organic and inorganic components. The biomass concentration in a biofilm can be ten times larger than its concentration in a liquid culture [2]. This reduces the volume of the equipment for a constant rate of elimination per unit volume.A theoretical analysis of the formation of biofilms was performed in [2], and in [3, 4] characterization and analysis of biofilm use in wastewater treatment were conducted. Experimental results on a pilot scale have demonstrated the efficiency of biofilm reactors for treating wastewater. For example, (1) that contained molasses as a source of carbon in different conditions of the influent [1]; (2) with a vertically moving biofilm reactor (VMBR) followed by a stratified sand filter [5]; and (3) in combination with treatment by activated sludge [6].The implementation of a pilot plant to study these procedures can be costly, which makes mathematical models and computer simulations essential to describe, predict, and control the complex interactions present in the reactor [7]. For several years, researchers have focused on developing mathematical methods that more realistically simulate the behavior of a biofilm reactor. For example, in [8] a biofilm reactor model was proposed based on physiological aspects. By time, the models began to involve more variables [9–11]. “However, these often turned out to be unsatisfactory due to certain intrinsic deficiencies. For example, in [8] the model only took account of the accumulation of active biomass and inactive biomass of biofilms and in [9–11] the model was mostly used to study the morphology and structure of biofilms, so mechanisms and kinetics processes of biofilm formation have not yet been well understood” [12].In [12] a model is proposed “which attempts to summarize the factors influencing biofilm formation and eliminating the deficiency of existing models and to study the kinetic mechanisms from a new perspective.” This model describes the formation of microbial populations in an aqueous medium inside a reactor and the mechanism under which bacteria can be suspended in the water or colonize the surface. The model in [12] considers extremely thin biofilms so it can be considered that all the microorganisms receive the substrate at the same time, therefore the model incorporates this condition by ignoring diffusion reactions. There are many applications where thin biofilms are desirable and detachment procedure is considered to keep such a condition [13–15].In this paper, the model in [12] is also considered. However, our bifurcational and dynamical analysis included not only the Monod kinetics but also the Haldane kinetics. Some bifurcations are found, including a transcritical bifurcation that does not observe the conditions in Sotomayor theorem.In [16] the model proposed in [12] is reformulated to study the behavior of the suspended microorganisms.


## 2. Materials and Methods

### 2.1. Mathematical Model

The continuous stirred-tank reactor described in [12] was chosen as the model reactor. The system consists of three ordinary nonlinear differential equations and models the formation of heterotrophic aerobic biofilms inside the reactor, taking into account both the microorganisms adhered to the biofilms and the microorganisms that are suspended. In addition to studying the influence of the various parameters in biofilm formation, in [12] the system is established and modeled using Monod kinetics, the equilibrium corresponding to the washing condition. In the present work, the system is also analyzed as modeled by Haldane kinetics and more equilibrium points are found and their stability is studied, in addition to characterizing the bifurcations that they represent. As bifurcation parameter of the system, the dilution coefficient D∈[0,1] was chosen, due to this parameter allows to control the residence time of wastewater in the reactor.The following assumptions are made with respect to the dynamics of the formation of biofilms. The three-dimensional structure of the biofilms is ignored because the biofilms are considered to be infinitely thin. There are a finite number of colonization sites available on the support surface, as well as maximum possible surface density of adhered microorganisms. The adhered cells are separated outside the fluid at a rate proportional to its density and the daughter cells of adhered microorganisms compete for space on the support surface—a fractionG of the daughter cells finds adhesion sites and a fraction 1-G does not. The function G decreases with the number of daughter cells (as the number of daughter cells increases, the number of available sites to occupy decreases), where G=(1-W)/(1.1-W), which is a function of probability [12].The reactor dynamics are described by the following equations:(1)dSdt=DS0-S-γ-1μSXu+Xw,dXudt=XuμS-D-k+βXw+XwμS1-GW-αXu1-W,dXwdt=XwμSGW-β-k+αXu1-W.The following parameters are used for the numerical simulation: (2)γ=0.5,k=0.01,Xw,m=5000,D=1,β=0.03,α=0.1,S0=500,μm=0.12,KS=80,Xu0=1000.These values are taken from [12].The growth of the suspended and adhered microorganisms in the reactor is a nonlinear natural phenomenon that can be analyzed using nonlinear dynamics, which provides a global perspective of the types of behavior in the system, such as equilibrium states, stability, and bifurcations.


### 2.2. Monod Kinetics

The suspended and adhered microorganisms inside the reactor consume the organic matter present in the wastewater to be treated. The specific rate of substrate consumption is the “speed” at which the organism consumes the substrate. Although the growth of microorganisms is a complex phenomenon, there are equations that can model this behavior. Monod’s equation is the simplest and most widely utilized equation for describing the kinetics of microbial growth. It is a function of the limiting substrate and is expressed asμ(S)=μmS/KS+S where μm is the maximum growth rate. KS is the concentration of the substrate corresponding to the half of the maximum growth rate of (S) and represents the affinity of the microorganisms for the substrate. If the organism has a great affinity for the limiting substrate, the value of KS is low [17].


#### 2.2.1. Equilibrium Points

A pointx0∈Rn is an equilibrium point or a critical point of x˙=f(x) if f(x0)=0. In this point, the system remains in a stationary state; that is, if ϕt:U→Rn is the flux of the system, then ϕt(x0)=x0∀t∈R [18, 19].The system has an equilibrium in the washing condition for any set of values of the parameters, in this equilibrium there aren’t microorganisms adhere or suspended, but there is substrate into the reactor. Exists other equilibrium physically feasible when the microorganisms adhere to the support surface and finally two equilibria in the case where suspended microorganisms do not adhere to the support surface.(i)
First equilibrium point, corresponding to the washing condition, isE1=(S0,0,0).(ii)
Ifα≠0, that is, there adherence of microorganisms to biofilms, then a nontrivial equilibrium point exists: E2=(S∗,Xu∗,Xw∗)≠(0,0,0). This equilibrium point cannot be found analytically.(iii)
Ifα=0, there exist two equilibrium points:(3)E3=S,Xu,Xw=D+kKSμm-D-k,γDS0-SD+k,0withμ(S)-D-k=0.This point is physically possible ifμm-D-k>0 and S0-S>0; that is, (4)D<S0μm-k-kKSS0+KS=0.093448.E4=(S,Xu,Xw), where (5)Xu=γS0-SμS-kμSXw=1.1β+k-μSβ+k-μSXwm.SubstitutingXu and Xw in the following equation, one can solve S. Consider(6)DSo-S-γ-1μSXu+Xw=0.


#### 2.2.2. Equilibrium Point Stability

An equilibrium point of a dynamical system is stable in Lyapunov sense if all trajectories with initial conditions near the equilibrium point remain in that vicinity. The equilibrium point is asymptotically stable, that is, the trajectories always tend toward the equilibrium point. Lyapunov’s indirect method is outlined in [20] and demonstrated in [21], and it provides a procedure to determine the local stability of a hyperbolic equilibrium point. Furthermore, the Hartman-Grobman theorem states that, near a hyperbolic equilibrium point, the nonlinear system presents a behavior qualitatively equivalent to that of the corresponding linear system [22, 23], where a hyperbolic equilibrium point is an equilibrium point in which its Jacobian matrix does not have eigenvalues with zero real part.In the following, the stability of the equilibrium points is analyzed by Lyapunov’s indirect method, linearizing the system and calculating its eigenvalues. To analyze the location in the plane of the eigenvalues, the Routh-Hurwitz criterion was utilized.At the equilibrium pointE1, the Jacobian matrix of the system is (7)-D-γ-1μS0-γ-1μS00μS0-D-k-αμS01-G0+β0αμS0G0-k-β.Let A be a submatrix, such that (8)A=μS0-D-k-αμS01-G0+βαμS0G0-k-β.The eigenvalues λ1, λ2, and λ3 of the matrix JE1 are the roots of the equation: (9)λ2-δλ+τλ+D=0,δ=trA,τ=detA,λ1,2=δ±δ2-4τ2,λ3=-D.It is readily apparent that (10)maxμ-D-k-α,μG-k-β<λ1.Substituting the values of the parameters, one obtains that μ(S0)G(0)-k-β>0. Therefore, E1 is unstable.For the equilibrium pointE2, numerically performing the stability analysis, we find that it is always stable.The algorithm is as follows: for a given value ofD the equilibrium point is found using the Newton-Raphson method, and the equilibrium curve is constructed using the same method taking the previous equilibrium as the initial condition. Every equilibrium point is evaluated in the Jacobian matrix, and the eigenvalues are calculated. If they are all negative, the equilibrium is stable. If either of them is negative, the equilibrium is unstable and if at least one has a zero real part, it is not possible to decide on the stability of the equilibrium.At the equilibrium pointE3, the Jacobian matrix of the system is (11)-D-γ-1Xuμ′-γ-1μ-γ-1μXuμ′0β+μ1-G000μG0-β-k=a11a12a13a210a2300a33and its characteristic polynomial is (12)λ3-a11+a33λ2+a11a33-a12a21λ+a12a21a33=0a11,a12<0,a21>0.We establish the following array: (13)λ31a11a33-a12a21λ2-a11+a33a12a21a33λ1a11a33-a11a12a21a11+a330λ0a12a21a33The Routh-Hurwitz stability criterion states that the number of roots of the characteristic polynomial with a positive real part is equal to the number of sign changes of the coefficients of the first column of the array.a
11
+
a
33 is negative. Indeed a11+a33=(DG(0)-D)+(kG(0)-k)-γ-1Xuμ′(S)-β<0 given that DG(0)-D<0 and kG(0)-k<0.Ifa33<0 then a11a33>0 and a11a33-a11a12a21/a11+a33>0.Ifa33<0 then D<β+k-kG(0)/G(0).Therefore, forD<0.034, a12a21a33>0.Ifa33>0 then D>β+k-kG(0)/G(0).Therefore, forD>0.034, a12a21a33<0.In summary, we have(14)D<0.0340.034<D<0.093448λ3+λ2+λ1+λ0+λ3+λ2+λ1xλ0-If D<0.034, there are no sign changes, and, thus, the real part of all the eigenvalues is negative, and the equilibrium point is stable. If 0.034<D<0.093448 there is always a sign change (x is unknown and it does not matter, because anyway there is a sign change), which means that there is an eigenvalue with a positive real part, and, thus, the equilibrium point is unstable.For the equilibrium pointE4, the characteristic polynomial is λ3+bλ2+cλ+d=0, where (15)
b
=
-
a
11
+
a
22
+
a
33
,
c
=
-
a
12
a
21
+
a
32
a
23
+
a
13
a
31
-
a
11
a
33
-
a
22
a
33
-
a
11
a
22
,
d
=
-
a
11
a
22
a
33
-
a
12
a
23
a
31
-
a
13
a
21
a
32
+
a
12
a
21
a
33
+
a
32
a
23
a
11
+
a
13
a
31
a
22
.By virtue of the Routh-Hurwitz criterion, one has (16)λ31cλ2bdλ1c1=bc-db0λ0dNumerically, it is found that, for values D<0.034, there is a sign change, and, thus, the equilibrium point is unstable, whereas for values D>0.034, there is no sign change, and, thus, the equilibrium point is stable.


#### 2.2.3. One-Parameter Bifurcations

The fact that the dynamics of a system can change drastically as one or more of its parameters are varied is well known. This qualitative or structural change is known as abifurcation. In general, bifurcation theory studies the structural changes that a dynamical system can undergo as its parameters are varied; further bifurcations only occur in nonhyperbolic equilibria [19, 24].In what follows, the transcritical bifurcation described in [25] is presented because it is precisely the one studied in this work.Letx=x0 be a nonhyperbolic equilibrium point. The structure of the system is characterized as follows.(i)
Two curves of equilibrium pointsϕ(t) and ψ(t) pass through x0.(ii)
Both curves exist on both sidesx0.(iii)
The stability of the curves is interchanged as they pass throughx0.The system presents a transcritical bifurcation in the equilibrium curvesE3 and E4 in fact.Whenα=0, for the point E3, an eigenvalue λ=0 exists if a12a21a33=0, given that a12<0 and a21>0; then a33 must be zero. That is, (17)μSG0-β-k=0withμS=D+k.Thus(18)D=β+k-kG0G0=β+k1-G0G0=0.03+0.011-10/1110/11=0.034.The equilibrium point in the bifurcation parameter is as follows: S,Xu,Xw=46.31,175.28,0.The equilibrium curveE3 is stable for D<0.034 and unstable for D>0.034. For E4 the equilibria are unstable for D<0.034 and stable for D>0.034. Figures 1 and 2 show how the equilibrium curves of the points E3 and E4 cross each other and interchange stability, such that a transcritical bifurcation appears.Figure 1
Transcritical bifurcation varyingD. Stable: blue. Unstable: green.Figure 2
Transcritical bifurcation varyingD. Stable: blue. Unstable: green.These results can be used at the moment to choose an efficient dilution coefficientD; that is, if the reactor starts its functioning without microorganisms adhered, the correct election for D is a value less than 0.034h-1, because in this interval for D the reactor is stable and this means that the reactor works well and the organic and inorganic material is being removed from the water. In the case that the reactor starts its functioning with any microorganisms adhered, the correct election for D is a value between 0.034h-1 and 1h-1, for the same reasons given before.On the other hand, the Sotomayor theorem establishes conditions for the presence of a transcritical bifurcation in a system. However, the system studied in this work presents a transcritical bifurcation but does not satisfy the Sotomayor theorem. So, the conditions of Sotomayor theorem are sufficient but not necessary and one example was found.The Sotomayor theorem is stated as follows.Consider a continuous time system that depends on a one-dimensional parameter:(19)x˙=fx,μ,x∈Rn,μ∈R,where f is smooth with respect to x and μ. Let x=x0 be an equilibrium point with a vanishing eigenvalue at μ=μ0. Assume that f(x0,μ0)=0 and the matrix A=Df(x0,μ0) has an eigenvalue λ=0 with eigenvector v and AT has an eigenvector w corresponding to the eigenvalue λ=0. Furthermore, assume that A has k eigenvalues with a negative real part and that n-k-1 with a positive real part. The system presents the following.A saddle-node bifurcation if
(1)
wTfμ(x0,μ0)≠0,(2)
wT[D2f(x0,μ0)(v,v)]≠0.A transcritical bifurcation if
(1)
wTfμ(x0,μ0)=0,(2)
wT[Dfμ(x0,μ0)v]≠0,(3)
wT[D2f(x0,μ0)(v,v)]≠0.The nature of the expressions above is as follows:(i)
fμ(x0,μ0)=∂f(x,μ)/∂μ(x0,μ0)(ii)
Dfμ(x0,μ0)→Jacobiano defμ(x0,μ0)(iii)
D2f(x0,μ0)(v,v)=∑i,jn‍∂2x/∂xixjvivjdondex=(x1,…,xn), v=(v1,…,vn).The modeled system presents a transcritical bifurcation, although the second condition of the Sotomayor theorem is not satisfied. In fact(20)
D
f
x
0
,
μ
0
=
-
D
-
γ
-
1
X
u
μ
′
-
γ
-
1
μ
-
γ
-
1
μ
X
u
μ
′
0
β
+
μ
11
0
0
0
=
a
11
a
12
a
13
a
21
0
a
23
0
0
0
,
w
=
0
0
1
,
v
=
v
1
v
2
v
3
=
v
3
-
a
23
a
21
a
11
a
23
-
a
21
a
13
a
21
1with v3 free constant.Condition  1. We have that(21)fμx0,μ0=S0-S-Xu0x0,D0,wTfμx0,μ0=0.Condition  2. We have that(22)Dfμx0,μ0=-1000-10000wTDfμx0,μ0v=wT-v1-v20=001-v1-v20=0.


#### 2.2.4. Bifurcations Diagram

Bifurcations are structural changes in the system. These types of changes can be appreciated by means of their correspondingbifurcation diagrams as it determines if the system converges to an n-cycle or if its behavior is random.Figure3 shows that for values D<0.034 the trajectories converge at the equilibrium point E3; that is, for a higher dilution rate, the system stabilizes at a larger substrate and suspended biomass concentrations. For this value, there is, however, no formation of a biofilm. For D>0.034 the trajectories converge at the equilibrium point E4, where there is biofilm formation. This coincides with the stability curves E3 and E4 shown in Figures 1 and 2.Figure 3
Bifurcation diagram varyingD with respect to S, Xu, and Xw. The adherence rate α is zero.Bifurcation Diagram Algorithm
Initial conditional =[40,170,20].
The parameter values used are the same data given at the beginning except forα, whose value is α=0.
Integration time =[0,6000].
Time step =0.05.
Integrator = ODE45.
Relative tolerance =1×10-5.
Absolute tolerance =1×10-6.
Interval partition =3000 points.Algorithm. The initial condition is taken. For each parameter we obtain the temporal evolution for the state variables using the same initial condition. We consider as Poincaré application the surface containing the local maxima in steady state that are plotted for each parameter value. To approximate the maximum is applied a simple linear interpolation scheme whenever a sign change is detected in the derivative (see equations (2.13) and (2.14) of [26]).


### 2.3. Haldane Kinetics

Haldane kinetics is used to describe the relation between the microorganism growth rate and the concentration in the inhibitory substrate. The model of the present work was used to study the possible inhibition of the substrate in the reactions taking place in the bioreactor. The growth rate of the microorganisms is determined byμ=μmS/KS+S+S2/KI, where KI is defined in an analogous way as KS; in this case, however, it represents the value of [S] when the inhibition is at its half maximum.Figure4 shows that each model of microbial growth presents a different behavior according to the values of the proposed parameters [27]. For the Monod curve, the speed of the biological growth process will tend asymptotically to the maximum value μm [28]. For the Haldane curve, the microbial growth rate is high for low concentrations of the substrate and low for high concentrations, which demonstrates the negative effect of a high substrate concentration. Both models exhibit a maximum value in the microbial growth rate [27, 29]. For the simulations KI=1500. The value of the constant is a characteristic typical of inhibitor. In this case when S=KI the degree of inhibition is about 50%.Figure 4
Growth rate for Monod and Haldane kinetics.


#### 2.3.1. Equilibrium Points

The system has an equilibrium at the washing condition for any value of the parameters, a physically feasible equilibrium when the microorganisms adhere to the support surface, and three equilibria when the suspended microorganisms do not adhere to the support surface.(i)
First equilibrium point, corresponding to the washing condition, is as follows:E1=(S0,0,0).(ii)
Ifα≠0 there is a nontrivial equilibrium point: E2=(S∗,Xu∗,Xw∗)≠(0,0,0). This point cannot be found analytically.(iii)
Ifα=0, there are three equilibrium points.LetXw=0 and μ(S)-D-k=0. From this, we obtain (23)S2+bS+c=0whereb=KID+kD+k-μmc=KSKI.S exists if Δ=b2-4c>0; that is, a1D2+b1D+c1>0 with (24)a1=KI-4KSb1=2ka1-2μmKIc1=k2a1-2μmKIk+KIμm2D=-b1±b12-4a1c12a1⟹D=KIμm-k+2KIKSKI-2KIKS=0.213D=KIμm-k-2KIKSKI+2KIKS=0.0721.It can then be seen thatS exists if D∈(0,0.0721)∪(0.213,1).WhenS exists, that is, when S is physically possible (the value of S is positive), two equilibrium points are obtained: (25)E3=S3,Xu,Xw=-b+b2-4c2,γDS0-SD+k,0E4=S4,Xu,Xw=-b-b2-4c2,γDS0-SD+k,0.Ifb>0, that is, D>μm-k=0.11, then S3,S4<0 given that 4c>0. Therefore, for D>0.11, the equilibria E3,E4 are not physically feasible.Ifb<0, that is, D<μm-k=0.11, then S3,S4>0.Xu>0 at the equilibrium point E3 if (26)D>KIμm-k+βkKI-β=0.0704withβ=-1S0S02+c.From the above,E3 has physical meaning if 0.0704<D and D<0.0721,E4 has physical meaning if 0<D<0.0721,ifXw≠0 there is an equilibrium point, E5=(S,Xu,Xw), where (27)Xu=γSo-SμS-kμSXw=1.1β+k-μSβ+k-μSXwm.SubstitutingXu and Xw in the following equation, we can solve S. Consider (28)DSo-S-γ-1μSXu+Xw=0.


#### 2.3.2. Equilibrium Points Stability

For equilibrium pointE1 the stability analysis is analogous to the case of Monod kinetics, and thus, it can be immediately concluded that E1 is also a saddle point.For equilibrium pointE2 the numerical stability analysis shows that it is always stable.For equilibrium pointsE3 and E4, we analytically evaluated the stability of the values of D where the equilibria are physically possible. Conducting a stability analysis analogous to the case of Monod kinetics, it is concluded that if D<0.034 then the equilibrium point E4 is stable, whereas for 0.034<D<0.0721 the equilibrium points E3 and E4 are unstable.For equilibrium pointE5 the process is analogous to point E4 of the previous section. It is observed that for values D<0.034 there is a sign change such that the equilibrium is unstable, whereas for values D>0.034 there are no sign changes, and the equilibrium is consequently stable.In the same way as in the Monod kinetics these results can be used at the moment to choose an efficient dilution coefficientD; that is, if the reactor starts its functioning without microorganisms adhered, the correct election for D is a value less than 0.034, because in this interval for D the reactor is stable. In the case that the reactor starts its functioning with any microorganisms adhered, the correct election for D is a value between 0.034 and 1, for the same reasons given before.


#### 2.3.3. One-Parameter Bifurcations

Transcritical Bifurcation. When α=0, for points E4 and E5 there is a transcritical bifurcation when (29)D=β+k1-G0G0=0.034.The analysis of the transcritical bifurcation using the Sotomayor theorem is analogous to that in the previous section. The second condition is not met.Saddle-Node Bifurcation. In what follows, the saddle-node bifurcation is described in [25] as it corresponds to that studied in this section.Letx=x0 be a nonhyperbolic equilibrium point, and let d be the bifurcation parameter. An eigenvalue λ=0 occurs at d=μ0. The structure of the system is characterized as follows.For valuesd<μ0 there are two equilibrium points.For valuesd=μ0 there is a unique nonhyperbolic equilibrium point.For valuesd>μ0 there are no equilibrium points.The system presents a saddle-node bifurcation in the equilibrium curves ofE3 and E4. Indeed, an eigenvalue λ=0 exists if Xuμ′(S)=0, from here μ′(S)=0.There are two values for the parameterD for which μ′(S)=0:(30)D1=KIμm-k-2KIKSKI+2KIKS=0.0721,D2=KIμm-k+2KIKSKI-2KIKS=0.213.The dynamics are as follows.ForD<0.0721 there are two physically sensible equilibrium points, E3 and E4.ForD=0.0721E3=E4.For0.0721<D<0.213 the equilibria E3 and E4 disappear.ForD=0.213E3=E4.ForD>0.213 the points E3 and E4 are restored, although they do not possess physical meaning.The system satisfied the Sotomayor theorem. Indeed,(31)
D
f
x
0
,
μ
0
=
-
D
-
γ
-
1
X
u
μ
′
-
γ
-
1
μ
-
γ
-
1
μ
0
0
β
+
μ
1
-
G
0
0
0
μ
G
0
-
β
-
k
=
a
11
a
12
a
13
0
0
a
23
0
0
a
33
w
=
w
1
w
2
w
3
=
w
3
0
-
a
33
a
23
1
w
1
=
0
a
n
d
w
3
free
.If w3=1 then w=0w21 and(32)
v
=
v
1
v
2
v
3
=
v
2
-
a
12
a
11
1
0
v
2
free
.If v2=1 then v=v110.Condition  1. We have that(33)fμx0,μ0=S0-S-Xu0x0,μ0wTfμx0,μ0=-w2Xu≠0.Condition  2. F=F(S,Xu,Xw). xi are state variables, x1=S, x2=Xu, and x3=Xw. Consider(34)D2fx0,μ0v,v=∂2F∂x1x1v1v1+∂2F∂x1x2v1+∂2F∂x2x1v1+∂2x∂x2x2=-γ-1Xuμ′′SXuμ′′S0,wTD2fx0,D0v,v=w2Xuμ′′S≠0.The system satisfied the Sotomayor theorem, meeting both conditions, and thus, it presents a saddle-node bifurcation at the two values of the parameters D, D1, and D2 as described above.


#### 2.3.4. Bifurcations Diagram

The bifurcation diagram of the system for varyingD when α=0 is analogous to the case of Monod kinetics. This is because the equilibrium points E3 and E4 are stable for D<0.034, and given that the initial condition [40,170,20] is closer to E4 the trajectories will tend to this equilibrium; that is, the suspended biomass concentration increases because there is no formation of biofilms. For D>0.034 the trajectories of the system tend to the stable equilibrium point E5.


## 2.1. Mathematical Model

The continuous stirred-tank reactor described in [12] was chosen as the model reactor. The system consists of three ordinary nonlinear differential equations and models the formation of heterotrophic aerobic biofilms inside the reactor, taking into account both the microorganisms adhered to the biofilms and the microorganisms that are suspended. In addition to studying the influence of the various parameters in biofilm formation, in [12] the system is established and modeled using Monod kinetics, the equilibrium corresponding to the washing condition. In the present work, the system is also analyzed as modeled by Haldane kinetics and more equilibrium points are found and their stability is studied, in addition to characterizing the bifurcations that they represent. As bifurcation parameter of the system, the dilution coefficient D∈[0,1] was chosen, due to this parameter allows to control the residence time of wastewater in the reactor.The following assumptions are made with respect to the dynamics of the formation of biofilms. The three-dimensional structure of the biofilms is ignored because the biofilms are considered to be infinitely thin. There are a finite number of colonization sites available on the support surface, as well as maximum possible surface density of adhered microorganisms. The adhered cells are separated outside the fluid at a rate proportional to its density and the daughter cells of adhered microorganisms compete for space on the support surface—a fractionG of the daughter cells finds adhesion sites and a fraction 1-G does not. The function G decreases with the number of daughter cells (as the number of daughter cells increases, the number of available sites to occupy decreases), where G=(1-W)/(1.1-W), which is a function of probability [12].The reactor dynamics are described by the following equations:(1)dSdt=DS0-S-γ-1μSXu+Xw,dXudt=XuμS-D-k+βXw+XwμS1-GW-αXu1-W,dXwdt=XwμSGW-β-k+αXu1-W.The following parameters are used for the numerical simulation: (2)γ=0.5,k=0.01,Xw,m=5000,D=1,β=0.03,α=0.1,S0=500,μm=0.12,KS=80,Xu0=1000.These values are taken from [12].The growth of the suspended and adhered microorganisms in the reactor is a nonlinear natural phenomenon that can be analyzed using nonlinear dynamics, which provides a global perspective of the types of behavior in the system, such as equilibrium states, stability, and bifurcations.


## 2.2. Monod Kinetics

The suspended and adhered microorganisms inside the reactor consume the organic matter present in the wastewater to be treated. The specific rate of substrate consumption is the “speed” at which the organism consumes the substrate. Although the growth of microorganisms is a complex phenomenon, there are equations that can model this behavior. Monod’s equation is the simplest and most widely utilized equation for describing the kinetics of microbial growth. It is a function of the limiting substrate and is expressed asμ(S)=μmS/KS+S where μm is the maximum growth rate. KS is the concentration of the substrate corresponding to the half of the maximum growth rate of (S) and represents the affinity of the microorganisms for the substrate. If the organism has a great affinity for the limiting substrate, the value of KS is low [17].


### 2.2.1. Equilibrium Points

A pointx0∈Rn is an equilibrium point or a critical point of x˙=f(x) if f(x0)=0. In this point, the system remains in a stationary state; that is, if ϕt:U→Rn is the flux of the system, then ϕt(x0)=x0∀t∈R [18, 19].The system has an equilibrium in the washing condition for any set of values of the parameters, in this equilibrium there aren’t microorganisms adhere or suspended, but there is substrate into the reactor. Exists other equilibrium physically feasible when the microorganisms adhere to the support surface and finally two equilibria in the case where suspended microorganisms do not adhere to the support surface.(i)
First equilibrium point, corresponding to the washing condition, isE1=(S0,0,0).(ii)
Ifα≠0, that is, there adherence of microorganisms to biofilms, then a nontrivial equilibrium point exists: E2=(S∗,Xu∗,Xw∗)≠(0,0,0). This equilibrium point cannot be found analytically.(iii)
Ifα=0, there exist two equilibrium points:(3)E3=S,Xu,Xw=D+kKSμm-D-k,γDS0-SD+k,0withμ(S)-D-k=0.This point is physically possible ifμm-D-k>0 and S0-S>0; that is, (4)D<S0μm-k-kKSS0+KS=0.093448.E4=(S,Xu,Xw), where (5)Xu=γS0-SμS-kμSXw=1.1β+k-μSβ+k-μSXwm.SubstitutingXu and Xw in the following equation, one can solve S. Consider(6)DSo-S-γ-1μSXu+Xw=0.


### 2.2.2. Equilibrium Point Stability

An equilibrium point of a dynamical system is stable in Lyapunov sense if all trajectories with initial conditions near the equilibrium point remain in that vicinity. The equilibrium point is asymptotically stable, that is, the trajectories always tend toward the equilibrium point. Lyapunov’s indirect method is outlined in [20] and demonstrated in [21], and it provides a procedure to determine the local stability of a hyperbolic equilibrium point. Furthermore, the Hartman-Grobman theorem states that, near a hyperbolic equilibrium point, the nonlinear system presents a behavior qualitatively equivalent to that of the corresponding linear system [22, 23], where a hyperbolic equilibrium point is an equilibrium point in which its Jacobian matrix does not have eigenvalues with zero real part.In the following, the stability of the equilibrium points is analyzed by Lyapunov’s indirect method, linearizing the system and calculating its eigenvalues. To analyze the location in the plane of the eigenvalues, the Routh-Hurwitz criterion was utilized.At the equilibrium pointE1, the Jacobian matrix of the system is (7)-D-γ-1μS0-γ-1μS00μS0-D-k-αμS01-G0+β0αμS0G0-k-β.Let A be a submatrix, such that (8)A=μS0-D-k-αμS01-G0+βαμS0G0-k-β.The eigenvalues λ1, λ2, and λ3 of the matrix JE1 are the roots of the equation: (9)λ2-δλ+τλ+D=0,δ=trA,τ=detA,λ1,2=δ±δ2-4τ2,λ3=-D.It is readily apparent that (10)maxμ-D-k-α,μG-k-β<λ1.Substituting the values of the parameters, one obtains that μ(S0)G(0)-k-β>0. Therefore, E1 is unstable.For the equilibrium pointE2, numerically performing the stability analysis, we find that it is always stable.The algorithm is as follows: for a given value ofD the equilibrium point is found using the Newton-Raphson method, and the equilibrium curve is constructed using the same method taking the previous equilibrium as the initial condition. Every equilibrium point is evaluated in the Jacobian matrix, and the eigenvalues are calculated. If they are all negative, the equilibrium is stable. If either of them is negative, the equilibrium is unstable and if at least one has a zero real part, it is not possible to decide on the stability of the equilibrium.At the equilibrium pointE3, the Jacobian matrix of the system is (11)-D-γ-1Xuμ′-γ-1μ-γ-1μXuμ′0β+μ1-G000μG0-β-k=a11a12a13a210a2300a33and its characteristic polynomial is (12)λ3-a11+a33λ2+a11a33-a12a21λ+a12a21a33=0a11,a12<0,a21>0.We establish the following array: (13)λ31a11a33-a12a21λ2-a11+a33a12a21a33λ1a11a33-a11a12a21a11+a330λ0a12a21a33The Routh-Hurwitz stability criterion states that the number of roots of the characteristic polynomial with a positive real part is equal to the number of sign changes of the coefficients of the first column of the array.a
11
+
a
33 is negative. Indeed a11+a33=(DG(0)-D)+(kG(0)-k)-γ-1Xuμ′(S)-β<0 given that DG(0)-D<0 and kG(0)-k<0.Ifa33<0 then a11a33>0 and a11a33-a11a12a21/a11+a33>0.Ifa33<0 then D<β+k-kG(0)/G(0).Therefore, forD<0.034, a12a21a33>0.Ifa33>0 then D>β+k-kG(0)/G(0).Therefore, forD>0.034, a12a21a33<0.In summary, we have(14)D<0.0340.034<D<0.093448λ3+λ2+λ1+λ0+λ3+λ2+λ1xλ0-If D<0.034, there are no sign changes, and, thus, the real part of all the eigenvalues is negative, and the equilibrium point is stable. If 0.034<D<0.093448 there is always a sign change (x is unknown and it does not matter, because anyway there is a sign change), which means that there is an eigenvalue with a positive real part, and, thus, the equilibrium point is unstable.For the equilibrium pointE4, the characteristic polynomial is λ3+bλ2+cλ+d=0, where (15)
b
=
-
a
11
+
a
22
+
a
33
,
c
=
-
a
12
a
21
+
a
32
a
23
+
a
13
a
31
-
a
11
a
33
-
a
22
a
33
-
a
11
a
22
,
d
=
-
a
11
a
22
a
33
-
a
12
a
23
a
31
-
a
13
a
21
a
32
+
a
12
a
21
a
33
+
a
32
a
23
a
11
+
a
13
a
31
a
22
.By virtue of the Routh-Hurwitz criterion, one has (16)λ31cλ2bdλ1c1=bc-db0λ0dNumerically, it is found that, for values D<0.034, there is a sign change, and, thus, the equilibrium point is unstable, whereas for values D>0.034, there is no sign change, and, thus, the equilibrium point is stable.


### 2.2.3. One-Parameter Bifurcations

The fact that the dynamics of a system can change drastically as one or more of its parameters are varied is well known. This qualitative or structural change is known as abifurcation. In general, bifurcation theory studies the structural changes that a dynamical system can undergo as its parameters are varied; further bifurcations only occur in nonhyperbolic equilibria [19, 24].In what follows, the transcritical bifurcation described in [25] is presented because it is precisely the one studied in this work.Letx=x0 be a nonhyperbolic equilibrium point. The structure of the system is characterized as follows.(i)
Two curves of equilibrium pointsϕ(t) and ψ(t) pass through x0.(ii)
Both curves exist on both sidesx0.(iii)
The stability of the curves is interchanged as they pass throughx0.The system presents a transcritical bifurcation in the equilibrium curvesE3 and E4 in fact.Whenα=0, for the point E3, an eigenvalue λ=0 exists if a12a21a33=0, given that a12<0 and a21>0; then a33 must be zero. That is, (17)μSG0-β-k=0withμS=D+k.Thus(18)D=β+k-kG0G0=β+k1-G0G0=0.03+0.011-10/1110/11=0.034.The equilibrium point in the bifurcation parameter is as follows: S,Xu,Xw=46.31,175.28,0.The equilibrium curveE3 is stable for D<0.034 and unstable for D>0.034. For E4 the equilibria are unstable for D<0.034 and stable for D>0.034. Figures 1 and 2 show how the equilibrium curves of the points E3 and E4 cross each other and interchange stability, such that a transcritical bifurcation appears.Figure 1
Transcritical bifurcation varyingD. Stable: blue. Unstable: green.Figure 2
Transcritical bifurcation varyingD. Stable: blue. Unstable: green.These results can be used at the moment to choose an efficient dilution coefficientD; that is, if the reactor starts its functioning without microorganisms adhered, the correct election for D is a value less than 0.034h-1, because in this interval for D the reactor is stable and this means that the reactor works well and the organic and inorganic material is being removed from the water. In the case that the reactor starts its functioning with any microorganisms adhered, the correct election for D is a value between 0.034h-1 and 1h-1, for the same reasons given before.On the other hand, the Sotomayor theorem establishes conditions for the presence of a transcritical bifurcation in a system. However, the system studied in this work presents a transcritical bifurcation but does not satisfy the Sotomayor theorem. So, the conditions of Sotomayor theorem are sufficient but not necessary and one example was found.The Sotomayor theorem is stated as follows.Consider a continuous time system that depends on a one-dimensional parameter:(19)x˙=fx,μ,x∈Rn,μ∈R,where f is smooth with respect to x and μ. Let x=x0 be an equilibrium point with a vanishing eigenvalue at μ=μ0. Assume that f(x0,μ0)=0 and the matrix A=Df(x0,μ0) has an eigenvalue λ=0 with eigenvector v and AT has an eigenvector w corresponding to the eigenvalue λ=0. Furthermore, assume that A has k eigenvalues with a negative real part and that n-k-1 with a positive real part. The system presents the following.A saddle-node bifurcation if
(1)
wTfμ(x0,μ0)≠0,(2)
wT[D2f(x0,μ0)(v,v)]≠0.A transcritical bifurcation if
(1)
wTfμ(x0,μ0)=0,(2)
wT[Dfμ(x0,μ0)v]≠0,(3)
wT[D2f(x0,μ0)(v,v)]≠0.The nature of the expressions above is as follows:(i)
fμ(x0,μ0)=∂f(x,μ)/∂μ(x0,μ0)(ii)
Dfμ(x0,μ0)→Jacobiano defμ(x0,μ0)(iii)
D2f(x0,μ0)(v,v)=∑i,jn‍∂2x/∂xixjvivjdondex=(x1,…,xn), v=(v1,…,vn).The modeled system presents a transcritical bifurcation, although the second condition of the Sotomayor theorem is not satisfied. In fact(20)
D
f
x
0
,
μ
0
=
-
D
-
γ
-
1
X
u
μ
′
-
γ
-
1
μ
-
γ
-
1
μ
X
u
μ
′
0
β
+
μ
11
0
0
0
=
a
11
a
12
a
13
a
21
0
a
23
0
0
0
,
w
=
0
0
1
,
v
=
v
1
v
2
v
3
=
v
3
-
a
23
a
21
a
11
a
23
-
a
21
a
13
a
21
1with v3 free constant.Condition  1. We have that(21)fμx0,μ0=S0-S-Xu0x0,D0,wTfμx0,μ0=0.Condition  2. We have that(22)Dfμx0,μ0=-1000-10000wTDfμx0,μ0v=wT-v1-v20=001-v1-v20=0.


### 2.2.4. Bifurcations Diagram

Bifurcations are structural changes in the system. These types of changes can be appreciated by means of their correspondingbifurcation diagrams as it determines if the system converges to an n-cycle or if its behavior is random.Figure3 shows that for values D<0.034 the trajectories converge at the equilibrium point E3; that is, for a higher dilution rate, the system stabilizes at a larger substrate and suspended biomass concentrations. For this value, there is, however, no formation of a biofilm. For D>0.034 the trajectories converge at the equilibrium point E4, where there is biofilm formation. This coincides with the stability curves E3 and E4 shown in Figures 1 and 2.Figure 3
Bifurcation diagram varyingD with respect to S, Xu, and Xw. The adherence rate α is zero.Bifurcation Diagram Algorithm
Initial conditional =[40,170,20].
The parameter values used are the same data given at the beginning except forα, whose value is α=0.
Integration time =[0,6000].
Time step =0.05.
Integrator = ODE45.
Relative tolerance =1×10-5.
Absolute tolerance =1×10-6.
Interval partition =3000 points.Algorithm. The initial condition is taken. For each parameter we obtain the temporal evolution for the state variables using the same initial condition. We consider as Poincaré application the surface containing the local maxima in steady state that are plotted for each parameter value. To approximate the maximum is applied a simple linear interpolation scheme whenever a sign change is detected in the derivative (see equations (2.13) and (2.14) of [26]).


## 2.2.1. Equilibrium Points

A pointx0∈Rn is an equilibrium point or a critical point of x˙=f(x) if f(x0)=0. In this point, the system remains in a stationary state; that is, if ϕt:U→Rn is the flux of the system, then ϕt(x0)=x0∀t∈R [18, 19].The system has an equilibrium in the washing condition for any set of values of the parameters, in this equilibrium there aren’t microorganisms adhere or suspended, but there is substrate into the reactor. Exists other equilibrium physically feasible when the microorganisms adhere to the support surface and finally two equilibria in the case where suspended microorganisms do not adhere to the support surface.(i)
First equilibrium point, corresponding to the washing condition, isE1=(S0,0,0).(ii)
Ifα≠0, that is, there adherence of microorganisms to biofilms, then a nontrivial equilibrium point exists: E2=(S∗,Xu∗,Xw∗)≠(0,0,0). This equilibrium point cannot be found analytically.(iii)
Ifα=0, there exist two equilibrium points:(3)E3=S,Xu,Xw=D+kKSμm-D-k,γDS0-SD+k,0withμ(S)-D-k=0.This point is physically possible ifμm-D-k>0 and S0-S>0; that is, (4)D<S0μm-k-kKSS0+KS=0.093448.E4=(S,Xu,Xw), where (5)Xu=γS0-SμS-kμSXw=1.1β+k-μSβ+k-μSXwm.SubstitutingXu and Xw in the following equation, one can solve S. Consider(6)DSo-S-γ-1μSXu+Xw=0.


## 2.2.2. Equilibrium Point Stability

An equilibrium point of a dynamical system is stable in Lyapunov sense if all trajectories with initial conditions near the equilibrium point remain in that vicinity. The equilibrium point is asymptotically stable, that is, the trajectories always tend toward the equilibrium point. Lyapunov’s indirect method is outlined in [20] and demonstrated in [21], and it provides a procedure to determine the local stability of a hyperbolic equilibrium point. Furthermore, the Hartman-Grobman theorem states that, near a hyperbolic equilibrium point, the nonlinear system presents a behavior qualitatively equivalent to that of the corresponding linear system [22, 23], where a hyperbolic equilibrium point is an equilibrium point in which its Jacobian matrix does not have eigenvalues with zero real part.In the following, the stability of the equilibrium points is analyzed by Lyapunov’s indirect method, linearizing the system and calculating its eigenvalues. To analyze the location in the plane of the eigenvalues, the Routh-Hurwitz criterion was utilized.At the equilibrium pointE1, the Jacobian matrix of the system is (7)-D-γ-1μS0-γ-1μS00μS0-D-k-αμS01-G0+β0αμS0G0-k-β.Let A be a submatrix, such that (8)A=μS0-D-k-αμS01-G0+βαμS0G0-k-β.The eigenvalues λ1, λ2, and λ3 of the matrix JE1 are the roots of the equation: (9)λ2-δλ+τλ+D=0,δ=trA,τ=detA,λ1,2=δ±δ2-4τ2,λ3=-D.It is readily apparent that (10)maxμ-D-k-α,μG-k-β<λ1.Substituting the values of the parameters, one obtains that μ(S0)G(0)-k-β>0. Therefore, E1 is unstable.For the equilibrium pointE2, numerically performing the stability analysis, we find that it is always stable.The algorithm is as follows: for a given value ofD the equilibrium point is found using the Newton-Raphson method, and the equilibrium curve is constructed using the same method taking the previous equilibrium as the initial condition. Every equilibrium point is evaluated in the Jacobian matrix, and the eigenvalues are calculated. If they are all negative, the equilibrium is stable. If either of them is negative, the equilibrium is unstable and if at least one has a zero real part, it is not possible to decide on the stability of the equilibrium.At the equilibrium pointE3, the Jacobian matrix of the system is (11)-D-γ-1Xuμ′-γ-1μ-γ-1μXuμ′0β+μ1-G000μG0-β-k=a11a12a13a210a2300a33and its characteristic polynomial is (12)λ3-a11+a33λ2+a11a33-a12a21λ+a12a21a33=0a11,a12<0,a21>0.We establish the following array: (13)λ31a11a33-a12a21λ2-a11+a33a12a21a33λ1a11a33-a11a12a21a11+a330λ0a12a21a33The Routh-Hurwitz stability criterion states that the number of roots of the characteristic polynomial with a positive real part is equal to the number of sign changes of the coefficients of the first column of the array.a
11
+
a
33 is negative. Indeed a11+a33=(DG(0)-D)+(kG(0)-k)-γ-1Xuμ′(S)-β<0 given that DG(0)-D<0 and kG(0)-k<0.Ifa33<0 then a11a33>0 and a11a33-a11a12a21/a11+a33>0.Ifa33<0 then D<β+k-kG(0)/G(0).Therefore, forD<0.034, a12a21a33>0.Ifa33>0 then D>β+k-kG(0)/G(0).Therefore, forD>0.034, a12a21a33<0.In summary, we have(14)D<0.0340.034<D<0.093448λ3+λ2+λ1+λ0+λ3+λ2+λ1xλ0-If D<0.034, there are no sign changes, and, thus, the real part of all the eigenvalues is negative, and the equilibrium point is stable. If 0.034<D<0.093448 there is always a sign change (x is unknown and it does not matter, because anyway there is a sign change), which means that there is an eigenvalue with a positive real part, and, thus, the equilibrium point is unstable.For the equilibrium pointE4, the characteristic polynomial is λ3+bλ2+cλ+d=0, where (15)
b
=
-
a
11
+
a
22
+
a
33
,
c
=
-
a
12
a
21
+
a
32
a
23
+
a
13
a
31
-
a
11
a
33
-
a
22
a
33
-
a
11
a
22
,
d
=
-
a
11
a
22
a
33
-
a
12
a
23
a
31
-
a
13
a
21
a
32
+
a
12
a
21
a
33
+
a
32
a
23
a
11
+
a
13
a
31
a
22
.By virtue of the Routh-Hurwitz criterion, one has (16)λ31cλ2bdλ1c1=bc-db0λ0dNumerically, it is found that, for values D<0.034, there is a sign change, and, thus, the equilibrium point is unstable, whereas for values D>0.034, there is no sign change, and, thus, the equilibrium point is stable.


## 2.2.3. One-Parameter Bifurcations

The fact that the dynamics of a system can change drastically as one or more of its parameters are varied is well known. This qualitative or structural change is known as abifurcation. In general, bifurcation theory studies the structural changes that a dynamical system can undergo as its parameters are varied; further bifurcations only occur in nonhyperbolic equilibria [19, 24].In what follows, the transcritical bifurcation described in [25] is presented because it is precisely the one studied in this work.Letx=x0 be a nonhyperbolic equilibrium point. The structure of the system is characterized as follows.(i)
Two curves of equilibrium pointsϕ(t) and ψ(t) pass through x0.(ii)
Both curves exist on both sidesx0.(iii)
The stability of the curves is interchanged as they pass throughx0.The system presents a transcritical bifurcation in the equilibrium curvesE3 and E4 in fact.Whenα=0, for the point E3, an eigenvalue λ=0 exists if a12a21a33=0, given that a12<0 and a21>0; then a33 must be zero. That is, (17)μSG0-β-k=0withμS=D+k.Thus(18)D=β+k-kG0G0=β+k1-G0G0=0.03+0.011-10/1110/11=0.034.The equilibrium point in the bifurcation parameter is as follows: S,Xu,Xw=46.31,175.28,0.The equilibrium curveE3 is stable for D<0.034 and unstable for D>0.034. For E4 the equilibria are unstable for D<0.034 and stable for D>0.034. Figures 1 and 2 show how the equilibrium curves of the points E3 and E4 cross each other and interchange stability, such that a transcritical bifurcation appears.Figure 1
Transcritical bifurcation varyingD. Stable: blue. Unstable: green.Figure 2
Transcritical bifurcation varyingD. Stable: blue. Unstable: green.These results can be used at the moment to choose an efficient dilution coefficientD; that is, if the reactor starts its functioning without microorganisms adhered, the correct election for D is a value less than 0.034h-1, because in this interval for D the reactor is stable and this means that the reactor works well and the organic and inorganic material is being removed from the water. In the case that the reactor starts its functioning with any microorganisms adhered, the correct election for D is a value between 0.034h-1 and 1h-1, for the same reasons given before.On the other hand, the Sotomayor theorem establishes conditions for the presence of a transcritical bifurcation in a system. However, the system studied in this work presents a transcritical bifurcation but does not satisfy the Sotomayor theorem. So, the conditions of Sotomayor theorem are sufficient but not necessary and one example was found.The Sotomayor theorem is stated as follows.Consider a continuous time system that depends on a one-dimensional parameter:(19)x˙=fx,μ,x∈Rn,μ∈R,where f is smooth with respect to x and μ. Let x=x0 be an equilibrium point with a vanishing eigenvalue at μ=μ0. Assume that f(x0,μ0)=0 and the matrix A=Df(x0,μ0) has an eigenvalue λ=0 with eigenvector v and AT has an eigenvector w corresponding to the eigenvalue λ=0. Furthermore, assume that A has k eigenvalues with a negative real part and that n-k-1 with a positive real part. The system presents the following.A saddle-node bifurcation if
(1)
wTfμ(x0,μ0)≠0,(2)
wT[D2f(x0,μ0)(v,v)]≠0.A transcritical bifurcation if
(1)
wTfμ(x0,μ0)=0,(2)
wT[Dfμ(x0,μ0)v]≠0,(3)
wT[D2f(x0,μ0)(v,v)]≠0.The nature of the expressions above is as follows:(i)
fμ(x0,μ0)=∂f(x,μ)/∂μ(x0,μ0)(ii)
Dfμ(x0,μ0)→Jacobiano defμ(x0,μ0)(iii)
D2f(x0,μ0)(v,v)=∑i,jn‍∂2x/∂xixjvivjdondex=(x1,…,xn), v=(v1,…,vn).The modeled system presents a transcritical bifurcation, although the second condition of the Sotomayor theorem is not satisfied. In fact(20)
D
f
x
0
,
μ
0
=
-
D
-
γ
-
1
X
u
μ
′
-
γ
-
1
μ
-
γ
-
1
μ
X
u
μ
′
0
β
+
μ
11
0
0
0
=
a
11
a
12
a
13
a
21
0
a
23
0
0
0
,
w
=
0
0
1
,
v
=
v
1
v
2
v
3
=
v
3
-
a
23
a
21
a
11
a
23
-
a
21
a
13
a
21
1with v3 free constant.Condition  1. We have that(21)fμx0,μ0=S0-S-Xu0x0,D0,wTfμx0,μ0=0.Condition  2. We have that(22)Dfμx0,μ0=-1000-10000wTDfμx0,μ0v=wT-v1-v20=001-v1-v20=0.


## 2.2.4. Bifurcations Diagram

Bifurcations are structural changes in the system. These types of changes can be appreciated by means of their correspondingbifurcation diagrams as it determines if the system converges to an n-cycle or if its behavior is random.Figure3 shows that for values D<0.034 the trajectories converge at the equilibrium point E3; that is, for a higher dilution rate, the system stabilizes at a larger substrate and suspended biomass concentrations. For this value, there is, however, no formation of a biofilm. For D>0.034 the trajectories converge at the equilibrium point E4, where there is biofilm formation. This coincides with the stability curves E3 and E4 shown in Figures 1 and 2.Figure 3
Bifurcation diagram varyingD with respect to S, Xu, and Xw. The adherence rate α is zero.Bifurcation Diagram Algorithm
Initial conditional =[40,170,20].
The parameter values used are the same data given at the beginning except forα, whose value is α=0.
Integration time =[0,6000].
Time step =0.05.
Integrator = ODE45.
Relative tolerance =1×10-5.
Absolute tolerance =1×10-6.
Interval partition =3000 points.Algorithm. The initial condition is taken. For each parameter we obtain the temporal evolution for the state variables using the same initial condition. We consider as Poincaré application the surface containing the local maxima in steady state that are plotted for each parameter value. To approximate the maximum is applied a simple linear interpolation scheme whenever a sign change is detected in the derivative (see equations (2.13) and (2.14) of [26]).


## 2.3. Haldane Kinetics

Haldane kinetics is used to describe the relation between the microorganism growth rate and the concentration in the inhibitory substrate. The model of the present work was used to study the possible inhibition of the substrate in the reactions taking place in the bioreactor. The growth rate of the microorganisms is determined byμ=μmS/KS+S+S2/KI, where KI is defined in an analogous way as KS; in this case, however, it represents the value of [S] when the inhibition is at its half maximum.Figure4 shows that each model of microbial growth presents a different behavior according to the values of the proposed parameters [27]. For the Monod curve, the speed of the biological growth process will tend asymptotically to the maximum value μm [28]. For the Haldane curve, the microbial growth rate is high for low concentrations of the substrate and low for high concentrations, which demonstrates the negative effect of a high substrate concentration. Both models exhibit a maximum value in the microbial growth rate [27, 29]. For the simulations KI=1500. The value of the constant is a characteristic typical of inhibitor. In this case when S=KI the degree of inhibition is about 50%.Figure 4
Growth rate for Monod and Haldane kinetics.


### 2.3.1. Equilibrium Points

The system has an equilibrium at the washing condition for any value of the parameters, a physically feasible equilibrium when the microorganisms adhere to the support surface, and three equilibria when the suspended microorganisms do not adhere to the support surface.(i)
First equilibrium point, corresponding to the washing condition, is as follows:E1=(S0,0,0).(ii)
Ifα≠0 there is a nontrivial equilibrium point: E2=(S∗,Xu∗,Xw∗)≠(0,0,0). This point cannot be found analytically.(iii)
Ifα=0, there are three equilibrium points.LetXw=0 and μ(S)-D-k=0. From this, we obtain (23)S2+bS+c=0whereb=KID+kD+k-μmc=KSKI.S exists if Δ=b2-4c>0; that is, a1D2+b1D+c1>0 with (24)a1=KI-4KSb1=2ka1-2μmKIc1=k2a1-2μmKIk+KIμm2D=-b1±b12-4a1c12a1⟹D=KIμm-k+2KIKSKI-2KIKS=0.213D=KIμm-k-2KIKSKI+2KIKS=0.0721.It can then be seen thatS exists if D∈(0,0.0721)∪(0.213,1).WhenS exists, that is, when S is physically possible (the value of S is positive), two equilibrium points are obtained: (25)E3=S3,Xu,Xw=-b+b2-4c2,γDS0-SD+k,0E4=S4,Xu,Xw=-b-b2-4c2,γDS0-SD+k,0.Ifb>0, that is, D>μm-k=0.11, then S3,S4<0 given that 4c>0. Therefore, for D>0.11, the equilibria E3,E4 are not physically feasible.Ifb<0, that is, D<μm-k=0.11, then S3,S4>0.Xu>0 at the equilibrium point E3 if (26)D>KIμm-k+βkKI-β=0.0704withβ=-1S0S02+c.From the above,E3 has physical meaning if 0.0704<D and D<0.0721,E4 has physical meaning if 0<D<0.0721,ifXw≠0 there is an equilibrium point, E5=(S,Xu,Xw), where (27)Xu=γSo-SμS-kμSXw=1.1β+k-μSβ+k-μSXwm.SubstitutingXu and Xw in the following equation, we can solve S. Consider (28)DSo-S-γ-1μSXu+Xw=0.


### 2.3.2. Equilibrium Points Stability

For equilibrium pointE1 the stability analysis is analogous to the case of Monod kinetics, and thus, it can be immediately concluded that E1 is also a saddle point.For equilibrium pointE2 the numerical stability analysis shows that it is always stable.For equilibrium pointsE3 and E4, we analytically evaluated the stability of the values of D where the equilibria are physically possible. Conducting a stability analysis analogous to the case of Monod kinetics, it is concluded that if D<0.034 then the equilibrium point E4 is stable, whereas for 0.034<D<0.0721 the equilibrium points E3 and E4 are unstable.For equilibrium pointE5 the process is analogous to point E4 of the previous section. It is observed that for values D<0.034 there is a sign change such that the equilibrium is unstable, whereas for values D>0.034 there are no sign changes, and the equilibrium is consequently stable.In the same way as in the Monod kinetics these results can be used at the moment to choose an efficient dilution coefficientD; that is, if the reactor starts its functioning without microorganisms adhered, the correct election for D is a value less than 0.034, because in this interval for D the reactor is stable. In the case that the reactor starts its functioning with any microorganisms adhered, the correct election for D is a value between 0.034 and 1, for the same reasons given before.


### 2.3.3. One-Parameter Bifurcations

Transcritical Bifurcation. When α=0, for points E4 and E5 there is a transcritical bifurcation when (29)D=β+k1-G0G0=0.034.The analysis of the transcritical bifurcation using the Sotomayor theorem is analogous to that in the previous section. The second condition is not met.Saddle-Node Bifurcation. In what follows, the saddle-node bifurcation is described in [25] as it corresponds to that studied in this section.Letx=x0 be a nonhyperbolic equilibrium point, and let d be the bifurcation parameter. An eigenvalue λ=0 occurs at d=μ0. The structure of the system is characterized as follows.For valuesd<μ0 there are two equilibrium points.For valuesd=μ0 there is a unique nonhyperbolic equilibrium point.For valuesd>μ0 there are no equilibrium points.The system presents a saddle-node bifurcation in the equilibrium curves ofE3 and E4. Indeed, an eigenvalue λ=0 exists if Xuμ′(S)=0, from here μ′(S)=0.There are two values for the parameterD for which μ′(S)=0:(30)D1=KIμm-k-2KIKSKI+2KIKS=0.0721,D2=KIμm-k+2KIKSKI-2KIKS=0.213.The dynamics are as follows.ForD<0.0721 there are two physically sensible equilibrium points, E3 and E4.ForD=0.0721E3=E4.For0.0721<D<0.213 the equilibria E3 and E4 disappear.ForD=0.213E3=E4.ForD>0.213 the points E3 and E4 are restored, although they do not possess physical meaning.The system satisfied the Sotomayor theorem. Indeed,(31)
D
f
x
0
,
μ
0
=
-
D
-
γ
-
1
X
u
μ
′
-
γ
-
1
μ
-
γ
-
1
μ
0
0
β
+
μ
1
-
G
0
0
0
μ
G
0
-
β
-
k
=
a
11
a
12
a
13
0
0
a
23
0
0
a
33
w
=
w
1
w
2
w
3
=
w
3
0
-
a
33
a
23
1
w
1
=
0
a
n
d
w
3
free
.If w3=1 then w=0w21 and(32)
v
=
v
1
v
2
v
3
=
v
2
-
a
12
a
11
1
0
v
2
free
.If v2=1 then v=v110.Condition  1. We have that(33)fμx0,μ0=S0-S-Xu0x0,μ0wTfμx0,μ0=-w2Xu≠0.Condition  2. F=F(S,Xu,Xw). xi are state variables, x1=S, x2=Xu, and x3=Xw. Consider(34)D2fx0,μ0v,v=∂2F∂x1x1v1v1+∂2F∂x1x2v1+∂2F∂x2x1v1+∂2x∂x2x2=-γ-1Xuμ′′SXuμ′′S0,wTD2fx0,D0v,v=w2Xuμ′′S≠0.The system satisfied the Sotomayor theorem, meeting both conditions, and thus, it presents a saddle-node bifurcation at the two values of the parameters D, D1, and D2 as described above.


### 2.3.4. Bifurcations Diagram

The bifurcation diagram of the system for varyingD when α=0 is analogous to the case of Monod kinetics. This is because the equilibrium points E3 and E4 are stable for D<0.034, and given that the initial condition [40,170,20] is closer to E4 the trajectories will tend to this equilibrium; that is, the suspended biomass concentration increases because there is no formation of biofilms. For D>0.034 the trajectories of the system tend to the stable equilibrium point E5.


## 2.3.1. Equilibrium Points

The system has an equilibrium at the washing condition for any value of the parameters, a physically feasible equilibrium when the microorganisms adhere to the support surface, and three equilibria when the suspended microorganisms do not adhere to the support surface.(i)
First equilibrium point, corresponding to the washing condition, is as follows:E1=(S0,0,0).(ii)
Ifα≠0 there is a nontrivial equilibrium point: E2=(S∗,Xu∗,Xw∗)≠(0,0,0). This point cannot be found analytically.(iii)
Ifα=0, there are three equilibrium points.LetXw=0 and μ(S)-D-k=0. From this, we obtain (23)S2+bS+c=0whereb=KID+kD+k-μmc=KSKI.S exists if Δ=b2-4c>0; that is, a1D2+b1D+c1>0 with (24)a1=KI-4KSb1=2ka1-2μmKIc1=k2a1-2μmKIk+KIμm2D=-b1±b12-4a1c12a1⟹D=KIμm-k+2KIKSKI-2KIKS=0.213D=KIμm-k-2KIKSKI+2KIKS=0.0721.It can then be seen thatS exists if D∈(0,0.0721)∪(0.213,1).WhenS exists, that is, when S is physically possible (the value of S is positive), two equilibrium points are obtained: (25)E3=S3,Xu,Xw=-b+b2-4c2,γDS0-SD+k,0E4=S4,Xu,Xw=-b-b2-4c2,γDS0-SD+k,0.Ifb>0, that is, D>μm-k=0.11, then S3,S4<0 given that 4c>0. Therefore, for D>0.11, the equilibria E3,E4 are not physically feasible.Ifb<0, that is, D<μm-k=0.11, then S3,S4>0.Xu>0 at the equilibrium point E3 if (26)D>KIμm-k+βkKI-β=0.0704withβ=-1S0S02+c.From the above,E3 has physical meaning if 0.0704<D and D<0.0721,E4 has physical meaning if 0<D<0.0721,ifXw≠0 there is an equilibrium point, E5=(S,Xu,Xw), where (27)Xu=γSo-SμS-kμSXw=1.1β+k-μSβ+k-μSXwm.SubstitutingXu and Xw in the following equation, we can solve S. Consider (28)DSo-S-γ-1μSXu+Xw=0.


## 2.3.2. Equilibrium Points Stability

For equilibrium pointE1 the stability analysis is analogous to the case of Monod kinetics, and thus, it can be immediately concluded that E1 is also a saddle point.For equilibrium pointE2 the numerical stability analysis shows that it is always stable.For equilibrium pointsE3 and E4, we analytically evaluated the stability of the values of D where the equilibria are physically possible. Conducting a stability analysis analogous to the case of Monod kinetics, it is concluded that if D<0.034 then the equilibrium point E4 is stable, whereas for 0.034<D<0.0721 the equilibrium points E3 and E4 are unstable.For equilibrium pointE5 the process is analogous to point E4 of the previous section. It is observed that for values D<0.034 there is a sign change such that the equilibrium is unstable, whereas for values D>0.034 there are no sign changes, and the equilibrium is consequently stable.In the same way as in the Monod kinetics these results can be used at the moment to choose an efficient dilution coefficientD; that is, if the reactor starts its functioning without microorganisms adhered, the correct election for D is a value less than 0.034, because in this interval for D the reactor is stable. In the case that the reactor starts its functioning with any microorganisms adhered, the correct election for D is a value between 0.034 and 1, for the same reasons given before.


## 2.3.3. One-Parameter Bifurcations

Transcritical Bifurcation. When α=0, for points E4 and E5 there is a transcritical bifurcation when (29)D=β+k1-G0G0=0.034.The analysis of the transcritical bifurcation using the Sotomayor theorem is analogous to that in the previous section. The second condition is not met.Saddle-Node Bifurcation. In what follows, the saddle-node bifurcation is described in [25] as it corresponds to that studied in this section.Letx=x0 be a nonhyperbolic equilibrium point, and let d be the bifurcation parameter. An eigenvalue λ=0 occurs at d=μ0. The structure of the system is characterized as follows.For valuesd<μ0 there are two equilibrium points.For valuesd=μ0 there is a unique nonhyperbolic equilibrium point.For valuesd>μ0 there are no equilibrium points.The system presents a saddle-node bifurcation in the equilibrium curves ofE3 and E4. Indeed, an eigenvalue λ=0 exists if Xuμ′(S)=0, from here μ′(S)=0.There are two values for the parameterD for which μ′(S)=0:(30)D1=KIμm-k-2KIKSKI+2KIKS=0.0721,D2=KIμm-k+2KIKSKI-2KIKS=0.213.The dynamics are as follows.ForD<0.0721 there are two physically sensible equilibrium points, E3 and E4.ForD=0.0721E3=E4.For0.0721<D<0.213 the equilibria E3 and E4 disappear.ForD=0.213E3=E4.ForD>0.213 the points E3 and E4 are restored, although they do not possess physical meaning.The system satisfied the Sotomayor theorem. Indeed,(31)
D
f
x
0
,
μ
0
=
-
D
-
γ
-
1
X
u
μ
′
-
γ
-
1
μ
-
γ
-
1
μ
0
0
β
+
μ
1
-
G
0
0
0
μ
G
0
-
β
-
k
=
a
11
a
12
a
13
0
0
a
23
0
0
a
33
w
=
w
1
w
2
w
3
=
w
3
0
-
a
33
a
23
1
w
1
=
0
a
n
d
w
3
free
.If w3=1 then w=0w21 and(32)
v
=
v
1
v
2
v
3
=
v
2
-
a
12
a
11
1
0
v
2
free
.If v2=1 then v=v110.Condition  1. We have that(33)fμx0,μ0=S0-S-Xu0x0,μ0wTfμx0,μ0=-w2Xu≠0.Condition  2. F=F(S,Xu,Xw). xi are state variables, x1=S, x2=Xu, and x3=Xw. Consider(34)D2fx0,μ0v,v=∂2F∂x1x1v1v1+∂2F∂x1x2v1+∂2F∂x2x1v1+∂2x∂x2x2=-γ-1Xuμ′′SXuμ′′S0,wTD2fx0,D0v,v=w2Xuμ′′S≠0.The system satisfied the Sotomayor theorem, meeting both conditions, and thus, it presents a saddle-node bifurcation at the two values of the parameters D, D1, and D2 as described above.


## 2.3.4. Bifurcations Diagram

The bifurcation diagram of the system for varyingD when α=0 is analogous to the case of Monod kinetics. This is because the equilibrium points E3 and E4 are stable for D<0.034, and given that the initial condition [40,170,20] is closer to E4 the trajectories will tend to this equilibrium; that is, the suspended biomass concentration increases because there is no formation of biofilms. For D>0.034 the trajectories of the system tend to the stable equilibrium point E5.


## 3. Conclusions

The reactor is stable when it works well, that is, when the organic and inorganic material is being removed from the water. The system analyzed with both Monod and Haldane kinetics is stable with the initial parameter values, that is, when the adherence rate of the microorganisms is nonzero.If the adherence rate (α) is zero (this may occur when the residence time of the wastewater inside the reactor is not long enough so that the suspended microorganisms adhere to biofilms) both systems (Monod and Haldane) present a transcritical bifurcation for the same value of D because the value of the dilution rate at the bifurcation does not depend on the type of kinetics being used to model the system, but instead it depends on the mortality and detachment rates of the adhered microorganisms. Also, the model studied in this paper is an example that the conditions of Sotomayor theorem are sufficient but not necessary, because this model presents a transcritical bifurcation but does not satisfy the Sotomayor theorem. On the other hand, these results can be used at the moment to choose an efficient dilution coefficient; in fact, for both systems we have the following: if the reactor starts its functioning without microorganisms adhered, the correct election for D is a value less than 0.034, because in this interval for D the reactor is stable; in the case that the reactor starts its functioning with any microorganisms adhered, the correct election for D is a value between 0.034 and 1, for the same reasons given before. In addition, the system modeled with Haldane kinetics presents a saddle-node bifurcation and is satisfies the Sotomayor theorem.The bifurcation diagrams show the same behavior for both models as the parameterD is varied. The saddle-node bifurcation predicted by the Haldane kinetics is not reflected in the bifurcation diagrams; this is because for values of the dilution rate for which the transcritical bifurcation appears the equilibria that generate the saddle-node bifurcation are unstable, and thus, the system trajectories will not tend toward them. For values of D before the bifurcation, the trajectories will tend toward one of the stable equilibrium points, depending on the initial conditions.It was also observed for both models that a higher dilution rate increases the formation of biofilms in the reactor. That is, the system stabilizes for a larger concentration of adhered microorganisms as the concentration of suspended microorganisms is reduced. The behavior for the system modeled by both Monod and Haldane kinetics is similar. After a certain period, the system stabilizes for approximately the same substrate concentrations and the same concentrations of suspended and adhered biomass. The difference lies in that, for the system modeled with Monod kinetics, the necessary time for biofilm formation is shorter.



---

*Source: 512404-2015-06-15.xml*