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

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

$$\frac{dX\left(t\right)}{dt}=\left(\mu (S_{ph},S_{cr})-D\right)X\left(t\right)$$

(1)

$$\frac{d{S}_{ph}\left(t\right)}{dt}=-k_{ph}\mu (S_{ph},S_{cr})X\left(t\right)+D(S_{ph}^0-S_{ph}\left(t\right))$$

(2)

$$\frac{d{S}_{cr}\left(t\right)}{dt}=-k_{cr}\mu (S_{ph},S_{cr})X\left(t\right)+D(S_{cr}^0-S_{cr}\left(t\right)),$$

(3)

where $\mu (S_{ph},S_{cr})$ is the specific growth rate, presented by

$$\mu (S_{ph},S_{cr})={\displaystyle \frac{{\mu }_{max\left(ph\right)}{S}_{ph}}{k_{s\left(ph\right)}+S_{ph}+\frac{S_{ph}^2}{k_{i\left(ph\right)}}+I_{cr/ph}{S}_{cr}}}+{\displaystyle \frac{{\mu }_{max\left(cr\right)}{S}_{cr}}{k_{s\left(cr\right)}+S_{cr}+\frac{S_{cr}^2}{k_{i\left(cr\right)}}+I_{ph/cr}{S}_{ph}}}.$$

(4)

The definition of the state variables X, $S_{ph}$ and $S_{cr}$ as well as of the model parameters is given in

Table 1. Model variables and parameters.

	Definitions	Values
X	biomass concentration [g/dm3]	-
$S_{ph}$	phenol concentration [g/dm3]	-
$S_{cr}$	p-cresol concentration [g/dm3]	-
D	dilution rate $\left[h^{-1}\right]$	-
$S_{ph}^0$	influent phenol concentration [g/dm3]	0.7
$S_{cr}^0$	influent p-cresol concentration [g/dm3]	0.3
$k_{ph}$	metabolic coefficient $[S_{ph}/X]$	11.7
$k_{cr}$	metabolic coefficient $[S_{cr}/X]$	5.8
$k_{i\left(ph\right)}$	inhibition constant for cell growth on phenol [g/dm3]	0.61
$k_{i\left(cr\right)}$	inhibition constant for cell growth on cresol [g/dm3]	0.45
$I_{ph/cr}$	interaction coefficient indicating the degree to which phenol affects the p-cresol biodegradation	0.3
$I_{cr/ph}$	interaction coefficient indicating the degree to which p-cresol affects the phenol biodegradation	8.6
${\mu }_{max\left(ph\right)}$	maximum specific growth rate on phenol as a single substrate $\left[h^{-1}\right]$	0.23
${\mu }_{max\left(cr\right)}$	maximum specific growth rate on p-cresol as a single substrate $\left[h^{-1}\right]$	0.17
$k_{s\left(ph\right)}$	saturation constant for cell growth on phenol [g/dm3]	0.11
$k_{s\left(cr\right)}$	saturation constant for cell growth on p-cresol [g/dm3]	0.35

. The numerical values in the last column are validated by laboratory experiments and given in [17].

The specific growth rate $\mu (S_{ph},S_{cr})$ represents a SKIP (sum kinetics with interaction parameters) model for biological degradation of the chemical compounds. The interaction parameters $I_{cr/ph}$ and $I_{ph/cr}$ 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_{cr/ph}$ (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 $\mu (S_{ph},S_{cr})$ also involves inhibition terms $\frac{S_{ph}^2}{k_{i\left(ph\right)}}$ and $\frac{S_{cr}^2}{k_{i\left(cr\right)}}$ for cell growth on phenol and p-cresol, respectively. Obviously, $\mu (S_{ph},0)$ and $\mu (0,S_{cr})$ are the well-known Andrews (or Haldane) model functions, which are unimodal and achieve their maximum at $S_{ph}=\sqrt{k_{s\left(ph\right)}{k}_{i\left(ph\right)}}$ and $S_{cr}=\sqrt{k_{s\left(cr\right)}{k}_{i\left(cr\right)}}$ respectively.

The influent concentrations $S_{ph}^0$, $S_{cr}^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_{ph}^0$ and $S_{cr}^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 $\mu (S_{ph},S_{cr})$ which does not involve the inhibition terms $\frac{S_{ph}^2}{k_{i\left(ph\right)}}$ and $\frac{S_{cr}^2}{k_{i\left(cr\right)}}$ 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

$$\left(\mu (S_{ph},S_{cr})-D\right)X=0$$

(5)

$$-k_{ph}\mu (S_{ph},S_{cr})X+D(S_{ph}^0-S_{ph})=0$$

(6)

$$-k_{cr}\mu (S_{ph},S_{cr})X+D(S_{cr}^0-S_{cr})=0.$$

(7)

Obviously, the point $E_0=(0,S_{ph}^0,S_{cr}^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\not\equiv 0$.

After multiplying Equation (6) by $-k_{cr}$, Equation (7) by $k_{ph}$ and summing the latter, we obtain

$$-k_{cr}(S_{ph}^0-S_{ph})+k_{ph}(S_{cr}^0-S_{cr})=0.$$

(8)

Let us express $S_{ph}$ from (8) as a function of $S_{cr}$. Denoting

$$K={\displaystyle \frac{k_{ph}}{k_{cr}}},S^0=S_{ph}^0-K{S}_{cr}^0,$$

(9)

We obtain

$$S_{ph}=S_{ph}^0-{\displaystyle \frac{k_{ph}}{k_{cr}}}\left(S_{cr}^0-S_{cr}\right)=S^0+K{S}_{cr}.$$

(10)

After replacing the latter presentation of $S_{ph}$ into the equation $\mu (S_{ph},S_{cr})=D$ from (5) we obtain an equation with respect to $S_{cr}$ of the form

$$\mu (S^0+K{S}_{cr},S_{cr})=D,$$

or equivalently

$$\begin{array}{c}{\displaystyle \frac{{\mu }_{max\left(ph\right)}\left(S^0+K{S}_{cr}\right)}{k_{s\left(ph\right)}+S^0+K{S}_{cr}+\frac{1}{k_{i\left(ph\right)}}{(S^0+K{S}_{cr})}^2+I_{cr/ph}{S}_{cr}}}\hfill \\ +{\displaystyle \frac{{\mu }_{max\left(cr\right)}{S}_{cr}}{k_{s\left(cr\right)}+S_{cr}+\frac{1}{k_{i\left(cr\right)}}{S}_{cr}^2+I_{ph/cr}(S^0+K{S}_{cr})}}=D.\hfill \end{array}$$

Straightforward calculations lead to a polynomial equation of the form

$$A_1{S}_{cr}^4+A_2{S}_{cr}^3+A_3{S}_{cr}^2+A_4{S}_{cr}+A_5=0,$$

(11)

where

$$\begin{array}{ccc}\hfill A_1& =& -D·\frac{1}{k_{i\left(cr\right)}}·\frac{1}{k_{i\left(ph\right)}};\hfill \\ \hfill A_2& =& {\mu }_{max\left(ph\right)}\frac{K}{k_{i\left(cr\right)}}+{\mu }_{max\left(cr\right)}\frac{1}{k_{i\left(ph\right)}}\hfill \\ & & -D\left[(1+I_{ph/cr}K)\frac{1}{k_{i\left(ph\right)}}+\frac{1}{k_{i\left(cr\right)}}\left(I_{cr/ph}+K+2K{S}^0\frac{1}{k_{i\left(ph\right)}}\right)\right];\hfill \\ \hfill A_3& =& {\mu }_{max\left(ph\right)}\left[S^0\frac{1}{k_{i\left(cr\right)}}+K(1+I_{ph/cr}K)\right]+{\mu }_{max\left(cr\right)}\left[K+I_{cr/ph}+2K{S}^0\frac{1}{k_{i\left(ph\right)}}\right]\hfill \\ & & -D\left[\frac{1}{k_{i\left(ph\right)}}(k_{s\left(cr\right)}+I_{ph/cr}{S}^0)+(1+I_{ph/cr}K)\left(K+I_{cr/ph}+2K{S}^0\frac{1}{k_{i\left(ph\right)}}\right)\right.\hfill \\ & & +\left.\frac{1}{k_{i\left(cr\right)}}\left(k_{s\left(ph\right)}+S^0+\frac{1}{k_{i\left(ph\right)}}{S^0}^2\right)\right];\hfill \\ \hfill A_4& =& {\mu }_{max\left(ph\right)}\left[S^0(1+I_{ph/cr}K)+K(k_{s\left(cr\right)}+I_{ph/cr}{S}^0)\right]\hfill \\ & & +{\mu }_{max\left(cr\right)}\left(k_{s\left(ph\right)}+S^0+\frac{1}{k_{i\left(ph\right)}}{S^0}^2\right)\hfill \\ & & -D\left[(k_{s\left(cr\right)}+I_{ph/cr}{S}^0)\left(K+I_{cr/ph}+2K{S}^0\frac{1}{k_{i\left(ph\right)}}\right)\right.\hfill \\ & & +\left.(1+I_{ph/cr}K)\left(k_{s\left(ph\right)}+S^0+\frac{1}{k_{i\left(ph\right)}}{S^0}^2\right)\right];\hfill \\ \hfill A_5& =& \left[{\mu }_{max\left(ph\right)}{S}^0-D\left(k_{s\left(ph\right)}+S^0+\frac{1}{k_{i\left(ph\right)}}{S^0}^2\right)\right](k_{s\left(cr\right)}+I_{ph/cr}{S}^0).\hfill \end{array}$$

All coefficients $A_i$, $i=1,2,\dots ,5$, depend on the parameter D.

Obviously, if $A_5=0$, then Equation (11) possesses a solution $S_{cr}=0$. We have

$$A_5=0⟺D=D_{cr}:=\frac{S^0{\mu }_{max\left(ph\right)}{k}_{i\left(ph\right)}}{{S^0}^2+k_{i\left(ph\right)}{S}^0+k_{i\left(ph\right)}{k}_{s\left(ph\right)}}=\mu (S^0,0).$$

(12)

The latter value of $D_{cr}$ 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_{ph}^0-K{S}_{cr}^0\approx 0.09483>0,$$

and so,

$$D_{cr}=\mu (S^0,0)\approx 0.09933.$$

This means that for $D=D_{cr}$ there exists an equilibrium point with $S_{cr}=0$. Further from (10) we compute the component of $S_{ph}=S^0$, and from (7) we get the corresponding component of $X={\displaystyle \frac{S_{cr}^0}{k_{cr}}}$. Thus, at $D=D_{cr}$ there exists a steady state

$$E_1=E_1\left(D_{cr}\right)=\left(\frac{S_{cr}^0}{k_{cr}},S^0,0\right)=(0.05172,0.09483,0).$$

(13)

Considering the cubic equation $A_1{S}_{cr}^3+A_2{S}_{cr}^2+A_3{S}_{cr}+A_4=0$ at $D=D_{cr}$ (i.e., with $A_5=0$), numerical computations produce the following roots of the latter equation

$$-4.484933737,0.2614282531\pm i0.2468184467,$$

so, the real root is negative and cannot serve as a component of the model equilibrium point.

If $D\ne D_{cr}$ then Equation (11) may possess up to 4 real positive solutions with respect to $S_{cr}$. If there exists at least one positive solution of (11), say $S_{cr}^{*}$, such that $S_{cr}^{*}<S_{cr}^0$ for some values of D, we shall have an interior (with positive components) equilibrium of the form

$$E^{*}=(X^{*},S_{ph}^{*},S_{cr}^{*}),S_{ph}^{*}=S^0+K{S}_{cr}^{*}<S_{ph}^0,X^{*}={\displaystyle \frac{S_{cr}^0-S_{cr}^{*}}{k_{cr}}}={\displaystyle \frac{S_{ph}^0-S_{ph}^{*}}{k_{ph}}}.$$

(14)

Remark 1. 

If we express $S_{cr}$ from (8) as a function of $S_{ph}$ and denote $\widehat{K}={\displaystyle \frac{k_{cr}}{k_{ph}}}={\displaystyle \frac{1}{K}}$, ${\widehat{S}}^0=S_{cr}^0-\widehat{K}{S}_{ph}^0$, then we shall have $S_{cr}={\widehat{S}}^0+\widehat{K}{S}_{ph}$. Similar calculations as above will produce a polynomial equation of the form ${\widehat{A}}_1{S}_{ph}^4+{\widehat{A}}_2{S}_{ph}^3+{\widehat{A}}_3{S}_{ph}^2+{\widehat{A}}_4{S}_{ph}+{\widehat{A}}_5=0$, where the coefficients ${\widehat{A}}_i$ are similar to $A_i$, $i=1,2,\dots ,5$, within ${\widehat{S}}^0$ and $\widehat{K}$ instead of $S^0$ and K, respectively. In this case we have

$${\widehat{A}}_5=\left({\mu }_{max,cr}-D\left(k_{s\left(cr\right)}+{\widehat{S}}^0+\frac{{\widehat{S^0}}^2}{k_{i\left(cr\right)}}\right)\right)(k_{s\left(ph\right)}+I_{cr/ph}{\widehat{S}}^0).$$

Obviously, ${\widehat{A}}_5=0$ at $\widehat{D}=\mu (0,{\widehat{S}}^0)$. But in this case ${\widehat{S}}^0=-\frac{1}{K}{S}^0\approx -0.047<0$, thus there is no value of D at which $S_{ph}=0$ is a root of the polynomial ${\sum }_{i=1}^5{\widehat{A}}_i{S}_{ph}^{5-i}=0$. As we shall see in the following, this is the case with the equilibrium component $S_{ph}$.

Numerical computations show that if $D>D_{cr}$ then there are no positive real roots of Equation (11). Therefore, we can expect interior (coexistence) equilibria of the form $E^{*}$ if $D\in (0,D_{cr})$, in case that the equilibrium components with respect to $S_{cr}$ satisfy the inequality $S_{cr}\le S_{cr}^0$. Further we obtain numerically the following results:

• There exists a value $D=D_{cr}^{\left(1\right)}\approx 0.0745599$, so that Equation (11) possesses a double root $S_{cr}\approx 0.04327$ for $D=D_{cr}^{\left(1\right)}$.
• If $D<D_{cr}^{\left(1\right)}$ then there are no positive roots of (11) which are less than or equal to $S_{cr}^0$.
• Denote $D_{cr}^{\left(2\right)}:=\mu (S_{ph}^0,S_{cr}^0)\approx 0.08651<D_{cr}$. If $D\in \left(D_{cr}^{\left(1\right)},D_{cr}^{\left(2\right)}\right)$ then there are two positive roots of (11) which are less than $S_{cr}^0$.
• If $D\in \left(D_{cr}^{\left(2\right)},D_{cr}\right)$, $D_{cr}=\mu (S^0,0)\approx 0.09933$, then there is only one positive root of (11) which is less than $S_{cr}^0$.

The left plot in Figure 1 shows the graph of the function $\mu (S^0+K{S}_{cr},S_{cr})$ for $S_{cr}\in [0,S_{cr}^0]=[0,0.3]$; the horizontal dash lines correspond to the values of $D_{cr}^{\left(1\right)}$, $D_{cr}^{\left(2\right)}$ and $D_{cr}$.

Therefore, the model (1)–(3) possesses two interior equilibrium points depending on the values of D. Denote them by

$$\begin{array}{ccc}\hfill E_2=E_2\left(D\right)& =& \left(X^{\left(2\right)},S_{ph}^{\left(2\right)},S_{cr}^{\left(2\right)}\right),D\in \left(D_{cr}^{\left(1\right)},D_{cr}\right);\hfill \\ \hfill E_3=E_3\left(D\right)& =& \left(X^{\left(3\right)},S_{ph}^{\left(3\right)},S_{cr}^{\left(3\right)}\right),D\in \left(D_{cr}^{\left(1\right)},D_{cr}^{\left(2\right)}\right),with{S}_{cr}^{\left(3\right)}>S_{cr}^{\left(2\right)}.\hfill \end{array}$$

Numerical computations also produce the following results:

$$\begin{array}{cc}{E}_2\left(D_{cr}\right)=E_1=(0.05172,0.09483,0),\hfill & D_{cr}=\mu (S^0,0)=0.09933;\hfill \\ E_2\left(D_{cr}^{\left(1\right)}\right)=E_3\left(D_{cr}^{\left(1\right)}\right)=(0.04426,0.18211,0.04327),\hfill & D_{cr}^{\left(1\right)}=0.0745599;\hfill \\ E_3\left(D_{cr}^{\left(2\right)}\right)=E_0=(0,S_{ph}^0,S_{cr}^0)=(0,0.7,0.3),\hfill & D_{cr}^{\left(2\right)}=\mu (S_{ph}^0,S_{cr}^0)=0.08651.\hfill \end{array}$$

Figure 1 (right plot) and Figure 2 visualize the components $S_{cr}$, $S_{ph}$ 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 $\mu$ instead of $\mu (S_{ph},S_{cr})$. The Jacobi matrix J related to the model Equations (1)–(3) has the form

$$J=\left(\begin{array}{ccc}\mu (S_{ph},S_{cr})-D\hfill & \frac{\partial \mu }{\partial S_{ph}}X\hfill & \frac{\partial \mu }{\partial S_{cr}}X\hfill \\ -k_{ph}\mu (S_{ph},S_{cr})\hfill & -k_{ph}\frac{\partial \mu }{\partial S_{ph}}X-D\hfill & -k_{ph}\frac{\partial \mu }{\partial S_{cr}}X\hfill \\ -k_{cr}\mu (S_{ph},S_{cr})\hfill & -k_{cr}\frac{\partial \mu }{\partial S_{ph}}X\hfill & -k_{cr}\frac{\partial \mu }{\partial S_{cr}}X-D\hfill \end{array}\right).$$

The characteristic polynomial corresponding to J is defined by $det(J-\lambda I_3)$, where $\lambda$ is any complex number and $I_3$ is the $(3\times 3)$–identity matrix

$$det(J-\lambda I_3)=\left|\begin{array}{ccc}\mu (S_{ph},S_{cr})-D-\lambda \hfill & \frac{\partial \mu }{\partial S_{ph}}X\hfill & \frac{\partial \mu }{\partial S_{cr}}X\hfill \\ -k_{ph}\mu (S_{ph},S_{cr})\hfill & -k_{ph}\frac{\partial \mu }{\partial S_{ph}}X-D-\lambda \hfill & -k_{ph}\frac{\partial \mu }{\partial S_{cr}}X\hfill \\ -k_{cr}\mu (S_{ph},S_{cr})\hfill & -k_{cr}\frac{\partial \mu }{\partial S_{ph}}X\hfill & -k_{cr}\frac{\partial \mu }{\partial S_{cr}}X-D-\lambda \hfill \end{array}\right|.$$

Multiplying the second row of the above determinant by $-{\displaystyle \frac{k_{cr}}{k_{ph}}}$ and adding the latter to the third row, we obtain

$$det(J-\lambda I_3)=\left|\begin{array}{ccc}\mu (S_{ph},S_{cr})-D-\lambda \hfill & \frac{\partial \mu }{\partial S_{ph}}X\hfill & \frac{\partial \mu }{\partial S_{cr}}X\hfill \\ -k_{ph}\mu (S_{ph},S_{cr})\hfill & -k_{ph}\frac{\partial \mu }{\partial S_{ph}}X-D-\lambda \hfill & -k_{ph}\frac{\partial \mu }{\partial S_{cr}}X\hfill \\ 0\hfill & \frac{k_{cr}}{k_{ph}}(D+\lambda )\hfill & -D-\lambda \hfill \end{array}\right|.$$

Straightforward calculations deliver the following characteristic polynomial

$$det(J-\lambda I_3)={(D+\lambda )}^2\left[\mu (S_{ph},S_{cr})-D-\lambda -X\left(k_{ph}\frac{\partial \mu }{\partial S_{ph}}+k_{cr}\frac{\partial \mu }{\partial S_{cr}}\right)\right].$$

(15)

Denote by $J\left(E_i\right)$ the Jacobian matrix evaluated at the equilibrium point $E_i$, $i=0,1,2,3$. It follows from (15) that ${\lambda }_{1,2}=-D<0$ are always eigenvalues of $J\left(E_i\right)$, $i=0,1,2,3$. The third eigenvalue ${\lambda }_3$ is determined from the second multiplier of (15).

Proposition 1. 

(i) 

If $D<D_{cr}^{\left(2\right)}=\mu (S_{ph}^0,S_{cr}^0)$ then the equilibrium point $E_0=\left(0,S_{ph}^0,S_{cr}^0\right)$ (with $X=0$) is locally asymptotically unstable (a saddle).

(ii) 

If $D>D_{cr}^{\left(2\right)}$ then $E_0$ is locally asymptotically stable (a stable node).

(iii) 

At $D=D_{cr}^{\left(2\right)}$ the equilibrium $E_0$ is neither stable, nor unstable: $J\left(E_0\right)$ possesses a zero eigenvalue, ${\lambda }_3=0$, thus $D_{cr}^{\left(2\right)}$ is a bifurcation parameter value.

(iv) 

The equilibrium point $E_1=E_1\left(D_{cr}\right)=\left({\displaystyle \frac{S_{cr}^0}{k_{cr}}},S^0,0\right)$, (see (13)), is locally asymptotically unstable.

Proof. 

(i)–(iii) We obtain from (15)

$$det(J\left(E_0\right)-\lambda I_3)={(D+\lambda )}^2(\mu (S_{ph}^0,S_{cr}^0)-D-\lambda ),$$

thus the third root ${\lambda }_3$ satisfies

$${\lambda }_3=\mu (S_{ph}^0,S_{cr}^0)-D\left\{\begin{array}{cc}>0,\hfill & ifD<D_{cr}^{\left(2\right)}=\mu (S_{ph}^0,S_{cr}^0),\hfill \\ <0,\hfill & ifD>D_{cr}^{\left(2\right)},\hfill \\ =0,\hfill & ifD=D_{cr}^{\left(2\right)}.\hfill \end{array}\right.$$

(iv) The characteristic polynomial corresponding to the equilibrium $E_1$ is presented by

$$det(J\left(E_1\right)-\lambda I_3)=-{(D_{cr}+\lambda )}^2\left(S_{cr}^0\left(K\frac{\partial \mu }{\partial S_{ph}}(S^0,0)+\frac{\partial \mu }{\partial S_{cr}}(S^0,0)\right)+\lambda \right).$$

The third root ${\lambda }_3$ of the latter polynomial is computed numerically and is equal to

$${\lambda }_3=-S_{cr}^0\left(K\frac{\partial \mu }{\partial S_{ph}}(S^0,0)+\frac{\partial \mu }{\partial S_{cr}}(S^0,0)\right)\approx -(-0.7574)>0,$$

which means that $E_1\left(D_{cr}\right)$ is a saddle equilibrium point. This proves the proposition. □

The equilibrium components $S_{ph}^{\left(i\right)}$ and $S_{cr}^{\left(i\right)}$ of the equilibria $E_i$, $i=2,3$, satisfy the equation $\mu (S_{ph},S_{cr})=D$, so that from (15) we obtain

$$\begin{array}{c}det(J\left(E_i\right)-\lambda I_3)=-{(D+\lambda )}^2\left[\lambda +X^{\left(i\right)}\left(k_{ph}{\displaystyle \frac{\partial \mu }{\partial S_{ph}}}(S_{ph}^{\left(i\right)},S_{cr}^{\left(i\right)})+k_{cr}{\displaystyle \frac{\partial \mu }{\partial S_{cr}}}(S_{ph}^{\left(i\right)},S_{cr}^{\left(i\right)})\right)\right],\hfill \\ i=2,3.\hfill \end{array}$$

The third root ${\lambda }_3^{\left(i\right)}=-X^{\left(i\right)}\left(k_{ph}\frac{\partial \mu }{\partial S_{ph}}(S_{ph}^{\left(i\right)},S_{cr}^{\left(i\right)})+k_{cr}\frac{\partial \mu }{\partial S_{cr}}(S_{ph}^{\left(i\right)},S_{cr}^{\left(i\right)})\right)$ is found numerically by computing the right-hand side expression on a discrete mesh of values for D, where $D\in (D_{cr}^{\left(1\right)},D_{cr})$ for $E_2$, and $D\in (D_{cr}^{\left(1\right)},D_{cr}^{\left(2\right)})$ for $E_3$. Figure 3 visualizes the three eigenvalues of $J\left(E_2\right)$ and $J\left(E_3\right)$. One can see that the eigenvalues of $J\left(E_3\right)$ are negative (right plot), and $J\left(E_2\right)$ possesses one real positive eigenvalue (left plot). Moreover, one eigenvalue of $J\left(E_2\right)$ approaches zero at $D=D_{cr}^{\left(1\right)}$, and one eigenvalue of $J\left(E_3\right)$ approaches zero at $D=D_{cr}^{\left(1\right)}$ and $D=D_{cr}^{\left(2\right)}$, thus $D_{cr}^{\left(1\right)}$ and $D_{cr}^{\left(2\right)}$ are bifurcation parameter values.

We summarize the above results in the next proposition.

Proposition 2. 

$\left(i\right)$

The equilibrium $E_2$, defined for $D\in (D_{cr}^{\left(1\right)},D_{cr})$, is locally asymptotically unstable (a saddle).

$\left(ii\right)$

The equilibrium $E_3$, defined for $D\in (D_{cr}^{\left(1\right)},D_{cr}^{\left(2\right)})$, is locally asymptotically stable (a stable node).

$\left(iii\right)$

At $D=D_{cr}^{\left(1\right)}$, the two interior equilibrium points, $E_2$ and $E_3$, are ’born’, thus $D_{cr}^{\left(1\right)}$ is a bifurcation value of the parameter D. At $D=D_{cr}^{\left(1\right)}$ the steady states $E_2$ and $E_3$ undergo a saddle-node bifurcation.

$\left(iv\right)$

At $D=D_{cr}^{\left(2\right)}$ the equilibrium points $E_3$ and $E_0$ coalesce and exchange stability for $D>D_{cr}^{\left(2\right)}$. Thus, at $D=D_{cr}^{\left(2\right)}$ 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_{cr}^{\left(1\right)}$, and $D=D_{cr}^{\left(2\right)}$ 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\left(0\right)\ge 0$, $S_{ph}\left(0\right)\ge 0$, $S_{cr}\left(0\right))\ge 0$.

(i) 

If $X\left(0\right)=0$ then all model solutions tend to the equilibrium point $E_0=(0,S_{ph}^0,S_{cr}^0)$.

(ii) 

If $X\left(0\right)>0$ then $X\left(t\right)>0$, $S_{ph}\left(t\right)>0$, $S_{cr}\left(t\right)>0$ for all $t>0$.

(iii) 

All solutions are uniformly bounded for all $t\ge 0$.

Proof. 

(i) Let $X\left(0\right)=0$ and $S_{ph}\left(0\right)\ge 0$, $S_{cr}\left(0\right)\ge 0$ be satisfied. It follows that $X\left(t\right)=0$ for all $t\ge 0$ due to uniqueness of solutions of the Cauchy problem. Then the model (1)–(3) reduces to

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d{S}_{ph}\left(t\right)}{dt}}& =& D(S_{ph}^0-S_{ph}\left(t\right))\hfill \\ \hfill {\displaystyle \frac{d{S}_{cr}\left(t\right)}{dt}}& =& D(S_{cr}^0-S_{cr}\left(t\right)).\hfill \end{array}$$

The latter equations imply that $S_{ph}\left(t\right)$ and $S_{cr}\left(t\right)$ converge exponentially to $S_{ph}^0$ and $S_{cr}^0$ respectively. The plane $X=0$ is invariant for the model.

$\left(ii\right)$–$\left(iii\right)$ Assume that $X\left(0\right)>0$, $S_{ph}\left(0\right)\ge 0$, $S_{cr}\left(0\right))\ge 0$. It follows from Equation (1) that

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle {\displaystyle \frac{dX}{X}}={\int }_0^t(\mu (S_{ph}\left(\tau \right),S_{cr}\left(\tau \right))-D)d\tau ,}\hfill \\ X\left(t\right)=X\left(0\right)e^{{\int }_0^t(\mu (S_{ph}\left(\tau \right),S_{cr}\left(\tau \right))-D)d\tau }>0 \mathrm{for} \mathrm{each} t\ge 0.\hfill \end{array}\end{array}$$

Denote ${\Sigma }_1\left(t\right)=S_{ph}\left(t\right)+k_{ph}X\left(t\right)-S_{ph}^0$. Then Equations (1) and (2) imply

$${\displaystyle \frac{d}{dt}}{\Sigma }_1\left(t\right)={\displaystyle \frac{d{S}_{ph}}{dt}}+k_{ph}{\displaystyle \frac{dX}{dt}}=D\left(S_{ph}^0-(S_{ph}+k_{ph}X)\right)=-D{\Sigma }_1\left(t\right),$$

which means that ${\Sigma }_1\left(t\right)=e^{-Dt}{\Sigma }_1\left(0\right)$, thus ${lim}_{t\to \infty }{\Sigma }_1\left(t\right)=0$, or equivalently

$$\underset{t\to \infty }{lim}\left(S_{ph}\left(t\right)+k_{ph}X\left(t\right)\right)=S_{ph}^0.$$

Since $X\left(t\right)>0$ for all $t>0$ this means $S_{ph}\left(t\right)>0$ for all $t>0$ as well. Moreover, $X\left(t\right)$ and $S_{ph}\left(t\right)$ are uniformly bounded.

Similarly, using Equations (1) and (3) and denoting ${\Sigma }_2\left(t\right)=S_{cr}\left(t\right)+k_{cr}X\left(t\right)-S_{cr}^0$ we obtain ${\Sigma }_2\left(t\right)=e^{-Dt}{\Sigma }_2\left(0\right)$, which means that

$$\underset{t\to \infty }{lim}\left(S_{cr}\left(t\right)+k_{cr}X\left(t\right)\right)=S_{cr}^0.$$

(16)

Therefore, $S_{cr}\left(t\right)>0$ for all $t>0$ and $S_{cr}\left(t\right)$ is uniformly bounded for $t\ge 0$. Hence, the model solutions $X\left(t\right)$, $S_{ph}\left(t\right)$, $S_{cr}\left(t\right)$ exist for all time $t\ge 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 $\left(D_{cr}^{\left(1\right)},D_{cr}^{\left(2\right)}\right)$, with $D_{cr}^{\left(2\right)}=\mu (S_{ph}^0,S_{cr}^0)$.

Similarly to the proof of Theorem 1, denote ${\Sigma }_3\left(t\right)=S_{ph}\left(t\right)-K{S}_{cr}\left(t\right)-S^0$, where K and $S^0$ are defined in (9). After multiplying Equation (3) by $-{\displaystyle \frac{k_{ph}}{k_{cr}}}$ and adding to Equation (2) we obtain

$$\begin{array}{ccc}\hfill \frac{d}{dt}{\Sigma }_3\left(t\right)& =& \frac{d}{dt}\left(S_{ph}\left(t\right)-K{S}_{cr}\left(t\right)\right)=D\left(S_{ph}^0-S_{ph}\left(t\right)-K{S}_{cr}^0+K{S}_{cr}\left(t\right)\right)\hfill \\ & =& D\left((S_{ph}^0-K{S}_{cr}^0)-(S_{ph}\left(t\right)-K{S}_{cr}\left(t\right))\right)\hfill \\ & =& -D\left(S^0-(S_{ph}\left(t\right)-K{S}_{cr}\left(t\right))\right)=-D{\Sigma }_3\left(t\right).\hfill \end{array}$$

This means that ${\Sigma }_3\left(t\right)=e^{-Dt}{\Sigma }_3\left(0\right)$, ${\Sigma }_3\left(0\right)\ge 0$, so ${lim}_{t\to \infty }{\Sigma }_3\left(t\right)=0$. Then system (1)–(3) may be written in the form

$$\begin{array}{ccc}\hfill \frac{d}{dt}{\Sigma }_3\left(t\right)& =& -D{\Sigma }_3\left(t\right)\hfill \\ \hfill \frac{d}{dt}X\left(t\right)& =& \left(\mu (S^0+K{S}_{cr}\left(t\right),S_{cr}\left(t\right))-D\right)X\left(t\right)\hfill \\ \hfill \frac{d}{dt}{S}_{cr}\left(t\right)& =& -k_{cr}\mu (S^0+K{S}_{cr}\left(t\right),S_{cr}\left(t\right))X\left(t\right)+D(S_{cr}^0-S_{cr}\left(t\right)).\hfill \end{array}$$

Since ${lim}_{t\to \infty }{\Sigma }_3\left(t\right)=0$, the positive $\omega$-limit set of any solution of system (1)–(3) is contained in the set

$${\Omega }_3=\left\{(X,S_{ph},S_{cr}):X>0,S_{ph}\ge 0,S_{cr}\ge 0,{\Sigma }_3=0\right\}.$$

Using the theory of the asymptotically autonomous systems (cf. [20,21]) it follows that all trajectories forming the $\omega$-limit set of any solution of (1)–(3) with initial conditions in ${\Omega }_3$ are solutions of the following limiting system

$$\begin{array}{ccc}{\displaystyle \frac{dX\left(t\right)}{dt}}\hfill & =& \left(\mu (S^0+K{S}_{cr}\left(t\right),S_{cr}\left(t\right))-D\right)X\left(t\right)\hfill \\ {\displaystyle \frac{d{S}_{cr}\left(t\right)}{dt}}\hfill & =& -k_{cr}\mu (S^0+K{S}_{cr}\left(t\right),S_{cr}\left(t\right))X\left(t\right)+D(S_{cr}^0-S_{cr}\left(t\right)).\hfill \end{array}$$

(17)

We consider Equation (17) on the set

$${\Omega }_2=\left\{(X,S_{cr}):X>0,S_{cr}\ge 0,S^0+K{S}_{cr}\ge 0\right\}.$$

Denote for simplicity ${\mu }_{cr}\left(S_{cr}\right)=\mu (S^0+K{S}_{cr},S_{cr})$. Then obviously ${\mu }_{cr}\left(S_{cr}^0\right)=\mu (S^0+K{S}_{cr}^0,S_{cr}^0)=\mu (S_{ph}^0,S_{cr}^0)=D_{cr}^{\left(2\right)}$ holds true.

Let us choose an arbitrary value $\overline{D}\in \left(D_{cr}^{\left(1\right)},D_{cr}^{\left(2\right)}\right)$, and consider the following system obtained from (17) after substituting $D=\overline{D}$ in the latter:

$${\displaystyle \frac{dX\left(t\right)}{dt}}=\left({\mu }_{cr}\left(S_{cr}\left(t\right)\right)-\overline{D}\right)X\left(t\right)$$

(18)

$${\displaystyle \frac{d{S}_{cr}\left(t\right)}{dt}}=-k_{cr}{\mu }_{cr}\left(S_{cr}\left(t\right)\right)X\left(t\right)+\overline{D}(S_{cr}^0-S_{cr}\left(t\right)).$$

(19)

Let us recall, that at $D=\overline{D}$ there are two interior equilibria of the model (1)–(3),

$$\begin{array}{c}{E}_2\left(\overline{D}\right)=(X^{\left(2\right)}\left(\overline{D}\right),S_{ph}^{\left(2\right)}\left(\overline{D}\right),S_{cr}^{\left(2\right)}\left(\overline{D}\right))and\hfill \\ E_3\left(\overline{D}\right)=(X^{\left(3\right)}\left(\overline{D}\right),S_{ph}^{\left(3\right)}\left(\overline{D}\right),S_{cr}^{\left(3\right)}\left(\overline{D}\right))with{S}_{cr}^{\left(2\right)}\left(\overline{D}\right)<S_{cr}^{\left(3\right)}\left(\overline{D}\right).\hfill \end{array}$$

Denote

$$\overline{E}=\left(\overline{X},{\overline{S}}_{cr}\right)=\left(X^{\left(3\right)}\left(\overline{D}\right),S_{cr}^{\left(3\right)}\left(\overline{D}\right)\right).$$

Obviously, $\overline{E}$ is an equilibrium point of (18) and (19).

We make the following assumption.

Assumption 1. 

There exist points $S_{cr}^{-}$ and $S_{cr}^{+}$ such that $0<S_{cr}^{-}<S_{cr}^{+}<S_{cr}^0$ and ${\mu }_{cr}\left(S_{cr}\right)$ is monotone increasing for all $S_{cr}\in \left(S_{cr}^{-},S_{cr}^{+}\right)$.

Assumption 1 identifies the equilibrium $\overline{E}$ with the projection of $E_3\left(\overline{D}\right)$ in the plane $S_{ph}-K{S}_{cr}=S^0$. If we choose for $S_{cr}^{-}$ the $S_{cr}$-component of the double root of Equation (11), i.e., $S_{cr}^{-}=S_{cr}^{\left(2\right)}\left(D_{cr}^{\left(1\right)}\right)=S_{cr}^{\left(3\right)}\left(D_{cr}^{\left(1\right)}\right)$, and $S_{cr}^{+}=S_{cr}^0$, then ${\mu }_{cr}\left(S_{cr}\right)$ is monotone increasing in $\left(S_{cr}^{-},S_{cr}^{+}\right)$, 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 $\overline{D}\in \left(D_{cr}^{\left(1\right)},D_{cr}^{\left(2\right)}\right)$. Then for any initial point $(X\left(0\right),S_{cr}\left(0\right))\in {\Omega }_2$ the corresponding solution of (18) and (19) converges asymptotically towards the equilibrium point $\overline{E}$.

Proof. 

Let us fix an arbitrary initial point $(X\left(0\right),S_{cr}\left(0\right))\in {\Omega }_2$.

First we shall show that there exists time $T>0$, such that $S_{cr}\left(t\right)<S_{cr}^0$ for all $t>T$. Assume that $S_{cr}\left(t\right)\ge S_{cr}^0$ holds true for each $t>0$. Then we have from (19) that

$${\displaystyle \frac{d{S}_{cr}\left(t\right)}{dt}}=-k_{cr}{\mu }_{cr}\left(S_{cr}\right)X\left(t\right)+\overline{D}(S_{cr}^0-S_{cr}\left(t\right))<0.$$

Barbălat’s Lemma [26] implies

$$0=\underset{t\to \infty }{lim}{\displaystyle \frac{d{S}_{cr}\left(t\right)}{dt}}=\underset{t\to \infty }{lim}\left(-k_{cr}{\mu }_{cr}\left(S_{cr}\right)X\left(t\right)+\overline{D}(S_{cr}^0-S_{cr}\left(t\right))\right),$$

which leads to $S_{cr}\left(t\right)\to S_{cr}^0$ and $X\left(t\right)\to 0$ as $t\to \infty$. Further we have that ${\mu }_{cr}\left({\overline{S}}_{cr}\right)=\overline{D}<D_{cr}^{\left(2\right)}={\mu }_{cr}\left(S_{cr}^0\right)$. The continuity of ${\mu }_{cr}(·)$ and the relation $S_{cr}\left(t\right)\to S_{cr}^0$ as $t\to \infty$ imply that there exists a number $\delta >0$ such that

$${\mu }_{cr}\left(S_{cr}\left(t\right)\right)-\overline{D}={\mu }_{cr}\left(S_{cr}\left(t\right)\right)-{\mu }_{cr}\left({\overline{S}}_{cr}\right)\ge \delta$$

for all sufficiently large t. Then it follows ${\displaystyle \frac{dX\left(t\right)}{dt}}=({\mu }_{cr}\left(S_{cr}\left(t\right)\right)-\overline{D})X\left(t\right)\ge \delta X\left(t\right)$ for all sufficiently large t, which contradicts the boundedness of $X\left(t\right)$. Hence, there exists a sufficiently large $T>0$ with $S_{cr}\left(T\right)\le S_{cr}^0$. If the equality $S_{cr}\left(T\right)=S_{cr}^0$ holds true, then we have

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d{S}_{cr}}{dt}}\left(T\right)& =& -k_{cr}{\mu }_{cr}\left(S_{cr}\left(T\right)\right)X\left(T\right)+\overline{D}(S_{cr}^0-S_{cr}\left(T\right))\hfill \\ & =& -k_{cr}{\mu }_{cr}\left(S_{cr}\left(T\right)\right)X\left(T\right)<0.\hfill \end{array}$$

The last inequality shows that $S_{cr}\left(t\right)<S_{cr}^0$ for each $t>T$.

Let us fix an arbitrary $\gamma \in \left(0,({\mu }_{cr}\left(S_{cr}^0\right)-{\mu }_{cr}\left({\overline{S}}_{cr}\right))/2\right)$. (Note that ${\mu }_{cr}\left(S_{cr}\right)$ is monotone increasing.) The continuity of ${\mu }_{cr}$ implies that there exists $\epsilon >0$ such that ${\mu }_{cr}\left({\overline{S}}_{cr}\right)+\gamma <{\mu }_{cr}\left(S_{cr}\right)$ for each $S_{cr}\in \left[S_{cr}^0-(1+k_{cr})\epsilon ,S_{cr}^0\right)$. It follows from (16) that there exists time $T_{\epsilon }>0$ so that $X\left(t\right)$ and $S_{cr}\left(t\right)$ satisfy

$$S_{cr}^0-\epsilon <S_{cr}\left(t\right)+k_{cr}X\left(t\right)<S_{cr}^0+\epsilon foreacht\ge T_{\epsilon }.$$

(20)

Assume now that $X\left(\overline{t}\right)\le \epsilon$ for some $\overline{t}\ge T_{\epsilon }$; then we obtain from (20)

$$S_{cr}^0>S_{cr}\left(\overline{t}\right)\ge S_{cr}^0-k_{cr}X\left(\overline{t}\right)-\epsilon \ge S_{cr}^0-(1+k_{cr})\epsilon ,$$

i.e., $S_{cr}\left(\overline{t}\right)\in \left[S_{cr}^0-(1+k_{cr})\epsilon ,S_{cr}^0\right)$. Hence,

$${\displaystyle \frac{d}{dt}}X\left(\overline{t}\right)=\left({\mu }_{cr}\left(S_{cr}\left(\overline{t}\right)\right)-\overline{D}\right)X\left(\overline{t}\right)=\left({\mu }_{cr}\left(S_{cr}\left(\overline{t}\right)\right)-{\mu }_{cr}\left({\overline{S}}_{cr}\right)\right)X\left(\overline{t}\right)\ge \gamma X\left(\overline{t}\right)>0.$$

It follows then that $X\left(t\right)\ge e^{(t-\overline{t})\gamma }X\left(\overline{t}\right)$. If there exists $t_1\ge \overline{t}$ such that $X\left(t_1\right)=\epsilon$, then at every time $t_2\ge t_1$ with $X\left(t_2\right)=\epsilon$ we have ${\displaystyle \frac{d}{dt}}X\left(t_2\right)=({\mu }_{cr}\left(S_{cr}\left(t_2\right)\right)-\overline{D})X\left(t_2\right)\ge \gamma \epsilon >0$. Hence there exists time $T_1>T$ such that $X\left(t\right)\ge \epsilon$ for each $t\ge T_1$.

The above considerations mean that the $\omega$-limit set of the corresponding trajectory of (18) and (19) lies in the set

$$\{(X,S_{cr}):X\ge \epsilon ,0\le S_{cr}\le S_{cr}^0\}.$$

For $X>0$ and $S_{cr}\in (0,S_{cr}^0)$ we define the following Lyapunov function

$$V=V(X,S_{cr})={\int }_{\overline{X}}^X{\displaystyle \frac{\eta -\overline{X}}{\eta }}d\eta +{\int }_{{\overline{S}}_{cr}}^{S_{cr}}{\displaystyle \frac{\overline{X}({\mu }_{cr}\left(\xi \right)-\overline{D})}{\overline{D}(S_{cr}^0-\xi )}}d\xi .$$

The derivative ${\displaystyle \frac{d}{dt}}V$ of V along the solutions of (18) and (19) is presented by

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d}{dt}}V& =& {\displaystyle \frac{X-\overline{X}}{X}}\left({\mu }_{cr}\left(S_{cr}\right)-\overline{D}\right)X+{\displaystyle \frac{\overline{X}({\mu }_{cr}\left(S_{cr}\right)-\overline{D})}{\overline{D}(S_{cr}^0-S_{cr})}}\left(-k_{cr}{\mu }_{cr}\left(S_{cr}\right)X+\overline{D}(S_{cr}^0-S_{cr})\right)\hfill \\ & =& X({\mu }_{cr}\left(S_{cr}\right)-\overline{D})\left(1-{\displaystyle \frac{\overline{X}{k}_{cr}{\mu }_{cr}\left(S_{cr}\right)}{\overline{D}(S_{cr}^0-S_{cr})}}\right)\hfill \\ & =& X({\mu }_{cr}\left(S_{cr}\right)-\overline{D})\left(1-{\displaystyle \frac{S_{cr}^0-{\overline{S}}_{cr}}{S_{cr}^0-S_{cr}}}·{\displaystyle \frac{{\mu }_{cr}\left(S_{cr}\right)}{\overline{D}}}\right)\le 0\hfill \end{array}$$

for each $S_{cr}\in (0,S_{cr}^0)$ and $X>0$. Applying LaSalle’s invariance principle it follows that each trajectory of (18) and (19) approaches the equilibrium point $\overline{E}$, i.e., $\overline{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_{cr}^{\left(2\right)}=\mu (S_{ph}^0,S_{cr}^0)$ then the model (1)–(3) possesses two equilibrium points—the wash-out equilibrium $E_0=(0,S_{ph}^0,S_{cr}^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 ${\overline{E}}_0=(0,S_{cr}^0)$ of the wash-out equilibrium $E_0$ is globally asymptotically stable if $D>D_{cr}^{\left(2\right)}={\mu }_{cr}\left(S_{cr}^0\right)$. 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_{cr}^{\left(2\right)}$ holds true. Then for any initial point $(X\left(0\right),S_{cr}\left(0\right))>0$ the corresponding solution of (17) converges asymptotically towards the equilibrium ${\overline{E}}_0=(0,S_{cr}^0)$.

Proof. 

Choose some ${\overline{D}}_0>D_{cr}^{\left(2\right)}$ and consider system (18) and (19), where $\overline{D}$ is replaced by ${\overline{D}}_0$. Assume that ${lim}_{t\to \infty }X\left(t\right)=X^{*}>0$. Then Barbălat’s Lemma [26] applied to Equation (18) implies $0={lim}_{t\to \infty }\frac{d}{dt}X\left(t\right)={lim}_{t\to \infty }({\mu }_{cr}\left(S_{cr}\left(t\right)\right)-{\overline{D}}_0)X^{*}$, which means that ${lim}_{t\to \infty }{\mu }_{cr}\left(S_{cr}\left(t\right)\right)={\overline{D}}_0>{\mu }_{cr}\left(S_{cr}^0\right)$. From the continuity of ${\mu }_{cr}(·)$ it follows that there exists time $T>0$ and a positive number $\delta$ such that ${\mu }_{cr}\left(S_{cr}\left(t\right)\right)-{\mu }_{cr}\left(S_{cr}^0\right)\ge \delta$ for all $t\ge T$. The latter inequality leads to $\frac{d}{dt}X\left(t\right)\ge \delta X\left(t\right)$, a contradiction with the boundedness of $X\left(t\right)$. Therefore, $X\left(t\right)\to 0$ as $t\to \infty$. 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 $\frac{d}{dt}{S}_{cr}\left(t\right)={\overline{D}}_0(S_{cr}^0-S_{cr}\left(t\right))$, which implies ${lim}_{t\to \infty }{S}_{cr}\left(t\right)=S_{cr}^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\in (D_{cr}^{\left(1\right)},D_{cr}^{\left(2\right)})$

In this case there exist two positive (coexistence) equilibrium points

$$E_2=(0.0491,0.126,0.0153)and{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_{cr})$, $(X,S_{ph})$ and $(S_{ph},S_{cr})$ 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_{cr}^{\left(2\right)}\approx 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_{cr})$, $(X,S_{ph})$ and $(S_{ph},S_{cr})$ 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 $\mu (S_{ph},S_{cr})$ 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=\left(0,S_{ph}^0,S_{cr}^0\right)$, 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\in \left(D_{cr}^{\left(1\right)},D_{cr}^{\left(2\right)}\right)$, and $E_2$ exists if $D\in \left(D_{cr}^{\left(1\right)},D_{cr}\right)$. 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_{cr}^{\left(2\right)}$. Two types of bifurcations of the equilibria occur, a saddle node bifurcation at $D=D_{cr}^{\left(1\right)}$ where $E_2$ and $E_3$ coalesce, and a transcritical bifurcation at $D=D_{cr}^{\left(2\right)}$, where $E_0$ coincides with $E_3$ and $E_3$ disappears for $D>D_{cr}^{\left(2\right)}$. Practically, the bifurcation values $D_{cr}^{\left(1\right)}$ and $D_{cr}^{\left(2\right)}$ 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 $\mu (·)$ 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\in \left(D_{cr}^{\left(1\right)},D_{cr}^{\left(2\right)}\right)$, $D_{cr}^{\left(2\right)}=\mu (S_{ph}^0,S_{cr}^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 $\left(D_{cr}^{\left(1\right)},D_{cr}^{\left(2\right)}\right)$ 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_{cr}^{\left(2\right)}=\mu (S_{ph}^0,S_{cr}^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_{ph}^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.