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

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.


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][10][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][10][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][14][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.

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  ∈ [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 surfacea fraction  of the daughter cells finds adhesion sites and a fraction 1 −  does not.The function  decreases with the number of daughter cells (as the number of daughter cells increases, the number of available sites to occupy decreases), where  = (1 − )/(1.1 − ), which is a function of probability [12].
The reactor dynamics are described by the following equations: The following parameters are used for the numerical simulation: ( 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.

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 () =   /(  + ) where   is the maximum growth rate.  is the concentration of the substrate corresponding to the half of the maximum growth rate of () and represents the affinity of the microorganisms for the substrate.If the organism has a great affinity for the limiting substrate, the value of   is low [17].

Equilibrium Points.
A point  0 ∈ R  is an equilibrium point or a critical point of ẋ = f() if f( 0 ) = 0.In this point, the system remains in a stationary state; that is, if   :  → R  is the flux of the system, then   ( 0 ) =  0 ∀ ∈ 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.

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 point  1 , the Jacobian matrix of the system is Let  be a submatrix, such that The eigenvalues  1 ,  2 , and  3 of the matrix   1 are the roots of the equation: It is readily apparent that Substituting the values of the parameters, one obtains that ( 0 )(0) −  −  > 0. Therefore,  1 is unstable.
For the equilibrium point  2 , numerically performing the stability analysis, we find that it is always stable.
The algorithm is as follows: for a given value of  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 point  3 , the Jacobian matrix of the system is and its characteristic polynomial is We establish the following array: If  < 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 <  < 0.093448 there is always a sign change ( 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 point  4 , the characteristic polynomial is  3 + By virtue of the Routh-Hurwitz criterion, one has Numerically, it is found that, for values  < 0.034, there is a sign change, and, thus, the equilibrium point is unstable, whereas for values  > 0.034, there is no sign change, and, thus, the equilibrium point is stable.

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 a bifurcation.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.
Let  =  0 be a nonhyperbolic equilibrium point.The structure of the system is characterized as follows.
(i) Two curves of equilibrium points () and () pass through  0 .(ii) Both curves exist on both sides  0 .(iii) The stability of the curves is interchanged as they pass through  0 .
The system presents a transcritical bifurcation in the equilibrium curves  3 and  4 in fact.
These results can be used at the moment to choose an efficient dilution coefficient ; that is, if the reactor starts its functioning without microorganisms adhered, the correct election for  is a value less than 0.034 h −1 , because in this interval for  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  is a value between 0.034 h −1 and 1 h −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: where f is smooth with respect to  and .Let  =  0 be an equilibrium point with a vanishing eigenvalue at  =  0 .Assume that f( 0 ,  0 ) = 0 and the matrix  = f( 0 ,  0 ) has an eigenvalue  = 0 with eigenvector k and   has an eigenvector w corresponding to the eigenvalue  = 0. Furthermore, assume that  has  eigenvalues with a negative real part and that −−1 with a positive real part.The system presents the following.
A saddle-node bifurcation if A transcritical bifurcation if The nature of the expressions above is as follows: The modeled system presents a transcritical bifurcation, although the second condition of the Sotomayor theorem is not satisfied.In fact with V 3 free constant.
Condition 1.We have that Condition 2. We have that 2.2.4.Bifurcations Diagram.Bifurcations are structural changes in the system.These types of changes can be appreciated by means of their corresponding bifurcation diagrams as it determines if the system converges to an -cycle or if its behavior is random.Figure 3 shows that for values  < 0.034 the trajectories converge at the equilibrium point  3 ; 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  > 0.034 the trajectories converge at the equilibrium point  4 , where there is biofilm formation.This coincides with the stability curves  3 and  4 shown in Figures 1 and 2. Interval partition = 3000 points.

Bifurcation Diagram Algorithm
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]).

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  =   /(  +  +  2 /  ), where   is defined in an analogous way as   ; in this case, however, it represents the value of [] when the inhibition is at its half maximum.
Figure 4 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   [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   = 1500.The value of the constant is a characteristic typical of inhibitor.In this case when  =   the degree of inhibition is about 50%.

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.
> 0 at the equilibrium point From the above,  3 has physical meaning if 0.0704 <  and  < 0.0721,  4 has physical meaning if 0 <  < 0.0721, if   ̸ = 0 there is an equilibrium point,  5 = (,   ,   ), where Substituting   and   in the following equation, we can solve .Consider

Equilibrium Points Stability.
For equilibrium point  1 the stability analysis is analogous to the case of Monod kinetics, and thus, it can be immediately concluded that  1 is also a saddle point.For equilibrium point  2 the numerical stability analysis shows that it is always stable.
For equilibrium points  3 and  4 , we analytically evaluated the stability of the values of  where the equilibria are physically possible.Conducting a stability analysis analogous to the case of Monod kinetics, it is concluded that if  < 0.034 then the equilibrium point  4 is stable, whereas for 0.034 <  < 0.0721 the equilibrium points  3 and  4 are unstable.
For equilibrium point  5 the process is analogous to point  4 of the previous section.It is observed that for values  < 0.034 there is a sign change such that the equilibrium is unstable, whereas for values  > 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 coefficient ; that is, if the reactor starts its functioning without microorganisms adhered, the correct election for  is a value less than 0.034, because in this interval for  the reactor is stable.In the case that the reactor starts its functioning with any microorganisms adhered, the correct election for  is a value between 0.034 and 1, for the same reasons given before.

One-Parameter Bifurcations
Transcritical Bifurcation.When  = 0, for points  4 and  5 there is a transcritical bifurcation when 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.
Let  =  0 be a nonhyperbolic equilibrium point, and let  be the bifurcation parameter.An eigenvalue  = 0 occurs at  =  0 .The structure of the system is characterized as follows.
For values  <  0 there are two equilibrium points.
For values  =  0 there is a unique nonhyperbolic equilibrium point.
For values  >  0 there are no equilibrium points.
There are two values for the parameter  for which   () = 0: The dynamics are as follows.
For  < 0.0721 there are two physically sensible equilibrium points,  3 and  4 .
For  > 0.213 the points  3 and  4 are restored, although they do not possess physical meaning.
Condition 1.We have that Condition 2.  = (,   ,   ).  are state variables,  1 = ,  2 =   , and  3 =   .Consider The system satisfied the Sotomayor theorem, meeting both conditions, and thus, it presents a saddle-node bifurcation at the two values of the parameters ,  1 , and  2 as described above.

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

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  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  is a value less than 0.034, because in this interval for  the reactor is stable; in the case that the reactor starts its functioning with any microorganisms adhered, the correct election for  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 parameter  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  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.
Parameters (these parameters are taken from [12]) : D i l u t i o n f a c t o r h −1 : M i c r o o r g a n i s m m o r t a l i t y r a t e h −1  0 : Feed substrate concentration mg/L : Y i e l d c o n s t a n t : Suspended microorganism adhesion rate h −1 : Detachment rate of the adhered microorganisms h −1  , : Maximal surface biomass density of the adhered microorganisms mg/L  =   /  : Surface occupied fraction.

Figure 3 :
Figure 3: Bifurcation diagram varying  with respect to ,   , and   .The adherence rate  is zero.

Figure 4 :
Figure 4: Growth rate for Monod and Haldane kinetics.