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

**Neli Dimitrova 1,\*,† and Plamena Zlateva 1,2,†**


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

**Keywords:** mathematical model; continuous bioreactor; biodegradation; phenol and *p*-cresol mixture; SKIP model; equilibrium points; stability analysis; global stabilizability; numerical simulation

#### **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–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

**Citation:** 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

Received: 18 November 2020 Accepted: 5 January 2021 Published: 8 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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 monoaromatic 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-methylphenol (*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(t)}{dt} \quad = \begin{pmatrix} \mu(S\_{ph}, S\_{cr}) - D \end{pmatrix} X(t) \tag{1}$$

$$\frac{dS\_{ph}(t)}{dt}\_{\dots} = -k\_{ph}\mu(S\_{ph}, S\_{cr})X(t) + D(S\_{ph}^0 - S\_{ph}(t))\tag{2}$$

$$\frac{dS\_{\mathcal{cr}}(t)}{dt} = -k\_{\mathcal{cr}} \,\mu(\mathcal{S}\_{\text{ph}}, \mathcal{S}\_{\text{cr}})X(t) + D(\mathcal{S}\_{\text{cr}}^0 - \mathcal{S}\_{\text{cr}}(t)),\tag{3}$$

where *μ*(*Sph*, *Scr*) is the specific growth rate, presented by

$$\mu(\mathcal{S}\_{\rm ph}, \mathcal{S}\_{\rm cr}) = \frac{\mu\_{\rm max(ph)} \mathcal{S}\_{\rm pl}}{k\_{\rm s(ph)} + \mathcal{S}\_{\rm pl} + \frac{\mathcal{S}\_{\rm pl}^2}{k\_{\rm i(ph)}} + I\_{\rm cr/ph} \mathcal{S}\_{\rm cr}} + \frac{\mu\_{\rm max(cr)} \mathcal{S}\_{\rm cr}}{k\_{\rm s(cr)} + \mathcal{S}\_{\rm cr} + \frac{\mathcal{S}\_{\rm cr}^2}{k\_{\rm i(cr)}} + I\_{\rm ph}/\sigma\_{\rm cr} \mathcal{S}\_{\rm ph}}}. \tag{4}$$

The definition of the state variables *X*, *Sph* and *Scr* 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 *μ*(*Sph*, *Scr*) represents a SKIP (sum kinetics with interaction parameters) model for biological degradation of the chemical compounds. The interaction parameters *Icr*/*ph* and *Iph*/*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 *Icr*/*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 *<sup>μ</sup>*(*Sph*, *Scr*) also involves inhibition terms *<sup>S</sup>*<sup>2</sup> *ph ki*(*ph*) and *<sup>S</sup>*<sup>2</sup> *cr ki*(*cr*) for cell growth on phenol and *p*-cresol, respectively. Obviously, *μ*(*Sph*, 0) and *μ*(0, *Scr*) are the well-known Andrews (or Haldane) model functions, which are unimodal and achieve their maximum at *Sph* <sup>=</sup> ( *ks*(*ph*)*ki*(*ph*) and *Scr* = ( *ks*(*cr*)*ki*(*cr*) respectively.

The influent concentrations *S*<sup>0</sup> *ph*, *<sup>S</sup>*<sup>0</sup> *cr* 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*<sup>0</sup> *ph* and *S*0 *cr* 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 *<sup>μ</sup>*(*Sph*, *Scr*) which does not involve the inhibition terms *<sup>S</sup>*<sup>2</sup> *ph ki*(*ph*) and *<sup>S</sup>*<sup>2</sup> *cr ki*(*cr*) 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.


**Table 1.** Model variables and parameters.

#### **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\_{pl\nu}, S\_{cr}) - D\right)X = 0\tag{5}$$

$$-k\_{\rm plr} \mu (S\_{\rm phr} S\_{\rm cr}) X + D (S\_{\rm pll}^0 - S\_{\rm plr}) = 0 \tag{6}$$

$$-k\_{\mathcal{cr}}\,\mu(S\_{\rm ph}, S\_{\rm cr})X + D(S\_{\rm cr}^0 - S\_{\rm cr}) = 0.\tag{7}$$

Obviously, the point *E*<sup>0</sup> = (0, *S*<sup>0</sup> *ph*, *<sup>S</sup>*<sup>0</sup> *cr*) (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 −*kcr*, Equation (7) by *kph* and summing the latter, we obtain

$$-k\_{cr}(S\_{ph}^0 - S\_{ph}) + k\_{ph}(S\_{cr}^0 - S\_{cr}) = 0.\tag{8}$$

Let us express *Sph* from (8) as a function of *Scr*. Denoting

$$K = \frac{k\_{ph}}{k\_{cr}}, \quad S^0 = S^0\_{ph} - K S^0\_{cr} \tag{9}$$

We obtain

$$S\_{ph} = S\_{ph}^0 - \frac{k\_{ph}}{k\_{cr}} \left( S\_{cr}^0 - S\_{cr} \right) = S^0 + KS\_{cr} \,. \tag{10}$$

After replacing the latter presentation of *Sph* into the equation *μ*(*Sph*, *Scr*) = *D* from (5) we obtain an equation with respect to *Scr* of the form

$$
\mu(S^0 + KS\_{cr\_\prime}S\_{cr}) = D\_{\prime\prime}
$$

or equivalently

$$\begin{split} \frac{\mu\_{\text{max}(ph)} \left(S^0 + KS\_{\text{cr}}\right)}{k\_{\text{s}(ph)} + S^0 + KS\_{\text{cr}} + \frac{1}{k\_{\text{i}(ph)}} (S^0 + KS\_{\text{cr}})^2 + I\_{\text{cr}/\text{ph}}S\_{\text{cr}}} \\ + \frac{\mu\_{\text{max}(\text{cr})}S\_{\text{cr}}}{k\_{\text{s}(\text{cr})} + S\_{\text{cr}} + \frac{1}{k\_{\text{i}(\text{cr})}}S\_{\text{cr}}^2 + I\_{\text{ph}/\text{cr}}(S^0 + KS\_{\text{cr}})} = D. \end{split}$$

Straightforward calculations lead to a polynomial equation of the form

$$A\_1 \mathcal{S}\_{cr}^4 + A\_2 \mathcal{S}\_{cr}^3 + A\_3 \mathcal{S}\_{cr}^2 + A\_4 \mathcal{S}\_{cr} + A\_5 = 0,\tag{11}$$

where

*<sup>A</sup>*<sup>1</sup> <sup>=</sup> <sup>−</sup>*<sup>D</sup>* · <sup>1</sup> *ki*(*cr*) · <sup>1</sup> *ki*(*ph*) ; *A*<sup>2</sup> = *μmax*(*ph*) *K ki*(*cr*) + *μmax*(*cr*) 1 *ki*(*ph*) − *D* (<sup>1</sup> <sup>+</sup> *Iph*/*crK*) <sup>1</sup> *ki*(*ph*) + 1 *ki*(*cr*) ) *Icr*/*ph* <sup>+</sup> *<sup>K</sup>* <sup>+</sup> <sup>2</sup>*KS*<sup>0</sup> <sup>1</sup> *ki*(*ph*) \*; *A*<sup>3</sup> = *μmax*(*ph*) *<sup>S</sup>*<sup>0</sup> <sup>1</sup> *ki*(*cr*) + *K*(1 + *Iph*/*crK*) + *μmax*(*cr*) *<sup>K</sup>* <sup>+</sup> *Icr*/*ph* <sup>+</sup> <sup>2</sup>*KS*<sup>0</sup> <sup>1</sup> *ki*(*ph*) − *D* 1 *ki*(*ph*) (*ks*(*cr*) + *Iph*/*crS*<sup>0</sup>)+(1 + *Iph*/*crK*) ) *<sup>K</sup>* <sup>+</sup> *Icr*/*ph* <sup>+</sup> <sup>2</sup>*KS*<sup>0</sup> <sup>1</sup> *ki*(*ph*) \* + 1 *ki*(*cr*) ) *ks*(*ph*) + *S*<sup>0</sup> + 1 *ki*(*ph*) *S*02 \*; *A*<sup>4</sup> = *μmax*(*ph*) *S*0(1 + *Iph*/*crK*) + *K*(*ks*(*cr*) + *Iph*/*crS*<sup>0</sup>) + *μmax*(*cr*) ) *ks*(*ph*) + *S*<sup>0</sup> + 1 *ki*(*ph*) *S*02 \* − *D* (*ks*(*cr*) + *Iph*/*crS*<sup>0</sup>) ) *<sup>K</sup>* <sup>+</sup> *Icr*/*ph* <sup>+</sup> <sup>2</sup>*KS*<sup>0</sup> <sup>1</sup> *ki*(*ph*) \* + (1 + *Iph*/*crK*) ) *ks*(*ph*) + *S*<sup>0</sup> + 1 *ki*(*ph*) *S*02 \*; *A*<sup>5</sup> = *<sup>μ</sup>max*(*ph*)*S*<sup>0</sup> <sup>−</sup> *<sup>D</sup>* ) *ks*(*ph*) + *S*<sup>0</sup> + 1 *ki*(*ph*) *S*02 \*(*ks*(*cr*) <sup>+</sup> *Iph*/*crS*<sup>0</sup>).

All coefficients *Ai*, *i* = 1, 2, . . . , 5, depend on the parameter *D*. Obviously, if *A*<sup>5</sup> = 0, then Equation (11) possesses a solution *Scr* = 0. We have

$$A\_5 = 0 \iff D = D\_{\mathcal{O}} := \frac{S^0 \mu\_{\max(ph)} k\_{i(ph)}}{S^{0^2} + k\_{i(ph)} S^0 + k\_{i(ph)} k\_{s(ph)}} = \mu(S^0, 0). \tag{12}$$

The latter value of *Dcr* is biologically reasonable only if *S*<sup>0</sup> > 0. Using the numerical values of the model coefficients in the last column of Table 1, we obtain

$$S^0 = S^0\_{\text{pl}} - KS^0\_{cr} \approx 0.09483 > 0,$$

and so,

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

This means that for *D* = *Dcr* there exists an equilibrium point with *Scr* = 0. Further from (10) we compute the component of *Sph* = *S*0, and from (7) we get the corresponding component of *<sup>X</sup>* <sup>=</sup> *<sup>S</sup>*<sup>0</sup> *cr kcr* . Thus, at *D* = *Dcr* there exists a steady state

$$E\_1 = E\_1(D\_{cr}) = \left(\frac{S\_{cr}^0}{k\_{cr}}, S^0, 0\right) = (0.05172, 0.09483, \, 0). \tag{13}$$

Considering the cubic equation *A*1*S*<sup>3</sup> *cr* + *A*2*S*<sup>2</sup> *cr* + *A*3*Scr* + *A*<sup>4</sup> = 0 at *D* = *Dcr* (i.e., with *A*<sup>5</sup> = 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* = *Dcr* then Equation (11) may possess up to 4 real positive solutions with respect to *Scr*. If there exists at least one positive solution of (11), say *S*∗ *cr*, such that *S*<sup>∗</sup> *cr* < *S*<sup>0</sup> *cr* for some values of *D*, we shall have an interior (with positive components) equilibrium of the form

$$E^\* = (X^\*, S^\*\_{ph'}, S^\*\_{cr}), \quad S^\*\_{ph} = S^0 + KS^\*\_{cr} < S^0\_{ph'} \quad X^\* = \frac{S^0\_{cr} - S^\*\_{cr}}{k\_{cr}} = \frac{S^0\_{ph} - S^\*\_{ph}}{k\_{ph}}.\tag{14}$$

**Remark 1.** *If we express Scr from (8) as a function of Sph and denote <sup>K</sup>*<sup>ˆ</sup> <sup>=</sup> *kcr kph* <sup>=</sup> <sup>1</sup> *<sup>K</sup>, <sup>S</sup>*ˆ0 <sup>=</sup> *S*0 *cr* <sup>−</sup> *KS*<sup>ˆ</sup> <sup>0</sup> *ph, then we shall have Scr* <sup>=</sup> *<sup>S</sup>*ˆ0 <sup>+</sup> *KS*<sup>ˆ</sup> *ph. Similar calculations as above will produce a polynomial equation of the form A*ˆ <sup>1</sup>*S*<sup>4</sup> *ph* <sup>+</sup> *<sup>A</sup>*<sup>ˆ</sup> <sup>2</sup>*S*<sup>3</sup> *ph* <sup>+</sup> *<sup>A</sup>*<sup>ˆ</sup> <sup>3</sup>*S*<sup>2</sup> *ph* <sup>+</sup> *<sup>A</sup>*<sup>ˆ</sup> <sup>4</sup>*Sph* <sup>+</sup> *<sup>A</sup>*<sup>ˆ</sup> <sup>5</sup> <sup>=</sup> <sup>0</sup>*, where the coefficients A*ˆ*<sup>i</sup> are similar to Ai, i* = 1, 2, ... , 5*, within S*ˆ0 *and K*ˆ *instead of S*<sup>0</sup> *and K, respectively. In this case we have*

$$\mathcal{A}\_5 = \left(\mu\_{\text{max},cr} - D\left(k\_{s(cr)} + \mathcal{S}^0 + \frac{\mathcal{S}^0}{k\_{\text{i}(cr)}}\right)\right) (k\_{s(ph)} + I\_{cr/ph}\mathcal{S}^0).$$

*Obviously, <sup>A</sup>*<sup>ˆ</sup> <sup>5</sup> <sup>=</sup> <sup>0</sup> *at <sup>D</sup>*<sup>ˆ</sup> <sup>=</sup> *<sup>μ</sup>*(0, *<sup>S</sup>*ˆ0)*. But in this case <sup>S</sup>*ˆ0 <sup>=</sup> <sup>−</sup> <sup>1</sup> *<sup>K</sup> <sup>S</sup>*<sup>0</sup> ≈ −0.047 <sup>&</sup>lt; <sup>0</sup>*, thus there is no value of D at which Sph* = 0 *is a root of the polynomial* ∑<sup>5</sup> *<sup>i</sup>*=<sup>1</sup> *<sup>A</sup>*ˆ*iS*<sup>5</sup>−*<sup>i</sup> ph* = 0*. As we shall see in the following, this is the case with the equilibrium component Sph.*

Numerical computations show that if *D* > *Dcr* then there are no positive real roots of Equation (11). Therefore, we can expect interior (coexistence) equilibria of the form *E*<sup>∗</sup> if *D* ∈ (0, *Dcr*), in case that the equilibrium components with respect to *Scr* satisfy the inequality *Scr* <sup>≤</sup> *<sup>S</sup>*<sup>0</sup> *cr*. Further we obtain numerically the following results:


The left plot in Figure <sup>1</sup> shows the graph of the function *<sup>μ</sup>*(*S*<sup>0</sup> <sup>+</sup> *KScr*, *Scr*) for *Scr* <sup>∈</sup> [0, *S*<sup>0</sup> *cr*]=[0, 0.3]; the horizontal dash lines correspond to the values of *<sup>D</sup>*(1) *cr* , *<sup>D</sup>*(2) *cr* and *Dcr*. Therefore, the model (1)–(3) possesses two interior equilibrium points depending on the values of *D*. Denote them by

$$\begin{aligned} E\_2 &= E\_2(D) &= \begin{pmatrix} X^{(2)}, S^{(2)}\_{ph}, S^{(2)}\_{cr} \end{pmatrix}, & D &\in \left( D^{(1)}\_{cr}, D\_{cr} \right); \\ E\_3 &= E\_3(D) &= \begin{pmatrix} X^{(3)}, S^{(3)}\_{ph}, S^{(3)}\_{cr} \end{pmatrix}, & D &\in \left( D^{(1)}\_{cr}, D^{(2)}\_{cr} \right), & \text{with } S^{(3)}\_{cr} > S^{(2)}\_{cr}. \end{aligned}$$

Numerical computations also produce the following results:

$$\begin{aligned} E\_2(D\_{\sigma7}) &= E\_1 = (0.05172, 0.09483, 0), & D\_{\sigma7} &= \mu(S^0, 0) = 0.09933; \\ E\_2(D\_{\sigma7}^{(1)}) &= E\_3(D\_{\sigma7}^{(1)}) = (0.04426, 0.18211, 0.04327), & D\_{\sigma7}^{(1)} &= 0.0745599; \\ E\_3(D\_{\sigma7}^{(2)}) &= E\_0 = (0, S\_{ph}^0, S\_{\sigma7}^0) = (0, 0.7, 0.3), & D\_{\sigma7}^{(2)} &= \mu(S\_{ph}^0, S\_{\sigma7}^0) = 0.08651. \end{aligned}$$

**Figure 1.** (**Left**): graph of the function *<sup>μ</sup>*(*S*<sup>0</sup> <sup>+</sup> *KScr*, *Scr*) for *Scr* <sup>∈</sup> [0, *<sup>S</sup>*<sup>0</sup> *cr*]. (**Right**): the equilibrium components *<sup>S</sup>*(2) *cr* (dash line) and *<sup>S</sup>*(3) *cr* (solid line), parameterized on *D*. The horizontal dash-dot&solid line passes trough *S*<sup>0</sup> *cr*. On the horizontal axis, the solid circle denotes *<sup>D</sup>*(1) *cr* , the solid box denotes *<sup>D</sup>*(2) *cr* , the diamond denotes *Dcr*. The vertical dot line passes through *<sup>D</sup>*(2) *cr* .

Figure 1 (right plot) and Figure 2 visualize the components *Scr*, *Sph* and *X* of the equilibria *E*0, *E*<sup>2</sup> and *E*3. In the three plots, the components of the equilibrium point *E*<sup>0</sup> are marked by horizontal dash-dot&solid lines, the components of *E*<sup>2</sup> are marked by dash lines and the ones of *E*<sup>3</sup> are shown by solid lines.

**Figure 2.** (**Left**): the equilibrium components *<sup>S</sup>*(2) *ph* (dash line) and *<sup>S</sup>*(3) *ph* (solid line), parameterized on *D*. The horizontal dashdot&solid line passes trough *S*<sup>0</sup> *ph*. (**Right**): the equilibrium components *<sup>X</sup>*(2) (dash line) and *<sup>X</sup>*(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 *<sup>D</sup>*(1) *cr* , the solid box denotes *<sup>D</sup>*(2) *cr* , the diamond denotes *Dcr*. The vertical dot line passes through *<sup>D</sup>*(2) *cr* .

#### **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 *μ*(*Sph*, *Scr*). The Jacobi matrix *J* related to the model Equations (1)–(3) has the form

$$J = \begin{pmatrix} \mu(S\_{ph}, S\_{cr}) - D & \frac{\partial \mu}{\partial S\_{ph}} X & \frac{\partial \mu}{\partial S\_{cr}} X \\\\ -k\_{phl} \mu(S\_{ph}, S\_{cr}) & -k\_{phl} \frac{\partial \mu}{\partial S\_{ph}} X - D & -k\_{phl} \frac{\partial \mu}{\partial S\_{cr}} X \\\\ -k\_{crl} \mu(S\_{ph}, S\_{cr}) & -k\_{crl} \frac{\partial \mu}{\partial S\_{ph}} X & -k\_{crl} \frac{\partial \mu}{\partial S\_{cr}} X - D \end{pmatrix}.$$

The characteristic polynomial corresponding to *J* is defined by *det*(*J* − *λI*3), where *λ* is any complex number and *I*<sup>3</sup> is the (3 × 3)–identity matrix

$$\det(I - \lambda I\_3) = \begin{vmatrix} \mu(\mathcal{S}\_{ph\prime}, \mathcal{S}\_{\mathcal{ST}}) - D - \lambda & \frac{\partial \mu}{\partial \mathcal{S}\_{ph}} X & \frac{\partial \mu}{\partial \mathcal{S}\_{\mathcal{ST}}} X \\\\ -k\_{phl} \mu(\mathcal{S}\_{ph\prime}, \mathcal{S}\_{\mathcal{ST}}) & -k\_{phl} \frac{\partial \mu}{\partial \mathcal{S}\_{ph}} X - D - \lambda & -k\_{phl} \frac{\partial \mu}{\partial \mathcal{S}\_{\mathcal{ST}}} X \\\\ -k\_{\mathcal{ST}} \mu(\mathcal{S}\_{ph\prime}, \mathcal{S}\_{\mathcal{ST}}) & -k\_{\mathcal{ST}} \frac{\partial \mu}{\partial \mathcal{S}\_{ph}} X & -k\_{\mathcal{ST}} \frac{\partial \mu}{\partial \mathcal{S}\_{\mathcal{T}}} X - D - \lambda \end{vmatrix}.$$

Multiplying the second row of the above determinant by <sup>−</sup> *kcr kph* and adding the latter to the third row, we obtain

$$\det(f - \lambda I\_3) = \begin{vmatrix} \mu(S\_{ph}, S\_{cr}) - D - \lambda & \frac{\partial \mu}{\partial S\_{ph}} X & \frac{\partial \mu}{\partial S\_{cr}} X \\\\ -k\_{ph} \mu(S\_{ph}, S\_{cr}) & -k\_{ph} \frac{\partial \mu}{\partial S\_{ph}} X - D - \lambda & -k\_{ph} \frac{\partial \mu}{\partial S\_{cr}} X \\\\ 0 & \frac{k\_{cr}}{k\_{ph}} (D + \lambda) & -D - \lambda \end{vmatrix}.$$

Straightforward calculations deliver the following characteristic polynomial

$$\det(I - \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]. \tag{15}$$

Denote by *J*(*Ei*) the Jacobian matrix evaluated at the equilibrium point *Ei*, *i* = 0, 1, 2, 3. It follows from (15) that *λ*1,2 = −*D* < 0 are always eigenvalues of *J*(*Ei*), *i* = 0, 1, 2, 3. The third eigenvalue *λ*<sup>3</sup> is determined from the second multiplier of (15).

#### **Proposition 1.**


**Proof.** (*i*)–(*iii*) We obtain from (15)

$$\det(f(E\_0) - \lambda I\_3) = (D + \lambda)^2 (\mu(S\_{\text{pl}\prime}^0 S\_{cr}^0) - D - \lambda)\_{\text{st}}$$

thus the third root *λ*<sup>3</sup> satisfies

$$\lambda\_3 = \mu(S^0\_{\text{pl}}, S^0\_{\text{cr}}) - D \left\{ \begin{aligned} &> 0, \quad \text{if } D < D^{(2)}\_{\text{cr}} = \mu(S^0\_{\text{pl}}, S^0\_{\text{cr}}), \\ &< 0, \quad \text{if } D > D^{(2)}\_{\text{cr}}, \\ &= 0, \quad \text{if } D = D^{(2)}\_{\text{cr}}. \end{aligned} \right.$$

(*iv*) The characteristic polynomial corresponding to the equilibrium *E*<sup>1</sup> is presented by

$$\det(I(E\_1) - \lambda I\_3) = -(D\_{cr} + \lambda)^2 \left( S\_{cr}^0 \left( K \frac{\partial \mu}{\partial S\_{pl}}(S^0, 0) + \frac{\partial \mu}{\partial S\_{cr}}(S^0, 0) \right) + \lambda \right).$$

The third root *λ*<sup>3</sup> of the latter polynomial is computed numerically and is equal to

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

which means that *E*1(*Dcr*) is a saddle equilibrium point. This proves the proposition.

The equilibrium components *<sup>S</sup>*(*i*) *ph* and *<sup>S</sup>*(*i*) *cr* of the equilibria *Ei*, *i* = 2, 3, satisfy the equation *μ*(*Sph*, *Scr*) = *D*, so that from (15) we obtain

$$\begin{split} \det(I(E\_i) - \lambda I\_3) &= -(D + \lambda)^2 \Big[ \lambda + X^{(i)} \Big( k\_{ph} \frac{\partial \mu}{\partial S\_{ph}} (S\_{ph'}^{(i)} S\_{cr}^{(i)}) + k\_{cr} \frac{\partial \mu}{\partial S\_{cr}} (S\_{ph'}^{(i)} S\_{cr}^{(i)}) \Big) \Big], \\ i = 2, 3. \end{split}$$

The third root *<sup>λ</sup>*(*i*) <sup>3</sup> <sup>=</sup> <sup>−</sup>*X*(*i*) \$ *kph ∂μ <sup>∂</sup>Sph* (*S*(*i*) *ph* , *<sup>S</sup>*(*i*) *cr* ) + *kcr ∂μ <sup>∂</sup>Scr* (*S*(*i*) *ph* , *<sup>S</sup>*(*i*) *cr* ) % is found numerically by computing the right-hand side expression on a discrete mesh of values for *D*, where *<sup>D</sup>* <sup>∈</sup> (*D*(1) *cr* , *Dcr*) for *<sup>E</sup>*2, and *<sup>D</sup>* <sup>∈</sup> (*D*(1) *cr* , *<sup>D</sup>*(2) *cr* ) for *<sup>E</sup>*3. Figure <sup>3</sup> 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 *<sup>J</sup>*(*E*2) approaches zero at *<sup>D</sup>* <sup>=</sup> *<sup>D</sup>*(1) *cr* , and one eigenvalue of *<sup>J</sup>*(*E*3) approaches zero at *<sup>D</sup>* <sup>=</sup> *<sup>D</sup>*(1) *cr* and *<sup>D</sup>* <sup>=</sup> *<sup>D</sup>*(2) *cr* , thus *<sup>D</sup>*(1) *cr* and *<sup>D</sup>*(2) *cr* are bifurcation parameter values.

**Figure 3.** Eigenvalues corresponding to the equilibrium points *E*<sup>2</sup> (**left**) and *E*<sup>3</sup> (**right**), parameterized on *D*. On the horizontal axis, the solid circle denotes *<sup>D</sup>*(1) *cr* , the solid box denotes *<sup>D</sup>*(2) *cr* , the diamond denotes *Dcr*.

We summarize the above results in the next proposition.

#### **Proposition 2.**


Figure 1 (right plot) and Figure 2 also visualize the stability of *E*0, *E*<sup>2</sup> 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 *<sup>D</sup>* <sup>=</sup> *<sup>D</sup>*(1) *cr* , and *<sup>D</sup>* <sup>=</sup> *<sup>D</sup>*(2) *cr* 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*, Sph*(0) ≥ 0*, Scr*(0)) ≥ 0*.*


**Proof.** (*i*) Let *X*(0) = 0 and *Sph*(0) ≥ 0, *Scr*(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

$$\begin{array}{rcl}\frac{dS\_{\mathrm{pl}}(t)}{dt} &=& D(S\_{\mathrm{pl}}^0 - S\_{\mathrm{pl}}(t))\\\frac{dS\_{\mathrm{cr}}(t)}{dt} &=& D(S\_{\mathrm{cr}}^0 - S\_{\mathrm{cr}}(t)).\end{array}$$

The latter equations imply that *Sph*(*t*) and *Scr*(*t*) converge exponentially to *S*<sup>0</sup> *ph* and *S*0 *cr* respectively. The plane *X* = 0 is invariant for the model.

(*ii*)–(*iii*) Assume that *X*(0) > 0, *Sph*(0) ≥ 0, *Scr*(0)) ≥ 0. It follows from Equation (1) that

$$\begin{aligned} \frac{dX}{X} &= \int\_0^t (\mu(S\_{ph}(\tau), S\_{cr}(\tau)) - D)d\tau, \\ X(t) &= X(0)e^{\int\_0^t (\mu(S\_{ph}(\tau), S\_{cr}(\tau)) - D)d\tau} > 0 \text{ for each } t \ge 0. \end{aligned}$$

Denote <sup>Σ</sup>1(*t*) = *Sph*(*t*) + *kphX*(*t*) <sup>−</sup> *<sup>S</sup>*<sup>0</sup> *ph*. Then Equations (1) and (2) imply

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

which means that <sup>Σ</sup>1(*t*) = *<sup>e</sup>*−*Dt*Σ1(0), thus lim*t*→<sup>∞</sup> <sup>Σ</sup>1(*t*) = 0, or equivalently

$$\lim\_{t \to \infty} \left( S\_{\text{ph}}(t) + k\_{\text{ph}} X(t) \right) = S\_{\text{ph}}^0 \dots$$

Since *X*(*t*) > 0 for all *t* > 0 this means *Sph*(*t*) > 0 for all *t* > 0 as well. Moreover, *X*(*t*) and *Sph*(*t*) are uniformly bounded.

Similarly, using Equations (1) and (3) and denoting <sup>Σ</sup>2(*t*) = *Scr*(*t*) + *kcrX*(*t*) <sup>−</sup> *<sup>S</sup>*<sup>0</sup> *cr* we obtain Σ2(*t*) = *e*−*Dt*Σ2(0), which means that

$$\lim\_{t \to \infty} (S\_{cr}(t) + k\_{cr}X(t)) = S\_{cr}^{0}. \tag{16}$$

Therefore, *Scr*(*t*) > 0 for all *t* > 0 and *Scr*(*t*) is uniformly bounded for *t* ≥ 0. Hence, the model solutions *X*(*t*), *Sph*(*t*), *Scr*(*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 *<sup>D</sup>* belongs to the interval \$ *<sup>D</sup>*(1) *cr* , *<sup>D</sup>*(2) *cr* % , with *<sup>D</sup>*(2) *cr* <sup>=</sup> *<sup>μ</sup>*(*S*<sup>0</sup> *ph*, *<sup>S</sup>*<sup>0</sup> *cr*). Similarly to the proof of Theorem 1, denote <sup>Σ</sup>3(*t*) = *Sph*(*t*) <sup>−</sup> *KScr*(*t*) <sup>−</sup> *<sup>S</sup>*0, where

*<sup>K</sup>* and *<sup>S</sup>*<sup>0</sup> are defined in (9). After multiplying Equation (3) by <sup>−</sup>*kph kcr* and adding to Equation (2) we obtain

$$\begin{aligned} \frac{d}{dt}\Sigma\_3(t) &= \quad \frac{d}{dt}\Big(S\_{ph}(t) - KS\_{cr}(t)\Big) = D\Big(S\_{ph}^0 - S\_{ph}(t) - KS\_{cr}^0 + KS\_{cr}(t)\Big) \\ &= \quad D\Big(\left(S\_{ph}^0 - KS\_{cr}^0\right) - \left(S\_{ph}(t) - KS\_{cr}(t)\right)\Big) \\ &= \quad -D\Big(S^0 - \left(S\_{ph}(t) - KS\_{cr}(t)\right)\Big) = -D\Sigma\_3(t). \end{aligned}$$

This means that <sup>Σ</sup>3(*t*) = *<sup>e</sup>*−*Dt*Σ3(0), <sup>Σ</sup>3(0) <sup>≥</sup> 0, so lim*t*→<sup>∞</sup> <sup>Σ</sup>3(*t*) = 0. Then system (1)–(3) may be written in the form

$$\begin{cases} \frac{d}{dt}\Sigma\_3(t) &=& -D\Sigma\_3(t) \\ \frac{d}{dt}X(t) &=& \left(\mu(S^0 + KS\_{cr}(t), S\_{cr}(t)) - D\right)X(t) \\ \frac{d}{dt}S\_{cr}(t) &=& -k\_{cr}\mu(S^0 + KS\_{cr}(t), S\_{cr}(t))X(t) + D(S\_{cr}^0 - S\_{cr}(t)). \end{cases}$$

Since lim*t*→<sup>∞</sup> <sup>Σ</sup>3(*t*) = 0, the positive *<sup>ω</sup>*-limit set of any solution of system (1)–(3) is contained in the set

$$\Omega\_3 = \left\{ (X, \mathcal{S}\_{\text{pl}}, \mathcal{S}\_{\text{cr}}) : X > 0, \mathcal{S}\_{\text{pl}} \ge 0, \mathcal{S}\_{\text{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 *ω*-limit set of any solution of (1)–(3) with initial conditions in Ω<sup>3</sup> are solutions of the following limiting system

$$\begin{array}{rcl}\frac{dX(t)}{dt} &=& \left(\mu(S^0 + KS\_{\rm cr}(t), S\_{\rm cr}(t)) - D\right)X(t) \\\\ \frac{dS\_{\rm cr}(t)}{dt} &=& -k\_{\rm cr}\,\mu(S^0 + KS\_{\rm cr}(t), S\_{\rm cr}(t))X(t) + D(S\_{\rm cr}^0 - S\_{\rm cr}(t)). \end{array} \tag{17}$$

We consider Equation (17) on the set

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

Denote for simplicity *μcr*(*Scr*) = *μ*(*S*<sup>0</sup> + *KScr*, *Scr*). Then obviously *μcr*(*S*<sup>0</sup> *cr*) = *μ*(*S*<sup>0</sup> + *KS*<sup>0</sup> *cr*, *S*<sup>0</sup> *cr*) = *μ*(*S*<sup>0</sup> *ph*, *<sup>S</sup>*<sup>0</sup> *cr*) = *<sup>D</sup>*(2) *cr* holds true.

Let us choose an arbitrary value *<sup>D</sup>*¯ <sup>∈</sup> \$ *<sup>D</sup>*(1) *cr* , *<sup>D</sup>*(2) *cr* % , and consider the following system obtained from (17) after substituting *D* = *D*¯ in the latter:

$$\frac{dX(t)}{dt}\_{\dots,\frac{dt}{dt}} = (\mu\_{cr}(S\_{cr}(t)) - D)X(t) \tag{18}$$

$$\frac{dS\_{\rm cr}(t)}{dt} = -k\_{\rm cr} \mu\_{\rm cr}(S\_{\rm cr}(t))X(t) + \vec{D}(S\_{\rm cr}^0 - S\_{\rm cr}(t)).\tag{19}$$

Let us recall, that at *D* = *D*¯ there are two interior equilibria of the model (1)–(3),

$$\begin{split} E\_2(\bar{D}) &= (X^{(2)}(\bar{D}), S^{(2)}\_{ph}(\bar{D}), S^{(2)}\_{cr}(\bar{D})) \quad \text{and} \\ E\_3(\bar{D}) &= (X^{(3)}(\bar{D}), S^{(3)}\_{ph}(\bar{D}), S^{(3)}\_{cr}(\bar{D})) \quad \text{with} \quad S^{(2)}\_{cr}(\bar{D}) < S^{(3)}\_{cr}(\bar{D}). \end{split}$$

Denote

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

Obviously, *E*¯ is an equilibrium point of (18) and (19). We make the following assumption.

**Assumption 1.** *There exist points S*− *cr and <sup>S</sup>*<sup>+</sup> *cr such that* 0 < *S*<sup>−</sup> *cr* < *<sup>S</sup>*<sup>+</sup> *cr* < *S*<sup>0</sup> *cr and μcr*(*Scr*) *is monotone increasing for all Scr* ∈ (*S*<sup>−</sup> *cr*, *<sup>S</sup>*<sup>+</sup> *cr*)*.*

Assumption 1 identifies the equilibrium *E*¯ with the projection of *E*3(*D*¯ ) in the plane *Sph* <sup>−</sup> *KScr* <sup>=</sup> *<sup>S</sup>*0. If we choose for *<sup>S</sup>*<sup>−</sup> *cr* the *Scr*-component of the double root of Equation (11), i.e., *S*− *cr* <sup>=</sup> *<sup>S</sup>*(2) *cr* (*D*(1) *cr* ) = *<sup>S</sup>*(3) *cr* (*D*(1) *cr* ), and *<sup>S</sup>*<sup>+</sup> *cr* = *S*<sup>0</sup> *cr*, then *μcr*(*Scr*) is monotone increasing in (*S*− *cr*, *<sup>S</sup>*<sup>+</sup> *cr*), 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–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 <sup>1</sup> be fulfilled. Assume that <sup>D</sup>*¯ <sup>∈</sup> \$ *<sup>D</sup>*(1) *cr* , *<sup>D</sup>*(2) *cr* % *. Then for any initial point* (*X*(0), *Scr*(0)) ∈ Ω<sup>2</sup> *the corresponding solution of (18) and (19) converges asymptotically towards the equilibrium point E.* ¯

**Proof.** Let us fix an arbitrary initial point (*X*(0), *Scr*(0)) ∈ Ω2.

First we shall show that there exists time *T* > 0, such that *Scr*(*t*) < *S*<sup>0</sup> *cr* for all *t* > *T*. Assume that *Scr*(*t*) <sup>≥</sup> *<sup>S</sup>*<sup>0</sup> *cr* holds true for each *t* > 0. Then we have from (19) that

$$\frac{d S\_{\mathcal{T}}(t)}{dt} = -k\_{\mathcal{T}} \,\mu\_{\mathcal{T}}(S\_{\mathcal{T}}) X(t) + \vec{D} (S\_{\mathcal{T}}^0 - S\_{\mathcal{T}}(t)) < 0.$$

Barb ˘alat's Lemma [26] implies

$$0 = \lim\_{t \to \infty} \frac{dS\_{cr}(t)}{dt} = \lim\_{t \to \infty} \left( -k\_{cr} \,\mu\_{cr}(S\_{cr}) X(t) + \mathcal{D}(S\_{cr}^0 - S\_{cr}(t)) \right),$$

which leads to *Scr*(*t*) <sup>→</sup> *<sup>S</sup>*<sup>0</sup> *cr* and *<sup>X</sup>*(*t*) <sup>→</sup> 0 as *<sup>t</sup>* <sup>→</sup> <sup>∞</sup>. Further we have that *<sup>μ</sup>cr*(*S*¯ *cr*) = *D*¯ < *<sup>D</sup>*(2) *cr* <sup>=</sup> *<sup>μ</sup>cr*(*S*<sup>0</sup> *cr*). The continuity of *<sup>μ</sup>cr*(·) and the relation *Scr*(*t*) <sup>→</sup> *<sup>S</sup>*<sup>0</sup> *cr* as *t* → ∞ imply that there exists a number *δ* > 0 such that

$$
\mu\_{cr}(\mathcal{S}\_{cr}(t)) - \mathcal{D} = \mu\_{cr}(\mathcal{S}\_{cr}(t)) - \mu\_{cr}(\mathcal{S}\_{cr}) \ge \delta
$$

for all sufficiently large *<sup>t</sup>*. Then it follows *dX*(*t*) *dt* = (*μcr*(*Scr*(*t*)) <sup>−</sup> *<sup>D</sup>*¯ )*X*(*t*) <sup>≥</sup> *<sup>δ</sup>X*(*t*) for all sufficiently large *t*, which contradicts the boundedness of *X*(*t*). Hence, there exists a sufficiently large *<sup>T</sup>* <sup>&</sup>gt; 0 with *Scr*(*T*) <sup>≤</sup> *<sup>S</sup>*<sup>0</sup> *cr*. If the equality *Scr*(*T*) = *S*<sup>0</sup> *cr* holds true, then we have

$$\begin{aligned} \frac{dS\_{cr}}{dt}(T) &= -k\_{cr} \mu\_{cr} (S\_{cr}(T)) X(T) + D(S\_{cr}^0 - S\_{cr}(T)),\\ &= -k\_{cr} \mu\_{cr} (S\_{cr}(T)) X(T) < 0. \end{aligned}$$

The last inequality shows that *Scr*(*t*) < *S*<sup>0</sup> *cr* for each *t* > *T*.

Let us fix an arbitrary *<sup>γ</sup>* <sup>∈</sup> - 0,(*μcr*(*S*<sup>0</sup> *cr*) <sup>−</sup> *<sup>μ</sup>cr*(*S*¯ *cr*))/2 . (Note that *μcr*(*Scr*) is monotone increasing.) The continuity of *μcr* implies that there exists *ε* > 0 such that *μcr*(*S*¯ *cr*) + *γ* < *<sup>μ</sup>cr*(*Scr*) for each *Scr* <sup>∈</sup> *S*0 *cr* <sup>−</sup> (<sup>1</sup> <sup>+</sup> *kcr*)*ε*, *<sup>S</sup>*<sup>0</sup> *cr* . It follows from (16) that there exists time *T<sup>ε</sup>* > 0 so that *X*(*t*) and *Scr*(*t*) satisfy

$$S\_{\mathcal{C}\mathcal{T}}^{0} - \varepsilon < S\_{\mathcal{C}\mathcal{T}}(t) + k\_{\mathcal{C}\mathcal{T}}X(t) < S\_{\mathcal{C}\mathcal{T}}^{0} + \varepsilon \text{ for each } t \ge T\_{\mathcal{C}}.\tag{20}$$

Assume now that *<sup>X</sup>*(¯*t*) ≤ *<sup>ε</sup>* for some ¯*<sup>t</sup>* ≥ *<sup>T</sup>ε*; then we obtain from (20)

$$S\_{cr}^{0} > S\_{cr}(\overline{\mathfrak{f}}) \ge S\_{cr}^{0} - k\_{cr}X(\overline{\mathfrak{f}}) - \varepsilon \ge S\_{cr}^{0} - (1 + k\_{cr})\varepsilon\_{rr}$$

i.e., *Scr*(¯*t*) <sup>∈</sup> *S*0 *cr* <sup>−</sup> (<sup>1</sup> <sup>+</sup> *kcr*)*ε*, *<sup>S</sup>*<sup>0</sup> *cr* . Hence,

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

It follows then that *<sup>X</sup>*(*t*) <sup>≥</sup> *<sup>e</sup>*(*t*−¯*t*)*γX*(¯*t*). If there exists *<sup>t</sup>*<sup>1</sup> <sup>≥</sup> ¯*<sup>t</sup>* such that *<sup>X</sup>*(*t*1) = *<sup>ε</sup>*, then at every time *<sup>t</sup>*<sup>2</sup> <sup>≥</sup> *<sup>t</sup>*<sup>1</sup> with *<sup>X</sup>*(*t*2) = *<sup>ε</sup>* we have *<sup>d</sup> dt <sup>X</sup>*(*t*2)=(*μcr*(*Scr*(*t*2)) <sup>−</sup> *<sup>D</sup>*¯ )*X*(*t*2) <sup>≥</sup> *γε* <sup>&</sup>gt; 0. Hence there exists time *T*<sup>1</sup> > *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, \mathcal{S}\_{cr}) : X \ge \varepsilon, \, 0 \le \mathcal{S}\_{cr} \le \mathcal{S}\_{cr}^{0} \}.$$

For *<sup>X</sup>* <sup>&</sup>gt; 0 and *Scr* <sup>∈</sup> (0, *<sup>S</sup>*<sup>0</sup> *cr*) we define the following Lyapunov function

$$V = V(X, S\_{cr}) = \int\_{\mathcal{X}}^{X} \frac{\eta - \mathcal{X}}{\eta} d\eta + \int\_{S\_{cr}}^{S\_{cr}} \frac{\mathcal{X}(\mu\_{cr}(\xi) - D)}{\vec{D}(S\_{cr}^{0} - \xi)} d\xi.$$

The derivative *<sup>d</sup> dt<sup>V</sup>* of *<sup>V</sup>* along the solutions of (18) and (19) is presented by

$$\begin{split} \frac{d}{dt}V &=& \frac{X - \bar{X}}{X} (\mu\_{cr}(S\_{cr}) - D)X + \frac{\bar{X}(\mu\_{cr}(S\_{cr}) - \bar{D})}{\bar{D}(S\_{cr}^{0} - S\_{cr})} \left(-k\_{cr}\mu\_{cr}(S\_{cr})X + D(S\_{cr}^{0} - S\_{cr})\right) \\ &=& X(\mu\_{cr}(S\_{cr}) - D) \left(1 - \frac{\bar{X}k\_{cr}\mu\_{cr}(S\_{cr})}{\bar{D}(S\_{cr}^{0} - S\_{cr})}\right) \\ &=& X(\mu\_{cr}(S\_{cr}) - D) \left(1 - \frac{S\_{cr}^{0} - S\_{cr}}{S\_{cr}^{0} - S\_{cr}} \cdot \frac{\mu\_{cr}(S\_{cr})}{\bar{D}}\right) \le 0 \end{split}$$

for each *Scr* <sup>∈</sup> (0, *<sup>S</sup>*<sup>0</sup> *cr*) 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 *<sup>D</sup>* <sup>&</sup>gt; *<sup>D</sup>*(2) *cr* <sup>=</sup> *<sup>μ</sup>*(*S*<sup>0</sup> *ph*, *<sup>S</sup>*<sup>0</sup> *cr*) then the model (1)–(3) possesses two equilibrium points—the wash-out equilibrium *E*<sup>0</sup> = (0, *S*<sup>0</sup> *ph*, *<sup>S</sup>*<sup>0</sup> *cr*) and the interior equilibrium *E*2, such that *E*<sup>0</sup> is locally asymptotically stable and *E*<sup>2</sup> is locally asymptotically unstable. Using the reduced model (17) it can be shown that the restriction *E*¯ <sup>0</sup> = (0, *S*<sup>0</sup> *cr*) of the wash-out equilibrium *<sup>E</sup>*<sup>0</sup> is globally asymptotically stable if *<sup>D</sup>* <sup>&</sup>gt; *<sup>D</sup>*(2) *cr* <sup>=</sup> *<sup>μ</sup>cr*(*S*<sup>0</sup> *cr*). 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 <sup>D</sup>* <sup>&</sup>gt; *<sup>D</sup>*(2) *cr holds true. Then for any initial point* (*X*(0), *Scr*(0)) <sup>&</sup>gt; <sup>0</sup> *the corresponding solution of (17) converges asymptotically towards the equilibrium E*¯ <sup>0</sup> = (0, *S*<sup>0</sup> *cr*)*.*

**Proof.** Choose some *<sup>D</sup>*¯ <sup>0</sup> <sup>&</sup>gt; *<sup>D</sup>*(2) *cr* and consider system (18) and (19), where *<sup>D</sup>*¯ is replaced by *<sup>D</sup>*¯ 0. Assume that lim*t*→<sup>∞</sup> *<sup>X</sup>*(*t*) = *<sup>X</sup>*<sup>∗</sup> <sup>&</sup>gt; 0. Then Barb˘alat's Lemma [26] applied to Equation (18) implies 0 <sup>=</sup> lim*t*→<sup>∞</sup> *<sup>d</sup> dt <sup>X</sup>*(*t*) = lim*t*→∞(*μcr*(*Scr*(*t*)) <sup>−</sup> *<sup>D</sup>*¯ <sup>0</sup>)*X*∗, which means that lim*t*→<sup>∞</sup> *<sup>μ</sup>cr*(*Scr*(*t*)) = *<sup>D</sup>*¯ <sup>0</sup> <sup>&</sup>gt; *<sup>μ</sup>cr*(*S*<sup>0</sup> *cr*). From the continuity of *μcr*(·) it follows that there exists time *<sup>T</sup>* <sup>&</sup>gt; 0 and a positive number *<sup>δ</sup>* such that *<sup>μ</sup>cr*(*Scr*(*t*)) <sup>−</sup> *<sup>μ</sup>cr*(*S*<sup>0</sup> *cr*) ≥ *δ* for all *<sup>t</sup>* <sup>≥</sup> *<sup>T</sup>*. The latter inequality leads to *<sup>d</sup> dt 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 *<sup>d</sup> dt Scr*(*t*) = *<sup>D</sup>*¯ <sup>0</sup>(*S*<sup>0</sup> *cr* <sup>−</sup> *Scr*(*t*)), which implies lim*t*→<sup>∞</sup> *Scr*(*t*) = *<sup>S</sup>*<sup>0</sup> *cr*, 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.** *<sup>D</sup>* <sup>=</sup> 0.08 <sup>∈</sup> (*D*(1) *cr* , *<sup>D</sup>*(2) *cr* )

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

*E*<sup>2</sup> = (0.0491, 0.126, 0.0153) and *E*<sup>3</sup> = (0.0318, 0.328, 0.116),

such that *E*<sup>2</sup> is locally asymptotically unstable, *E*<sup>3</sup> is the globally asymptotically stable equilibrium point according to Theorem 2. The wash-out equilibrium *E*<sup>0</sup> = (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*<sup>3</sup> 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*, *Scr*), (*X*, *Sph*) and (*Sph*, *Scr*) 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.

**Figure 4.** *D* = 0.08. (**Left**): time evolution of solutions *X* (solid line), *Sph* (dash-dot line) and *Scr* (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*, *Scr*)-phase plane with three different initial points, denoted by circles. The corresponding equilibrium components of *E*<sup>3</sup> are marked by a solid box, of *E*<sup>2</sup> are denoted by a box. The dash line presents the a projection of invariant plane *Scr* + *kcrX* = *S*<sup>0</sup> *cr*.

**Example 2.** *<sup>D</sup>* <sup>=</sup> 0.095 <sup>&</sup>gt; *<sup>D</sup>*(2) *cr* <sup>≈</sup> 0.0865

In this case there exists only one interior equilibrium point *E*<sup>2</sup> = (0.0514, 0.0987, 0.00191) which is locally asymptotically unstable. The wash-out equilibrium *E*<sup>0</sup> = (0, 0.7, 0.3) is the globally asymptotically stable steady state according to Theorem 3.

The global stability of *E*<sup>0</sup> 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*<sup>0</sup> 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*, *Scr*), (*X*, *Sph*) and (*Sph*, *Scr*) 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.

**Figure 5.** *D* = 0.08. Projections of the trajectories in the (*X*, *Sph*)-phase plane (**left**) and in the (*Sph*, *Scr*)-phase plane (**right**) with three different initial points, denoted by circles. The corresponding equilibrium components of *E*<sup>3</sup> are marked by solid boxes, of *E*<sup>2</sup> are denoted by boxes. The dash lines present projections of the invariant planes *Sph* + *kphX* = *S*<sup>0</sup> *ph* (left) and *Sph* <sup>−</sup> *KScr* <sup>=</sup> *<sup>S</sup>*<sup>0</sup> (right).

**Figure 6.** *D* = 0.095. (**Left**): time evolution of solutions *X* (solid line), *Sph* (dash-dot line) and *Scr* (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*, *Scr*)-phase plane with three different initial points, denoted by circles. The corresponding equilibrium components of *E*<sup>0</sup> are marked by a solid box, of *E*<sup>2</sup> are denoted by a box. The dash line presents a projection of the invariant plane *Scr* + *kcrX* = *S*<sup>0</sup> *cr*.

**Figure 7.** *D* = 0.095. Projections of the trajectories in the (*X*, *Sph*)-phase plane (**left**) and in the (*Sph*, *Scr*)-phase plane (**right**) with three different initial points, denoted by circles. The corresponding components of *E*<sup>0</sup> are marked by solid boxes, of *E*<sup>2</sup> are denoted by boxes. The dash lines present projections of the invariant planes *Sph* + *kphX* = *S*<sup>0</sup> *ph* (left) and *Sph* <sup>−</sup> *KScr* <sup>=</sup> *<sup>S</sup>*<sup>0</sup> (**right**).

#### **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 *μ*(*Sph*, *Scr*) 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*<sup>0</sup> = \$ 0, *S*<sup>0</sup> *ph*, *<sup>S</sup>*<sup>0</sup> *cr*% , 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*<sup>2</sup> and *E*<sup>3</sup> do exist: *E*<sup>3</sup> is defined for *D* ∈ \$ *<sup>D</sup>*(1) *cr* , *<sup>D</sup>*(2) *cr* % , and *E*<sup>2</sup> exists if *D* ∈ \$ *<sup>D</sup>*(1) *cr* , *Dcr*% . Local stability analysis shows that *E*<sup>2</sup> is locally asymptotically unstable and *E*<sup>3</sup> is locally asymptotically stable where they exist, *<sup>E</sup>*<sup>0</sup> is locally asymptotically stable for *<sup>D</sup>* <sup>&</sup>gt; *<sup>D</sup>*(2) *cr* . Two types of bifurcations of the equilibria occur, a saddle node bifurcation at *<sup>D</sup>* <sup>=</sup> *<sup>D</sup>*(1) *cr* where *<sup>E</sup>*<sup>2</sup> and *<sup>E</sup>*<sup>3</sup> coalesce, and a transcritical bifurcation at *<sup>D</sup>* <sup>=</sup> *<sup>D</sup>*(2) *cr* , where *<sup>E</sup>*<sup>0</sup> coincides with *<sup>E</sup>*<sup>3</sup> and *<sup>E</sup>*<sup>3</sup> disappears for *<sup>D</sup>* <sup>&</sup>gt; *<sup>D</sup>*(2) *cr* . Practically, the bifurcation values *<sup>D</sup>*(1) *cr* and *<sup>D</sup>*(2) *cr* of *<sup>D</sup>* 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*<sup>3</sup> when *D* takes values within certain bounds, *D* ∈ \$ *<sup>D</sup>*(1) *cr* , *<sup>D</sup>*(2) *cr* % , *<sup>D</sup>*(2) *cr* <sup>=</sup> *<sup>μ</sup>*(*S*<sup>0</sup> *ph*, *<sup>S</sup>*<sup>0</sup> *cr*). 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 *<sup>D</sup>* [27]. Choosing *<sup>D</sup>* in the interval \$ *<sup>D</sup>*(1) *cr* , *<sup>D</sup>*(2) *cr* % ensures practically long-term sustainability of the bioremediation process in the bioreactor. On the other hand, large values of the dilution rate *<sup>D</sup>*, *<sup>D</sup>* <sup>&</sup>gt; *<sup>D</sup>*(2) *cr* <sup>=</sup> *<sup>μ</sup>*(*S*<sup>0</sup> *ph*, *<sup>S</sup>*<sup>0</sup> *cr*), 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*<sup>0</sup> = (0, *S*<sup>0</sup> *ph*, *Scr*<sup>0</sup> ) 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. Both 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**

