Article pubs.acs.org/IECR
Method for Regulating Oscillatory Dynamic Behavior in a Zymomonas mobiliz Continuous Fermentation Process Hangzhou Wang, Nan Zhang, Tong Qiu, Jinsong Zhao, and Bingzhen Chen* Department of Chemical Engineering, Tsinghua University, Beijing 100084, China S Supporting Information *
ABSTRACT: Periodically oscillating concentrations of product, biomass, and substrate have been observed in a continuous fermentations process with Zymomonas mobiliz both in experiments and simulations. This paper describes the characteristics of oscillatory phenomena based on Hopf bifurcation analysis and proposes a method to regulate the oscillatory dynamic behavior. In this method, these dynamic behaviors of oscillations in concentrations of product, biomass, and substrate were simulated, the relationship between amplitude/period and operational condition was modeled. This model could be used as a guide to adjust the operating conditions to attenuate unfavorable oscillations and mitigate corresponding operability problems when the oscillation itself is unavoidable.
■
Outlet ethanol concentration fluctuations thus significantly affect the subsequent ethanol distillation step. Moreover, more sugars could be discharged with the effluent under oscillatory conditions, and industrial ethanol yields calculated based on the sugars fed into fermentation systems, without deduction of unfermented sugars, would be compromised.20 Therefore, these oscillations risk increasing the average residual sugar, and corresponding decreases in the ethanol yield,4 so oscillations must be attenuated. This motivates studies of mechanisms that lead to this oscillatory behavior. A considerable amount of work has identified the causes of oscillatory behavior in fermentation processes, both experimentally and numerically. Lee et al.7,8 demonstrated that Z. mobiliz is a very promising organism for ethanol production because it has higher specific rates of growth and ethanol production, as well as a higher tolerance to ethanol. After studying continuous culture with glucose media, it was concluded that ethanol inhibition of growth is the cause of oscillatory phenomena. Jobses et al.9 also assumed that oscillations were mainly caused by product inhibition of enzymatic reactions and experimental studies revealed that high temperatures may enhance such inhibitory effects. In addition, oscillations at high product concentrations were difficult to avoid. Ghommidh et al.10 introduced a structuration of the cell population into a mathematical model, allowing any situation, from completely stable to completely unstable continuous cultures, to be described. It was further found that oscillatory steady states were achieved at low dilution rates. Bruce et al. 11 used extractive fermentation to improve the performance of Z. mobiliz in continuous culture during the conversion of concentrated substrates to ethanol, and used this to eliminate the oscillatory behavior which arises routinely in conventional fermentations of Z. mobiliz.
INTRODUCTION Bioethanol is an important renewable and sustainable alternative fuel source,1 both directly in the form of fuel ethanol and when blended with gasoline. The yeast Saccharomyces cerevisiae2,3 and the bacterium Zymomonas mobiliz2−4 are currently the most important microorganisms for industrial production of ethanol. Z. mobiliz produces less biomass than S. cerevisiae, while more carbon is directed to ethanol synthesis. It was reported that the ethanol yield of Z. mobiliz could be as high as 97% of the theoretical yield of ethanol from glucose,5 whereas only 90−93% can be achieved by S. cerevisiae.4 In addition, Z. mobiliz maintains a higher glucose metabolic flux that guarantees its higher ethanol productivity.5 Thus, Z. mobiliz has been promoted as a more promising microorganism than yeast for the industrial production of ethanol. Compared with batch operation, continuous fermentation can improve productivity by saving labor and maintenance costs and reducing capital investment on production facilities, thus it has been practiced for large scale production of fuel ethanol in industry.6 However, continuous ethanol fermentation with Z. mobiliz tends to be rich in dynamic instability, in which sustainable or damped oscillations have routinely been observed.7−25 The concentrations of product, biomass and substrate vary periodically in these oscillatory phenomena, especially at low growth rate and high ethanol concentrations,10 and these are characterized by long oscillatory periods and large oscillation amplitudes. Under certain operating conditions, the oscillation periods are sustained for 50 h for glucose, ethanol and biomass, and the amplitudes of oscillation are about 20− 100 kg/m3 for glucose and 50−90 kg/m3 for ethanol, as reported for continuous ethanol fermentation with Z. mobiliz.16 These oscillations negatively affect ethanol fermentation performance and result in increased residual sugar and decreased ethanol yields in the ethanol fermentation industry.26 The process oscillation causes downstream process fluctuations, particularly ethanol distillation, which requires relatively constant ethanol concentrations in the fermentation broth.21 © 2014 American Chemical Society
Received: Revised: Accepted: Published: 12399
March 4, 2014 June 20, 2014 July 10, 2014 July 10, 2014 dx.doi.org/10.1021/ie5009234 | Ind. Eng. Chem. Res. 2014, 53, 12399−12410
Industrial & Engineering Chemistry Research
Article
Jarzebski12 developed a tricompartment cell population model with kinetics, predicting a decrease in product inhibition with decrease in substrate concentration, and this model was able to portray all dynamic situations observed in continuous ethanol fermentations, depending on substrate concentration, from sustained oscillations to completely stable cultures. Hobley et al.13 proposed that the oscillatory mechanism is driven by the ethanol concentration history but Li et al.14 observed that the effect of ethanol concentration history on the fermentation capacity of Z. mobiliz is negligible, while an upward ethanol concentration change rate had an intense inhibitory effect on Z. mobiliz. Daugulis et al.15 and McLellan et al. 16 showed some numerical simulations of mathematical models that perfectly matched the results in their experiments. Gandhian et al.17−19 reported that ethanol yield and the production rate would be favored by minimizing ethanol inhibition, and explored the conditions for increasing the productivity and yield of ethanol by fermentation using mathematical modeling and nonlinear analysis. Membranes were positioned within a fermentor to continuously remove ethanol to increase productivity and stabilize the system by eliminating oscillations. In addition, experimental and theoretical studies illustrated that variables such as dilution rate and feed substrate concentration have an effect on static and dynamic behaviors in a continuous stirred tank fermentor. Bai et al.20,21 showed parameter variation could cause oscillations in a very high-gravity medium continuous ethanol fermentation, and the experimental results revealed that only Intalox ceramic saddles and wood chips could attenuate the oscillations, improving the ethanol production performance of the system. Diehl et al.22 analyzed the nonlinear dynamic behavior in a continuous ethanol fermentation process with Z. mobiliz, and presented a multivariable control strategy with a high yield solution branch inside the steady-state multiplicity region, that could maintain process stability. Li et al.23 developed three models to study the inhibitory effects that cause oscillatory behavior in alcohol fermentation via Z. mobiliz. The results showed that ethanol concentration alone did not cause periodic behavior during ethanol production but computer simulations indicated that oscillations were caused under certain circumstances. Sridhar24 reported that oscillatory behavior was caused by the existence of Hopf bifurcations, and demonstrated simple strategies with very minor changes in the input conditions that could eliminate the Hopf bifurcation points, to avoid oscillation. Our previous work25 analyzed and calculated Hopf points distribution in a Z. mobiliz continuous fermentation process to enhance the stability of the fermentation process and maintain process stability. As stated above, there are complex nonlinear dynamic behaviors,27 especially unfavorable oscillatory phenomena, in continuous ethanol fermentation processes with Z. mobiliz. Operating parameters, such as dilution rate and feed substrate concentration, influence both static and dynamic behaviors. Previous researchers, either by experimentation or numerical simulation, showed that optimal operating points and regions can be identified to avoid the oscillatory behavior. Additionally, with proper selection of operating conditions, these Hopf points can be eliminated, to avoid oscillations, directly resulting in a more stable and productive fermentation process. However, in some situations the oscillatory phenomena are unavoidable, especially when the Hopf points are very close to the maximum product concentration operating point, and a supercritical Hopf bifurcation is more likely to generate sustainable oscillations. As a result, it is necessary to attenuate
the oscillations, as a next-best or suboptimum solution, to improve ethanol production. Consequently, the relationship between characteristics under the dynamic behavior of oscillatory phenomena and operating conditions is helpful in a very practical sense. However, no quantitative study of such a relationship has yet been reported. In this paper, the oscillatory phenomena were thoroughly studied, and quantitative relationships between the amplitude/period and operating conditions were formulated. Experimental research suffers some disadvantages, such as (I) being unable to locate all points corresponding to certain dynamic behaviors because of the slow fermentation process dynamics, and (II) only a limited number of experiments can be performed, which may lead to some dynamic behavior being missed or neglected. In contrast, numerical analysis methods can enable complete and thorough analysis of dynamic systems. In our research, we numerically simulated the oscillations in a continuous fermentation process to study the relationships mentioned above. This paper proceeds as follows. First, we establish the basic procedure for regulation of oscillatory phenomena. Following this, the numerical model for continuous Z. mobiliz fermentation process is described. We then present simulations of oscillatory phenomena with varied operating conditions to obtain a series of data. By analyzing the collected data, quantitative relationships between amplitude/period and operating conditions are formulated, then further validated by simulation. Finally, we provide discussion and conclusions.
■
METHOD FOR REGULATING OSCILLATORY PHENOMENON
The basic procedure for regulating oscillatory phenomena in a dynamic system is summarized in Figure 1. Starting with a model of the fermentation process, this method was successively progressed by calculating the steady state solution of the dynamic system along with varying parameters, then identifying the steady state solutions to determine first whether or not there were Hopf points, and then locating them if they existed. After simulating any oscillatory phenomena occurring near to the Hopf singularity point, the relationship between oscillation parameters and operating conditions were formulated. The formulated model was useful for guiding practical operation. If the simulated result did not agree with the simulation or experimental results, more experimental data were gathered and coefficients of the dynamic model were regressed with the data. After updating the dynamic system model for the chemical process, the above steps were repeated until the calculated result agreed with simulated/experimental results.
■
PACKAGES FOR NUMERICAL ANALYSIS METHODS
Bifurcation theory suggests that limit cycles, which represent oscillatory phenomena, are intrinsically concomitant with Hopf singularity points. Obtaining the periodic solution of Hopf singularity point is a tedious numerical analysis problem but, fortunately, some software packages, such as Auto,28−30 content,31 matcont,32,33 DsTool,34 PyDSTool,35 XPPAUT36 and some work25,37−45 by our research group, can be used to implement numerical calculations. 12400
dx.doi.org/10.1021/ie5009234 | Ind. Eng. Chem. Res. 2014, 53, 12399−12410
Industrial & Engineering Chemistry Research
Article
Figure 2. Fermentation process.25
The mathematical model analyzed is based on the following simplifying assumptions:46 (1) The fermentor is a well-mixed system and is not controlled by diffusion. (2) The feed contains only substrate. (3) The specific growth rate follows a Monodtype equation. (4) The maximum specific growth rate is a function of the key component concentration. (5) For some of the parameters, k1, k2, k3, Ks, ms, mp, Ysx and Ypx, for which the temperature dependence had a negligible effect on the results, constant values were used. The equation of substrate consumption is as follows: ⎛ μ ⎞ + ms⎟Cx rs = ⎜ ⎝ Ysx ⎠
The product generation rate equation is defined as:
⎛ μ ⎞ + mp⎟⎟Cx rp = ⎜⎜ ⎝ Ypx ⎠ The growth rate of biomass rx is expressed in the form:
rx = μCx The rate of the key compound re is formulated as:
Figure 1. Diagram of the regulation of oscillatory phenomena
■
re = f (Cs)f (Cp)Ce
FERMENTATION PROCESS FLOWCHART The dynamic mathematical model for Z. mobiliz fermentation process has been developed over several years, and there are several reported mathematical models.9,10,12,15,16,24,27 A recently developed model27 incorporating a membrane to separate product is a more practical approach, and the calculations in this article are based on this model. Ethanol production by Z. mobiliz is subject to end-product inhibition. The product alters the cell membrane composition and inhibits enzymatic reactions such as carrier-mediated transport processes and metabolic syntheses. Figure 2 25 shows a flow sheet of the Z. mobiliz continuous fermentation process, and a membrane was placed at the bottom of the fermentor to remove ethanol produced during fermentation. The membrane has the selectivity toward ethanol within the reaction environment, with a gas/liquid sweep stream on the permeate side to remove product from the membrane surface. With the removal of ethanol, inhibitory effects are reduced and cell activity is maintained at a high level, making it possible to use a higher glucose feed concentration, and further resulting in increased fermentor cell density. The increased cell activity and cell density contribute to higher volumetric productivity and a smaller fermentor volume. Additionally, the use of a concentrated glucose feed reduces the amount of process water required.27
where f(Cs) is given by the Monod equation: f (Cs) =
Cs K s + Cs
and the function f(Cp) is defined in polynomial form: f (Cp) = k1 − k 2Cp + k 3Cp2
Thus, the equation for re can be substituted into the following expression: re =
k1 − k 2Cp + k 3Cp2 K s + Cs
CsCe
The exit fermentor dilution rate can be obtained as follows: A P Dout = Din − M m (Cp − CpM) VFρ and the corresponding exit membrane dilution rate becomes: A P DMout = DMin + M m (Cp − CpM) VMρ The dynamic model for the continuous fermentation process for substrate (s), microorganism (x), key compound (e), product in fermentor (p) and product on the membrane 12401
dx.doi.org/10.1021/ie5009234 | Ind. Eng. Chem. Res. 2014, 53, 12399−12410
Industrial & Engineering Chemistry Research
Article
Lyapunov coefficient l1 equals to −0.00153. The negative sign indicates that this point corresponds to a supercritical Hopf bifurcation and the limit cycle is stable; that is, operating points near this Hopf singularity point are prone to generating sustained oscillations. In Figure 3, the abscissa denotes the feedstock dilution rate, ranging from zero to a reasonable maximum. The dilution rate indicates the degree of feed stock volume diffusing into the reactor volume and has a relatively small value in practice. As mentioned above, the Hopf point was very close to the maximum conversion operation point, and the first Lyapunov coefficient of this point indicated that the Hopf point was prone to cause sustained oscillatory phenomena. These oscillations associated with limit cycles are described in Figure 4.25 Each
permeate side (pM) are described by the following set of ordinary differential equations: dCx dt
= μCx + DinCx0 − Dout Cx
dCs dt
⎛ μ ⎞ = −⎜ + ms⎟Cx + DinCs0 − Dout Cs ⎝ Ysx ⎠
dCe dt
=
dCp dt
dCpM dt
k1 − k 2Cp + k 3Cp2 K s + Cs
CsCe + DinCe0 − Dout Ce
⎛ μ ⎞ = ⎜⎜ + mp⎟⎟Cx + DinCp0 − Dout Cp ⎝ Ypx ⎠ AM Pm − (Cp − CpM) VF =
AM Pm (Cp − CpM) + DMinCpM0 − DMout CpM VF
The model parameters used in the current investigation were taken from the published literature and are given in Table S1.27
■
HOPF POINTS IN THE SYSTEM AND THE RELATED LIMIT CYCLES After calculating the steady state points of this dynamic system model for continuous fermentation process with Z. mobiliz, Figure 3 25 shows the product concentration of steady state
Figure 4. Homology limit cycles adjoining the Hopf point.25
limit cycle represents a sustained amplitude and fixed-period oscillation, in which the concentration of production, biomass and substrate were varied periodically. To regulate the oscillation, a quantitative study of the relationship between oscillatory characteristics and operating conditions was required.
■
DYNAMIC BIFURCATION ANALYSIS OF OSCILLATORY PHENOMENA According to bifurcation theory, whether experimentally observed or mathematically simulated, oscillations in a dynamic system must start at certain critical points in the system. From this singularity point, with varied initial parameter values, the process will suffer sustained oscillation. It has been proven that the singular points and limit cycles found in Z. mobiliz continuous fermentation processes are fairly common. As analyzed, these oscillations are produced from Hopf singularity points. The values for this Hopf point are listed in Table 1, and this point is referred to as HP1. In these oscillatory phenomena, the concentrations of product, biomass and substrate vary periodically with time and the characteristics of the oscillations are affected by operating conditions. As a result, it is important to investigate how these detailed oscillatory phenomena characteristics are influenced by operating conditions. Here, we numerically simulated the oscillations to obtain data and analyzed these to
Figure 3. Changes in product concentration with dilution rate.25
operating points with changes in dilution rate. In the case of Cs0 = 150.3 kg/m3, a Hopf point was detected at Din = 0.052 h−1, which coincides well with reported literature values.18,24 This Hopf point corresponds to the situation where the system Jacobian has one pair of complex conjugate eigenvalues crossing the imaginary axis and at the point of crossing, a limit cycle is born. Here, at the Hopf point, the product concentration is Cp = 58.050 kg/m3, which is quite close to the maximum point of product yield, Din = 0.046 h−1, Cp = 58.103 kg/m3. We can see that the Hopf point lies close to the maximum conversion point. At the same time, the Hopf point has a first 12402
dx.doi.org/10.1021/ie5009234 | Ind. Eng. Chem. Res. 2014, 53, 12399−12410
Industrial & Engineering Chemistry Research
Article
new operating condition. We then maintained this condition to study the oscillatory phenomena. We observed, first, that the oscillation was generated near Hopf singularity points and, second, that when operational conditions were varied slightly, the oscillation remained but its characteristics changed. We continued to raise the dilution rate, to 0.2 h−1, and then performed the dynamic simulation. The dynamic behavior curves are shown in Figure 6.
Table 1. Values of Hopf Point, HP1 parameter
value
unit
Din Ce Cp CpM Cs0 Cs Cx
0.05 0.06 58.05 57.98 150.30 6.10 1.77
h−1 kg/m3 kg/m3 kg/m3 kg/m3 kg/m3 kg/m3
formulate the relationship model between oscillatory amplitude/period and operating parameter values. We carried out simulations with gradually increased dilution rate and recorded the oscillatory characteristics data, finally obtaining the relationship between amplitude/period and operating conditions.
■
SIMULATION OF OSCILLATORY DYNAMIC BEHAVIOR Simulation of Oscillatory Dynamic Behavior When Din = 0.1 h−1. Because unfavorable oscillations were generated from the Hopf singularity point, we simulated the dynamic behavior near the Hopf point with slight changes in dilution rate to 0.1 h−1. The simulated curve is given in Figure 5.
Figure 6. Simulated dynamic behavior when Din =0.2 h−1.
Simulation of Oscillatory Dynamic Behavior when Din = 0.2 h−1. From Figure 6, we can see that, with the new dilution rate and other values at the Hopf point, the system dynamics behaved as a sustained oscillatory curve. After a short time, the system performed a stable periodic oscillation in concentration of product, for which the amplitude was calculated. From the dynamic curve, or corresponding data, the characteristics of the oscillatory phenomena can be found. The amplitude of the product concentration oscillation, Amp,Cpm, was 17.42 kg/m3, and the period time, T, was 9.46 h. Using the same simulation method, we gradually changed dilution rate, together with other unchanged values of Hopf point as initial values, to follow the dynamic behavior of the mathematical model and obtain the amplitude and period for every oscillation. These data are presented in Table 2, where the dilution rate varied from 0.1 to 0.9 h−1 in steps of 0.1 h−1.
Figure 5. Numerical run, illustrating limit cycle behavior: Din = 0.1 h−1, Cs0 = 150.3(kg/m3), with other conditions given in Table S1.
From Figure 5, we can see that when we took the Hopf points as the initial point in the mathematical model for the continuous fermentation process with Z. mobiliz, the system dynamic behaved as expected. There was an oscillation in the system, sustained variations in product concentration occurred periodically with time, and similar behaviors were observed in the concentrations of biomass and substrate. From Figure 5 the characteristics can be extracted from the dynamic result. The oscillation amplitude of product concentration, Cpm, was 12.19 kg/m3, while the period, T, was 15.12 h. Because the unfavorable oscillation was generated near to the Hopf singularity point, we took the values of the Hopf point as the initial value for the ordinary differential equation (ODE) model to perform the dynamic simulation, and changed the dilution rate from 0.05 h−1 at the Hopf point to 0.1 h−1 as a
Table 2. Oscillatory Amplitude and Period Change with Dilution Rate
12403
Din (h−1)
Amp,Cpm (kg/m3)
T (h)
0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90
12.19 17.42 20.14 22.20 23.97 25.55 26.99 28.32 29.57
15.12 9.46 7.09 5.76 4.91 4.32 3.88 3.53 3.26
dx.doi.org/10.1021/ie5009234 | Ind. Eng. Chem. Res. 2014, 53, 12399−12410
Industrial & Engineering Chemistry Research
Article
was obtained using the dynamic curve from Figure 8, with a corresponding amplitude 20.06 kg/m3 and period of 7.1 h. The
From Table 2, we can see that detailed oscillation characteristics, the amplitude and period, changed with dilution rate. From these data, the relationship between amplitude/ period and operating conditions can be developed. To obtain a more accurate model, we repeatedly simulated the oscillations to obtain dynamic behavior data. After analyzing more simulated results, the relationship between dilution and oscillatory phenomena characteristics, the amplitude and period time, was formulated as follows. The relationship between the product concentration oscillation amplitude and operating condition parameter, dilution rate is Din = 0.000037A mp,Cpm 3 − 0.00012A mp,Cpm 2 − 0.00028 A mp,Cpm + 0.054
with a polynomial fitting correlation coefficient of R2 = 0.9996. The relationship between oscillation period and dilution rate is T = 3.03Din−0.7
with a correlation coefficient R2 = 0.9998.
■
Figure 8. Predicted simulated dynamic behavior curve for Amp = 20.00 kg/m3.
PREDICTION To illustrate the capability of our formulated model, we predicted the amplitudes of product concentration in the oscillations generated by this Hopf singularity point, with given dilution rates. Predicted Oscillation Amplitude 10.00 kg/m3. We found that when the oscillation amplitude was 10.00 kg/m3, the dilution rate was 0.076 h−1 and the period was 18.42 h. In comparison with the simulation results, the corresponding predictions were obtained from the dynamic curve in Figure 7 with an amplitude of 9.3 kg/m3 and corresponding period of 17.8 h. These predicted results were judged to be sufficiently close to the simulated results. Predicted Oscillation Amplitude 20.00 kg/m3. Similarly, we found that when oscillation amplitude was 20.00 kg/ m3, the dilution rate was 0.297 h−1 and the period was 7.10 h. In comparison with previously simulated values, this prediction
predicted results of amplitude and period thus coincided well with the simulation results for dynamic behavior.
■
DAMPED OSCILLATORY PHENOMENA From the relationship model, it can be concluded that when the dilution rate was less than 0.054 h−1, the oscillation amplitude was close to zero, which means the oscillations tend to vanish at this dilution rate. However, from the simulated result, it can be seen when the dilution rate value ranges from 0.050 to 0.054 h−1, there would be damped oscillatory phenomena. When Din was 0.052 and 0.050 h−1, there were damped oscillations, as depicted in Figures 9 and 10. When Din was further decreased, the oscillatory phenomena rapidly vanished. Although there was a supercritical bifurcation, which is prone to generating sustained oscillations near the detected Hopf singularity point, the characteristics of the oscillatory phenomena, such as amplitude and period, can be determined
Figure 7. Predicted simulated dynamic behavior curve for Amp = 10.00 kg/m3.
Figure 9. Dynamic behavior when Din = 0.052 h−1. 12404
dx.doi.org/10.1021/ie5009234 | Ind. Eng. Chem. Res. 2014, 53, 12399−12410
Industrial & Engineering Chemistry Research
Article
It should be noted that the regulation strategy presented herein is based on parametric variations because of their ease of manipulation. Feedback control design is also possible for controlling the amplitude/period of limit cycles, but these must be implemented in closed loop systems. During the investigation reported in this paper, a number of detailed problems that require further discussion were identified, as follows: Maximum and Minimum Product Concentrations under Oscillatory Phenomena. With the given Hopf singularity point, varying the dilution rate changes the amplitude of product concentration variation. The result is shown in Figure 11. As we can see, with increasing dilution rate, the amplitude of product concentration gradually increases.
Figure 10. Dynamic behavior when Din = 0.05 h−1.
by the dilution rate. In other words, with strategic selection of certain operating conditions, the initial presence of oscillation can be gradually attenuated, or in some certain situations, be eliminated.
■
DISCUSSION In this paper, we investigated the quantitative relationship between the characteristics of oscillatory phenomena and operating conditions near Hopf singularity points. Because the singular points and limit cycles found in Z. mobiliz continuous fermentation processes are fairly common, the relationship between oscillation characteristic parameters and the operating conditions is useful for guiding the regulation of oscillations, especially for cases in which the oscillations cannot be avoided. The procedure is as follows. (1) Model the process in ODE mathematical formulations. (2) Calculate the steady state solutions, so these solutions can be selected as operating points. (3) Determine and locate any particular steady state operating point corresponding to a Hopf singularity point. (4) If there are no Hopf points or these singularity points play very little roles in the process, the dynamic oscillation behavior can be neglected, but in practice, oscillations are often generated from Hopf singularity points near optimal operating conditions. As a result, (5) the dynamic behavior must be studied, so we simulate the oscillations and identify the relationship, (6) formulate the model between oscillation amplitude and operating conditions, and find the oscillation period, then (7) use the model to predict the oscillation amplitude under certain conditions, using this result as a guide to attenuate unfavorable oscillations in the process. (8) If the predicted result does not agree with the simulation or experimental result, the model must be adjusted, by obtaining more simulated results to construct the relationship model, and collection of more data to validate the dynamic process system model may be needed in some cases. The above steps must then be repeated until we have a better predictive model to regulate oscillations. It can be seen that in this procedure, the Hopf singularity is a crucial factor that triggers dynamic oscillations, and the oscillation amplitude and frequency can be adjusted by changing the operating conditions. In this way, we can regulate oscillatory phenomena.
Figure 11. Product concentration changes with dilution rate.
The maximum and the minimum concentrations of product re shown in Figure 11. From this, we can see that (1) the amplitude increases with increasing dilution rate; (2) the maximum product concentration is confined by an upper bound, which means there exists a limitation in this fermentation process; (3) the minimum concentration decreases with increasing dilution rate; and (4) the average product concentration decreases with increasing dilution rate. As a result, oscillatory phenomena should be avoided because oscillations have a negative impact on total product quantity because of the decrease in average product concentration. Maximum and Minimum Substrate Concentrations under Oscillatory Phenomena. With the same Hopf singularity point, how would other state variables, such as substrate concentration change under these oscillatory phenomena? The result of substrate concentration amplitude changes with dilution rate is given in Figure 12. The relationship between substrate concentration oscillation amplitude and the operating condition parameter, dilution rate, was Din = 0.000031A mp,Cs 3 − 0.00038A mp,Cs 2 + 0.00024A mp,Cs + 0.054
with a regression coefficient R2 = 0.9996. The periodic time relationship could be analyzed through product oscillations because these dynamic behaviors in oscillation were varied simultaneously. That is 12405
dx.doi.org/10.1021/ie5009234 | Ind. Eng. Chem. Res. 2014, 53, 12399−12410
Industrial & Engineering Chemistry Research
Article
because of the simultaneously varying dynamic behavior of the system. Oscillations at Other Hopf Points. The above-mentioned results were obtained from a single Hopf singularity point. As mentioned in our previous research work, there exists a series of Hopf points in the operating conditions panel and every Hopf point will give birth to rich dynamic behavior in this fermentation system. Here, an attempt was made to calculate and show all the Hopf points in the practically feasible operating panel. It can be seen in Figure 14, from the operation panel constructed by dilution rates and substrate concentrations,
Figure 12. Substrate concentration changes with dilution rate.
T = 3.03Din−0.70 From Figure 12, we can see that, in these oscillatory phenomena, the amplitude of substrate concentration increased as the dilution rate increases. Compared with Figure 11, when product concentration was at a maximum, the substrate concentration was at its minimum, because at this point, more substrate was fermented to form product. On the other hand, when the substrate was at its maximum, there was less substrate fermented to form product. As a result, the production reached its minimum point. Figure 14. All Hopf points in the parameter panel.
there is a resultant curve composed of Hopf points. The Hopf point that was analyzed in detail above is marked in red. This “Hopf curve” indicates that there exists a series of Hopf points and corresponding limit cycles with quite different values of parameters (Cpm, Cs) and oscillatory phenomena. For example, we can select another Hopf point from the Hopf curve, where the dilution rate was 0.2 h−1. Other corresponding values for this point are given in Table 3 below. This point is referred to as HP2. Table 3. Values of Hopf Point, HP2
Figure 13. Microorganism concentration changes with dilution rate.
From Figure 13, the relationship between oscillation amplitude in microorganism concentration and the operating condition parameter dilution rate was found to be Din = 0.20A mp,Cx 3 − 0.041A mp,Cx 2 − 0.0070A mp,Cx
parameter
value
unit
Ce Din Cp Cs0 Cs CpM Cx
0.21 0.20 55.59 136.15 12.46 55.52 2.40
kg/m3 h−1 kg/m3 kg/m3 kg/m3 kg/m3 kg/m3
From HP2, the dynamic behaviors are illustrated in Figure 15. The relationship between period and Din can be formulated from the simulated data in the form shown below.
+ 0.054
with a regression coefficient of R2 = 0.9995. The periodical time relationship was related to product oscillation by
T = 2.84Din−0.65
T = 3.03Din−0.7
with a correlation coefficient of R2 = 0.9996. 12406
dx.doi.org/10.1021/ie5009234 | Ind. Eng. Chem. Res. 2014, 53, 12399−12410
Industrial & Engineering Chemistry Research
Article
In addition, the concentration of biomass can be simulated in these oscillations, and the relationship between amplitude of the biomass and the dilution rate is shown in Figure 17.
Figure 15. Oscillatory product concentration changes with dilution rate.
Using the same method, the relationship between the amplitude of product concentration and the dilution rate, has the following form:
Figure 17. Oscillatory biomass concentration changes with dilution rate.
Din = 0.0000077A mp,Cpm 3 + 0.0015A mp,Cpm 2 + 0.0042
A quantitative relationship in mathematical terms between amplitude of biomass concentration and the dilution rate was constructed as follows:
A mp,Cpm + 0.20
with a correlation coefficient of R2 = 0.9999. Another important substance, the concentration of substrate, was similarly influenced by dilution rate. The relationship between amplitude in substrate concentration and the dilution rate was simulated in these oscillations and the results are given in Figure 16.
Din = 0.051A mp,Cx 3 + 0.45A mp,Cx 2 + 0.078A mp,Cx + 0.20
with a correlation coefficient of R2 = 0.9999. These results were derived from another Hopf point, where the dilution rate was 0.2 h−1 and the initial substrate concentration started at 136.15 kg/m3. These models can be used as guidelines for attenuating the oscillations that arose from this Hopf singularity point. Oscillations with Other Parameter Influences. By combining other parameters, an alternative Hopf point curve can be obtained. Consequently, with Hopf points displaying larger differences, such as Hopf points lying in different curves, the relationship between the oscillation amplitude/period and operating conditions would vary dramatically. In this continuous fermentation process with Z. mobiliz, the fermentor dilution rate, Din, the membrane permeate-side dilution rate, and the initial substrate concentration of feed stock, Cs0, could be selected as bifurcation parameters, not only because of their decisive role as operating variables in this system but also because of their ease of manipulation and control in practice. We calculated Hopf points using an operational condition panel spanned by the dilution rate and the initial substrate concentration, with different degrees of membrane permeateside dilution rate of 0.1, 0.5, and 0.9 h−1. These curves collected three comprehensive sets of Hopf points, shown in Figure 18. The point HP1, studied at the beginning of this paper, is marked in red in the Hopf curve for DMin = 0.1 h−1. It is understandable that the oscillation characteristics derived from different Hopf curves should differ. For example, we selected another Hopf point from the curve calculated from DMin = 0.9 h−1 and Din = 0.6 h−1, and the values for this point are listed in Table 4 below. This point is referred to as HP3. Calculating the oscillations generated from HP3, and then analyzing these data, we obtained a mathematical model for the
Figure 16. Oscillatory substrate concentration changes with dilution rate.
The quantitative model was regressed from these oscillation data, and the equation was thus determined to be as follows Din = 0.00000078A mp,Cs 3 + 0.00027A mp,Cs 2 + 0.0018 A mp,Cs + 0.20
with a correlation coefficient of R2 = 0.9999. 12407
dx.doi.org/10.1021/ie5009234 | Ind. Eng. Chem. Res. 2014, 53, 12399−12410
Industrial & Engineering Chemistry Research
Article
Figure 18. All Hopf points in the parameter panel.
Table 4. Values of Hopf Point, HP3 parameter
value
unit
Cs0 Din Cs Cp Ce CpM Cx
160.34 0.60 35.34 51.11 0.60 50.64 2.84
kg/m3 h−1 kg/m3 kg/m3 kg/m3 kg/m3 kg/m3
relationship between the oscillatory period and dilution rate, with the result formulated as below: T = 2.82Din−0.51
with a correlation coefficient of R2 = 0.9953. At the selected dilution rate, the oscillations varied simultaneously, with the same period as the concentration variations, while the amplitudes of the different component concentrations varied to different extents. From the simulation results, concentrations of product, substrate and biomass were observed to change with dilution rate under oscillations, as shown in Figure 19. The relationship between product concentration and dilution rate is shown in Figure 19(a). The oscillation generated from the Hopf point tends to be more vigorous, with gradually increased amplitude as the dilution rate increases. After analyzing these oscillatory data, a relationship model between the product concentration amplitude and dilution rate was formulated as below: Din = 0.000025A mp,Cpm 3 − 0.00038A mp,Cpm 2 + 0.0055 Figure 19. Oscillatory changes in product, substrate and biomass concentrations with dilution rate.
A mp,Cpm + 0.60
with a correlation coefficient of R2 = 0.9998. Applying the same method, using simulations of oscillatory phenomena generated from the Hopf singularity points, results for the substrate concentration are shown in Figure 19(b). The quantitative relation between amplitude of the substrate and dilution rate was then constructed. The result was formulated as follows:
Din = 0.0000020A mp,Cs 3 − 0.000084A mp,Cs 2 + 0.0026 A mp,Cs + 0.60
with a correlation coefficient of R2 = 0.9998. 12408
dx.doi.org/10.1021/ie5009234 | Ind. Eng. Chem. Res. 2014, 53, 12399−12410
Industrial & Engineering Chemistry Research
Article
Ce0
Finally, the relationship between the concentration of biomass and the dilution rate was also studied, with the result shown in Figure 19(c). The mathematical model for the relationship between amplitude in biomass concentration and dilution rate had the following form:
Cp0 CpM0
Din = 0.13A mp,Cx 3 − 0.13A mp,Cx 2 + 0.10A mp,Cx + 0.60
Cs0
2
with a correlation coefficient of R = 0.9998. Having analyzed the above correlations, we can see that quantitative models relating the characteristics of oscillatory phenomena and operating conditions can be formulated. From the points of view of process and product design, quantitative models such as those presented here are a useful indicator of how and to what extent we can attenuate oscillation behaviors in cases where oscillations literally cannot be avoided.
Cx0 Ce Cp CpM
■
Cs Cx Din Dout DMin
CONCLUSIONS In this paper, we investigated the oscillatory phenomena in a continuous fermentation process with Z. mobiliz. In these oscillations, the concentrations of product, biomass, and substrate sustainable varied periodically. It was reported that Hopf singularity points in the process caused these oscillations but previous research work did not give a quantitative method for regulating these oscillations. In this paper, we studied the process model and calculated the steady state solution, located the Hopf singularity points, and simulated oscillations around Hopf points. A general explicit formula for the relationship between concentration amplitudes/period and operating conditions was derived as a regulation approach. As a result, oscillations could be quantitatively attenuated, as expected, where they cannot be avoided. Simulation results were presented to confirm the analytical predictions. Also the result can be validated with experimental approach, in which the exact state variable values and initial operating condition values are required to perform the oscillatory phenomenon.
■
DMout Pm T VM VF Ypx Ysx k1 k2 k3 mp ms re rp rs rx ρ μ
ASSOCIATED CONTENT
S Supporting Information *
The model parameters used in the current investigation given in Table S1.27 This material is available free of charge via the Internet at http://pubs.acs.org.
■
■
AUTHOR INFORMATION
concentration of key component in fermentor feed, kg/m3 concentration of product in fermentor feed, kg/ m3 concentration of product in membrane permeate side feed, kg/m3 concentration of product in fermentor feed, kg/ m3 concentration of biomass in fermentor feed, kg/m3 concentration of key component in fermentor, kg/m3 concentration of product in fermentor, kg/m3 concentration of product in membrane permeate side, kg/m3 concentration of substrate in fermentor, kg/m3 concentration of biomass in fermentor, kg/m3 inlet dilution rate of fermentor, h−1 outlet dilution rate of fermentor, h−1 inlet dilution rate in membrane permeate side, h−1 outlet dilution rate in membrane permeate side, h−1 permeability of the membrane for ethanol, m/h period of oscillation, h volume of membrane permeate side, m3 volume of fermentor, m3 yield factors of biomass on product, kg/kg yield factors of biomass on substrate, kg/kg empirical constants, h−1 empirical constants, m3/kgh empirical constants, m6/kg2h maintenance factors based on product formation, kg/kgh maintenance factors based on substrate requirements, kg/kgh rate of the key compound, kg/m3h product formation rate, kg/m3h substrate consumption rate, kg/m3h growth rate of biomass, kg/m3h density of mixture, kg/m3 specific growth rate, kg/m3h
REFERENCES
(1) Abashar, M. E. E.; Elnashaie, S. S. E. H. Multistablity, bistability and bubbles phenomena in a periodically forced ethanol fermentor. Chem. Eng. Sci. 2011, 66 (23), 6146−6158. (2) Shen, Y.; Zhao, X. Q.; Ge, X. M.; Bai, F. W. Metabolic flux and cell cycle analysis indicating new mechanism underlying process oscillation in continuous ethanol fermentation with Saccharomyces cerevisiae under VHG conditions. Biotechnol. Adv. 2009, 27 (6), 1118−1123. (3) Paz Astudillo, I. C.; Cardona Alzate, C. A. Importance of stability study of continuous systems for ethanol production. J. Biotechnol. 2011, 151 (1), 43−55. (4) Bai, F. W.; Anderson, W. A.; Moo-Young, M. Ethanol fermentation technologies from sugar and starch feedstocks. Biotechnol. Adv. 2008, 26 (1), 89−105. (5) Sprenger, G. A. Carbohydrate metabolism in Zymomonas mobilis: a catabolic highway with some scenic routes. FEMS Microbiol. Lett. 1996, 145 (3), 301−307. (6) Wang, L.; Zhao, X.-Q.; Xue, C.; Bai, F.-W. Impact of osmotic stress and ethanol inhibition in yeast cells on process oscillation associated with continuous very-high-gravity ethanol fermentation. Biotechnol. Biofuels 2013, 6 (1), 133.
Corresponding Author
*Tel.: +86 10 62784572. Fax: +86 10 62770304. E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS The authors gratefully acknowledge financial support from the NSFC (No. 21306100) and the National Basic Research Program of China (No. 2012CB720500). Nomenclature AM permeation area, m2 Amp,Cpm amplitude of product concentration in oscillation, kg/m3 Amp,Cs amplitude of substrate concentration in oscillation, kg/m3 Amp,Cx amplitude of biomass concentration in oscillation, kg/m3 12409
dx.doi.org/10.1021/ie5009234 | Ind. Eng. Chem. Res. 2014, 53, 12399−12410
Industrial & Engineering Chemistry Research
Article
(7) Lee, K. J.; Tribe, D.; Rogers, P. Ethanol production byZymomonas mobilis in continuous culture at high glucose concentrations. Biotechnol. Lett. 1979, 1 (10), 421−426. (8) Lee, K.; Skotnicki, M.; Tribe, D.; Rogers, P. Kinetic studies on a highly productive strain of Zymomonas mobilis. Biotechnol. Lett. 1980, 2 (8), 339−344. (9) Jöbses, I.; Egberts, G.; Luyben, K.; Roels, J. Fermentation kinetics of Zymomonas mobilis at high ethanol concentrations: oscillations in continuous cultures. Biotechnol. Bioeng. 1986, 28 (6), 868−877. (10) Ghommidh, C.; Vaija, J.; Bolarinwa, S.; Navarro, J. M. Oscillatory behaviour of Zymomonas in continuous cultures: A simple stochastic model. Biotechnol. Lett. 1989, 11 (9), 659−664. (11) Bruce, L.; Axford, D.; Ciszek, B.; Daugulis, A. Extractive fermentation by Zymomonas mobilis and the control of oscillatory behavior. Biotechnol. Lett. 1991, 13 (4), 291−296. (12) Jarzebski, A. B. Modeling of oscillatory behavior in continuous ethanol fermentation. Biotechnol. Lett. 1992, 14 (2), 137−142. (13) Hobley, T.; Pamment, N. Differences in response of Zymomonas mobilis and Saccharomyces cerevisiae to change in extracellular ethanol concentration. Biotechnol. Bioeng. 1994, 43 (2), 155−158. (14) Li, J.; McLellan, P. J.; Daugulis, A. J. Inhibition effects of ethanol concentration history and ethanol concentration change rate onZymomonas mobilis. Biotechnol. Lett. 1995, 17 (3), 321−326. (15) Daugulis, A. J.; McLellan, P. J.; Li, J. H. Experimental investigation and modeling of oscillatory behavior in the continuous culture of Zymomonas mobilis. Biotechnol. Bioeng. 1997, 56 (1), 99− 105. (16) McLellan, P. J.; Daugulis, A. J.; Li, J. H. The incidence of oscillatory behavior in the continuous fermentation of Zymomonas mobilis. Biotechnol. Prog. 1999, 15 (4), 667−680. (17) Garhyan, P.; Elnashaie, S. S. E. H.; Al-Haddad, S. M.; Ibrahim, G.; Elshishini, S. S. Exploration and exploitation of bifurcation/chaotic behavior of a continuous fermentor for the production of ethanol. Chem. Eng. Sci. 2003, 58 (8), 1479−1496. (18) Garhyan, P.; Elnashaie, S. Static/dynamic bifurcation and chaotic behavior of an ethanol fermentor. Ind. Eng. Chem. Res. 2004, 43 (5), 1260−1273. (19) Garhyan, P.; Elnashaie, S. S. E. H. Experimental Investigation and Confirmation of Static/Dynamic Bifurcation Behavior in a Continuous Ethanol Fermentor. Practical Relevance of Bifurcation and the Contribution of Harmon Ray. Ind. Eng. Chem. Res. 2004, 44 (8), 2525−2531. (20) Bai, F. W.; Chen, L. J.; Anderson, W. A.; Moo-Young, M. Parameter oscillations in a very high gravity medium continuous ethanol fermentation and their attenuation on a multistage packed column bioreactor system. Biotechnol. Bioeng. 2004, 88 (5), 558−566. (21) Bai, F. W.; Ge, X. M.; Anderson, W. A.; Moo-Young, M. Parameter Oscillation Attenuation and Mechanism Exploration for Continuous VHG Ethanol Fermentation. Biotechnol. Bioeng. 2009, 102 (1), 113−121. (22) Diehl, F. C.; Trierweiler, J. O., Control Strategy for a Zymomonas mobilis Bioreactor Used in Ethanol Production. In Computer-Aided Chemical Engineering, de Brito Alves, R. M. d. N., Evaristo Chalbaud, B. C. A. O., Eds.; Elsevier: Amsterdam, 2009; Vol. 27, pp 1605−1610. (23) Li, C. C. Mathematical models of ethanol inhibition effects during alcohol fermentation. Nonlinear Anal. 2009, 71 (12), 1608− 1619. (24) Sridhar, L. N. Elimination of Oscillations in Fermentation Processes. AIChE J. 2011, 57 (9), 2397−2405. (25) Wang, H. Z.; Zhang, N.; Qiu, T.; Zhao, J. S.; He, X. R.; Chen, B. Z. Analysis of Hopf Points for a Zymomonas mobilis Continuous Fermentation Process Producing Ethanol. Ind. Eng. Chem. Res. 2012, 52 (4), 1645−1655. (26) Bai, F.; Ge, X.; Anderson, W.; Moo-Young, M. Parameter oscillation attenuation and mechanism exploration for continuous VHG ethanol fermentation. Biotechnol. Bioeng. 2009, 102 (1), 113− 121.
(27) Mahecha-Botero, A.; Garhyan, P.; Elnashaie, S. Non-linear characteristics of a membrane fermentor for ethanol production and their implications. Nonlinear Anal.: Real World Appl. 2006, 7 (3), 432− 457. (28) Doedel, E. J. AUTO: A program for the automatic bifurcation analysis of autonomous systems. Congr. Numer 1981, 30, 265−284. (29) Doedel, E.; Kernevez, J. P., AUTO: Software for Continuation and Bifurcation Problems in Ordinary Differential Equations: Including the AUTO 86 User Manual. In Applied Mathematics, California Institute of Technology: 1986. (30) Doedel, E.; Champneys, A.; Fairgrieve, T.; Kuznetsov, Y.; Sandstede, B.; Wang, X., AUTO 97: Continuation and Bifurcation Software for Ordinary Differential Equations, User’s Manual. Center for Research on Parallel Computing, California Institute of Technology, Pasadena 1997. (31) Kuznetsov, Y.; Levitin, V., CONTENT: dynamical system software. http://www.computeralgebra.nl/systemsoverview/special/ diffeqns/content/gcbody.html 1997. (32) Dhooge, A.; Govaerts, W.; Kuznetsov, Y. A. MATCONT: a MATLAB package for numerical bifurcation analysis of ODEs. ACM Trans. Math. Software (TOMS) 2003, 29 (2), 141−164. (33) Dhooge, A.; Govaerts, W.; Kuznetsov, Y. A. MATCONT: a Matlab package for numerical bifurcation analysis of ODEs. ACM SIGSAM Bull. 2004, 38 (1), 21−22. (34) Back, A.; Guckenheimer, J.; Myers, M.; Wicklin, F.; Worfolk, P. DsTool: Computer assisted exploration of dynamical systems. Notices Am. Math. Soc 1992, 39 (4), 303−309. (35) Clewley, R.; Sherwood, W.; LaMar, M.; Guckenheimer, J., PyDSTool, a software environment for dynamical systems modeling. http://pydstool.sourceforge.net 2007. (36) Ermentrout, B., Simulating, analyzing, and animating dynamical systems: a guide to XPPAUT for researchers and students. Society for Industrial Mathematics: 2002; Vol. 14. (37) Wang, H. Z.; Chen, B. Z.; He, X. R.; Zhao, J. S.; Qiu, T. Numerical Analysis Tool for Obtaining Steady-State Solutions and Analyzing Their Stability Characteristics for Nonlinear Dynamic Systems. J. Chem. Eng. Jpn. 2010, 43 (4), 394−400. (38) Yuan, Z. H.; Wang, H. Z.; Chen, B. Z.; Zhao, J. S. Operating zone segregation of chemical reaction systems based on stability and non-minimum phase behavior analysis. Chem. Eng. J. 2009, 155 (1−2), 304−311. (39) Wang, H. Z.; Chen, B. Z.; He, X. R.; Qiu, T.; Zhao, J. S. Singularity Theory Based Stability Analysis of Reacting Systems. Comput.-Aided Chem. Eng. 2009, 27, 645−650. (40) Wang, H. Z.; Chen, B. Z.; He, X. R.; Zhao, J. S.; Qiu, T. Modeling, simulation and analysis of the liquid-phase catalytic oxidation of toluene. Chem. Eng. J. 2010, 158 (2), 220−224. (41) Wang, H. Z.; Yuan, Z. H.; Chen, B. Z.; He, X. R.; Zhao, J. S.; Qiu, T. Analysis of the stability and controllability of chemical processes. Comput. Chem. Eng. 2011, 35 (6), 1101−1109. (42) Wang, H. Z.; Chen, B. Z.; Qiu, T.; He, X. R.; Zhao, J. S., An integrated quantitative index of stable steady state points in chemical process design. In 11th International Symposium on Process Systems Engineering-PSE2012, Singapore, 2012. (43) Wang, H. Z.; Chen, B. Z.; Qiu, T.; He, X. R.; Zhao, J. S., An Approach Considering Both Operation Stability and System’s Hopf Bifurcations to Chemical Process Design. In 2012 AIChE Annual Meeting, Pittsburg, PA., 2012. (44) Wang, H. Z; Zhang, N.; Qiu, T.; Zhao, J. S.; He, X. R.; Chen, B. Z. A process design framework for considering the stability of steady state operating points and Hopf singularity points in chemical processes. Chem. Eng. Sci. 2013, 99, 252−264. (45) Wang, H. Z.; Zhang, N.; Qiu, T.; Zhao, J. S.; He, X. R.; Chen, B. Z. Optimization of a continuous fermentation process producing 1,3propane diol with Hopf singularity and unstable operating points as constraints. Chem. Eng. Sci. 2014, 116, 668−681. (46) Abashar, M.; Elnashaie, S. Dynamic and chaotic behavior of periodically forced fermentors for bioethanol production. Chem. Eng. Sci. 2010, 65 (16), 4894−4905. 12410
dx.doi.org/10.1021/ie5009234 | Ind. Eng. Chem. Res. 2014, 53, 12399−12410
