Global stability of a continuous bioreactor model under persistent variation of the dilution rate

1.   Introduction

Continuous tank bioreactors are useful for producing chemical compounds, cultivating biomass, and for biological wastewater treatment ^[1,2,3]. In these bioreactors, the dilution rate and inlet substrate concentration have an important effect on the behavior of the biomass and substrate concentrations. Also, an overlarge dilution rate or equivalently low hydraulic retention time leads to biomass washout in most cases ^[4,5,6,7]. Adequate operation conditions of continuous bioreactors with avoidance of biomass washout can be studied through global and local stability analysis ^[3,7,8]. It has been found that in the chemostat model with non-monotonic growth function, two stable biomass equilibria may occur for certain values of the dilution rate, corresponding to biomass washout and to positive biomass. The initial condition determines which equilibrium is reached ^[3,7]. Some examples are the model for phenol and p-cresol mixture degradation in a continuously stirred tank bioreactor ^[9]; the microalgae Droop model ^[10]; the model for microalgae culture in presence of nitrifying bacteria ^[11]; the model of anaerobic digester ^[12]. In addition, if the yield coefficient is a function of the substrate concentration, a stable limit cycle exists for a certain range of the dilution rate, what is consistent with experimental results. Bioreactor operation with these limit cycles is usually avoided in practice ^[13,14,15].

Local stability can be studied by the indirect Lyapunov method, the main results include equilibrium points related to biomass existence and biomass washout, and the effect of dilution rate and inlet substrate concentration on the existence and stability of equilibria ^[7,16,17]. In contrast, global stability for constant dilution rate can be assessed by Lyapunov function theory, and the main results include the determination of global stability of the equilibrium associated with positive biomass concentration ^[9,11,18,19], stability of the equilibrium associated to biomass extinction ^[9,19,20], effect of kinetic parameters on the stability of equilibria ^[20], convergence of the additive combination of biomass, substrate and product concentrations ^[9,19]; domain of attraction of stable equilibria (see ^[11]); range of the input substrate concentration for global stability of the positive equilibrium point ^[21].

Some of these works are discussed at what follows. In ^[19], the global stability of a chemostat with overflow metabolism is studied. The state variables of the model are the concentrations of biomass, substrate, and by-product. The by-product secretion has an inhibitory effect on the growth of microorganisms, whereas the by-product excretion rate function is non-smooth. The existence and global stability of the equilibrium corresponding to positive biomass are proved, and the parameter conditions related to the existence of by-product are determined. In ^[9], the local and global stability of a model for phenol and cresol mixture degradation in a continuously stirred bioreactor is studied. The state variables are the concentrations of biomass, phenol and cresol. The biomass-specific growth rate involves inhibition caused by the substrate concentrations. The local asymptotic stability of the equilibrium point and the existence of local bifurcation in dependence on the dilution rate ($D$) are studied. It was found that two positive biomass equilibria exist, and the interval of $D$ for the existence of each is determined. Global asymptotic stability of one equilibrium is established. In ^[18], a stability analysis is performed for a simplified model of anaerobic digestion. The model corresponds to continuously stirred bioreactor, with two biomasses (acidogenic and methanogenic microorganisms) and three substrates (simple substrate, acetic acid, and ammonia). The specific growth rate of the methanogenic microorganisms involves inhibition dependence on acetic acid and ammonia. There are four equilibria, some of which involve the existence of one biomass. Their local stability analysis is proved for different cases of the input concentration of simple substrates (${s}^{o}$) whereas global stability is proved for some cases of low ${s}^{o}$. In ^[11], a stability analysis is performed for a microalgal pond with nitrification, whose model corresponds to that of a continuously stirred bioreactor with two biomass species (nitrifying biomass and microalgal biomass) and two substrates (ammonium concentration and nitrate concentration). Based on the Lyapunov function method, the global asymptotic stability of the equilibria is proved. There are different equilibrium points, corresponding to the following cases: washout of both biomasses; existence of one biomass; and existence of two biomasses. The conditions for the existence of these equilibrium points are determined, and their local and global stability are determined. In ^[22] a global stability analysis is performed for a wastewater treatment process whose model corresponds to continuously stirred bioreactor with three substrates (organic nitrogen, NH₄⁺-N, and NO₂^-+NO₃^-). The specific reaction rates are Monod type, and the dilution rate is considered constant, whereas the inlet concentration of organic nitrogen is considered time-varying and bounded, so that the system converges to a compact set instead of an equilibrium point. Only normal operation is considered, with no biomass washout. The convergence sets are determined, and global convergence is proved. To this end, dead-zone Lyapunov functions are used.

In cases that the dilution rate ($D$) is neither fixed nor controlled, the state variables converge to a region instead of an equilibrium point. This case is addressed in ^[3], where the stability of a simple chemostat model is studied, considering the concentrations of biomass and substrate as state variables, non-monotonic specific growth rate, and random disturbances on the input flow rate. The perturbations on the input flow rate are modeled through bounded random fluctuations. Positive constant lower and upper bounds are considered for the dilution rate. Biomass concentration can exhibit either extinction, weak persistence or strong persistence, depending on the bounds of the dilution rate constraint. The occurrence of these biomass behaviors is proved, and the respective dilution rate range is given. In addition, the convergence regions of substrate and biomass concentrations are determined for the case of biomass persistence. To this end, a comparison of solutions of scalar differential equations, and theory of asymptotically autonomous dynamical systems are used.

In this work, the global stability of a model of continuous bioreactor is studied. The concentrations of biomass and substrate are the state variables, whereas non-monotonic and monotonic functions of substrate concentration are used for the specific growth rate, and the inlet substrate concentration is considered constant. Also, the dilution rate is time varying but bounded, thus leading to the convergence of biomass and substrate concentrations to a compact set instead of an equilibrium point. Based on Lyapunov function theory with dead-zone modification, global stability is studied and proved. Finally, these results are illustrated by numerical simulation of a continuous bioreactor, considering non-monotonic and a monotonic specific growth rates. The main contributions with respect to closely related studies are:

●     Contribution CA. The convergence regions of the substrate and biomass concentrations are determined as a function of the upper bound of the dilution rate ($D$), and the global asymptotic convergence to these compact sets is proved, considering the persistent variation of $D$. In contrast, in stability studies of open loop continuous bioreactor model based on Lyapunov function, $D$ is commonly considered constant and substrate concentration converges to a point, not a compact set (for instance ^{[9,11,18,19,23]}. In addition, the working volume is not assumed constant, which is in contrast to other studies based on the Lyapunov function, for instance ^[3,11,18,22].

●     Contribution CB. Several improvements are proposed in the stability analysis, including the definition of a new dead zone Lyapunov function and the properties of its gradient. These improvements allow proving global convergence of substrate and biomass concentrations to their compact sets, while tackling: a) The interwoven and nonlinear nature of the dynamics of biomass and substrate concentrations, b) the non-monotonic nature of the specific growth rate; c) the time-varying nature of the dilution rate. These improvements are a basis for further global stability analysis of bioreactor models exhibiting convergence to a compact set instead of convergence to an equilibrium point. In contrast: i) Challenges a and c are not considered in ^[22], which is based on dead-zone Lyapunov function theory; ii) other global stability studies based on dead-zone Lyapunov function (for instance ^{[24,25,26,27]}) are applied to closed-loop systems and are not straightforwardly applicable to the considered bioreactor; iii) global stability studies based on Lyapunov function prove state convergence to compact sets, including the finite-time Lyapunov theory ^{[28,29,30,31,32,33]} and the ultimate bound approach ^{[34,35,36,37,38]} but they are not straightforwardly applicable to the considered bioreactor.

The work is organized as follows. The bioreactor model is presented in Section 2. The main results are presented in Section 3. The discussion is presented in Section 4. Numerical simulations are presented in Section 5. The conclusions are drawn in Section 6.

2.   Bioreactor model

Consider a continuous flow bioreactor described by model ^[3,8,39,40]:

where $x$ and $s$ are the concentrations of biomass and substrate, respectively; $\mu$ is the specific growth rate; ${y}_{s}$ is the substrate to biomass yield, ${y}_{s}\mu x$ is the specific substrate consumption rate; ${s}_{in}$ is the input substrate concentration; ${Q}_{in}$ is the input flow rate of medium, ${Q}_{out}$ is the output flow rate of medium, and $v$ is the volume of the medium in the bioreactor vessel; and $D\left(t\right) = {Q}_{in}/v$ is the dilution rate. If the bioreactor works as an original chemostat, the input and output flow rates are the same (${Q}_{in}$ = ${Q}_{out}$) so that the working volume $v$ is constant ^[3]. The $t$ argument is hereafter dropped for the sake of simplicity.

The following assumptions are considered:

Assumption 1. The biomass and substrate concentrations feature $x\in \left(0, \infty \right)$, $s\in \left(0, {S}_{in}\right)$ ^[33].

Assumption 2. ${y}_{s}$ and ${s}_{in}$ are positive and constant ^[3,39,40].

Assumption 3. The specific growth rate $\mu$ is a continuous function of the substrate concentration ($s$): it is non-negative, it is either monotonic or non-monotonic, it may involve growth threshold, and it satisfies the following features:

●     Non-monotonic case:

where ${S}_{\mu mx} = \underset{s\in (0, {s}_{\mathit{in}})}{\mathit{arg}\mathit{sup}}\mu$, and $\underline{s}$ is a constant value in the range $(0, {S}_{\mu mx})$.

●     Monotonic case:

where $\underline{s}$ is a constant value in the range $(0, {S}_{in})$.

Remark 2.1. The features of the non-monotonic growth rate (2) are based on ^[3,40,41], considering the growth threshold of ^[42]. Nonmonotone functions describe the case that substrate limits growth at low concentrations, but it is inhibitory at high concentrations ^[43,44,45]. As some examples of substrate induced inhibition, Candida utilis is inhibited by butanol ^[46]. A. succinogenes is inhibited by sugar ^[47], Nitrobacter is inhibited by nitrite and Nitrosomas is inhibited by ammonia ^[48]. A common nonmonotone function is the Andrews model ^[3,48,49].

where ${k}_{s}$, ${k}_{i}$ are saturation and inhibition constants, and $s$ is the substrate concentration. The coordinates of the maximum are:

where ${\mu }_{max}$ is the maximum specific growth rate attainable in presence of inhibition, and ${S}_{\mu mx}$ is the substrate concentration at maximum specific growth rate. Growth is inhibited and ${\mu }_{s}$ is decreasing for ${s > S}_{\mu mx}$ ^[48].

Remark 2.2. The features of the monotonic growth rate (3) are based on ^[17,39,41], considering the growth threshold of ^[42]. A basic example of monotonic growth rate (3) with growth threshold is obtained by using the Monod growth rate function:

where ${\mu }_{mx} > {k}_{o}$, whereas ${\mu }_{mx}$, ${k}_{s}$, ${k}_{o}$ are positive.

Remark 2.3. Growth ceasing for low values of substrate concentration has been verified for some strains. This is represented by using zero reaction rate for substrate concentrations lower than a threshold $\underline{s}$ ^[42,50].

Remark 2.4. In the case that $s = {s}_{in}$ for $t = {t}_{o}$, Eq (1b) leads to $ds/dt = -{y}_{s}\mu x < 0$, what implies $s < {s}_{to}$ for $t > {t}_{o}$, where ${s}_{to}$ is $s$ at initial time.

3.   Analysis of convergence to compact set

In this section, the global stability of model (1) with convergence to a compact set is studied, considering non-monotonic and monotonic growth functions:

ⅰ)    The dynamics of transformed states of biomass and substrate concentrations and the summary of the stability analysis procedure is given in subsection 3.1. The dynamics of the transformed states is a necessity for the Lyapunov-based stability analysis.

ⅱ)     The convergence regions of the substrate and biomass concentrations are given in terms of the upper bound of the dilution rate (D), and the global asymptotic convergence to these compact sets is stated in subsection 3.2; the non-monotonic growth function (2) is considered in Theorem 3.1, and the monotonic growth function (3) is considered in Theorem 3.2.

ⅲ)     An overview of the proposed improvements of the stability analysis necessary for proving convergence of the substrate and biomass concentrations to their compact sets is given in subsection 3.3, including the proposed dead-zone Lyapunov function and the defined properties of the gradient $W$.

ⅳ)    The relationship of results stated in Theorems 3.1 and 3.2 with equilibrium points is given in subsection 3.4.

3.1.   Dynamics of the transformed states and overview of the stability analysis procedure

Consider the bioreactor model (1), subject to assumptions 1 to 3. Let $z = {y}_{s}x+s$. Differentiating with respect to time and using the time derivatives (1a), (1b) yields

Let $\overline{z} = {y}_{s}x+s-{s}_{in}$. Differentiating with respect to time and arranging, yields

The error of substrate concentration is defined as $e = s-{s}^{*}$. The constant ${s}^{*}$ is:

●     The lowest constant value that satisfies $\mu = {D}_{max}$ at $s = {s}^{*}$ for non-monotonic growth function (2). In this case, two $s$ values satisfy $\mu = {D}_{max}$, and ${s}^{*}$ is the lowest of them.

●     The constant substrate concentration value that satisfies $\mu = {D}_{max}$ at $s = {s}^{*}$ for monotonic growth function (3). In this case, only one $s$ value satisfies $\mu = {D}_{max}$.

Differentiating $e$ with respect to time, yields $\dot{e} = \dot{s}$. Substituting the time derivative of $s$ (1b) and arranging, yields

Arranging, yields

The overall steps of the procedure used in the proofs of Theorems 3.1 and 3.2 are: definition of the error $e$ as the difference between the substrate concentration and the upper limit of its convergence region; determination of the time derivative $de/dt$; definition of the new state $z$, consisting on the weighted addition of the biomass and substrate concentrations; determination of the time derivative $dz/dt$; definition of the Lyapunov function for $z\left({V}_{z}\right)$ and determination of its time derivative $\left(d{V}_{z}/dt\right)$; definition of the features of the Lyapunov function for $e\left({V}_{e}\right)$ and the features of its gradient $\left(W\right)$, using dead-zone modification; determination of the $\left(d{V}_{e}/dt\right)$ expression in terms of the gradient $W$; definition of the dead zone $W$ function; determination of the ${V}_{e}$ function using the defined $W$ function; determination of the $d\left({V}_{e}+{V}_{z}\right)/dt$ expression in terms of $W$; arrangement of the $d\left({V}_{z}+{V}_{e}\right)/dt$ expression in terms of a non-positive function of $W$, by applying the properties of the defined dead zone $W$ function and determining the properties of its interaction with the $\mu$, $\overline{z}$ terms; determination of the convergence of the substrate concentration; determination of the convergence of biomass concentration.

This procedure aims at identifying the convergence region of the substrate concentration that leads to non-positive nature of the $d\left({V}_{e}+{V}_{z}\right)/dt$ expression, what in turn allows to prove the convergence. The dead-zone modification is used in the Lyapunov function ${V}_{e}$ to facilitate the study of convergence to a compact set instead of an equilibrium point. Dead-zone Lyapunov functions have been mainly used for robust control design: recent studies are presented in ^{[24,25,26,27]} whereas early studies are presented in ^[51,52,53]. In addition, an application to open loop systems is presented in ^[22].

3.2.   Convergence of the substrate and biomass concentrations

Theorem 3.1. Consider the bioreactor model (1) subject to assumptions 1 to 3, with the non-monotonic growth function $\mu$ satisfying conditions (2), and the dilution rate $D$ is time-varying and subject to saturation:

where ${\mu }_{sin}$ is the $\mu$ value at $s = {s}_{in}$. Then:

ⅰ) the substrate concentration $s$ converges asymptotically to ${{\mathrm{\Omega }}}_{s} = \left\{s:0 < s\le {s}^{*}\right\}$ where ${s}^{*}$ is the lowest constant value that satisfies $\mu = {D}_{max}$ at $s = {s}^{*}$;

ⅱ) the biomass concentration $\left(x\right)$ converges asymptotically to ${{\mathrm{\Omega }}}_{x} = \left\{x:\frac{1}{{y}_{s}}\left({s}_{in}-{s}^{*}\right)\le x < \frac{1}{{y}_{s}}{s}_{in}\right\}$.

Proof. Expression (8) involves both $s$ and $\overline{z}$ with nonlinear dependence on $s$, and time varying dilution rate (D). Then, choosing a Lyapunov function that allows proving the convergence of $e$ is overly complex. To this end, the Lyapunov function ${V}_{e}$ is defined in terms of its gradient $W$, and the properties of ${V}_{e}$, $W$ are also defined at this point, but the ${V}_{e}$, $W$ expressions are defined at the last part of the stability analysis. The main features of ${V}_{e}$ are:

The time derivative of ${V}_{e}$ can be expressed as

Where the gradient of ${V}_{e}$ is:

The properties of $W$ in terms of $e$ are:

The ${V}_{e}$ function can be determined from the integrated form of Eq (13):

The resulting ${V}_{e}$ function fulfills properties (11). Combining Eq (12) with the $de/dt$ expression (8), yields

Which can be rewritten as

Where ${\varepsilon }_{a}$, ${\varepsilon }_{b}$ are positive constants fulfilling

From the definition of ${s}^{*}$, the dependence of $\mu$ on $s$ (2) and the limitation of ${D}_{max}$ (10) it follows that

Hence, $\mu -{D}_{max} > 0fore > 0$. Considering the limitation of $D$ (9), we get

Combining with the $W$ properties (14), we have $W\left(\mu -D\right) = 0$ for $s\le {s}^{*}$ and $W\left(\mu -D\right) > 0$ for $s > {s}^{*}$. Therefore,

Hence, a preliminary $W$ function is chosen to be

Which does not satisfy the last of $W$ properties (14), because of the non-monotonic nature of $\mu$. From Eq (17) and the fact that $\mu = {\mu }_{sin}$ for $s = {s}_{in}$, it follows that $\mu$ is higher than the line that connects the points $\left({s}^{*}, {D}_{max}\right)$ and $\left({s}_{in}, {\mu }_{sin}\right)$:

Therefore, the $W$ function is chosen to be

Which satisfies $W$ properties (14). In addition, the resulting Lyapunov function ${V}_{e}$ can be obtained by using Eqs (15) and (20):

Equations (14), (18), (19) and (20) imply

The proof is given in Appendix A1. From this property and the positiveness of $x$ stated in assumption 1, it follows that

Substituting in Eq (16), yields

In view of the effect of $\overline{z}$, we consider the Lyapunov function for it:

which exhibits the following properties:

Differentiating ${V}_{z}$ with respect to time, using Eq (7), yields

Combining with Eq (23), yields

Where

Substituting into Eq (26), yields

Considering assumption 1, yields

Arranging and integrating, yields

where ${V}_{eto}$ is ${V}_{e}$ at initial time, whereas ${V}_{zto}$ is ${V}_{z}$ at initial time. From (28) it follows that ${V}_{e}\in {\mathcal{L}}_{\mathcal{\infty }}, {V}_{z}\in {\mathcal{L}}_{\mathcal{\infty }}, {W}_{a}^{2}\in {\mathcal{L}}_{1}$. Applying the Barbalat's lemma ^[54], yields

And consequently,

Accounting for the $W$ properties (14) and the substrate positiveness constraint stated in assumption 1, we deduce that $s$ converges asymptotically to ${{\mathrm{\Omega }}}_{s} = \left\{s:0 < s\le {s}^{*}\right\}$.

From Eq (7) it follows that $\overline{z} = {y}_{s}x+s-{s}_{in}$ converges exponentially to zero. This result and the asymptotic convergence of $s$ to ${{\mathrm{\Omega }}}_{s}$ imply that $x$ converges asymptotically to

This completes the proof.

Remark 3.1. The value $\mu = {D}_{max}$ holds true for two values of substrate concentration ($s$), for non-monotonic specific growth rate $\mu$, and ${s}^{*}$ is the lowest of them.

Remark 3.2. The fulfillment of the $W$ properties (14) leads to the fulfillment of the ${V}_{e}$ properties (11). In addition, ${V}_{z}$ (24) exhibits the properties (25).

Remark 3.3. In Theorem 3.1, the conditions $0 < {D}_{max} < {\mu }_{sin}$, ${D}_{max} < \underset{s\in \mathit{ }(0, {s}_{\mathit{in}})}{\mathit{sup}}\mu$ imply

because ${\mu }_{sin}\le \underset{s\in \mathit{ }(0, {s}_{\mathit{in}})}{\mathit{sup}}\mu$.

Remark 3.4. Condition $0 < {D}_{max} < {\mu }_{sin}$ implies that only one $s$ value $\left({s}^{*}\right)$ satisfies both ${\left.\mu \right|}_{s = {s}^{*}} = {D}_{max}$ and $s < {s}_{in}$.

Theorem 3.2. Consider the bioreactor model (1) subject to assumptions 1 to 3, with the monotonic growth function $\mu$, satisfying conditions (3), and the dilution rate $D$ is time varying and subject to saturation:

Where ${\mu }_{sin}$ is the $\mu$ value at $s = {s}_{in}$. Then:

ⅰ)    the substrate concentration converges asymptotically to ${{\mathrm{\Omega }}}_{s} = \left\{s:0 < s\le {s}^{*}\right\}$ where ${s}^{*}$ is a constant substrate concentration value that satisfies $\mu = {D}_{max}$ at $s = {s}^{*}$;

ⅱ)   the biomass concentration $\left(x\right)$ converges asymptotically to ${{\mathrm{\Omega }}}_{x} = \left\{x:\frac{1}{{y}_{s}}\left({s}_{in}-{s}^{*}\right)\le x < \frac{1}{{y}_{s}}{s}_{in}\right\}$.

Proof. The transformed states $e$, $\overline{z}$ are defined as

The time derivative of $e$ is given by Eq (8), whereas the time derivative of $\overline{z}$ is given by Eq (7). Consider the Lyapunov function

Where

and the gradient function $W$ satisfies the properties in terms of $e$:

Differentiating $V$ (31) with respect to time, using the time derivatives of $e$ and $\overline{z}$ given by Eqs (8), (7), yields

From the definition of ${s}^{*}$, the dependence of $\mu$ on $s$ (3), and the limitation of ${D}_{max}$ (30), it follows that $\mu > {D}_{max}$ for $s > {s}^{*}$ and $\mu = {D}_{max}$ for $s = {s}^{*}$. Hence, $\mu -{D}_{max} > 0$ for $e > 0$ and $\mu -{D}_{max} = 0$ for $e = 0$. Considering the limitation of $D$ (29), we have

Then, considering the $W$ properties (32), the exact $W$ function is chosen as

Which satisfies properties (32). In addition, the resulting Lyapunov function ${V}_{e}$ can be obtained from Eq (15):

Where ${\mu }_{e}$ is the $\mu$ function with the state transformation $s = e+{s}^{*}$. From the $W$ definition (35) and properties (34), (32) it follows that

Combining with the second subequation of (32), yields

From the first subequation of (32), it follows that

Combining with (37), yields

Using this result in Eq (33), yields

Considering assumption 1, yields

Arranging and integrating, yields

Hence, ${V}_{e}\in {\mathcal{L}}_{\mathcal{\infty }}, {V}_{z}\in {\mathcal{L}}_{\mathcal{\infty }}, {W}^{2}\in {\mathcal{L}}_{1}$. Applying the Barbalat's lemma ^[54], yields

And consequently,

Accounting for the $W$ definition (35) and the substrate positiveness constraint stated in assumption 1, we deduce that $s$ converges asymptotically to ${{\mathrm{\Omega }}}_{s} = \left\{s:0 < s\le {s}^{*}\right\}$.

From Eq (7) it follows that $\overline{z} = {y}_{s}x+s-{s}_{in}$ converges exponentially to zero. This result and the asymptotic convergence of $s$ to ${{\mathrm{\Omega }}}_{s}$ imply that $x$ converges asymptotically to

This completes the proof.

Remark 3.5. Theorem 3.1 considers nonmonotonic growth rate function $\mu$ whereas Theorem 3.2 considers monotonic growth rate function.

Remark 3.6. Any time-varying $D$ behavior satisfying $D\in \left[0, {D}_{max}\right]$, ${0 < D}_{max} < {\mu }_{sin}$ is allowed in Theorems 3.1 and 3.2, including the sinusoidal, random type or a combination of them. Also, constant $D$ behavior is allowed.

Remark 3.7. For the nonmonotonic case, the geometric properties of ${V}_{e}$ stated in Eq (11) arise as a consequence of the properties of $W$ stated in Eq (14). Indeed, properties (11) are straightforwardly verified for the Lyapunov function (21). For the monotonic case, the geometric properties of ${V}_{e}$ arise as a consequence of the properties of $W$ stated in Eq (32).

3.3.   Summary of the proposed improvements of the stability analysis and dead zone Lyapunov function

Proposed improvements of the stability analysis. The improvements proposed in the stability analysis are detailed in the proofs of Theorems 3.1 and 3.2, most of them related to the proposed dead zone Lyapunov function for the tracking error $e = s-{s}^{*}$, ${V}_{e}$, and its gradient $W$. These improvements can be summarized as follows:

●     Before defining ${V}_{e}$, the $d{V}_{e}/dt$ equation is expressed in terms of the gradient $W$ (see Eqs (12) and (16)) and the properties of $W$ relevant for proving convergence of the tracking error $e$ are defined using dead zone modification, see Eq (14).

●     Then, the Lyapunov function ${V}_{e}$ is expressed as an integral in terms of its gradient $W$, see Eq (15).

●     A particular dead zone $W$ function is proposed, see Eqs (20) and (36), fulfilling the stated $W$ properties. As a result, a new expression is obtained for ${V}_{e}$, with dead zone modification, see Eqs (21) and (36).

●     The $d\left({V}_{z}+{V}_{e}\right)/dt$ expression is arranged in terms of a non-positive function of $W$, by applying the properties of the $W$ stated in Eq (14) and determining the properties of the interaction of $W$ with the $\mu$, $\overline{z}$ terms, see Eq (22) and subsequent equations.

Proposed dead zone Lyapunov function and its gradient. The overall Lyapunov function is

Where ${k}_{z} = {D}_{max}/\left(4{\varepsilon }_{b}{y}_{s}{x}_{min}\right)$, $\overline{z} = {y}_{s}x+s-{s}_{in}$, and ${V}_{e}$ is function of the error $e = s-{s}^{*}$, and it is defined in terms of its gradient:

The properties of $W$ in terms of the error $e = s-{s}^{*}$ are:

A particular $W$ function that satisfies these properties, for the case of non monotonic specific growth rate (2), is:

For which the Lyapunov function ${V}_{e}$ is found by using Eq (38):

A particular $W$ function that satisfies properties (39), for the case of monotonic specific growth rate (3), is:

For which the Lyapunov function ${V}_{e}$ is:

Where ${\mu }_{e}$ is the $\mu$ function with the state transformation $s = e+{s}^{*}$. Both $W$ functions (40) and (42) fulfill

In the case of Monod growth function with growth threshold (6), the Lyapunov expression ${V}_{e}$ (43) gives:

Remark 3.8. The results stated in Theorems 3.1 and 3.2 and their proofs do not allow stating the input to state stability (ISS) of system (1), in the context of the ISS notion stated in ^[55]. One way to prove that a system is ISS is by proving that the existence of an ISS Lyapunov function. An ISS Lyapunov function satisfies

Where $y$ is the state vector, $u$ is the input; ${\alpha }_{1}$, ${\alpha }_{2}$, ${\alpha }_{3}$, ${\alpha }_{4}$ are ${\mathcal{K}}_{\mathcal{\infty }}$ functions, according to ^[55]. In Eqs (26) and (27):

ⅰ)    The term $-{\varepsilon }_{a}{y}_{s}{W}^{2}{x}_{min}$ is not ${\mathcal{K}}_{\mathcal{\infty }}$, because ${W}^{2}$ is not strictly increasing, according to Eq (14). Indeed, it is not increasing for $e\le 0$.

ⅱ)     The Lyapunov function $V = {V}_{e}+{V}_{z}$ does not satisfy condition ($45$) because ${V}_{e}$ is not strictly increasing.

Then, ISS property cannot be concluded. However, Theorems 3.1 and 3.2 and their proofs allow stating the boundedness for bounded input $D$. Equation (28) implies that $e, z$ are bounded, and consequently $x, s$ are bounded, provided bounded $D$.

3.4.   Relationship with equilibrium points

Notice that the limits of the convergence regions of $x, s$ stated in Theorem 3.1 are related to one of the equilibrium points of model (1). Model (1) with nonmonotonic growth rate (4) and constant dilution rate $D$ exhibits the equilibrium points:

The upper limit of the convergence region of $s$, that is, ${s}^{*} = \frac{-\left({D}_{max}-{\mu }_{o}\right)-\sqrt{{\left({D}_{max}-{\mu }_{o}\right)}^{2}-4\left(\frac{{k}_{s}}{{k}_{i}}\right){D}_{max}^{2}}}{2{D}_{max}\;/{k}_{i}}$, obtained from condition $\mu = {D}_{max}$ at $s = {s}^{*}$, that is, ${\left.{\mu = \mu }_{o}\frac{s}{{k}_{s}+s+\frac{{s}^{2}}{{k}_{i}}}\right|}_{s = {s}^{*}} = {D}_{max}$ corresponds to equilibrium point ${E}_{1}$ with constant ${D = D}_{max}$. The limits of the convergence region of $x$, that is $\frac{1}{{y}_{s}}\left({s}_{in}-{s}^{*}\right) = \frac{1}{{y}_{s}}\left\{{s}_{in}+\frac{\left({D}_{max}-{\mu }_{o}\right)+\sqrt{{\left({D}_{max}-{\mu }_{o}\right)}^{2}-4\left({k}_{s}/{k}_{i}\right){D}_{max}^{2}}}{2{D}_{max}\;/{k}_{i}}\right\}$ and $\frac{1}{{y}_{s}}{s}_{in}$, correspond to the equilibrium point ${E}_{1}$, with constant ${D = D}_{max}$ and constant $D =$0, respectively.

Also, notice that the limits of the convergence regions of $x, s$ stated in Theorem 3.2 are related to one of the equilibrium points of model (1). Model (1) with monotonic growth rate ${\mu = \mu }_{mx}\frac{s}{{k}_{s}+s}$ and constant dilution rate $D$ exhibits the equilibrium points

The upper limit of the convergence region of s, obtained from condition $\mu = {D}_{max}$ at $s = {s}^{*}$, is ${s}^{*} = \frac{{k}_{s}{D}_{max}}{{\mu }_{mx}-{D}_{max}}$, and it corresponds to equilibrium point ${E}_{1}$ with constant ${D = D}_{max} < {\mu }_{mx}$. The limits of the convergence region of $x$, that is, $\frac{1}{{y}_{s}}\left({s}_{in}-{s}^{*}\right) = \frac{1}{{y}_{s}}\left({s}_{in}-\frac{{k}_{s}{D}_{max}}{{\mu }_{mx}-{D}_{max}}\right)$ and $\frac{1}{{y}_{s}}{s}_{in}$, correspond to the equilibrium point ${E}_{1}$ with constant ${D = D}_{max}$ and constant $D = 0$, respectively. In summary, the limits of the convergence regions of $x, s$ stated in Theorems 3.1 and 3.2 correspond to the equilibrium point ${E}_{1}$ with constant ${D = D}_{max}$ and constant $D = 0$.

Remark 3.9. Consider the case that model (1) with constant dilution rate

and Haldane kinetics (4), so that

This $D$ value does not fulfill conditions (9) and (10), so that neither Theorem 3.1 nor 3.2 apply, and the constant nature of $D$ implies that the states converge to equilibrium point, not to compact set, so that dead-zone modification of the Lyapunov function is not necessary. In this case, the $D$ value implies that there is only two equilibrium points, and Eqs (46)–(48) lead to:

In addition, equilibrium ${E}_{1}$ is a pitchfork bifurcation point with respect to parameter $D$, so that its local stability cannot be stated through eigenvalues, and other stability tools for bifurcations are necessary, see for instance ^[56,57].

4.   Discussion

4.1.   General discussion of results

The stability analysis summarized in Theorems 3.1 and 3.2 determines global convergence of biomass and substrate concentrations from the model (1) towards a compact set. Assumptions 1 to 3 are considered, the specific growth rate function is assumed to be a nonlinear function of substrate concentration, either non-monotonic (2) or monotonic (3), with constant inflow substrate concentration. The dilution rate $\left(D\right)$ is varying within the compact set ${\Omega }_{D} = \left\{D: 0\le D\le {D}_{max}\right\}$, implying that the states $x, s$ converge to the compact sets ${{\mathrm{\Omega }}}_{s}$, ${{\mathrm{\Omega }}}_{x}$ instead of an equilibrium point. Common global stability of open loop continuous bioreactor model based on the Lyapunov function (for instance ^{[9,11,18,19,23]}) cannot be used for this case, because they consider constant dilution rate $D$ and convergence of substrate and biomass concentrations towards a point, not a compact set.

To this end, methods that consider state convergence to compact set include the comparison of solutions of scalar differential equations and theory of asymptotically autonomous dynamical systems ^[3]; dead zone Lyapunov functions ^{[24,25,26,27]}, and ultimate bound approach ^{[34,35,36,37,38]}. The dead zone Lyapunov method is commonly used for control design, so that its application for bioreactor systems as the model (1) is hampered by the nonlinear $\mu$ term and the complex nonlinear connection between the dynamics of $x, s$. The application of the ultimate bound approach is also hampered by these facts.

Results related to Contribution CA. Theorems 3.1 and 3.2 indicate the global convergence of substrate concentration $s$ to the compact set ${{\mathrm{\Omega }}}_{s} = \left\{s:0 < s\le {s}^{*}\right\}$ and the global convergence of biomass concentration $x$ to the compact set ${{\mathrm{\Omega }}}_{x} = \left\{x:\left(1/{y}_{s}\right)\left({s}_{in}-{s}^{*}\right)\le x < \left(1/{y}_{s}\right){s}_{in}\right\}$, and also provide a persistent condition $\left(0 < {D}_{max} < {\mu }_{sin}\right)$ that leads to biomass persistence (avoidance of washout) for the case of nonmonotonic growth rate function, where ${s}^{*}$ is the solution of $\mu = {D}_{max}$ for $s = {s}^{*}$. The limit value ${s}^{*}$ is a function of ${D}_{max}$, so that the widths of ${\Omega }_{s}$, ${\Omega }_{x}$ depend on the width of ${\Omega }_{D}$.

Results related to Contribution CB. The stability proofs of Theorems 3.1 and 3.2 in subsection 3.2 use the Lyapunov function method with dead zone modification. However, the stability analysis is significantly improved, and a new dead-zone Lyapunov function $\left({V}_{e}\right)$ is proposed, in order to prove global convergence of substrate and biomass concentrations to their compact sets ${\Omega }_{s}$, ${\Omega }_{x}$ while tackling the non-linear nature of the specific growth rate $\mu$, the interwined dynamics of the concentrations of substrate and biomass, and the variation of the dilution rate $\left(D\right)$ within the compact set ${\Omega }_{D}$. Indeed, the properties of ${V}_{e}$ and of its gradient $W$ are defined and used for the stability analysis but the exact $W$ and ${V}_{e}$ expressions are provided at the last part of the stability analysis. Some remarkable features are:

●     The properties of $W$ in Eq (39) indicate that $W$ is continuous, piecewise, dependent on $e = s-{s}^{*}$, increasing for $e > 0$ and involving a vanishing region for $e\le 0$.

●     In the case of non-monotonic $\mu$ (2), $\mu$ is only increasing for a certain range of $e > 0$, so that it does not lead to increasing $W$. Thus, $W$ is defined as a piecewise linear function of $e$, including no section of the $\mu$ curve, see Eq (40).

●     In the case of monotonic $\mu$ (3), $\mu$ leads to increasing $W$ function. Thus, the $W$ definition includes a section of the $\mu$ curve, see Eq (42).

●     The equation of ${V}_{e}$ as an integral in terms of $W$ (38) allows obtaining the ${V}_{e}$ functions (41) and (43).

In addition, the improved stability analysis (Theorems 3.1 and 3.2 and their proofs) is useful as the procedure for further global stability analysis of bioreactor models exhibiting convergence to a compact set instead of convergence to an equilibrium point, allowing the determination of:

●     R1: range of initial values of the state variables and parameter values leading to global convergence to the non-washout convergence set (washout avoidance) in the case in the case of nonmonotonic growth rate function.

●     R2: upper bound of the transient values of the state variables as a function of their initial values and model parameters

●     R3: exponential, asymptotic, monotonic or nonmonotonic character of the global convergence of the state variables to the compact set

In agreement with R1, the global stability of continuous bioreactors through Lyapunov function allows checkinh domains of attraction of normal operation equilibrium point, under a constant dilution rate ^[8,18,58].

From the analysis of the relationship between the results stated in Theorems 3.1 and 3.2 with the equilibrium points (in subsection 3.4) it would seem that local stability (equilibrium points and eigenvalues) is enough to study the convergence of the state variables to compact sets. However, local stability does not allow us to obtain results R1, R2 and R3, because it is limited to a close neighborhood of the equilibrium point, and it disregards the transient behavior in regions far from the equilibrium point, as can be deduced from ^[12,56].

4.2.   Future work

An important application field of the stability analysis based on the Lyapunov function presented in Theorems 3.1 and 3.2 and their proofs, is the anaerobic digestion process ^[18]. Anaerobic digestion is useful for wastewater treatment with biogas production. It involves a complex process with hydrolysis, acidogenesis, acetogenesis and methanogenesis steps. It is prone to biomass washout, which consists of microorganism extinction, it occurs for overlarge values of dilution rate and is related to high substrate concentration ^[8,12,59]. The stability of equilibrium points related to normal operation point and washout has been studied through equilibrium points, diagrams of equilibrium versus dilution rate, bifurcation diagrams with dilution rate and input substrate concentration as bifurcation parameters ^[12,56,59]. In contrast, global stability based on the Lyapunov function gives as result the conditions and the state space region that implies global stability of the normal operation equilibrium point, considering constant dilution rate ^[18,58].

Also, the developed stability analysis can be extended to more complex bioreactor models, and modifications of the model (1), including:

●     Model of anaerobic upflow fixed-bed digester with partial biomass attachment on support:

where $\alpha$ is a positive constant that represents the model heterogeneity ^[7,49].

●     Model with maintenance term $mx$:

where $m$ is a postive constant ^[42].

●     Model with several biomass species:

where ${x}_{i}$ is the mass density of species i ^[60].

These models, combined with either time varying dilution rate $\left(D\right)$ or time varying inflow substrate concentration $\left({s}_{in}\right)$ leads to different convergence sets for the states $x, s$. In the case of varying ${s}_{in}$, the model for $\overline{z} = {y}_{s}x+s-{s}_{in}$ is modified. The limits of the convergence regions depend on the limits of $D$, ${s}_{in}$.

Some biological models exhibit limit cycles under constant parameters (see ^[61,62]). It has been suggested that their oscillations are a consequence of competitive interaction and low nutrient supply ^[63,64]. Simple models with two microbial populations and Monod uptake terms are capable of describing the states of the system with stable oscillations ^[64]. As an example, in nitrification at the soil with contaminant ammonium plume, bacterial populations exhibit boom and boost dynamics ^[62]. In these cases, the Lyapunov function is a possible tool for examining the global convergence of the states towards the limit cycle, but other tools are also necessary for proving the existence and stability of the limit cycles (see ^[65,66,67]). Indeed, some biological systems are described in predator-prey type models, whereas these models may be taken to polar coordinates so that the differential equation of the radial coordinate indicates that the radius is stable, and a Lyapunov function can be defined in terms of its gradient, see ^[51,68]. In this case, the dead zone Lyapunov function can provide a more complete stability result.

5.   Numerical simulation

In this section, the convergence results stated in theorems 3.1 and 3.2 are illustrated through simulation, by showing the convergence regions of biomass and substrate concentrations $\left(x, s\right)$ and the asymptotic convergence of their trajectories towards to convergence regions. Bioreactor model (1) subject to assumptions 1 to 3 is considered, Theorem 3.1 is illustrated in the first simulation example, using a Haldane growth rate function, whereas theorem 3.2 is illustrated in the second simulation example, using a Monod growth rate function.

First simulation example. The simulation conditions correspond to chemostat with bacterial species and Haldane growth rate function (4):

The parameter values are given in Table 1.

For the growth rate function (51), the coordinates of the maximum are given by Eq (5). The ${D}_{max}$ value satisfies condition (10), and ${s}^{*}$ is the lowest solution of $\mu = {D}_{max}$:

what leads to

An example of ${D}_{max}$, ${s}^{*}$ values is given in Table 1. The growth function (51), its main points, the gradient $W$ (40) and Lyapunov function ${V}_{e}$ (41) are illustrated in Figure 1.

The time varying dilution rate is:

where ${a}_{1} = 0.93, {\tau }_{d1} = 45$ ${\tau }_{d2} = 45/10$ and ${{\mathrm{\Omega }}}_{t1} = \left[100,230\right]$ for ${D}_{max} = 1.65$ and for ${D}_{max} = 1$.

Simulations of substrate and biomass concentrations show that (Figure 2):

ⅰ)    For each ${D}_{max}$ value, the convergence region of the simulated $s$ trajectory is inside the region bounded by the computed limit ${s}^{*}$, thus verifying the validity of theorem 3.1, and the residual set ${{\mathrm{\Omega }}}_{s} = \left\{s:0 < s\le {s}^{*}\right\}.$

ⅱ)   the substrate concentration $s$ is outside the residual set ${\Omega }_{s}$ at initial time, it gets inside at 5.22 for ${D}_{max}$ = 1.65, and at 3.73 for ${D}_{max} = 1$, and it remains inside afterwards.

ⅲ)    For each ${D}_{max}$ value, the convergence region of the simulated $x$ trajectory is inside the region bounded by the computed limits $\left(1/{y}_{s}\right)\left({s}_{in}-{s}^{*}\right)$, $\left(1/{y}_{s}\right){s}_{in}$, thus verifying the validity of theorem 3.1, and the residual set ${{\mathrm{\Omega }}}_{x} = \left\{x:\left(1/{y}_{s}\right)\left({s}_{in}-{s}^{*}\right)\le x < \left(1/{y}_{s}\right){s}_{in}\right\}.$

Second simulation example. Continuous cultivation of Gluconacetobacter diazotrophicus is simulated, using conditions and data of batch cultivation with cane molasses as carbon source (see ^[69]), considering the Monod growth rate function:

Where ${\mu }_{mx} > {k}_{o}$.

The parameters of model (1) with this growth rate function were fitted to the batch data and are given in Table 1. For the growth rate function (52), the ${D}_{max}$ value satisfies condition (30) and ${s}^{*}$ is obtained from $\mu = {D}_{max}$:

what leads to

An example of ${D}_{max}$, ${s}^{*}$ values is given in Table 1. The growth function (52), its main points, the gradient $W$ (42) and Lyapunov function ${V}_{e}$ (44) are illustrated in Figure 3.

The time-varying dilution rate is:

Where ${a}_{1} = 0.9, {\tau }_{d1} = 120$ h, ${\tau }_{d2} = 120/9$ h and ${{\mathrm{\Omega }}}_{t1} = \left[300,550\right]$ h for ${D}_{max} = 0.0656$ h^-1; whereas ${a}_{1} = 0.91,$ ${\tau }_{d1} = 190$ h, ${\tau }_{d2} = 190/9$ h and ${{\mathrm{\Omega }}}_{t1} = \left[350,600\right]$ h for ${D}_{max} = 0.0291$ h^-1.

Simulations of substrate and biomass concentrations show that (Figure 4):

ⅰ)     for each ${D}_{max}$ value, the convergence region of the simulated $s$ trajectory is inside the region bounded by the computed limit ${s}^{*}$, thus verifying the validity of theorem 3.2, and the residual set ${{\mathrm{\Omega }}}_{s} = \left\{s:0 < s\le {s}^{*}\right\}$.

ⅱ)    the substrate concentration $s$ is inside the residual set ${\Omega }_{s}$ at the initial time, and it is outside only during [40.7,131.1] h for ${D}_{max} = 0.0656$ h^-1 and during [16.2,138.7] h for ${D}_{max} = 0.0291$ h^-1.

ⅲ)    for each ${D}_{max}$ value, the convergence region of the simulated $x$ trajectory is inside the region bounded by the computed limits $\left(1/{y}_{s}\right)\left({s}_{in}-{s}^{*}\right)$, $\left(1/{y}_{s}\right){s}_{in}$, thus verifying the validity of Theorem 3.2, and the residual set ${{\mathrm{\Omega }}}_{x} = \left\{x:\left(1/{y}_{s}\right)\left({s}_{in}-{s}^{*}\right)\le x < \left(1/{y}_{s}\right){s}_{in}\right\}$.

General discussion on the numerical simulations. In summary, the first and second simulation examples show that the considered dilution rate constraint (9) and (29) with ${D}_{max}$ conditions (10) and (30) leads to asymptotic convergence of the substrate concentration to the determined compact set ${{\mathrm{\Omega }}}_{s} = \left\{s:0 < s\le {s}^{*}\right\}$ and asymptotic convergence of biomass concentration to the determined compact set ${{\mathrm{\Omega }}}_{x} = \left\{x:\left(1/{y}_{s}\right)\left({s}_{in}-{s}^{*}\right)\le x < \left(1/{y}_{s}\right){s}_{in}\right\}$, in presence of varying dilution rate within the defined set ${\Omega }_{D}$. This illustrates Theorems 3.1 and 3.2, and verifies the validity of the residual sets ${\Omega }_{s} = \left\{s:0 < s\le {s}^{*}\right\}$, ${{\mathrm{\Omega }}}_{x} = \left\{x:\left(1/{y}_{s}\right)\left({s}_{in}-{s}^{*}\right)\le x < \left(1/{y}_{s}\right){s}_{in}\right\}$, where ${s}^{*}$ is the solution of $\mu = {D}_{max}$ for $s = {s}^{*}$.

6.   Conclusions

In this work, the global stability of a continuous bioreactor model is studied, with the concentrations of biomass and substrate as state variables. Non-monotonic and monotonic functions of substrate concentration are considered for the specific growth rate, whereas the dilution rate is time-varying but bounded, thus leading to state convergence to a compact set instead of an equilibrium point. Based on the Lyapunov function theory with dead-zone modification, the asymptotic convergence of the substrate and biomass concentrations is studied.

The main contributions with respect to closely related studies are: i) the convergence region of the substrate and biomass concentrations are determined as a function of the upper bound of the dilution rate ($D$), and the global asymptotic convergence to these compact sets is proved, considering monotonic and non-monotonic growth functions separately; ii) several improvements are proposed in the stability analysis, including the definition of a new dead zone Lyapunov function and the properties of its gradient. These improvements allow proving convergence of substrate and biomass concentrations to their compact sets, while tackling the interwoven and nonlinear nature of the dynamics of biomass and substrate concentrations, the non-monotonic nature of the specific growth rate, and the time-varying nature of the dilution rate. The proposed modifications are a basis for further global stability analysis of bioreactor models exhibiting convergence to a compact set instead of an equilibrium point. Finally, the convergence of the states and the biomass persistence under varying dilution rates are illustrated by simulation.

Acknowledgements

Acknowledgments: The work of A.R. and G.R. was supported by Universidad Católica de Manizales. The work of F.E.H. was supported by Universidad Nacional de Colombia Sede Medellín.

Conflict of interest

The authors declare that they have no conflict of interests.

Appendix

Proof of Eq (22).

From the $W$ properties (14) or alternatively from the definition of $W$(20) it follows that

From property (19) and the $W$ definition (20) it follows that

Combining with properties (18), (14), yields

That is,

Combining with result (A1), yields

From the $W$ definition (20) and properties (14) it follows that

Combining with Eq (A2), yields

This completes the proof.