Stability and Performance Analysis of Bioethanol Production with Delay and Growth Inhibition in a Continuous Bioreactor with Recycle

: This paper proposes and analyzes a mathematical model for the production of bioethanol in a continuous bioreactor with recycling. The kinetics correspond to the use of Saccharomyces bayanus for the fermentation of sugars found in wastewater from soft drinks. The proposed model considers product growth latency, which was experimentally found in batch studies of ethanol production. Fur-thermore, the inhibition effect of ethanol is expressed by a modified version of the classical Andrew’s model for substrate inhibition. The proposed model consists of only three ordinary differential equations containing a minimal number of operating parameters, which include the bioreactor residence time, glucose feed concentration, recycle ratio and the fraction of biomass removed from the reactor by the flow. The positivity and the boundedness of solutions of the model were confirmed under reasonable restrictions of parameters. The stability analysis showed that there is a value of residence time at which an exchange of stability occurs between the trivial washout and non-washout solutions. This critical value depends only on the substrate feed concentration, biomass death rate, recycle ratio and purge fraction. Dynamic simulations of the model were carried out for substrate concentration in the range of 100–250 g/L, commonly used for the production of ethanol. An inverse response due to the inhibition effects of ethanol was observed in the time evolution of substrate and biomass concentrations. Parametric studies showed that ethanol concentration increases with the recycle ratio, with the inverse of residence time and with the inverse of purge fraction. The effect of ethanol latency has, on the other hand, a substantial effect on ethanol concentration. Despite its unstructured nature and the fact that some parameters such as temperature and acidity were not taken into consideration, the proposed model managed to provide useful results on the bioreactor-settler stability and the effect of key parameters on its dynamic behavior, which could pave the way for future optimization studies.


Introduction
The growing demand for petrifaction and fossil fuels due to the swift revolution of automotive industries and modern societies has been inspired the study of developing alternative fuels for combustion engines for decades. Environmental pollution issues are also pushing scientists and industrialists to search for alternatives to fossil fuels [1]. Sincere attention has been given to inventing/discovering a sustainable, clean and cheap alternative to satisfy modern smart environmental needs [2,3].
Bioethanol is one of the best alternatives to oil-derived fossil fuels. Ethanol extracted from suitable renewable sources is carbon friendly and an attractive green energy to control environmental contamination and reduce dependence on petrifaction and fossilized carbon

Mathematical Model
In this section, we present the model equations that describe the dynamics of a wellstirred continuous bioreactor with cell recycling shown schematically in the diagram of Figure 1. The mass balances of substrate (S), biomass (X) and ethanol (E) yield the following ordinary differential equations V dE dt = −Fβ 2 E + Y e/x µ max · M(S, E)e −(pS) q · X · V + VγX, V is the rector volume, F the volumetric flow rate, S f the concentration of the substrate feed, µ max the maximum specific growth rate, Y x/s the biomass yield coefficient, Y e/x the ethanol/biomass yield coefficient, b H the biomass death rate coefficient and γ is the kinetic constant of ethanol production by maintenance. The specific growth rate is denoted by M(S, E) and is a function of both substrate and ethanol concentrations, since it is known that the growth rate is inhibited by ethanol concentration [2,3]. A number of growth rate models were proposed in the literature, but here, we adopt the following expression proposed by [3] M(S, E) = S K s + S + K e E 2 , where K s is the saturation constant and K e the inhibition constant by ethanol. The form of Equation (4) is inspired by the expression used in the case of pure substrate inhibition S K s +S+K e S 2 , which is known as Andrew's growth model [17]. As was mentioned in the introduction, the choice of the expression in Equation (4) was guided by the experimental work carried out in [3] on the growth of Saccharomyces bayanus on wastewater from soft drinks. Moreover, to account for latency in ethanol production, a term dependent only on sugar concentration [3] was included in the ethanol mass balance (Equation (3)). The reactor residence time is defined by In Equation (2), the term CX represents the biomass concentration in the flow leaving the separating unit. The value of the concentrating factor (C) depends upon the design and operation of the settling unit. It is also highly dependent on recycle properties such as settling, and compressibility behavior. The term β 1 , on the other hand, represents the fraction of the biomass that leaves the reactor (purge fraction). A mass balance around the settling unit shows that the maximum value of the concentrating factor is given by Thus, the maximum value of the product R(C − 1) is The coefficient β 2 in Equation (3) represents, on the other hand, the fraction of ethanol in the feed. For practical reasons, this fraction is negligeable or nonexistent, and therefore a zero value is assumed for β 2 . The model is supplemented with an initial profile satisfying It can be seen that, for particular laboratory environmental conditions, the kinetic parameters K e , K s , Y x/s , Y e/x , b H , γ, p,q and µ max are all fixed. The operating parameters that one may vary are S 0 , R, β 1 and τ.
The values of the kinetic parameters are shown in Table 1.

The Dimensionless Model
The model is rendered dimensionless using the following variables: [S * = S/K s ], [X * = X/Y x/s K s ], [E * = E/K s ] and [t * = µ max t]. The system of ODE (1)-(3) can be written in a dimensionless form by All the dimensionless parameters here are non-negative. A short tabular description of the model parameters is given in Nomenclature section. Looking at the dimensionless parameter R * = R(C − 1) that represents the effective recycle ratio, it can be seen that the maximum value of R * is 1. On the other hand, the choice of β 1 = 1 gives a continuous flow reactor with recycle. Therefore, the cases R * = 1 and 0 < R * < 1 with β 1 = 1 define a flow reactor with ideal and non-ideal recycle respectively. The values of the model dimensionless parameters (computed using the dimensional values of Table 1) are given in Table 2. Besides the fixed kinetic parameters, the nominal values of operating parameters are also shown in the table. Specifically, the glucose feed concentration S 0 has a nominal value of 100 g/L, which corresponds to S * 0 = 1.53. It should be noted that the range of sugar concentrations commonly used for the production of ethanol is about 100-250 g/L.

Mathematical Analysis
Uniqueness, uniform boundedness and invariance of solutions play an important role in mathematical modeling, since these stated properties reflect the physically meaningful features of the problem. To that end, we begin here by studying the uniqueness properties of solutions followed by the boundedness and invariance of model solutions over the following physically meaningful domain To facilitate our study a variable Z ∈ R 3 is defined by Using the new set of redefined variables, the model (9)-(11) is written as where F : R 3 → R 3 , and F i : R 3 → R, i = 1, 2, 3 are given by Here, the modified function It is obvious that 0 ≤ M(Z) ≤ 1, ∀ Z ∈ Ω. The time evolution of the growth rate of biomass (Z 1 ), consumption of substrate (Z 2 ), and ethanol production (Z 3 ) are described by (12). Now, it is our aim to investigate if the solutions of (12) are positive. This step is necessary to see if the model is physically meaningful. To facilitate the analysis, we enforce a few physically reasonable restrictions H 1 : Parameters used here belong to R + .
The assumptions H 1 , H 2 , H 3 and H 4 represent physically meaningful restrictions on the model parameters.

Theorem 1.
Assume that H 1 -H 4 hold, then there exists a bounded unique solution to the initial value problem (12) over J.
Proof. Here, the gradient of F(Z) exists for all Z ∈ Ω. We can see that |M(Z)| < 1 for all Z ∈ Ω. Applying H 1 -H 4 leads to (12) has a unique solution which is also bounded [19,20].

Proof. Let us define
Then, using the fact that where α = |M(Z)e −(p * Z 2 ) q | << 1 which completes the proof. Now, we focus on showing the invariance of integral Z(t) of (12) in Ω.

Theorem 3. Ω is a positively invariant set for the integral Z if H 1 holds.
Proof. Consider the domains on the boundary and the edges of Ω. Then for the Set 1. Thus, Z 1 is inside Ω. For Set 2, we have and this guarantees that Z 2 resides inside the domain Ω. Now, for Set 3, we have which also shows that Z 3 ∈ Ω. We can show that for Set 4. Thus, all the cases confirm that Z ∈ Ω. Hence, the solution domain Ω is invariant [21]. Now, we analyze the parameter dependence of the solution trajectories of the dimensionless model (9)-(11). Let Now, let recall that We need the following lemma to analyze parameter estimation.

Lemma 1.
Let Z(t) be the trajectory of the model with initial value Z 0 , and γ, C ≥ 0 be the constants such that the following linear growth condition holds [22] then the inequality Theorem 4. Let Y(t) and Z(t) be two integrals of the system of differential Equations (9)- (11) with parameter sets S1 and S2, R * ≥ B * 2 and initial profiles Y 0 and Z 0 , respectively; then holds where C * ≥ 0 and C 1 ≥ 0 are constants.
Proof. Let Y(t) and Z(t) be two integrals for the system (9)-(11), then where C 1 and C 2 depend on the parameters. Similarly, using the fact that 0 ≤ M where C 3 depends on the parameters. Additionally, and A combination of all the above inequalities yields which completes the proof.

Equilibria
The objective in this section is to find steady-state solutions of the nonlinear system It can be noted that Z * 0 = (S * 0 , 0, 0) is always a washout solution of the nonlinear system F(Z) = 0 in Ω. This solution is not desirable since it indicates that no biological process has occurred. We now look for a non-trivial solution (non-washout) to the nonlinear system of equations F(z) = 0 in Ω, considering some physically meaningful restrictions on parameter values. Let where Now, adding F 1 (Z) = 0 and F 2 (Z) = 0 yields Substituting the value of X * (Equation (15)) in F 1 (Z) = 0 yields From equation (F 1 (Z) = 0), we have Thus, substituting Equation (15) in Equation (17) results in Note that the right side of Equation (18) must be positive since the left side is E * 2 . Thus, the expression inside the square root in Equation (16) is positive. Now, substituting Equations (15), (17) and (18) where

Stability of the Equilibria
Stability analysis is at the heart of dynamical analysis. Only stable solutions can be noticed experimentally. Thus, it is very important to analyze the stability of an ODE system. Mostly researchers rely on the eigenvalues of the Jacobian matrix relevant to the nonlinear dynamical system. The Jacobian of the system of ODE (9)-(11) can be presented by The eigenvalues of J(S * , X * , E * ) evaluated at the washout steady-state can be presented by Therefore, it follows from the sign of eigenvalues that the washout branch is always stable if the following restriction holds Performing similar results as shown in [23], we can see that the washout solution is globally stable when the following restriction of parameters holds b * H ≥ we conclude from Equation (20) that the washout steady-state is always stable provided that It can be noted that the washout solution is never stable when R * = β 1 . Therefore, we end up with the following result. The Jacobian matrix relevant to the non-washout solution branch can be presented by The term J(2, 2) was evaluated using the following relation Here the solutions are physically meaningful if A i > 0, i = 1, ..., 9, i = 6, 7. Thus the characteristic polynomial of J(S * , X * , E * ) is given by, Thus, it is sufficient (using the well-known Routh-Hurwitz criteria [24]) to establish that c 1 > 0, c 3 > 0 and c 1 c 2 − c 3 > 0, to confirm that all the roots of (λ) have negative real parts.
Clearly c 1 > 0. To prove that c 3 > 0, we have For positive solutions, we have A i > 0, i = 1, ..., 9, i = 6, 7, Hence, c 3 > 0. We now examine the expression c 1 c 2 − c 3 . After some manipulations, we find that Thus, we conclude, using Routh-Hurwitz criteria, that all the roots of the characteristic polynomial (λ) have negative real parts. Therefore, we have the following result: Lemma 3. Whenever the no-washout branch X * = (S * , X * , E * ) is positive, it is locally asymptotically stable.

Numerical Discussion
In this section, we move on to study the time evolution of the model solutions of Equations (9)- (11). The initial value problem was integrated using the Runge-Kutta based ODE45 solver of MATLAB. We start by showing the bioreactor dynamics for the case of β 1 = 1 and different values of the recycle ratio (R * = 0, 0.5 and 1). This corresponds to a continuous flow with no recycle (R * = 0), with a non-ideal recycle (R * = 0.5) and with an ideal recycle (R * = 1). The rest of the model's dimensionless parameters are shown in Table 2.
The simulations are shown in Figure 2 for startup conditions S = 100 g/L, X = 4.32 g/L and E = 0. This corresponds to (S * = 1.53, X * = 1, E * = 0). The figures show a similar trend for all values of R * . The substrate concentration decreases quickly for a short time, reaches a minimum and then increases to reach a constant value. This is a typical example of an inverse response. This behavior is also shown in the biomass profile as the biomass concentration reaches a maximum and then decreases due to the inhibition effects of ethanol. The profile of ethanol concentration shows, on the other hand, a steady increase until an asymptotic value is reached.
The effect of recycle ratio is shown in the diagram. It can be seen that bioethanol concentration substantially increases as the recycling ratio is increased. In fact the increase in bioethanol concentration is around 43% when the bioreactor is operated with no recycle (R * = 0) compared to R * = 0.5.
The effect of latency of ethanol production is investigated by numerically varying the parameter p * for the same residence time τ = 12 and start-up conditions as Figure 2 but with an ideal recycle R * = 1, since it yields the highest ethanol concentration. Although p * is not an operating parameter, it is important to study its effect on the performance of the process. Figure 3 shows the effect of different values of p * (p * = 0, 0.5, 1 and 10). First, it can be seen that the inverse response, which is the feature of the inhibition effect of ethanol, is also maintained in this figure; however, it is attenuated as the delay is increased. As expected, the ethanol concentration decreases as the delay is pronounced. For example, the decrease in ethanol concentration is around 23.8% when the delay is increased from p * = 0 to p * = 1.  The effect of varying the bioreactor residence time τ * is shown in Figure 4. τ * was varied from 1, 6, 12, to 20. A large ethanol concentration is obtained for small residence times (i.e., larger dilution rates). The increase in the residence time rate almost linearly affects the decrease in ethanol concentration. An increase in residence time from τ * = 1 to τ * = 6 decreases the asymptotic value of ethanol concentration by 10.5%.    Table 2.
Finally, we examine the effect of the purge fraction β 1 . In order to avoid a washout solution, as predicted by Equation (22) for R * = β 1 , and to allow for the study of variations of β 1 , the effective recycle ratio R * was fixed at 0.5 and β 1 was allowed to vary from 0.5, 0.75 to 1. It can be seen from Figure 5 that the highest ethanol concentration occurs at the smallest values of the purge fraction. However, a careful choice of values of β 1 and R * should be made to avoid washout conditions.

Conclusions
In this study, a mathematical reactor model with Andrew's growth rate and latency of ethanol production was developed to study the dynamical behaviour of continuous fermentation in reactor with recycle using parameter values from batch experimental data [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. The positivity and the boundedness of model solutions were confirmed under reasonable restrictions on parameters. The stability analysis of the steady state solutions yielded an analytical expression for the critical residence time at which a stability exchange between washout and non-trivial solution occurs. We observed that this critical value of residence time is dependent on the feed substrate concentration, purge fraction and the recycle ratio. Therefore, a careful choice of these operating parameters can avoid washout conditions. Numerical simulations provided useful qualitative trends of the dynamics of the process. An inverse response was observed in the profile of substrate and biomass concentrations. This dynamic behavior has practical implications when feedback control is desired to maintain the process at certain set points [12], since such systems are more difficult to control. It was found that the maximum ethanol concentration is obtained at the maximum allowable recycle ratio, the smallest purge fraction and largest dilution rate. However, care should be taken to optimize the selection of such parameters to avoid washout conditions. A final note should be made about the proposed model's strengths and limitations. Since the model kinetics were extracted from a validated batch experiment, and since the operating parameters were numerically varied around realistic values, this provided the obtained results with a fair degree of credibility. However, the proposed model could be improved by taking some important parameters that affect the fermentation process into consideration, such as medium temperature and acidity. The analysis carried out in the paper was limited to steady state behavior, but could be used in the future to address the complex phenomena associated with the occurrence of oscillatory behavior in continuous ethanol fermentation as result of dynamic bifurcations [11].