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.

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 [911]. “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 [911] 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 [1315].

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 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 fraction of the daughter cells finds adhesion sites and a fraction 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 , 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.

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 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].

2.2.1. Equilibrium Points

A point is an equilibrium point or a critical point of if . In this point, the system remains in a stationary state; that is, if is the flux of the system, then [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, is .(ii)If , that is, there adherence of microorganisms to biofilms, then a nontrivial equilibrium point exists: . This equilibrium point cannot be found analytically.(iii)If , there exist two equilibrium points:with .This point is physically possible if and ; that is, , where Substituting and in the following equation, one can solve . Consider

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 point , the Jacobian matrix of the system is Let be a submatrix, such that The eigenvalues , , and of the matrix are the roots of the equation: It is readily apparent that Substituting the values of the parameters, one obtains that . Therefore, is unstable.

For the equilibrium point , 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 , the Jacobian matrix of the system is and its characteristic polynomial is We establish the following array: The 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.

is negative. Indeed given that and .

If then and .

If then .

Therefore, for , .

If then .

Therefore, for , .

In summary, we have If , there are no sign changes, and, thus, the real part of all the eigenvalues is negative, and the equilibrium point is stable. If 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 , the characteristic polynomial is , where By virtue of the Routh-Hurwitz criterion, one has Numerically, it is found that, for values , there is a sign change, and, thus, the equilibrium point is unstable, whereas for values , 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 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 be a nonhyperbolic equilibrium point. The structure of the system is characterized as follows.(i)Two curves of equilibrium points and pass through .(ii)Both curves exist on both sides .(iii)The stability of the curves is interchanged as they pass through .

The system presents a transcritical bifurcation in the equilibrium curves and in fact.

When , for the point , an eigenvalue exists if , given that and ; then must be zero. That is,

Thus The equilibrium point in the bifurcation parameter is as follows: .

The equilibrium curve is stable for and unstable for . For the equilibria are unstable for and stable for . Figures 1 and 2 show how the equilibrium curves of the points and cross each other and interchange stability, such that a transcritical bifurcation appears.

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 , 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 and , 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 is smooth with respect to and . Let be an equilibrium point with a vanishing eigenvalue at . Assume that and the matrix has an eigenvalue with eigenvector and has an eigenvector corresponding to the eigenvalue . Furthermore, assume that has eigenvalues with a negative real part and that with a positive real part. The system presents the following.A saddle-node bifurcation if(1)  ,(2)  .A transcritical bifurcation if(1)  ,(2)  ,(3)  .

The nature of the expressions above is as follows:(i)(ii) Jacobiano de (iii) donde , .

The modeled system presents a transcritical bifurcation, although the second condition of the Sotomayor theorem is not satisfied. In factwith 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 the trajectories converge at the equilibrium point ; 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 the trajectories converge at the equilibrium point , where there is biofilm formation. This coincides with the stability curves and shown in Figures 1 and 2.

Bifurcation Diagram AlgorithmInitial conditional = .The parameter values used are the same data given at the beginning except for , whose value is .Integration time = .Time step = .Integrator = ODE45.Relative tolerance = .Absolute tolerance = .Interval partition = 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 , 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 . The value of the constant is a characteristic typical of inhibitor. In this case when the degree of inhibition is about 50%.

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: .(ii)If there is a nontrivial equilibrium point: . This point cannot be found analytically.(iii)If , there are three equilibrium points.Let and . From this, we obtain exists if ; that is, with It can then be seen that exists if .When exists, that is, when is physically possible (the value of is positive), two equilibrium points are obtained: If , that is, , then given that . Therefore, for , the equilibria are not physically feasible.If , that is, , then . at the equilibrium point if

From the above, has physical meaning if and , has physical meaning if ,if there is an equilibrium point, , where Substituting and in the following equation, we can solve . Consider

2.3.2. Equilibrium Points Stability

For equilibrium point the stability analysis is analogous to the case of Monod kinetics, and thus, it can be immediately concluded that is also a saddle point.

For equilibrium point the numerical stability analysis shows that it is always stable.

For equilibrium points and , 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 then the equilibrium point is stable, whereas for the equilibrium points and are unstable.

For equilibrium point the process is analogous to point of the previous section. It is observed that for values there is a sign change such that the equilibrium is unstable, whereas for values 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 , 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 and , for the same reasons given before.

2.3.3. One-Parameter Bifurcations

Transcritical Bifurcation. When , for points and 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 be a nonhyperbolic equilibrium point, and let be the bifurcation parameter. An eigenvalue occurs at . The structure of the system is characterized as follows.

For values there are two equilibrium points.

For values there is a unique nonhyperbolic equilibrium point.

For values there are no equilibrium points.

The system presents a saddle-node bifurcation in the equilibrium curves of and . Indeed, an eigenvalue exists if , from here .

There are two values for the parameter for which :The dynamics are as follows.For there are two physically sensible equilibrium points, and .For .For the equilibria and disappear.For .For the points and are restored, although they do not possess physical meaning.

The system satisfied the Sotomayor theorem. Indeed, If then andIf then .

Condition  1. We have that

Condition  2. . are state variables, , , and . ConsiderThe system satisfied the Sotomayor theorem, meeting both conditions, and thus, it presents a saddle-node bifurcation at the two values of the parameters , , and as described above.

2.3.4. Bifurcations Diagram

The bifurcation diagram of the system for varying when is analogous to the case of Monod kinetics. This is because the equilibrium points and are stable for , and given that the initial condition is closer to the trajectories will tend to this equilibrium; that is, the suspended biomass concentration increases because there is no formation of biofilms. For the trajectories of the system tend to the stable equilibrium point .

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 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 , 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 and , 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.

Nomenclature

Variables
:Substrate concentration mg/L
:Suspended microorganism biomass concentration mg/L
:Adhered microorganism biomass concentration .
Parameters (these parameters are taken from [12])
:Dilution factor
:Microorganism mortality rate
:Feed substrate concentration mg/L
:Yield constant
:Suspended microorganism adhesion rate
:Detachment rate of the adhered microorganisms
:Maximal surface biomass density of the adhered microorganisms mg/L
:Surface occupied fraction.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This research was funded through the project “Dynamic Analysis of a Continuous-Stirred Tank Reactor (CSTR) with Biofilm Formation for the Treatment of Residual Wastewater” from the Program of Support for Postgraduate Theses DIMA 2012 of the Universidad Nacional de Colombia Sede Manizales. The authors thank professor Juan Carlos Higuita at the Universidad Nacional de Colombia Sede Manizales, for his collaboration in the phenomenological explanation of reactor.