Next Article in Journal
Biodiesel from Crude Tall Oil and Its NOx and Aldehydes Emissions in a Diesel Engine Fueled by Biodiesel-Diesel Blends with Water Emulsions
Next Article in Special Issue
FPGA-Based Implementation of an Optimization Algorithm to Maximize the Productivity of a Microbial Electrolysis Cell
Previous Article in Journal
Special Issue on “Recent Advances in Population Balance Modeling”
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Global Stability Analysis of a Bioreactor Model for Phenol and Cresol Mixture Degradation

by
Neli Dimitrova
1,*,† and
Plamena Zlateva
1,2,†
1
Institute of Mathematics and Informatics, Bulgarian Academy of Sciences, Acad. G. Bonchev Str. Block 8, 1113 Sofia, Bulgaria
2
Institute of Robotics, Bulgarian Academy of Sciences, Acad. G. Bonchev Str. Block 2, 1113 Sofia, Bulgaria
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Processes 2021, 9(1), 124; https://doi.org/10.3390/pr9010124
Submission received: 18 November 2020 / Revised: 4 January 2021 / Accepted: 5 January 2021 / Published: 8 January 2021
(This article belongs to the Special Issue Mathematical Modeling and Control of Bioprocesses)

Abstract

:
We propose a mathematical model for phenol and p-cresol mixture degradation in a continuously stirred bioreactor. The model is described by three nonlinear ordinary differential equations. The novel idea in the model design is the biomass specific growth rate, known as sum kinetics with interaction parameters (SKIP) and involving inhibition effects. We determine the equilibrium points of the model and study their local asymptotic stability and bifurcations with respect to a practically important parameter. Existence and uniqueness of positive solutions are proved. Global stabilizability of the model dynamics towards equilibrium points is established. The dynamic behavior of the solutions is demonstrated on some numerical examples.

1. Introduction

Organic chemical mixtures are among the most persistent environmental pollutants. Different aromatic compounds such as phenol, cresols, nitrophenols, benzene, etc. coexist as complex mixtures in wastewaters from petroleum refineries, coal mining and other industrial chemical sources [1]. Biological degradation has recently become a viable technology for remediation of organic pollutants as an alternative to the traditional physical and chemical methods that can be costly and produce hazardous products. Most of the current research has been directed to the isolation and study of microbial species with high-degradation activity and capabilities of degrading chemical compounds. The review paper [2] reports on hundreds of isolated bacteria capable of degrading aromatic compounds, among them different strains of Aspergillus awamori, Arthrobacter, Burkholderia, Mycobacterium, Pseudomonas, Rhodococcus, Staphylococcus, Trametes hirsute etc. The biodegradation of one or all chemical components depends on the composition of the particular mixture and the used microorganisms [3,4,5]. The adequate analysis of interactions between the compounds and their influence on microbial growth is very important for understanding the simultaneous metabolism of phenolic mixtures [6].
Most research on microbial potentials to degrade chemical pollutants has been performed on a laboratory scale. Based on batch processes various mathematical biodegradation kinetic models have been recently developed and widely used. Among them are Monod’s, Haldane’s (known also as Andrews), sum kinetic models, sum kinetics with interaction parameter (SKIP) models, etc. [7,8]. It is known that Monod’s and Haldane’s models are appropriate for single substrate utilization. The Monod model describes the biodegradation rate in dependence of the biomass concentration. When a substrate inhibits its own degradation then Haldane’s model is more appropriate. In [9] the Haldane equation modified with a Monod-like switching function is proposed and applied to the biological removal of mixtures of phenolic compounds in sequential batch bioreactors. In [10] the aerobic biodegradability of phenol, resorcinol and 5-methylresorcinol and their different two-component mixtures is investigated and various kinetic models are tested to obtain the best curve fit.
In the case when a mixture of two or more substrates occurs, the sum kinetic and SKIP models predict better the outcome of biodegradation experiments. The latter have been proposed for the first time in [11] and widely used by many researchers. The (no-interaction) sum kinetics model for cell growth is usually represented as a sum of the specific growth rates on each substrate, e.g., as a sum of Monod- and/or Haldane-type specific growth rates. These models were evaluated in [12,13] for biodegradation of benzene, toluene and phenol mixtures using Pseudomonas putida F1 and Burkholderia sp. strain JS150 and found that the interactions between these substrates could not be described by sum kinetics models. On the contrary, the SKIP model predicts better the outcome of the mixed-culture experiments. This is due the fact that the SKIP models extend the sum kinetics models by incorporating interaction parameters to describe more accurately the biodegradation of the chemical mixture.
The biodegradation of benzene, toluene and phenol is studied in [14] by adaptation of Pseudomonas putida F1 ATCC 700007. For the substrate mixtures, a SKIP model is used. The latter provides an excellent prediction of the growth kinetics and the interactions between these substrates.
In [15] biodegradation kinetics of different multiple substrate mixtures of mono-aromatic volatile organic carbon (VOCs) such as toluene, ethyl benzene and o-xylene are studied. A general mixed-substrate biodegradation model is developed which can describe the biodegradation kinetics of common industrial VOCs when present as a mixture, incorporating parameters for interaction effects.
The paper [16] examines biodegradation kinetics of styrene and ethylbenzene, independently and as binary mixtures, using a series of aerobic batch degradation. The SKIP model and the purely competitive enzyme kinetics model are employed to evaluate any interactions. The SKIP model is found to more accurately describe the interactions.
Here, we propose a mathematical model for biodegradation of phenol and 4-methylph- enol (p-cresol) in a continuously stirred tank bioreactor, in which the biodegradation kinetics is described by a SKIP model. The bioreactor model presents an extension of the growth kinetic model proposed in [17]. There, the growth behavior and degradation capacity of Aspergillus awamori NRRL 3112 microbial strain on the binary mixture phenol/p-cresol are investigated. Based on laboratory experiments, the growth kinetic model is first evaluated by a sum kinetic model involving Haldane’s specific growth rate. An alternative model is then formulated by adding interaction parameters into the sum kinetics model to produce the SKIP model. It is shown that the SKIP model describes better the degradation patterns in the biological system.
The paper is organized as follows. The next Section 2 presents a short description of the proposed mathematical model. Section 3 includes steady states computations. Local stability analysis and bifurcations of the equilibrium points are presented in Section 4. Section 5 reports on general and important properties of the model solutions and provides results on the global stabilizability of the system towards an interior equilibrium point. The last Section 6 presents numerical examples as illustration of the theoretical studies on the model dynamics.

2. Model Description

We consider the following mathematical model for phenol and p-cresol mixture degradation in a continuously stirred bioreactor
d X ( t ) d t = μ ( S p h , S c r ) D X ( t )
d S p h ( t ) d t = k p h μ ( S p h , S c r ) X ( t ) + D ( S p h 0 S p h ( t ) )
d S c r ( t ) d t = k c r μ ( S p h , S c r ) X ( t ) + D ( S c r 0 S c r ( t ) ) ,
where μ ( S p h , S c r ) is the specific growth rate, presented by
μ ( S p h , S c r ) = μ m a x ( p h ) S p h k s ( p h ) + S p h + S p h 2 k i ( p h ) + I c r / p h S c r + μ m a x ( c r ) S c r k s ( c r ) + S c r + S c r 2 k i ( c r ) + I p h / c r S p h .
The definition of the state variables X, S p h and S c r as well as of the model parameters is given in Table 1. The numerical values in the last column are validated by laboratory experiments and given in [17].
The specific growth rate μ ( S p h , S c r ) represents a SKIP (sum kinetics with interaction parameters) model for biological degradation of the chemical compounds. The interaction parameters I c r / p h and I p h / c r indicate the degree to which substrate p-cresol affects the biodegradation of substrate phenol, and substrate phenol affects the biodegradation of substrate p-cresol, respectively. The larger value of I c r / p h (see Table 1) indicates that p-cresol inhibits the utilization of phenol much more than phenol inhibits the utilization of p-cresol. The kinetic function μ ( S p h , S c r ) also involves inhibition terms S p h 2 k i ( p h ) and S c r 2 k i ( c r ) for cell growth on phenol and p-cresol, respectively. Obviously, μ ( S p h , 0 ) and μ ( 0 , S c r ) are the well-known Andrews (or Haldane) model functions, which are unimodal and achieve their maximum at S p h = k s ( p h ) k i ( p h ) and S c r = k s ( c r ) k i ( c r ) respectively.
The influent concentrations S p h 0 , S c r 0 and the dilution rate D are the parameters that can be manipulated by the operator of the bioreactor. In our analysis we assume that S p h 0 and S c r 0 are constant and consider the dilution rate D as a varying control parameter. Clearly, D > 0 should be fulfilled.
The same model (1)–(3) has been considered in [18] using a more simple specific growth rate function μ ( S p h , S c r ) which does not involve the inhibition terms S p h 2 k i ( p h ) and S c r 2 k i ( c r ) for cell growth on phenol and p-cresol. Adding these terms makes the dynamics (1)–(3) more complicated, but as shown in [17], see also [5], the SKIP model (4) describes the trend of experimental data much better than other kinetic models.

3. Existence of Equilibrium Points

We shall investigate existence of the model equilibrium points in dependence of the control parameter D.
The equilibrium points of (1)–(3) are solutions of the following system of algebraic equations
μ ( S p h , S c r ) D X = 0
k p h μ ( S p h , S c r ) X + D ( S p h 0 S p h ) = 0
k c r μ ( S p h , S c r ) X + D ( S c r 0 S c r ) = 0 .
Obviously, the point E 0 = ( 0 , S p h 0 , S c r 0 ) (with X = 0 ) is an equilibrium point of the model for all D > 0 .
We are looking now for solutions of (5)–(7) assuming that X 0 .
After multiplying Equation (6) by k c r , Equation (7) by k p h and summing the latter, we obtain
k c r ( S p h 0 S p h ) + k p h ( S c r 0 S c r ) = 0 .
Let us express S p h from (8) as a function of S c r . Denoting
K = k p h k c r , S 0 = S p h 0 K S c r 0 ,
We obtain
S p h = S p h 0 k p h k c r S c r 0 S c r = S 0 + K S c r .
After replacing the latter presentation of S p h into the equation μ ( S p h , S c r ) = D from (5) we obtain an equation with respect to S c r of the form
μ ( S 0 + K S c r , S c r ) = D ,
or equivalently
μ m a x ( p h ) S 0 + K S c r k s ( p h ) + S 0 + K S c r + 1 k i ( p h ) ( S 0 + K S c r ) 2 + I c r / p h S c r + μ m a x ( c r ) S c r k s ( c r ) + S c r + 1 k i ( c r ) S c r 2 + I p h / c r ( S 0 + K S c r ) = D .
Straightforward calculations lead to a polynomial equation of the form
A 1 S c r 4 + A 2 S c r 3 + A 3 S c r 2 + A 4 S c r + A 5 = 0 ,
where
A 1 = D · 1 k i ( c r ) · 1 k i ( p h ) ; A 2 = μ m a x ( p h ) K k i ( c r ) + μ m a x ( c r ) 1 k i ( p h ) D ( 1 + I p h / c r K ) 1 k i ( p h ) + 1 k i ( c r ) I c r / p h + K + 2 K S 0 1 k i ( p h ) ; A 3 = μ m a x ( p h ) S 0 1 k i ( c r ) + K ( 1 + I p h / c r K ) + μ m a x ( c r ) K + I c r / p h + 2 K S 0 1 k i ( p h ) D 1 k i ( p h ) ( k s ( c r ) + I p h / c r S 0 ) + ( 1 + I p h / c r K ) K + I c r / p h + 2 K S 0 1 k i ( p h ) + 1 k i ( c r ) k s ( p h ) + S 0 + 1 k i ( p h ) S 0 2 ; A 4 = μ m a x ( p h ) S 0 ( 1 + I p h / c r K ) + K ( k s ( c r ) + I p h / c r S 0 ) + μ m a x ( c r ) k s ( p h ) + S 0 + 1 k i ( p h ) S 0 2 D ( k s ( c r ) + I p h / c r S 0 ) K + I c r / p h + 2 K S 0 1 k i ( p h ) + ( 1 + I p h / c r K ) k s ( p h ) + S 0 + 1 k i ( p h ) S 0 2 ; A 5 = μ m a x ( p h ) S 0 D k s ( p h ) + S 0 + 1 k i ( p h ) S 0 2 ( k s ( c r ) + I p h / c r S 0 ) .
All coefficients A i , i = 1 , 2 , , 5 , depend on the parameter D.
Obviously, if A 5 = 0 , then Equation (11) possesses a solution S c r = 0 . We have
A 5 = 0 D = D c r : = S 0 μ m a x ( p h ) k i ( p h ) S 0 2 + k i ( p h ) S 0 + k i ( p h ) k s ( p h ) = μ ( S 0 , 0 ) .
The latter value of D c r is biologically reasonable only if S 0 > 0 . Using the numerical values of the model coefficients in the last column of Table 1, we obtain
S 0 = S p h 0 K S c r 0 0.09483 > 0 ,
and so,
D c r = μ ( S 0 , 0 ) 0.09933 .
This means that for D = D c r there exists an equilibrium point with S c r = 0 . Further from (10) we compute the component of S p h = S 0 , and from (7) we get the corresponding component of X = S c r 0 k c r . Thus, at D = D c r there exists a steady state
E 1 = E 1 ( D c r ) = S c r 0 k c r , S 0 , 0 = ( 0.05172 , 0.09483 , 0 ) .
Considering the cubic equation A 1 S c r 3 + A 2 S c r 2 + A 3 S c r + A 4 = 0 at D = D c r (i.e., with A 5 = 0 ), numerical computations produce the following roots of the latter equation
4.484933737 , 0.2614282531 ± i 0.2468184467 ,
so, the real root is negative and cannot serve as a component of the model equilibrium point.
If D D c r then Equation (11) may possess up to 4 real positive solutions with respect to S c r . If there exists at least one positive solution of (11), say S c r * , such that S c r * < S c r 0 for some values of D, we shall have an interior (with positive components) equilibrium of the form
E * = ( X * , S p h * , S c r * ) , S p h * = S 0 + K S c r * < S p h 0 , X * = S c r 0 S c r * k c r = S p h 0 S p h * k p h .
Remark 1. 
If we express S c r from (8) as a function of S p h and denote K ^ = k c r k p h = 1 K , S ^ 0 = S c r 0 K ^ S p h 0 , then we shall have S c r = S ^ 0 + K ^ S p h . Similar calculations as above will produce a polynomial equation of the form A ^ 1 S p h 4 + A ^ 2 S p h 3 + A ^ 3 S p h 2 + A ^ 4 S p h + A ^ 5 = 0 , where the coefficients A ^ i are similar to A i , i = 1 , 2 , , 5 , within S ^ 0 and K ^ instead of S 0 and K, respectively. In this case we have
A ^ 5 = μ m a x , c r D k s ( c r ) + S ^ 0 + S 0 ^ 2 k i ( c r ) ( k s ( p h ) + I c r / p h S ^ 0 ) .
Obviously, A ^ 5 = 0 at D ^ = μ ( 0 , S ^ 0 ) . But in this case S ^ 0 = 1 K S 0 0.047 < 0 , thus there is no value of D at which S p h = 0 is a root of the polynomial i = 1 5 A ^ i S p h 5 i = 0 . As we shall see in the following, this is the case with the equilibrium component S p h .
Numerical computations show that if D > D c r then there are no positive real roots of Equation (11). Therefore, we can expect interior (coexistence) equilibria of the form E * if D ( 0 , D c r ) , in case that the equilibrium components with respect to S c r satisfy the inequality S c r S c r 0 . Further we obtain numerically the following results:
  • There exists a value D = D c r ( 1 ) 0.0745599 , so that Equation (11) possesses a double root S c r 0.04327 for D = D c r ( 1 ) .
  • If D < D c r ( 1 ) then there are no positive roots of (11) which are less than or equal to S c r 0 .
  • Denote D c r ( 2 ) : = μ ( S p h 0 , S c r 0 ) 0.08651 < D c r . If D D c r ( 1 ) , D c r ( 2 ) then there are two positive roots of (11) which are less than S c r 0 .
  • If D D c r ( 2 ) , D c r , D c r = μ ( S 0 , 0 ) 0.09933 , then there is only one positive root of (11) which is less than S c r 0 .
The left plot in Figure 1 shows the graph of the function μ ( S 0 + K S c r , S c r ) for S c r [ 0 , S c r 0 ] = [ 0 , 0.3 ] ; the horizontal dash lines correspond to the values of D c r ( 1 ) , D c r ( 2 ) and D c r .
Therefore, the model (1)–(3) possesses two interior equilibrium points depending on the values of D. Denote them by
E 2 = E 2 ( D ) = X ( 2 ) , S p h ( 2 ) , S c r ( 2 ) , D D c r ( 1 ) , D c r ; E 3 = E 3 ( D ) = X ( 3 ) , S p h ( 3 ) , S c r ( 3 ) , D D c r ( 1 ) , D c r ( 2 ) , w i t h S c r ( 3 ) > S c r ( 2 ) .
Numerical computations also produce the following results:
E 2 ( D c r ) = E 1 = ( 0.05172 , 0.09483 , 0 ) , D c r = μ ( S 0 , 0 ) = 0.09933 ; E 2 ( D c r ( 1 ) ) = E 3 ( D c r ( 1 ) ) = ( 0.04426 , 0.18211 , 0.04327 ) , D c r ( 1 ) = 0.0745599 ; E 3 ( D c r ( 2 ) ) = E 0 = ( 0 , S p h 0 , S c r 0 ) = ( 0 , 0.7 , 0.3 ) , D c r ( 2 ) = μ ( S p h 0 , S c r 0 ) = 0.08651 .
Figure 1 (right plot) and Figure 2 visualize the components S c r , S p h and X of the equilibria E 0 , E 2 and E 3 . In the three plots, the components of the equilibrium point E 0 are marked by horizontal dash-dot&solid lines, the components of E 2 are marked by dash lines and the ones of E 3 are shown by solid lines.

4. Local Stability of the Equilibrium Points

In this section we shall study the conditions for local asymptotic stability of the model equilibrium points.
It is well known that an equilibrium point is locally asymptotically stable, if all eigenvalues of the Jacobi matrix evaluated at this equilibrium have negative real parts, cf. e.g., [19]. The eigenvalues of the Jacobi matrix coincide with the roots of the corresponding characteristic polynomial.
To simplify notations, in the following we shall sometimes write μ instead of μ ( S p h , S c r ) . The Jacobi matrix J related to the model Equations (1)–(3) has the form
J = μ ( S p h , S c r ) D μ S p h X μ S c r X k p h μ ( S p h , S c r ) k p h μ S p h X D k p h μ S c r X k c r μ ( S p h , S c r ) k c r μ S p h X k c r μ S c r X D .
The characteristic polynomial corresponding to J is defined by d e t ( J λ I 3 ) , where λ is any complex number and I 3 is the ( 3 × 3 ) –identity matrix
d e t ( J λ I 3 ) = μ ( S p h , S c r ) D λ μ S p h X μ S c r X k p h μ ( S p h , S c r ) k p h μ S p h X D λ k p h μ S c r X k c r μ ( S p h , S c r ) k c r μ S p h X k c r μ S c r X D λ .
Multiplying the second row of the above determinant by k c r k p h and adding the latter to the third row, we obtain
d e t ( J λ I 3 ) = μ ( S p h , S c r ) D λ μ S p h X μ S c r X k p h μ ( S p h , S c r ) k p h μ S p h X D λ k p h μ S c r X 0 k c r k p h ( D + λ ) D λ .
Straightforward calculations deliver the following characteristic polynomial
d e t ( J λ I 3 ) = ( D + λ ) 2 μ ( S p h , S c r ) D λ X k p h μ S p h + k c r μ S c r .
Denote by J ( E i ) the Jacobian matrix evaluated at the equilibrium point E i , i = 0 , 1 , 2 , 3 . It follows from (15) that λ 1 , 2 = D < 0 are always eigenvalues of J ( E i ) , i = 0 , 1 , 2 , 3 . The third eigenvalue λ 3 is determined from the second multiplier of (15).
Proposition 1. 
(i) 
If D < D c r ( 2 ) = μ ( S p h 0 , S c r 0 ) then the equilibrium point E 0 = 0 , S p h 0 , S c r 0 (with X = 0 ) is locally asymptotically unstable (a saddle).
(ii) 
If D > D c r ( 2 ) then E 0 is locally asymptotically stable (a stable node).
(iii) 
At D = D c r ( 2 ) the equilibrium E 0 is neither stable, nor unstable: J ( E 0 ) possesses a zero eigenvalue, λ 3 = 0 , thus D c r ( 2 ) is a bifurcation parameter value.
(iv) 
The equilibrium point E 1 = E 1 ( D c r ) = S c r 0 k c r , S 0 , 0 , (see (13)), is locally asymptotically unstable.
Proof. 
(i)(iii) We obtain from (15)
d e t ( J ( E 0 ) λ I 3 ) = ( D + λ ) 2 ( μ ( S p h 0 , S c r 0 ) D λ ) ,
thus the third root λ 3 satisfies
λ 3 = μ ( S p h 0 , S c r 0 ) D > 0 , i f D < D c r ( 2 ) = μ ( S p h 0 , S c r 0 ) , < 0 , i f D > D c r ( 2 ) , = 0 , i f D = D c r ( 2 ) .
(iv) The characteristic polynomial corresponding to the equilibrium E 1 is presented by
d e t ( J ( E 1 ) λ I 3 ) = ( D c r + λ ) 2 S c r 0 K μ S p h ( S 0 , 0 ) + μ S c r ( S 0 , 0 ) + λ .
The third root λ 3 of the latter polynomial is computed numerically and is equal to
λ 3 = S c r 0 K μ S p h ( S 0 , 0 ) + μ S c r ( S 0 , 0 ) ( 0.7574 ) > 0 ,
which means that E 1 ( D c r ) is a saddle equilibrium point. This proves the proposition. □
The equilibrium components S p h ( i ) and S c r ( i ) of the equilibria E i , i = 2 , 3 , satisfy the equation μ ( S p h , S c r ) = D , so that from (15) we obtain
d e t ( J ( E i ) λ I 3 ) = ( D + λ ) 2 λ + X ( i ) k p h μ S p h ( S p h ( i ) , S c r ( i ) ) + k c r μ S c r ( S p h ( i ) , S c r ( i ) ) , i = 2 , 3 .
The third root λ 3 ( i ) = X ( i ) k p h μ S p h ( S p h ( i ) , S c r ( i ) ) + k c r μ S c r ( S p h ( i ) , S c r ( i ) ) is found numerically by computing the right-hand side expression on a discrete mesh of values for D, where D ( D c r ( 1 ) , D c r ) for E 2 , and D ( D c r ( 1 ) , D c r ( 2 ) ) for E 3 . Figure 3 visualizes the three eigenvalues of J ( E 2 ) and J ( E 3 ) . One can see that the eigenvalues of J ( E 3 ) are negative (right plot), and J ( E 2 ) possesses one real positive eigenvalue (left plot). Moreover, one eigenvalue of J ( E 2 ) approaches zero at D = D c r ( 1 ) , and one eigenvalue of J ( E 3 ) approaches zero at D = D c r ( 1 ) and D = D c r ( 2 ) , thus D c r ( 1 ) and D c r ( 2 ) are bifurcation parameter values.
We summarize the above results in the next proposition.
Proposition 2. 
( i )
The equilibrium E 2 , defined for D ( D c r ( 1 ) , D c r ) , is locally asymptotically unstable (a saddle).
( i i )
The equilibrium E 3 , defined for D ( D c r ( 1 ) , D c r ( 2 ) ) , is locally asymptotically stable (a stable node).
( i i i )
At D = D c r ( 1 ) , the two interior equilibrium points, E 2 and E 3 , are ’born’, thus D c r ( 1 ) is a bifurcation value of the parameter D. At D = D c r ( 1 ) the steady states E 2 and E 3 undergo a saddle-node bifurcation.
( i v )
At D = D c r ( 2 ) the equilibrium points E 3 and E 0 coalesce and exchange stability for D > D c r ( 2 ) . Thus, at D = D c r ( 2 ) the steady states E 3 and E 0 undergo a transcritical bifurcation.
Figure 1 (right plot) and Figure 2 also visualize the stability of E 0 , E 2 and E 3 : the solid lines correspond to the components of the stable equilibria, the dash and the dash-dot lines mark the components of the unstable equilibria. Therefore, these three plots can also be considered as bifurcation diagrams: a saddle node bifurcation occurs at the parameter value D = D c r ( 1 ) , and D = D c r ( 2 ) serves as a transcritical bifurcation point.

5. Global Stabilizability of the Model Dynamics

First we prove that the model (1)–(3) exhibits the standard properties that we would expect from a bioreactor model, namely uniqueness and positiveness of solutions for non-negative initial conditions.
Theorem 1. 
Consider the model (1)–(3) and assume that X ( 0 ) 0 , S p h ( 0 ) 0 , S c r ( 0 ) ) 0 .
(i) 
If X ( 0 ) = 0 then all model solutions tend to the equilibrium point E 0 = ( 0 , S p h 0 , S c r 0 ) .
(ii) 
If X ( 0 ) > 0 then X ( t ) > 0 , S p h ( t ) > 0 , S c r ( t ) > 0 for all t > 0 .
(iii) 
All solutions are uniformly bounded for all t 0 .
Proof. 
(i) Let X ( 0 ) = 0 and S p h ( 0 ) 0 , S c r ( 0 ) 0 be satisfied. It follows that X ( t ) = 0 for all t 0 due to uniqueness of solutions of the Cauchy problem. Then the model (1)–(3) reduces to
d S p h ( t ) d t = D ( S p h 0 S p h ( t ) ) d S c r ( t ) d t = D ( S c r 0 S c r ( t ) ) .
The latter equations imply that S p h ( t ) and S c r ( t ) converge exponentially to S p h 0 and S c r 0 respectively. The plane X = 0 is invariant for the model.
( i i ) ( i i i ) Assume that X ( 0 ) > 0 , S p h ( 0 ) 0 , S c r ( 0 ) ) 0 . It follows from Equation (1) that
d X X = 0 t ( μ ( S p h ( τ ) , S c r ( τ ) ) D ) d τ , X ( t ) = X ( 0 ) e 0 t ( μ ( S p h ( τ ) , S c r ( τ ) ) D ) d τ > 0   for   each   t 0 .
Denote Σ 1 ( t ) = S p h ( t ) + k p h X ( t ) S p h 0 . Then Equations (1) and (2) imply
d d t Σ 1 ( t ) = d S p h d t + k p h d X d t = D S p h 0 ( S p h + k p h X ) = D Σ 1 ( t ) ,
which means that Σ 1 ( t ) = e D t Σ 1 ( 0 ) , thus lim t Σ 1 ( t ) = 0 , or equivalently
lim t S p h ( t ) + k p h X ( t ) = S p h 0 .
Since X ( t ) > 0 for all t > 0 this means S p h ( t ) > 0 for all t > 0 as well. Moreover, X ( t ) and S p h ( t ) are uniformly bounded.
Similarly, using Equations (1) and (3) and denoting Σ 2 ( t ) = S c r ( t ) + k c r X ( t ) S c r 0 we obtain Σ 2 ( t ) = e D t Σ 2 ( 0 ) , which means that
lim t S c r ( t ) + k c r X ( t ) = S c r 0 .
Therefore, S c r ( t ) > 0 for all t > 0 and S c r ( t ) is uniformly bounded for t 0 . Hence, the model solutions X ( t ) , S p h ( t ) , S c r ( t ) exist for all time t 0 . This completes the proof of Theorem 1. □
In the following we shall prove the global asymptotic stabilizability of system (1)–(3) when the control parameter D belongs to the interval D c r ( 1 ) , D c r ( 2 ) , with D c r ( 2 ) = μ ( S p h 0 , S c r 0 ) .
Similarly to the proof of Theorem 1, denote Σ 3 ( t ) = S p h ( t ) K S c r ( t ) S 0 , where K and S 0 are defined in (9). After multiplying Equation (3) by k p h k c r and adding to Equation (2) we obtain
d d t Σ 3 ( t ) = d d t S p h ( t ) K S c r ( t ) = D S p h 0 S p h ( t ) K S c r 0 + K S c r ( t ) = D ( S p h 0 K S c r 0 ) ( S p h ( t ) K S c r ( t ) ) = D S 0 ( S p h ( t ) K S c r ( t ) ) = D Σ 3 ( t ) .
This means that Σ 3 ( t ) = e D t Σ 3 ( 0 ) , Σ 3 ( 0 ) 0 , so lim t Σ 3 ( t ) = 0 . Then system (1)–(3) may be written in the form
d d t Σ 3 ( t ) = D Σ 3 ( t ) d d t X ( t ) = μ ( S 0 + K S c r ( t ) , S c r ( t ) ) D X ( t ) d d t S c r ( t ) = k c r μ ( S 0 + K S c r ( t ) , S c r ( t ) ) X ( t ) + D ( S c r 0 S c r ( t ) ) .
Since lim t Σ 3 ( t ) = 0 , the positive ω -limit set of any solution of system (1)–(3) is contained in the set
Ω 3 = ( X , S p h , S c r ) : X > 0 , S p h 0 , S c r 0 , Σ 3 = 0 .
Using the theory of the asymptotically autonomous systems (cf. [20,21]) it follows that all trajectories forming the ω -limit set of any solution of (1)–(3) with initial conditions in Ω 3 are solutions of the following limiting system
d X ( t ) d t = μ ( S 0 + K S c r ( t ) , S c r ( t ) ) D X ( t ) d S c r ( t ) d t = k c r μ ( S 0 + K S c r ( t ) , S c r ( t ) ) X ( t ) + D ( S c r 0 S c r ( t ) ) .
We consider Equation (17) on the set
Ω 2 = ( X , S c r ) : X > 0 , S c r 0 , S 0 + K S c r 0 .
Denote for simplicity μ c r ( S c r ) = μ ( S 0 + K S c r , S c r ) . Then obviously μ c r ( S c r 0 ) = μ ( S 0 + K S c r 0 , S c r 0 ) = μ ( S p h 0 , S c r 0 ) = D c r ( 2 ) holds true.
Let us choose an arbitrary value D ¯ D c r ( 1 ) , D c r ( 2 ) , and consider the following system obtained from (17) after substituting D = D ¯ in the latter:
d X ( t ) d t = μ c r ( S c r ( t ) ) D ¯ X ( t )
d S c r ( t ) d t = k c r μ c r ( S c r ( t ) ) X ( t ) + D ¯ ( S c r 0 S c r ( t ) ) .
Let us recall, that at D = D ¯ there are two interior equilibria of the model (1)–(3),
E 2 ( D ¯ ) = ( X ( 2 ) ( D ¯ ) , S p h ( 2 ) ( D ¯ ) , S c r ( 2 ) ( D ¯ ) ) a n d E 3 ( D ¯ ) = ( X ( 3 ) ( D ¯ ) , S p h ( 3 ) ( D ¯ ) , S c r ( 3 ) ( D ¯ ) ) w i t h S c r ( 2 ) ( D ¯ ) < S c r ( 3 ) ( D ¯ ) .
Denote
E ¯ = X ¯ , S ¯ c r = X ( 3 ) ( D ¯ ) , S c r ( 3 ) ( D ¯ ) .
Obviously, E ¯ is an equilibrium point of (18) and (19).
We make the following assumption.
Assumption 1. 
There exist points S c r and S c r + such that 0 < S c r < S c r + < S c r 0 and μ c r ( S c r ) is monotone increasing for all S c r S c r , S c r + .
Assumption 1 identifies the equilibrium E ¯ with the projection of E 3 ( D ¯ ) in the plane S p h K S c r = S 0 . If we choose for S c r the S c r -component of the double root of Equation (11), i.e., S c r = S c r ( 2 ) ( D c r ( 1 ) ) = S c r ( 3 ) ( D c r ( 1 ) ) , and S c r + = S c r 0 , then μ c r ( S c r ) is monotone increasing in S c r , S c r + , see the left plot in Figure 1.
Based on the above considerations, the problem for global stabilizability of the model (1)–(3) is reduced to proving the global stabilizability of the well known basic bioreactor (chemostat) model (17), which is well studied in the literature, see e.g., [20,22,23,24] and the references therein. The next Theorem 2 is also a corollary from Theorem 2.1 in [25]. We present the proof here for reader’s convenience.
Theorem 2. 
Let Assumption 1 be fulfilled. Assume that D ¯ D c r ( 1 ) , D c r ( 2 ) . Then for any initial point ( X ( 0 ) , S c r ( 0 ) ) Ω 2 the corresponding solution of (18) and (19) converges asymptotically towards the equilibrium point E ¯ .
Proof. 
Let us fix an arbitrary initial point ( X ( 0 ) , S c r ( 0 ) ) Ω 2 .
First we shall show that there exists time T > 0 , such that S c r ( t ) < S c r 0 for all t > T . Assume that S c r ( t ) S c r 0 holds true for each t > 0 . Then we have from (19) that
d S c r ( t ) d t = k c r μ c r ( S c r ) X ( t ) + D ¯ ( S c r 0 S c r ( t ) ) < 0 .
Barbălat’s Lemma [26] implies
0 = lim t d S c r ( t ) d t = lim t k c r μ c r ( S c r ) X ( t ) + D ¯ ( S c r 0 S c r ( t ) ) ,
which leads to S c r ( t ) S c r 0 and X ( t ) 0 as t . Further we have that μ c r ( S ¯ c r ) = D ¯ < D c r ( 2 ) = μ c r ( S c r 0 ) . The continuity of μ c r ( · ) and the relation S c r ( t ) S c r 0 as t imply that there exists a number δ > 0 such that
μ c r ( S c r ( t ) ) D ¯ = μ c r ( S c r ( t ) ) μ c r ( S ¯ c r ) δ
for all sufficiently large t. Then it follows d X ( t ) d t = ( μ c r ( S c r ( t ) ) D ¯ ) X ( t ) δ X ( t ) for all sufficiently large t, which contradicts the boundedness of X ( t ) . Hence, there exists a sufficiently large T > 0 with S c r ( T ) S c r 0 . If the equality S c r ( T ) = S c r 0 holds true, then we have
d S c r d t ( T ) = k c r μ c r ( S c r ( T ) ) X ( T ) + D ¯ ( S c r 0 S c r ( T ) ) = k c r μ c r ( S c r ( T ) ) X ( T ) < 0 .
The last inequality shows that S c r ( t ) < S c r 0 for each t > T .
Let us fix an arbitrary γ 0 , ( μ c r ( S c r 0 ) μ c r ( S ¯ c r ) ) / 2 . (Note that μ c r ( S c r ) is monotone increasing.) The continuity of μ c r implies that there exists ε > 0 such that μ c r ( S ¯ c r ) + γ < μ c r ( S c r ) for each S c r S c r 0 ( 1 + k c r ) ε , S c r 0 . It follows from (16) that there exists time T ε > 0 so that X ( t ) and S c r ( t ) satisfy
S c r 0 ε < S c r ( t ) + k c r X ( t ) < S c r 0 + ε f o r e a c h t T ε .
Assume now that X ( t ¯ ) ε for some t ¯ T ε ; then we obtain from (20)
S c r 0 > S c r ( t ¯ ) S c r 0 k c r X ( t ¯ ) ε S c r 0 ( 1 + k c r ) ε ,
i.e., S c r ( t ¯ ) S c r 0 ( 1 + k c r ) ε , S c r 0 . Hence,
d d t X ( t ¯ ) = μ c r ( S c r ( t ¯ ) ) D ¯ X ( t ¯ ) = μ c r ( S c r ( t ¯ ) ) μ c r ( S ¯ c r ) X ( t ¯ ) γ X ( t ¯ ) > 0 .
It follows then that X ( t ) e ( t t ¯ ) γ X ( t ¯ ) . If there exists t 1 t ¯ such that X ( t 1 ) = ε , then at every time t 2 t 1 with X ( t 2 ) = ε we have d d t X ( t 2 ) = ( μ c r ( S c r ( t 2 ) ) D ¯ ) X ( t 2 ) γ ε > 0 . Hence there exists time T 1 > T such that X ( t ) ε for each t T 1 .
The above considerations mean that the ω -limit set of the corresponding trajectory of (18) and (19) lies in the set
{ ( X , S c r ) : X ε , 0 S c r S c r 0 } .
For X > 0 and S c r ( 0 , S c r 0 ) we define the following Lyapunov function
V = V ( X , S c r ) = X ¯ X η X ¯ η d η + S ¯ c r S c r X ¯ ( μ c r ( ξ ) D ¯ ) D ¯ ( S c r 0 ξ ) d ξ .
The derivative d d t V of V along the solutions of (18) and (19) is presented by
d d t V = X X ¯ X μ c r ( S c r ) D ¯ X + X ¯ ( μ c r ( S c r ) D ¯ ) D ¯ ( S c r 0 S c r ) k c r μ c r ( S c r ) X + D ¯ ( S c r 0 S c r ) = X ( μ c r ( S c r ) D ¯ ) 1 X ¯ k c r μ c r ( S c r ) D ¯ ( S c r 0 S c r ) = X ( μ c r ( S c r ) D ¯ ) 1 S c r 0 S ¯ c r S c r 0 S c r · μ c r ( S c r ) D ¯ 0
for each S c r ( 0 , S c r 0 ) and X > 0 . Applying LaSalle’s invariance principle it follows that each trajectory of (18) and (19) approaches the equilibrium point E ¯ , i.e., E ¯ is globally asymptotically stable. This proves the theorem. □
It follows from Propositions 1 and 2 that when the control input D takes values D > D c r ( 2 ) = μ ( S p h 0 , S c r 0 ) then the model (1)–(3) possesses two equilibrium points—the wash-out equilibrium E 0 = ( 0 , S p h 0 , S c r 0 ) and the interior equilibrium E 2 , such that E 0 is locally asymptotically stable and E 2 is locally asymptotically unstable. Using the reduced model (17) it can be shown that the restriction E ¯ 0 = ( 0 , S c r 0 ) of the wash-out equilibrium E 0 is globally asymptotically stable if D > D c r ( 2 ) = μ c r ( S c r 0 ) . Although the proof can be extracted from the more general Lemma 2.2 in [24], we present it below for completeness.
Theorem 3. 
Assume that D > D c r ( 2 ) holds true. Then for any initial point ( X ( 0 ) , S c r ( 0 ) ) > 0 the corresponding solution of (17) converges asymptotically towards the equilibrium E ¯ 0 = ( 0 , S c r 0 ) .
Proof. 
Choose some D ¯ 0 > D c r ( 2 ) and consider system (18) and (19), where D ¯ is replaced by D ¯ 0 . Assume that lim t X ( t ) = X * > 0 . Then Barbălat’s Lemma [26] applied to Equation (18) implies 0 = lim t d d t X ( t ) = lim t ( μ c r ( S c r ( t ) ) D ¯ 0 ) X * , which means that lim t μ c r ( S c r ( t ) ) = D ¯ 0 > μ c r ( S c r 0 ) . From the continuity of μ c r ( · ) it follows that there exists time T > 0 and a positive number δ such that μ c r ( S c r ( t ) ) μ c r ( S c r 0 ) δ for all t T . The latter inequality leads to d d t X ( t ) δ X ( t ) , a contradiction with the boundedness of X ( t ) . Therefore, X ( t ) 0 as t . From the theory of the asymptotically autonomous systems (cf. [20,21]) it follows that the dynamics (18) and (19) can be reduced to the limiting equation d d t S c r ( t ) = D ¯ 0 ( S c r 0 S c r ( t ) ) , which implies lim t S c r ( t ) = S c r 0 , and this completes the proof. □

6. Dynamic Behavior of the Model Solutions: Numerical Simulation

In this section we present two numerical examples that illustrate the dynamic behavior of the model solutions.
Example 1. 
D = 0.08 ( D c r ( 1 ) , D c r ( 2 ) )
In this case there exist two positive (coexistence) equilibrium points
E 2 = ( 0.0491 , 0.126 , 0.0153 ) a n d E 3 = ( 0.0318 , 0.328 , 0.116 ) ,
such that E 2 is locally asymptotically unstable, E 3 is the globally asymptotically stable equilibrium point according to Theorem 2. The wash-out equilibrium E 0 = ( 0 , 0.7 , 0.3 ) is locally asymptotically unstable.
The left plot in Figure 4 visualizes the convergence of the solutions towards the corresponding equilibrium components of E 3 using two different starting points. The right plot of Figure 4 as well as Figure 5 visualize projections of the trajectories in the phase planes ( X , S c r ) , ( X , S p h ) and ( S p h , S c r ) respectively with three different initial points, denoted by circles. The corresponding projections of the invariant planes are marked by dash lines in the three plots.
Example 2. 
D = 0.095 > D c r ( 2 ) 0.0865
In this case there exists only one interior equilibrium point E 2 = ( 0.0514 , 0.0987 , 0.00191 ) which is locally asymptotically unstable. The wash-out equilibrium E 0 = ( 0 , 0.7 , 0.3 ) is the globally asymptotically stable steady state according to Theorem 3.
The global stability of E 0 for large values of the control parameter D means total wash-out of the biomass X and thus no detoxification of the bioreactor medium.
The left plot in Figure 6 visualizes the convergence of the solutions towards the corresponding components of E 0 using two different starting points. The right plot of Figure 6 as well as Figure 7 visualize projections of the trajectories in the phase planes ( X , S c r ) , ( X , S p h ) and ( S p h , S c r ) respectively with three different initial points, marked by circles. The latter three plots also visualize the projections of the invariant planes in the corresponding phase planes.

7. Conclusions

We perform a mathematical analysis of a dynamic model, describing phenol and 4-methylphenol (p-cresol) biodegradation in a continuously stirred tank bioreactor. The model is described by three nonlinear ordinary differential equations and presents an extension of the batch growth model given in [17] to perform the ability of Aspergillus awamori strain to degrade the mixture of phenol and p-cresol. The novel idea is the usage of sum kinetic interaction parameters in the analytic expression of the microorganisms specific growth rate μ ( S p h , S c r ) in the medium, as well as inhibition terms with respect to both phenol and p-cresol concentrations. The advantages of using such kind of specific growth rates is validated by practical laboratory experiments [17], see also [5]. To our knowledge, such kind of dynamic models, describing biodegradation in continuous biorectors (chemostats), are not studied in the literature until now.
We compute the equilibrium points of the model and investigate their local asymptotic stability as well as existence of bifurcations in dependence of the input control parameter D, the dilution rate. It is shown that an equilibrium E 0 = 0 , S p h 0 , S c r 0 , corresponding to total wash-out of the biomass in the bioreactor, exists for all D > 0 . We find values of D such that two interior (coexistence) equilibria E 2 and E 3 do exist: E 3 is defined for D D c r ( 1 ) , D c r ( 2 ) , and E 2 exists if D D c r ( 1 ) , D c r . Local stability analysis shows that E 2 is locally asymptotically unstable and E 3 is locally asymptotically stable where they exist, E 0 is locally asymptotically stable for D > D c r ( 2 ) . Two types of bifurcations of the equilibria occur, a saddle node bifurcation at D = D c r ( 1 ) where E 2 and E 3 coalesce, and a transcritical bifurcation at D = D c r ( 2 ) , where E 0 coincides with E 3 and E 3 disappears for D > D c r ( 2 ) . Practically, the bifurcation values D c r ( 1 ) and D c r ( 2 ) of D should be carefully avoided, because small nearby perturbations may cause destabilization of the process, leading to total wash-out of the biomass. Most of the computations are carried out numerically due to the complicated expression of the model function μ ( · ) and the large number of model parameters. The computations are performed in the computer algebra system Maple. The most important property of the model solutions—existence, uniqueness and uniform boundedness—is established theoretically in Theorem 1. We also prove (Theorem 2) the global asymptotic stability of the interior equilibrium point E 3 when D takes values within certain bounds, D D c r ( 1 ) , D c r ( 2 ) , D c r ( 2 ) = μ ( S p h 0 , S c r 0 ) . The existence of these bounds for D is not restrictive in practical applications, since the dilution rate D is proportional to the speed of the pumping mechanism which feeds the bioreactor, thus there always exist a lower and an upper bound for D [27]. Choosing D in the interval D c r ( 1 ) , D c r ( 2 ) ensures practically long-term sustainability of the bioremediation process in the bioreactor. On the other hand, large values of the dilution rate D, D > D c r ( 2 ) = μ ( S p h 0 , S c r 0 ) , may cause total wash-out of the biomass in the reactor and may lead to process breakdown. This is due to the fact that the wash-out equilibrium E 0 = ( 0 , S p h 0 , S c r 0 ) is the global attractor of the dynamics (Theorem 3). The dynamic behavior of the model solutions is illustrated by some numerical examples for different values of the dilution rate.

Author Contributions

Conceptualization, N.D. and P.Z.; methodology, N.D.; software, N.D.; theoretical investigation, N.D.; validation, N.D. and P.Z.; data curation, P.Z.; writing—original draft preparation, N.D.; writing—review and editing, N.D. and P.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

This work has been partially supported by the National Scientific Program “Information and Communication Technologies for a Single Digital Market in Science, Education and Security (ICTinSES)”, contract No DO1–205/23.11.2018, financed by the Ministry of Education and Science in Bulgaria. The work of the first author has been partially supported by grant No BG05M2OP001-1.001-0003, financed by the Science and Education for Smart Growth Operational Program (2014–2020) in Bulgaria and co-financed by the European Union through the European Structural and Investment Funds.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Singh, P.; Jain, R.; Srivastava, N.; Borthakur, A.; Pal, D.B.; Singh, R.; Madhav, S.; Srivastava, P.; Tiwary, D.; Mishra, P.K. Current and emerging trends in bioremediation of petrochemical waste: A review. Crit. Rev. Environ. Sci. Technol. 2017, 47, 155–201. [Google Scholar] [CrossRef]
  2. Seo, J.-S.; Keum, Y.-S.; Li, Q.X. Bacterial degradation of aromatic compounds. Int. J. Environ. Res. Public Health 2009, 6, 278–309. [Google Scholar] [CrossRef] [PubMed]
  3. Sharma, N.K.; Philip, L.; Bhallamudi, S.M. Aerobic degradation of phenolics and aromatic hydrocarbons in presence of cyanide. Bioresour. Technol. 2012, 121, 263–273. [Google Scholar] [CrossRef] [PubMed]
  4. Tomei, M.C.; Annesini, M.C. Biodegradation of phenolic mixtures in a sequencing batch reactor: A kinetic study. Environ. Sci. Pollut. Res. 2008, 15, 188–195. [Google Scholar] [CrossRef] [PubMed]
  5. Yemendzhiev, H.; Zlateva, P.; Alexieva, Z. Comparison of the biodegradation capacity of two fungal strains toward a mixture of phenol and cresol by mathematical modeling. Biotechnol. Biotechnol. Equip. 2012, 26, 3278–3281. [Google Scholar] [CrossRef] [Green Version]
  6. Kietkwanboot, A.; Chaiprapat, S.; Müller, R.; Suttinun, O. Biodegradation of phenolic compounds present in palm oil mill effluent as single and mixed substrates by Trameteshirsuta AK04. J. Environ. Sci. Heal. Part A Toxic/Hazard. Subst. Environ. Eng. 2020, 55, 989–1002. [Google Scholar] [CrossRef]
  7. Liu, J.; Jia, X.; Wen, J.; Zhou, Z. Substrate interactions and kinetics study of phenolic compounds biodegradation by Pseudomonas sp. cbp1-3. Biochem. Eng. J. 2012, 67, 156–166. [Google Scholar] [CrossRef]
  8. Kumar, S.; Arya, D.; Malhotra, A.; Kumar, S.; Kumar, B. Biodegradation of dual phenolic substrates in simulated wastewater by Gliomastixindicus MTCC 3869. J. Environ. Chem. Eng. 2013, 1, 865–874. [Google Scholar] [CrossRef]
  9. Angelucci, D.M.; Annesini, M.C.; Tomei, M.C. Modelling of biodegradation kinetics for binary mixtures of substituted phenols in sequential bioreactors. Chem. Eng. Trans. 2013, 32, 1081–1086. [Google Scholar] [CrossRef]
  10. Lepik, R.; Tenno, T. Biodegradability of phenol, resorcinol and 5-methylresorcinol as single and mixed substrates by activated sludge. Oil Shale 2011, 28, 425–446. [Google Scholar] [CrossRef] [Green Version]
  11. Yoon, H.; Klinzing, G.; Blanch, H.W. Competition for mixed substrate by microbial populations. Biotechnol. Bioeng. 1977, 19, 1193–1210. [Google Scholar] [CrossRef] [PubMed]
  12. Reardon, K.F.; Mosteller, D.C.; Rogers, J.D.B. Biodegradation kinetics of benzene, toluene, and phenol as single and mixed substrates for Pseudomonas Putida F1. Biotechnol. Bioeng. 2000, 69, 385–400. [Google Scholar] [CrossRef]
  13. Reardon, K.F.; Mosteller, D.C.; Rogers, J.B.; DuTeau, N.M.; Kim, K.-H. Biodegradation kinetics of aromatic hydrocarbon mixtures by pure and mixed bacterial cultures. Environ. Health 2002, 110, 1005–1011. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Abuhamed, T.; Bayraktar, E.; Mehmetoǧlu, T.; Mehmetoǧlu, Ü. Kinetics model for growth of Pseudomonas Putida F1 Benzene, Toluene Phenol Biodegrad. Process Biochem. 2004, 39, 983–988. [Google Scholar] [CrossRef]
  15. Datta, A.; Philip, L.; Bhallamudi, S.M. Modeling the biodegradation kinetics of aromatic and aliphatic volatile pollutant mixture in liquid phase. Chem. Eng. J. 2014, 241, 288–300. [Google Scholar] [CrossRef]
  16. Hazrati, H.; Shayegan, J.; Seyedi, S.M. Biodegradation kinetics and interactions of styrene and ethylbenzene as single and dual substrates for a mixed bacterial culture. J. Environ. Health Sci. Eng. 2015, 13, 72. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  17. Yemendzhiev, H.; Gerginova, M.; Zlateva, P.; Stoilova, I.; Krastanov, A.; Alexieva, Z. Phenol and cresol mixture degradation by Aspergillus awamori strain: Biochemical and kinetic substrate interactions. Proc. ECOpole 2008, 2, 153–159. Available online: https://ecesociety.com/proceedings-of-ecopole-peco (accessed on 7 January 2021).
  18. Dimitrova, N.; Zlateva, P. Stability Analysis of a Model for Phenol and Cresol Mixture Degradation. In Proceedings of the IOP Conference Series: Earth and Environmental Science (EES), Xiamen, China, 7–9 June 2019; Volume 356, p. 012209. [Google Scholar] [CrossRef]
  19. Wiggins, S. Introduction to Applied Nonlinear Dynamical Systems and Chaos; Springer: New York, NY, USA, 1990. [Google Scholar] [CrossRef]
  20. Smith, H.L.; Waltman, P. The Theory of the Chemostat: Dynamics of Microbial Competition; Cambridge University Press: Cambridge, UK, 1995. [Google Scholar] [CrossRef]
  21. Thieme, H.R. Convergence results and a Poincaré–Bendixon trichotomy for asymptotically autonomous differential equations. J. Math. Biol. 1992, 30, 755–763. [Google Scholar] [CrossRef]
  22. Harmand, J.; Lobry, C.; Rapaport, A.; Sari, T. The Chemostat: Mathematical Theory of Microorganism Cultures; Volume 1 of Chemical Engineering Series; Wiley: Hoboken, NJ, USA, 2017. [Google Scholar] [CrossRef]
  23. Hsu, S.B. Limiting behavior of competing species. SIAM J. Appl. Math. 1978, 34, 760–763. [Google Scholar] [CrossRef] [Green Version]
  24. Wolkowicz, G.S.K.; Lu, Z. Global dynamics of a mathematical model of competition in the chemostat: General response function and differential death rates. SIAM J. Appl. Math. 1992, 52, 222–233. [Google Scholar] [CrossRef]
  25. Dimitrova, N.S.; Krastanov, M.I. Model-based optimization of biogas production in an anaerobic biodegradation process. Comput. Math. Appl. 2014, 68, 986–993. [Google Scholar] [CrossRef]
  26. Gopalsamy, K. Stability and Oscillations in Delay Differential Equations of Population Dynamics; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1992. [Google Scholar] [CrossRef]
  27. Grognard, F.; Bernard, O. Stability analysis of a wastewater treatment plant with saturated control. Water Sci. Technol. 2006, 53, 149–157. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Figure 1. (Left): graph of the function μ ( S 0 + K S c r , S c r ) for S c r [ 0 , S c r 0 ] . (Right): the equilibrium components S c r ( 2 ) (dash line) and S c r ( 3 ) (solid line), parameterized on D. The horizontal dash-dot&solid line passes trough S c r 0 . On the horizontal axis, the solid circle denotes D c r ( 1 ) , the solid box denotes D c r ( 2 ) , the diamond denotes D c r . The vertical dot line passes through D c r ( 2 ) .
Figure 1. (Left): graph of the function μ ( S 0 + K S c r , S c r ) for S c r [ 0 , S c r 0 ] . (Right): the equilibrium components S c r ( 2 ) (dash line) and S c r ( 3 ) (solid line), parameterized on D. The horizontal dash-dot&solid line passes trough S c r 0 . On the horizontal axis, the solid circle denotes D c r ( 1 ) , the solid box denotes D c r ( 2 ) , the diamond denotes D c r . The vertical dot line passes through D c r ( 2 ) .
Processes 09 00124 g001
Figure 2. (Left): the equilibrium components S p h ( 2 ) (dash line) and S p h ( 3 ) (solid line), parameterized on D. The horizontal dash-dot&solid line passes trough S p h 0 . (Right): the equilibrium components X ( 2 ) (dash line) and X ( 3 ) (solid line), parameterized on D. The horizontal dash-dot&solid line passes trough 0. On the horizontal axis (left and right plot), the solid circle denotes D c r ( 1 ) , the solid box denotes D c r ( 2 ) , the diamond denotes D c r . The vertical dot line passes through D c r ( 2 ) .
Figure 2. (Left): the equilibrium components S p h ( 2 ) (dash line) and S p h ( 3 ) (solid line), parameterized on D. The horizontal dash-dot&solid line passes trough S p h 0 . (Right): the equilibrium components X ( 2 ) (dash line) and X ( 3 ) (solid line), parameterized on D. The horizontal dash-dot&solid line passes trough 0. On the horizontal axis (left and right plot), the solid circle denotes D c r ( 1 ) , the solid box denotes D c r ( 2 ) , the diamond denotes D c r . The vertical dot line passes through D c r ( 2 ) .
Processes 09 00124 g002
Figure 3. Eigenvalues corresponding to the equilibrium points E 2 (left) and E 3 (right), parameterized on D. On the horizontal axis, the solid circle denotes D c r ( 1 ) , the solid box denotes D c r ( 2 ) , the diamond denotes D c r .
Figure 3. Eigenvalues corresponding to the equilibrium points E 2 (left) and E 3 (right), parameterized on D. On the horizontal axis, the solid circle denotes D c r ( 1 ) , the solid box denotes D c r ( 2 ) , the diamond denotes D c r .
Processes 09 00124 g003
Figure 4. D = 0.08 . (Left): time evolution of solutions X (solid line), S p h (dash-dot line) and S c r (dash line); the horizontal dot lines pass through the corresponding components of the equilibrium point E 3 . (Right): projections of the trajectories in the ( X , S c r ) -phase plane with three different initial points, denoted by circles. The corresponding equilibrium components of E 3 are marked by a solid box, of E 2 are denoted by a box. The dash line presents the a projection of invariant plane S c r + k c r X = S c r 0 .
Figure 4. D = 0.08 . (Left): time evolution of solutions X (solid line), S p h (dash-dot line) and S c r (dash line); the horizontal dot lines pass through the corresponding components of the equilibrium point E 3 . (Right): projections of the trajectories in the ( X , S c r ) -phase plane with three different initial points, denoted by circles. The corresponding equilibrium components of E 3 are marked by a solid box, of E 2 are denoted by a box. The dash line presents the a projection of invariant plane S c r + k c r X = S c r 0 .
Processes 09 00124 g004
Figure 5. D = 0.08 . Projections of the trajectories in the ( X , S p h ) -phase plane (left) and in the ( S p h , S c r ) -phase plane (right) with three different initial points, denoted by circles. The corresponding equilibrium components of E 3 are marked by solid boxes, of E 2 are denoted by boxes. The dash lines present projections of the invariant planes S p h + k p h X = S p h 0 (left) and S p h K S c r = S 0 (right).
Figure 5. D = 0.08 . Projections of the trajectories in the ( X , S p h ) -phase plane (left) and in the ( S p h , S c r ) -phase plane (right) with three different initial points, denoted by circles. The corresponding equilibrium components of E 3 are marked by solid boxes, of E 2 are denoted by boxes. The dash lines present projections of the invariant planes S p h + k p h X = S p h 0 (left) and S p h K S c r = S 0 (right).
Processes 09 00124 g005
Figure 6. D = 0.095 . (Left): time evolution of solutions X (solid line), S p h (dash-dot line) and S c r (dash line); the horizontal dot lines pass through the corresponding components of the equilibrium point E 0 . (Right): projections of the trajectories in the ( X , S c r ) -phase plane with three different initial points, denoted by circles. The corresponding equilibrium components of E 0 are marked by a solid box, of E 2 are denoted by a box. The dash line presents a projection of the invariant plane S c r + k c r X = S c r 0 .
Figure 6. D = 0.095 . (Left): time evolution of solutions X (solid line), S p h (dash-dot line) and S c r (dash line); the horizontal dot lines pass through the corresponding components of the equilibrium point E 0 . (Right): projections of the trajectories in the ( X , S c r ) -phase plane with three different initial points, denoted by circles. The corresponding equilibrium components of E 0 are marked by a solid box, of E 2 are denoted by a box. The dash line presents a projection of the invariant plane S c r + k c r X = S c r 0 .
Processes 09 00124 g006
Figure 7. D = 0.095 . Projections of the trajectories in the ( X , S p h ) -phase plane (left) and in the ( S p h , S c r ) -phase plane (right) with three different initial points, denoted by circles. The corresponding components of E 0 are marked by solid boxes, of E 2 are denoted by boxes. The dash lines present projections of the invariant planes S p h + k p h X = S p h 0 (left) and S p h K S c r = S 0 (right).
Figure 7. D = 0.095 . Projections of the trajectories in the ( X , S p h ) -phase plane (left) and in the ( S p h , S c r ) -phase plane (right) with three different initial points, denoted by circles. The corresponding components of E 0 are marked by solid boxes, of E 2 are denoted by boxes. The dash lines present projections of the invariant planes S p h + k p h X = S p h 0 (left) and S p h K S c r = S 0 (right).
Processes 09 00124 g007
Table 1. Model variables and parameters.
Table 1. Model variables and parameters.
DefinitionsValues
Xbiomass concentration [g/dm3]-
S p h phenol concentration [g/dm3]-
S c r p-cresol concentration [g/dm3]-
Ddilution rate [ h 1 ] -
S p h 0 influent phenol concentration [g/dm3]0.7
S c r 0 influent p-cresol concentration [g/dm3]0.3
k p h metabolic coefficient [ S p h / X ] 11.7
k c r metabolic coefficient [ S c r / X ] 5.8
k i ( p h ) inhibition constant for cell growth on phenol [g/dm3]0.61
k i ( c r ) inhibition constant for cell growth on cresol [g/dm3]0.45
I p h / c r interaction coefficient indicating the degree to which phenol affects the p-cresol biodegradation0.3
I c r / p h interaction coefficient indicating the degree to which p-cresol affects the phenol biodegradation8.6
μ m a x ( p h ) maximum specific growth rate on phenol as a single substrate [ h 1 ] 0.23
μ m a x ( c r ) maximum specific growth rate on p-cresol as a single substrate [ h 1 ] 0.17
k s ( p h ) saturation constant for cell growth on phenol [g/dm3]0.11
k s ( c r ) saturation constant for cell growth on p-cresol [g/dm3]0.35
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Dimitrova, N.; Zlateva, P. Global Stability Analysis of a Bioreactor Model for Phenol and Cresol Mixture Degradation. Processes 2021, 9, 124. https://doi.org/10.3390/pr9010124

AMA Style

Dimitrova N, Zlateva P. Global Stability Analysis of a Bioreactor Model for Phenol and Cresol Mixture Degradation. Processes. 2021; 9(1):124. https://doi.org/10.3390/pr9010124

Chicago/Turabian Style

Dimitrova, Neli, and Plamena Zlateva. 2021. "Global Stability Analysis of a Bioreactor Model for Phenol and Cresol Mixture Degradation" Processes 9, no. 1: 124. https://doi.org/10.3390/pr9010124

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop