Mechanistic Modeling of Biodiesel Production Using a Liquid Lipase Formulation Jason Price, Bj€orn Hofmann, Vanessa T. L. Silva, Mathias Nordblad, John M. Woodley, and Jakob K. Huusom Dept. of Chemical and Biochemical Engineering, Technical University of Denmark, DK-2800 Kgs. Lyngby, Denmark DOI 10.1002/btpr.1985 Published online 00 Month 2014 in Wiley Online Library (wileyonlinelibrary.com)
In this article, a kinetic model for the enzymatic transesterification of rapeseed oil with methanol using CalleraTM Trans L (a liquid formulation of a modified Thermomyces lanuginosus lipase) was developed from first principles. We base the model formulation on a PingPong Bi-Bi mechanism. Methanol inhibition, along with the interfacial and bulk concentrations of the enzyme was also modeled. The model was developed to describe the effect of different oil compositions, as well as different water, enzyme, and methanol concentrations, which are relevant conditions needed for process evaluation, with respect to the industrial production of biodiesel. The developed kinetic model, coupled with a mass balance of the system, was fitted to and validated on experimental results for the fed-batch transesterification of rapeseed oil. The confidence intervals of the parameter estimates, along with the identifiability of the model parameters were presented. The predictive capability of the model was tested for a case using 0.5% (wt. Enzyme/wt. Oil), 0.5% (wt. Water /wt. Oil) and feeding 1.5 times the stoichiometric amount of methanol in total over 24 h. For this case, an optimized methanol feeding profile that constrains the amount of methanol in the reactor was computed and the predictions experimentally validated. Monte-Carlo simulations were then used to characterize the effect of the parameter uncertainty on the model outputs, giving a C 2014 American Institute biodiesel yield, based on the mass of oil, of 90.8 6 0.55 mass %. V of Chemical Engineers Biotechnol. Prog., 000:000–000, 2014 Keywords: enzymatic biodiesel, kinetic modeling, transesterification, uncertainty analysis
Introduction For enzymatic biodiesel production, the conventional biocatalyst is usually immobilized in order to improve enzyme recovery, which combined with increased stability, allows for reuse.1 However, the immobilization carrier, as well as the immobilization process significantly contributes to the price of the biocatalyst, which further necessitate the reuse of the enzymes for the process to be competitive.1–3 Within the last few years, the use of liquid lipase formulations for enzymatic biodiesel production has resulted in a significant reduction in the biocatalyst cost.4–6 Furthermore, compared to the conventional alkali catalysts used to produce biodiesel on an industrial scale, the use of an enzymatic catalyst has the advantage that low-quality feedstock and waste oils that have a high free fatty acid (FFA) content can be treated. This is due to the fact that lipases are able to esterify the FFA contained in waste oils to esters as well as transesterify the acyl-glycerides in the oil.7 This further reduces operating costs of the enzyme-catalysed biodiesel process. Nevertheless, when developing an industrial enzymatic biodiesel process a few issues related to the biocatalyst need to be addressed; methanol inhibition, deactivation at high methanol concentrations, the limited lifespan of the lipase and the relatively high cost of the enzyme per kilogram of oil treated
Correspondence concerning this article should be addressed to J. K. Huusom at [email protected]. C 2014 American Institute of Chemical Engineers V
compared to traditional chemical catalysts (e.g., sodium hydroxide).8–10 All of the issues outlined can be mitigated by operating the enzyme-catalysed biodiesel process in an optimal manner. In terms of developing and optimizing the enzymecatalyzed biodiesel process efficiently, process modeling is a valuable tool to help focus the experimental work needed for process understanding and to support further process development.11 For example, a mechanistic process model of the system can be used to simulate the process performance over a wide range of conditions and to optimize how the process is operated. Hence, a process model can further be used to develop control strategies to mitigate enzyme inactivation and improve enzyme stability. Likewise, modelling can help in the formulation of innovative process designs and configurations.12 Integral to the development of a mechanistic model of any process, is the availability of reliable kinetic models. Over the years, various kinetic models for the enzymatic transesterification of vegetable oils have been proposed.2,13–18 However, only the model by Lv et al. deal with model development for a liquid lipase formulation.2 Likewise, in terms of engineering design purposes such as reactor configuration, optimization, and control, the gap we noticed in the literature for kinetic models for enzymatic biodiesel production is that there was a lack of emphasis placed on: 1. Use of the proposed kinetic models in the combined act of process optimization and experimental validation of the process optimization. This would enable the modeler to 1
2
Biotechnol. Prog., 2014, Vol. 00, No. 00
ascertain how well the model could be extended outside the model fitting. 2. Statistical analysis of the working bounds of the model. This remains a weak point in the credibility of these kinetic models and hence their applicability for engineering design purposes. In this contribution we address the two aforementioned points so that that we can perform model-based engineering design. In using the model for engineering design, it is desired that the kinetic model describing enzymatic biodiesel production, can predict the concentration of all the major species in the reaction. It is also essential to be able to characterize how the process responds to changes in the process conditions over the entire course of the reaction for changes in: 1. Alcohol/oil molar ratios: This is important given the need to balance the amount of methanol needed to shift the final equilibrium conversion of fatty acid methyl esters (FAME, biodiesel) while minimizing the effects of methanol inhibition and deactivation.1,15,19 2. Concentration of reactants in the reactor: The oil composition from batch to batch can vary. The oil composition then needs to be characterized so as to ascertain when the reaction has reached within specification. 3. Enzyme loading: The amount of biocatalyst added will affect reaction time, enzyme efficiency, and potential reuse of the enzyme. 4. Area of the oil–water interface: The liquid lipase CalleraTM Trans L (a liquid formulation of a modified Thermomyces lanuginosus lipase) being used is interfacially activated and hence this phenomena needs to be taken into account when modeling the system.20,21 None of the aforementioned models take into consideration all of the process conditions outlined for a liquid lipase formulation (As mentioned previously, the model by Lv et al. deal with model development for a liquid lipase formulation.2 However, their model formulation does not enable the evaluation of changes in enzyme concentration and the area of the oil-water interface). Hence in this work we: Develop and validate a dynamic mechanistic model from first principles for the transesterification of rapeseed oil with methanol using CalleraTM Trans L, which takes into consideration the effects of the process conditions outlined. Test and experimentally validate the predictive capabilities of the model by optimizing the methanol feed profile. Evaluate the model structure for identifiability of the kinetic parameters (Identifiability Analysis-Correlation matrix and Collinearity index) and characterize the uncertainty in the model outputs due to the uncertainty in the parameter estimates (Uncertainty Analysis-Monte-Carlo Simulations). The article is organized as follows: The proposed reaction mechanism is presented, followed by the experimental and numerical methods used. Subsequently, the results from parameter estimation and identifiability analysis are discussed. Finally, the model is used to predict an optimal methanol feeding profile and the uncertainty analysis is used to characterize the uncertainty in the model outputs.
Reaction Mechanism Model formulation The mathematical model describing the transesterification reaction in the biphasic oil–water system with a liquid lipase,
Figure 1.
Diagrammatic representation of the enzyme at the oil– water interface. The polar phase contains water, methanol, glycerol, and the free enzyme (Ebulk). The nonpolar phase contains the oil components along with the biodiesel formed. At the interface is the penetrated enzyme (E) and the acyl enzyme complex (EX).
CalleraTM Trans L was formulated on the basis of the following assumptions: 1. The reaction proceeds via a Ping-Pong Bi-Bi mechanism.13,15,18,22 2. Alcohol inhibition is competitive.15,18 3. Deactivation due to the alcohol could be ignored at low methanol concentrations.15,18 4. The interfacial and bulk concentrations of the substrate and products are the same (mass transfer from the bulk to the interface is instantaneous).15,18 5. Acyl migration is neglected.18 6. All reaction steps are reversible.18,23 7. The reactor mixture is homogenous (ideal mixing) and the density of the mixture is constant.18 Of the various kinetic models reported,2,13–18 the mathematical formulation of the kinetic model presented by Fedosov and coworkers fulfil the first three of the process conditions outlined in the introduction section. Since their model was developed for an immobilized enzyme, it does not describe the behavior of the liquid lipase at the oil-water interface. Hence, we extend their work by modeling the oil– water interfacial area and present a fully mechanistic formulation for the transesterification reaction using a liquid lipase. In their formulation, various pseudo-components were introduced to imitate the phase boundaries of the system, for the transesterification of rapeseed oil using Novozym 435 (immobilized Candida Antarctica lipase B).18 Key to how we describe the reaction mechanistically for the liquid lipase, is the interaction of the enzyme at the oil–water interface. A schematic illustration of the oil water interface along with the enzymes and its complexes is presented in Figure 1. The triglycerides (T), diglycerides (D), monoglycerides (M), biodiesel (BD), and FFA occupy the nonpolar phase while the bulk enzyme (Ebulk), Water (W), Methanol (CH), and glycerol (G) occupy the polar phase. The lipase used in this study exhibits a pronounced interfacial activation and the reaction is assumed to proceed exclusively at the interface.24 By including the interfacial enzyme concentration (E), the reaction scheme proceeds as shown in Table 1 and Figure 2. The specific interfacial area per unit volume of the oil water mixture (aT [m2/m3]) can be represented as:
Biotechnol. Prog., 2014, Vol. 00, No. 00
3
Table 1. Kinetic Mechanism for the Enzymatic Transesterification i
Reactions
1 2 3 4 5 6 7
Ebulk 1 Af $ E T 1 E $ E.T E.T $ EX 1 D D 1 E $ E.D E.D $ EX 1 M M 1 E $ E.M E.M $ EX 1 G
8 9 10
EX 1 W $ FA 1 E EX 1 CH $ BD 1 E CH 1 E $ E.CH
Figure 2.
Rate of Reaction (ri)
Enzyme in bulk absorbed at the interface In reactions 2, 4, and 6 the penetrated enzyme can react with the substrate to form an enzyme substrate complex E.T, E.D, or E.M (Ping) In reactions 3, 5, and 7 the enzyme substrate complex forms the Acyl enzyme complex and releases the first product D, M or G (Pong) Intermediate steps for the reactions were grouped together given interest is in the overall rate. The acyl enzyme complex can then react with water or methanol (Pong) and then release the second product FA or BD (Ping). Reversible competitive methanol inhibition.
k1 ½Ebulk ½Af 2k21 ½E k2 ½T ½E2k22 ½ET k3 ½ET2k23 ½EX ½D k4 ½D ½E2k24 ½ED k5 ½ED2k25 ½EX ½M k6 ½M ½E2k26 ½EM k7 ½EM2k27 ½EX ½G k8 ½EX ½W2k28 ½FA ½E k9 ½EX ½CH2k29 ½BD ½E k10 ½CH ½E2k210 ½ECH
Conceptual scheme for the overall reaction mechanism describing the enzymatic transesterification.
aT 5
6 Vp : ds V
(1)
where ds is the Sauter mean diameter of the droplets in the system [m], Vp is the polar volume [m3] and V is the bulk volume [m3].25,26 An enzyme coverage Ae, of 2.98 3 107 m2/mol, was used in this study, which denotes the interfacial area that is required for the adsorption of 1 mol enzyme to the interface.20 It is assumed that all forms of the adsorbed enzyme molecules, meaning both the free enzyme as well as all
forms of the enzyme complexes and the inhibition complex, occupy the same area at the interface.20 Given Ae, it is possible to calculate the free specific interfacial area, af [m2/m3] as shown in Eq. 2. af 5aT 2Ae :ðE1EX1ET1ED1EM1ECHÞ
(2)
The free specific interfacial area can then be expressed as a volumetric concentration (Af [mol/m3]), by using the enzyme coverage to estimate a theoretical upper limit of the moles of enzyme molecules that can occupy the free interfacial area.
4
Biotechnol. Prog., 2014, Vol. 00, No. 00 Af 5af =Ae
(3)
In the model formulation, all concentrations used are lumped concentrations (concentration of the component in the entire reaction volume). The mass balance for the system can then be combined with the kinetics to give the system of ordinary differential equations presented below. dð½TVÞ
. dt
dð½DVÞ
52Vðr2 Þ
.
dt
dð½MVÞ
dð½FAVÞ
dð½GVÞ
.
dt
.
dt
.
dt
dð½WVÞ
.
5Vðr1 1r8 1r9 2r2 2r4 2r6 2r10 Þ
dt
dð½EXVÞ
52Vðr9 1r10 Þ
dt
. .
.
5Vðr2 2r3 Þ
dt
dð½E:DVÞ
.
dt
dð½E:MVÞ
5Vðr5 2r6 Þ
.
dt
dð½E:CHVÞ
dð½Ebulk VÞ
dðVp Þ
dðVÞ
.
dt
dt
dt
dt
.
5Vðr6 2r7 Þ
.
.
(4)
5Vðr3 1r5 1r7 2r8 2r9 Þ
dt
dð½E:TVÞ
Fed-batch experiments
5Vðr8 Þ
52Vðr8 Þ
dt
dð½EVÞ
5Vðr9 Þ
5Vðr7 Þ
.
dð½CHVÞ
CalleraTM Trans L with a hydrolytic activity of 1 3 105 LU/g was kindly donated by Novozymes A/S (Bagsværd, Denmark). One LU is defined as the activity required to produce 1 lmol butyric acid in the hydrolysis of tributyrin under standard conditions (pH 7.5, 0.2M substrate).24
5Vðr5 2r6 Þ
dt
dð½BDVÞ
Biocatalyst
5Vðr3 2r4 Þ
.
5Vðr10 Þ
52Vðr1 Þ
For the 13 experiments (Table 3), the water and enzyme content were varied from 3 to 7 and 0.1 to 0.5 wt % oil, respectively. In all the experiments 1.5 equivalents (Eq.) of methanol was reacted with the rapeseed oil. One equivalents corresponds to the stoichiometric amount of alcohol needed to convert all fatty acid residues in the oil to biodiesel (i.e., 1 mol oil:3 mol alcohol). The reaction was carried out in a 0.25 L glass reactor with a tank diameter of 55 mm (T) and 2 baffles, each 0.18 3 T wide. The reactor was immersed in a water bath with temperature control (Julabo Labor-technik GmbH, Seelbach, Germany) maintained at 35 C. A rushton turbine (impeller diameter 0.44 T), spinning at 1400 rpm provided the mixing. Initially 0.2 Eq methanol was charged with the oil in the reactor. When the reaction mixture reached the reaction temperature, the amount of water and enzyme to be used in the experiment (Table 2), was then added to the reactor and methanol feeding started. Methanol feeding was provided by a KNF STEPDOS .03 pump (KNF Neuberger AB, Stockholm, Sweden), calibrated prior to each experiment. Sample preparation
5RG 1RW
5ðFa Þ
The measurement vector ym,i is then shown in Eq. 5 where ym is the measurement matrix [mass %], x are the state variables [mol/L], V is the bulk volume and rmm is the relative molecular mass of component i.
xi V rmmi ym;i 5 X5 3100; x V rmmi i51 i
(Brïndby, Denmark). Absolute methanol (99.8%, technical grade) was purchased from VWR Bie & Berntsen A/S (Herlev, Denmark). n-Heptane (99%), acetic acid (99%), isopropanol (99%), and tert-butyl methyl ether (99.8%) for HPLCAnalysis were obtained from Sigma-Aldrich A/S (Brïndby, Denmark).
0 TAG DAG ym 5 MAG FAME FFA
0 T D and x5 M BD FA (5)
Experimental Materials and Methods Chemicals Rapeseed oil was obtained from Emmelev A/S (Otterup, Denmark) and oleic acid was purchased from Sigma-Aldrich
Fifty-microliter samples were taken from the reactor and mixed with 500 lL solvent A (acetic acid and n-heptane 4:1000 v/v – mobile phase). Samples were then centrifuged at 14,500 rpm for 5 min and 10 lL of the supernatant was mixed with 990 lL of solvent A prior to the HPLC analysis. HPLC analysis Forty-microliters of the prepared sample was injected in the HPLC (Ultimate 3000, Dionex A/S, Hvidovre, Denmark) for analysis of triglycerides (TAGs), diglycerides (DAGs), monoglycerides (MAGs), FFAs and fatty acid methyl esters (FAME). The separation of the different compounds was carried out with a cyanopropyl column (0.25 3 0.004 m) (DisR , Cyano, Sigma Aldrich A/S, Brïndby, Denmark), covery V U3000 auto-sampler, TCC-3000SD column oven, U3400A R Charged Aerosol quaternary pump modules and a CoronaV Detector (Thermo Scientific Dionex, Chelmsford, MA). A binary gradient program was employed for the separation of the different compounds using Solvent A, Solvent B (99.6% v/v tert-butyl methyl ether and 0.4% v/v acetic acid) and iso-propanol as Solvent C.27,28 The detection of the different compounds after separation with the column was carried out R Charged Aerosol Detector from Thermo Sciby a CoronaV entific Dionex (Chelmsford, MA) with nitrogen gas at a
Biotechnol. Prog., 2014, Vol. 00, No. 00
5
Table 2. Experiments for the Data Fitting, Validation, and Optimization
Parameter Fitting
Validation
Optimization
Exp.
Methanol Feed Rate [Eq./h]
Initial Dose Methanol [Eq]
Water [wt.% oil ]
Enzyme [wt.% oil ]
1 2 3 4 5 6 7 8 9 10 11 12 13
0.06 0.06 0.06 0.06 0.06 0.1 2hrs. 0.06 2hrs. 0.06 2hrs. 0.06 0 0 0 3hrs. 0.02
0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.4 0.4 0.4 0.525
3 3 3 5 5 5 5 5 5 5 7 7 5
0.1 0.2 0.3 0.2 0.5 0.3 0.2 0.5 0.3 0.5 0.2 0.5 0.5
0.185 first 0.185 first 0.185 first
0.152 first
thereafter thereafter thereafter
thereafter
pressure of 241 kPa. The composition of the reaction samples was reported on a mass percentage basis, relative to the sum of quantified mass of the five analyzed components (TAG, DAG, MAG, FFA, and FAME). From previous experiments the standard deviation of the measurements, for TAG, DAG, MAG, FFA, and FAME were found to be 0.40, 0.75, 0.18, 0.28, and 0.26 mass %, respectively. Given that the oil composition can vary from different batches of oil, the initial concentration of the oil at the start of each experiment is measured.
Numerical Methods Simulation environment R The model was implemented and simulated in MatlabV (The Mathworks, Natick, MA). All the methods for performing the identifiability, sensitivity, and uncertainty analysis were built on the toolbox based on the work by Sin et al.29 The following sections give further details of the methods used.
Parameter estimation, confidence intervals, and identifiability analysis Parameter Estimation and Confidence Intervals. The 20 unknown kinetic constants (k12k10, k212k210), were estimated by fitting the model equations with full time course data using experiments 1–8, which covered the span of our operating conditions (Table 3). To judge the quality of the fitting, the validation data sets were chosen so that the operating conditions fall within the parameter fitting dataset (Experiments 9) and also to evaluate how the model predicts the initial rates of the components (Experiments 10–12). The differential equations were solved using a stiff variable order solver based on numerical differentiation formulas (ode15s). For the parameter fitting, the squared-sum of the relative errors between the simulated and experimental values for TAG, DAG, MAG, FAME, and FFA were minimized using fminsearch (based on a simplex search algorithm).30 To quickly assess the quality of the data fitting, the histogram of residuals were used to examine the underlying statistical assumptions of the residuals having zero mean and being normally distribution. Scott’s method, was used to determine the number of bins and is based upon the sample standard deviation and number of data points.31 A robust method to calculate realistic confidence intervals for the parameter estimates h^ was made by using a bootstrap method.32,33 To perform the bootstrap analysis, the residuals from the parameter estimates were used to generate synthetic data, which was subsequently used in a Monte–Carlo method
10,000 times to generate simulated data (Increasing the sample size from 5,000 to 10,000 did not change parameter distribution significantly). The simulated data was then used to generate a set of re-estimated parameters (matrix size [10,000 3 20]), where the mean of the distribution was used as the mean parameter estimate. The 95th and 5th percentiles of the re-estimated parameters were then used as the upper and lower bounds of the confidence intervals for the parameters estimates, respectively.34 Identifiability Analysis. The identifiability method based on the work by Brun et al. was used to ascertain, which parameters could actually be identified from the available experimental data given the model structure.35 The method follows three main steps. (1) Calculation of the sensitivity matrix, (2) scaling of the sensitivity matrix, and (3) calculation of the collinearity index for the subset of parameters. Step 1: Calculate sensitivity matrix S: The sensitivities of the model outputs (T, D, M, BD, and FA) to the parameters were calculated by the direct differential method.36 S has dimensions N 3 m (N is the number of experimental data points and m is the number of parameters). The sensitivities of the model outputs were placed in the sensitivity matrix, S: @xi S5 @hj h^
(6)
where i is the total number of observations for the five model outputs for the 8 experiments, and j 5 1:m for each of the parameters. Step 2: Scaling of the sensitivity matrix, S: The nondimensional sensitivity matrix si,j and the normalised sensitivity matrix S^i;j were then computed: Si;j 5Si;j
Dhj si;j and S^i;j 5 sci ksj k
(7)
Where the mean estimate of hj is used for Dhj, the mean of the experimental observations for each model output (i 5 1:5) is used for sci and kSi;j k is the Euclidean norm of the jth column of si,j. Step 3: Calculation of the collinearity index for the different subset of parameters: The collinearity index, ck for the subset of parameters k (k 5 2:M) is then: 1 ck 5 pffiffiffiffiffi kk
(8)
where kk is the smallest eigenvalue of ^s kT^s k with ^s k being a N 3 k submatrix of S^i;j whose columns correspond to the
FAME
0.00 0.90 3.09 6.19 8.32 14.61 18.79 22.66 68.00 70.00 71.00 71.50
0.00 10.75 16.97 22.12 26.00 30.50 33.77 37.36 73.44 75.50 76.92
0.00 7.73 16.85 24.45 29.85 36.58 39.35 44.54 49.00 52.20 77.11 79.30 81.10
Time(min)
0 30 60 90 120 180 240 300 1260 1380 1440 1500
0 30 60 120 180 240 300 360 1320 1380 1440
0 30 60 90 120 150 180 240 300 360 1200 1320 1440
Exp. 1
Exp. 4
Exp. 7
5.75 10.84 12.48 12.16 11.11 9.79 8.79 7.60 6.31 5.48 2.91 2.99 2.96
8.10 14.97 16.80 16.50 15.00 14.00 12.50 11.00 4.51 4.40 4.32
5.73 9.86 11.74 14.44 15.69 17.89 18.91 18.50 5.90 5.64 5.76 5.71
FFA
84.52 50.69 29.50 18.71 12.36 8.68 6.46 4.44 3.49 2.85 1.09 1.19 1.23
85.49 48.58 31.71 18.51 13.00 10.78 9.33 7.92 2.12 2.18 2.23
92.07 77.76 62.90 51.76 42.38 30.55 23.53 19.12 1.17 1.05 0.98 0.88
TAG
Mass %
9.58 26.66 31.24 31.34 30.27 26.33 24.70 22.18 19.82 18.03 7.86 7.32 7.08
5.97 22.08 28.02 32.32 33.00 31.35 29.88 29.11 11.66 11.00 10.58
2.18 11.16 21.28 25.76 30.77 33.63 34.46 34.32 14.00 12.80 11.75 11.88
DAG
0.16 4.08 9.94 13.33 16.41 18.62 20.71 21.24 21.39 21.43 11.02 9.19 7.63
0.44 3.62 6.50 10.95 14.56 14.90 15.10 15.03 8.00 7.00 5.95
0.02 0.32 0.99 1.86 2.84 3.33 4.31 5.40 12.67 12.13 11.50 11.14
MAG
Exp.8
Exp.5
Exp. 2
0 15 30 45 60 120 180 240 300 360 420 480 540
0 30 60 120 180 240 300 360 1320 1380 1440
0 30 60 90 120 180 240 300 1260 1320 1380 1440 1500
Time(min)
Table 3. Experimental Data Used in the Parameter Fitting and Model Validation
0.00 7.10 9.68 11.86 14.71 26.75 32.00 40.33 46.71 51.97 55.10 60.33 65.75
0.00 9.91 16.99 22.53 29.59 34.00 38.59 41.00 85.46 88.51 89.50
0.00 2.94 6.33 10.04 12.80 20.11 24.91 31.74 76.50 78.32 80.20 81.00 81.20
FAME
0.74 9.00 10.84 11.77 11.72 11.17 9.69 8.38 7.07 5.72 5.06 4.22 3.81
7.55 16.02 16.91 16.90 15.50 14.00 13.00 12.46 5.50 5.00 4.60
5.73 11.47 15.25 17.60 18.69 19.59 17.73 16.39 5.06 4.80 4.65 4.65 4.64
FFA
97.80 48.72 38.75 34.23 27.87 14.91 9.93 7.46 6.13 5.44 4.68 4.24 3.99
84.36 39.92 24.11 14.90 11.34 11.00 10.12 9.42 1.48 1.26 1.22
92.07 67.70 51.62 39.99 32.54 23.46 18.01 14.79 1.30 1.13 1.01 1.05 1.03
TAG
Mass %
1.40 29.07 32.71 32.01 34.09 29.70 31.06 26.00 23.31 21.43 20.89 19.03 16.11
6.26 28.80 32.65 34.62 33.00 31.00 29.50 28.22 4.87 3.63 3.09
2.18 17.21 25.23 29.75 32.56 33.30 34.82 32.71 11.42 10.67 10.20 9.01 9.00
DAG
0.06 6.11 8.01 10.14 11.61 17.47 17.31 17.83 16.78 15.44 14.28 12.18 10.35
1.83 5.35 9.34 12.36 12.32 11.50 10.82 10.07 1.98 1.60 1.55
0.02 0.67 1.57 2.62 3.41 3.54 4.54 4.37 5.50 5.08 4.28 4.20 4.10
MAG
Exp. 13
Exp. 6
Exp. 3
0 15 30 45 60 120 180 240 300 360 420 480 540
0 30 60 90 120 180 240 300 360 1320 1380 1440
0 30 60 90 120 180 240 300 1260 1320 1380 1440 1500
Time(min)
0.00 16.30 20.34 23.54 25.68 34.44 44.36 51.97 55.04 60.50 64.29 67.35 71.08
0.00 10.91 20.22 23.06 26.00 31.00 37.00 42.87 48.00 86.85 88.00 89.00
0.00 2.87 6.54 11.79 16.12 20.10 29.22 35.01 86.00 86.22 87.26 87.91 88.62
FAME
0.22 8.31 9.48 9.79 10.13 9.97 8.59 7.29 6.07 5.69 5.53 5.38 5.08
8.03 15.70 17.05 16.49 16.18 13.56 10.81 9.11 7.72 3.22 3.10 3.00
5.73 11.78 16.18 19.15 20.04 21.00 18.26 15.96 4.40 4.35 4.30 4.29 4.20
FFA
98.10 34.23 27.67 21.94 18.09 9.71 5.84 3.85 2.99 2.92 3.20 3.19 3.21
86.59 44.50 27.58 21.29 16.08 11.71 8.23 6.71 5.10 1.16 1.27 1.26
92.07 63.61 48.14 37.10 30.06 22.25 18.08 15.14 0.99 1.03 0.97 0.90 1.00
TAG
Mass %
1.33 30.79 30.09 31.20 30.35 26.18 21.46 17.61 16.49 14.55 13.62 12.92 11.62
5.24 24.73 28.45 30.80 32.08 31.25 29.66 27.00 24.34 5.50 4.96 4.29
2.18 20.90 27.07 29.27 30.96 32.00 31.10 30.19 7.00 5.89 5.43 4.97 4.53
DAG
0.35 10.37 12.43 13.54 15.75 19.69 19.75 19.28 19.41 16.34 13.36 11.16 9.02
0.14 4.16 6.70 8.36 10.51 14.18 15.69 17.28 15.63 4.00 3.84 3.27
0.02 0.85 2.08 2.69 3.00 4.00 3.80 3.69 1.90 1.90 1.89 1.93 1.38
MAG
6 Biotechnol. Prog., 2014, Vol. 00, No. 00
0.44 1.26 2.70 4.83 0.37 1.03 2.73 3.78 5.72 13.33 17.60 24.39 87.26 76.44 64.61 49.86
9.43 8.01 6.71 5.82 4.72 3.81 3.23 2.66 2.20 1.70 1.31 1.10 0.89 0.89 0.67 15.24 12.03 11.85 10.54 8.46 7.34 6.65 5.10 3.99 3.36 2.68 2.58 1.77 1.77 1.31
DAG TAG
3.76 3.33 2.30 2.60 2.36 1.76 0.98 1.14 1.14 0.73 0.88 0.63 0.71 0.71 0.31 600 660 720 780 840 900 960 1020 1080 1140 1200 1260 1320 1380 1440
MAG
Exp. 12
0 5 10 20
0.00 5.44 9.37 18.00
5.85 8.47 8.83 11.31
87.44 69.67 57.30 39.80 6.65 7.37 8.97 10.57 0.00 1.82 6.10 11.39 0 5 10 20
6.27 15.17 21.80 26.06
7.95 7.04 5.67 4.84 4.12 3.51 3.25 2.77 2.36 2.13 1.92 1.66 1.54 1.38 1.27 4.91 4.95 4.65 4.64 4.20 4.07 3.40 4.08 3.71 3.57 3.15 2.81 2.72 3.02 2.71
FFA FAME
72.80 74.40 77.46 79.27 81.50 83.42 85.37 85.96 87.96 89.07 90.50 91.76 92.41 92.54 93.57 600 660 720 780 840 900 960 1020 1080 1140 1200 1260 1320 1380 1440
Time(min)
3.02 3.02 2.94 2.74 2.48 2.35 1.70 1.90 1.66 1.50 1.06 0.84 0.73 0.75 0.53 3.69 3.32 2.91 2.96 3.16 3.05 2.28 2.46 2.13 1.72 2.72 2.11 2.67 2.67 2.14
Results and Discussion
TAG FFA
Uncertainty analysis: Monte–Carlo simulations In this work, the Monte–Carlo technique was used for the evaluation of uncertainty in the kinetic parameters, on the model outputs. This method offers global results due to the large number of model evaluations performed using the randomly sampled parameters, to obtain the distribution of the model outputs. The method samples from the input parameter space, h and generates model outputs, y. The Monte– Carlo analysis of uncertainty involved three steps: Step 1: Specifying parameter input uncertainty: The confidence intervals of the parameters estimates, h^ were used as the upper and lower bounds of a given kinetic parameter. Step 2: Sampling input uncertainty: Latin hypercube sampling with correlation control was used to sample within the input parameter space given the dense stratification over the range of each sampled variable.37 Samples were selected from the input parameter space, where each sample, hi contained one value for each input parameter creating a [m 3 Nlh] matrix. Where m stands for the total number of model parameters, and Nlh is the total number of Latin-Hypercube samples. Step 3: Simulating the model using the sampling matrix: Nlh dynamic simulations were then performed using the [m 3 Nlh] sampled input matrix. Each simulation result was stored in a [Nt 3 Nu 3 Nlh] size array where, Nt is the length of the discrete time series and Nu is the number of model outputs. The complete Monte–Carlo results provide a cumulative distribution function for each output variable at each time instant. The uncertainty of the model outputs were then represented using mean and percentile calculations.
Parameter estimates and confidence intervals
0.89 0.74 1.11 2.21 4.95 10.92 14.87 21.51 84.38 76.76 70.71 59.93 0.00 1.47 2.58 4.09 0 5 10 20 Exp. 10
9.78 10.12 10.73 12.26
0.19 5.77 11.21 14.20 16.90 18.22 18.52 19.55 17.66 16.45 3.12 2.30 1.72 8.26 27.09 31.07 30.66 27.57 24.62 23.94 20.87 19.18 18.10 4.67 3.50 2.89
DAG TAG
84.61 44.78 25.19 15.48 10.79 8.11 6.26 4.42 4.35 3.92 1.05 0.86 0.71 0.00 8.14 17.38 25.58 32.17 37.60 41.33 47.10 51.00 54.19 87.01 89.02 90.36 0 30 60 90 120 150 180 240 300 360 1200 1320 1440
FFA FAME Time(min)
Exp. 9
6.93 14.23 15.16 14.07 12.58 11.44 9.95 8.05 7.81 7.34 4.15 4.32 4.31
MAG
Exp. 11
Time(min)
Mass % Mass % TABLE 3. Continued
parameters in k. Brun et al. determined an empirical threshold of kk being below 15, which is used in this study.35
67.88 73.31 76.22 78.07 81.30 84.03 86.86 88.65 90.55 92.49 92.41 93.57 93.96 93.96 95.56
DAG
7
FAME
Mass %
11.32 10.59 9.27 8.51 7.69 6.65 6.28 5.29 4.31 3.74 3.38 2.93 2.60 2.32 1.93
MAG
Biotechnol. Prog., 2014, Vol. 00, No. 00
The histogram of the residuals for the fitting of Exp. 1–8 were used to assess the quality of the model fitting. For the proposed model, it can be seen in Figure 3 that the histogram is not skewed which gives an indication that the complexity and choice of model is appropriate. Also the residuals have a mean of 20.05 mass % (approximately zero mean) and standard deviation of 2.64 mass %. This signifies that 95% of the residuals are within 20.05 6 5.28 mass %. This result is reasonable, given that a mass balance on the acyl groups for the experimental data close within 3 mass %. The parameter estimates for the proposed model are shown in Table 4 along with the confidence intervals and correlation matrix. Generally, the narrower the confidence interval, the higher the quality of the parameter estimate. The confidence intervals for the kinetic parameters k22, k3, and k26, (the reverse kinetic constants for formation of the TAG enzyme substrate complex (E.T), the rate of DAG production from E.T and the reverse kinetic constants for formation of the MAG enzyme substrate complex (E.M), respectively) deviate more than 40% from the mean estimates signifying a low sensitivity of the model outputs to those parameters. This may be due to the data set not having sufficient information given that the intermediate enzyme substrate complexes were not measured.
1.00 20.26 1.00 20.53 0.61 1.00 0.41 0.01 20.51 1.00 20.33 0.17 0.69 20.07 1.00 20.64 0.08 0.28 20.45 0.04 1.00 0.69 20.03 20.43 0.50 20.25 20.81 1.00 1.00 0.19 20.31 20.23 20.06 20.16 0.07 0.00 1.00 20.01 20.34 0.11 0.12 20.26 20.03 0.36 20.24 1.00 0.18 20.19 20.61 0.18 0.72 20.36 0.78 0.46 20.65 1.00 20.17 0.05 20.16 20.01 20.25 20.35 0.10 20.33 0.08 0.04 1.00 0.02 0.29 0.18 0.17 20.27 20.02 0.30 20.44 0.14 0.36 20.35 1.00 20.04 20.10 20.55 20.23 0.19 0.51 20.13 20.57 0.49 20.48 20.33 0.38 1.00 20.03 20.20 0.50 0.26 20.11 20.26 20.10 20.21 0.09 20.06 0.18 0.01 20.13 1.00 20.03 0.46 20.46 0.14 20.52 20.30 0.06 0.44 20.18 20.61 0.70 20.29 20.48 0.54 1.00 20.20 20.06 20.21 0.25 20.17 0.20 0.12 0.09 20.18 0.07 0.31 20.45 0.11 0.27 20.26 1.00 0.51 20.77 0.18 20.32 0.40 0.01 0.47 0.20 0.02 20.30 0.00 0.42 20.58 0.18 0.40 20.49 1.00 0.47 0.16 20.37 0.46 20.52 0.24 20.08 0.80 20.04 20.08 20.40 20.09 0.63 20.40 0.72 0.15 20.44 1.00 20.32 20.07 20.15 20.07 20.32 0.01 20.01 20.02 20.15 0.19 0.01 20.17 20.03 20.22 0.02 20.25 0.18 0.00
k10 k29 k9 k28 k8 k27 k7 k26 k6 k25 k5 k24 k4 k23 k3 k22 k2 k21 k1
1.00 0.60 20.81 20.38 20.13 20.21 20.40 0.30 20.12 0.08 20.70 0.14 0.10 0.21 20.12 20.54 0.04 20.68 0.03 0.26 – 冑 – 冑 – 2 冑 – – 冑 – – – 冑 冑 冑 冑 冑 冑 104 106 106 104 107 107 106 107 107 105 105 105 106 100 104 106 104 105 102 1024 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
Upper
5.88 7.78 1.96 1.63 3.32 2.51 4.15 1.59 1.79 2.21 1.01 8.41 9.41 5.86 2.84 4.84 2.85 2.44 4.08 5.08 104 100 106 103 107 107 106 107 107 105 104 105 106 100 104 106 104 105 102 1024
Lower
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
Correlation Matrix
Identifiable parameter Subset Confidence Intervals
3.30 4.26 1.45 8.78 1.53 1.75 2.68 1.06 1.36 1.35 8.31 4.14 4.84 4.45 2.00 2.66 2.08 1.66 2.69 3.05 104 100 106 104 107 106 106 107 107 105 104 107 106 100 104 106 104 105 102 1024 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
Having generated parameter estimates, key to the analysis is to find out how specific the parameter estimates actually are, given the experimental data used. One simple method is to look at the correlation between the parameters (Table 4). In this study 60.75 was used to signify highly correlated parameters. This value is chosen based on the correlation value of the inhibition constant (20.81) given the inhibition
Parameter
Correlated and identifiable parameters
Mean Estimate
For Exp. 9–12 the histogram of the residuals (Figure 4) for the validation data set is slightly skewed to the left indicating the model tends to underestimate the concentrations. The performance of the parameter estimates over the entire time course of the reaction is illustrated in Figure 5, using validation data set Exp. 9. The proposed model captures the dynamics for the five components over the entire course of the reaction, although the prediction for FFA and MAG shows some deviation from the experimental data. We also investigated the performance of the model in the prediction of the concentrations of TAG, DAG, MAG, FAME, FFA, and CH for the first 20 min of the reaction (Figure 6), for various enzyme and water concentrations (Note: the water concentrations for Exp. 11 and 12 are outside the range used for the model fitting but are included to investigate the extrapolability of the mechanistic model). The model follows the expected trends of the experimental data. For example, the FAME production increases, with increasing water concentration. The available interfacial area is larger, which means there is an increased chance for substrates to react and hence an increase in the rate of FAME production.38 Likewise, as the enzyme concentration increases from 0.2 to 0.5 wt % enzyme, the production of FAME increases as expected. The same reasoning extends for the other components plotted. However, at the higher water concentrations, the model tends to over predict the amount of FFA produced, reducing the FAME production. The model mismatch observed may be due to process phenomena not taken into account. For example, the viscosity of the reaction media changes one order of magnitude over the 24 h.39 Hence the parameter estimates are average values of the rate constants over the entire course of the reaction. Likewise, the uncertainty in the parameter estimates plays a part given some of the parameter estimates are strongly correlated.
Table 4. Parameter Estimates, Related Confidence Intervals and Correlation of the Parameters for 10000 Bootstrap Samples. Highlighted Values Represent Highly Correlated Values
Figure 3. Histogram of residuals for the fitting of the proposed model for Exp. 1–8. The distribution has a mean of 20.05 mass % and a standard deviation of 2.64 mass %.
4.95 6.60 1.69 1.11 2.07 2.20 3.41 1.33 1.55 1.81 9.13 5.43 7.06 4.93 2.36 3.51 2.54 2.05 3.23 4.39
k210
Biotechnol. Prog., 2014, Vol. 00, No. 00
k1 [L/mol.s] k21 [1/s] k2 [L/mol.s] k22 [1/s] k3 [1/s] k23 [L/mol s] k4 [L/mol s] k24 [1/s] k5 [1/s] k25 [L/mol s] k6 [L/mol s] k26 [1/s] k7 [L/mol s] k27 [L/mol s] k8 [L/mol s] k28 [L/mol s] k9 [L/mol s] k29 [L/mol s] k10 [L/mol s] k210 [1/s]
8
Biotechnol. Prog., 2014, Vol. 00, No. 00
constants are usually strongly correlated.40 For two highly correlated parameters, the change in model output due to changing one of the correlated parameters, can be compensated for by an appropriate change of the other parameter value, preventing a unique estimate of the parameter value. This can be due to the model structure or the similarity of the parameters of the underlying biological system.41 What this may signify is that enzyme activities can be modified by changing one of the correlated parameters. For example k10 is highly correlated with k210. The inhibition constant has a negative coefficient for the correlation value. If one parameter value increases, the other decreases. Hence a modification the enzyme structure that affects inhibition can potentially have a twofold effect. If the forward rate where the dead end complex E.CH is reduced, then the rate of disassociation of E.CH will increase. The correlation seen amongst some of
9
the parameters should be expected given the complex parallel and sequential reactions occurring simultaneously. The plot of the Collinearity index in Figure 7 shows that the Collinearity index increases with the number of parameters and that only a maximum of 10 of the 20 parameters can be identified with the available experimental data. In Table 4 one potential subset was identified, that takes into account the parameters that are correlated. The ticked (冑) parameters were the ones estimated and the others were fixed. Hence, the procedure is iterative, although in this case it only gives a reduction in the squared-sum of the relative errors between the simulated and experimental of 0.01%. It should be noted that fixing parameters, while estimating others, results in reasonable parameter values rather than “true parameter values.” 35 Given the uncertainty in the parameter estimates, we then look at how the model can be used for engineering purposes. Uncertainty analysis: Monte–Carlo simulations
Figure 4.
Histogram of residuals for the fitting of the proposed model for Exp. 9–12. The distribution has a mean of 20.005 mass % and a standard deviation of 2.33 mass %.
The uncertainty in the model outputs for the typically measured variables (TAG, DAG, MAG, FAME, and FFA) during the transesterification reaction are seen in the spaghetti plots of the 500 Latin hypercube samples in Figure 8. When the experimental values are overlaid on the Monte– Carlo simulations, we can see that over the complete course of the reaction, most of the experimental values fall within the bounds of the spaghetti plots. The narrow prediction bands for FAME and TAG reflect the robustness of the predictions for those model outputs over the entire course of the reaction, while the wide bands observed for FFA and MAG show the need for a more accurate estimate of the parameters in order to obtain more certain model predictions. Using the cumulative frequency distribution plots (Figure 9) it is possible to put bounds on the model predictions, which can give the modeler some insight into the reliability of the model to make predictions. Take for example the FAME predictions. At the end of the reaction the model output has a mean value of 87.9 mass % with a standard deviation of 0.64 mass %.
C O L O R Figure 5. Comparison of the model fitting (2) to the Validation data set Exp 9 (*).
10
Figure 6.
Biotechnol. Prog., 2014, Vol. 00, No. 00
Comparison of the model fitting for the first 20 min of the reaction for experiments 10 (* experimental, - simulated), 11 ( experimental,- simulated) and 12 (䉮 experimental, . . . simulated).
Given the uncertainty in the model parameters, the model gives excellent predictions of the FAME and TAG values and shows the deficiencies in the FFA and DAG predictions during the course of the reaction. We hence see this method as a valuable tool to gauge the robustness of a model to parameter uncertainty. Engineering application of the model given the parameter uncertainty For process development, a reliable kinetic model can be used, for example, to simulate and evaluate variations in feed composition, alternative reactor configurations and feeding strategies to mitigate methanol inhibition to name a few. The results from the simulations can then be tested experimentally. Below we investigate a methanol feeding strategy to mitigate inhibition and the uncertainty in the model outputs due to the parameter estimates. Control of Methanol Feeding. One hurdle to industrial implementation of enzymatic biodiesel production is inhibition and deactivation of the biocatalyst by the alcohol substrate. Simulations can be used to devise an optimal feeding policy. We followed a methanol feeding strategy similar to the one proposed by Samukawa et al.42 They found that they could increase the reuse of the immobilized enzyme (a clear indication of a reduction in enzyme deactivation), by keeping the methanol content in the reactor below the concentration that gave the highest initial rate of FAME production (CHcriti42 In our previous work we found CHcritical to be 0.525 cal). 43 Eq. We then combined the model of the system (Eq. 4) with the objective function in Eq. 9. By minimizing the objective function in Eq. 9 we ensure that the methanol concentration in the reactor never goes above the critical value CHcritical at each time step ti by manipulating the methanol feed Fa. To simplify the experimental procedure only two step changes in the methanol feed rate were used. min JEq 5ðfCHEq gti 2fCHcritical gti Þ2
(9)
Fa
The optimization results are now compared to the experimental results as illustrated in Figure 10. The optimization
Figure 7. Plot of how the collinearity index varies with all parameters (Top plot) and how the collinearity index varies with the identifiable parameters for a threshold value of 15 (bottom plot).
objective to constrain the amount of methanol in the reactor was validated experimentally (Exp. 13) and compared to the Exp. 9 where the methanol concentration is over 0.6 Eq. for the last 7 h of the reaction (Figure 8). The optimised case respects the constraint, with the methanol concentration never going above the critical value of 0.525 Eq. for the entire reaction. However, for the optimized case, there was a 2% reduction in the final biodiesel concentration compared to Exp. 8, which had the highest FAME conversion. The FAME equilibrium concentrations at the end of the reaction can be increased by increasing the methanol concentration at the end of the reaction. However, we are then exposing the enzyme to higher concentrations of methanol, which potentially reduces the number of times the enzyme could be reused. What is interesting is the trade-off between downstream processing to bring the final biodiesel concentration within specifications and potential increase in enzyme reuse. However, this analysis is beyond the scope of this article.
Biotechnol. Prog., 2014, Vol. 00, No. 00
11
Figure 8. Uncertainty analysis for the validation experiment (Exp 9). The experimental values (*) are overlaid on the 500 Monte– Carlo simulations (2).
Figure 9. Cumulative distribution function for the 500 Monte–Carlo simulations (Exp 9) at time 24 h.
Characterization of the Model Uncertainty for the Optimization. The same parameter input uncertainty used for the validation data set was used for the uncertainty analysis given we are still operating within the calibrated range of the model. In Figure 10 the narrow prediction bands for TAG, DAG, and FAME reflects the robustness of the predictions. When the experimental values are overlaid on the 500 Monte–Carlo simulations, we can see over the course of the reaction that most of the model outputs fall within the bounds of the spaghetti plots except most notably for the FFA prediction. The concentrations for the FFA predictions are on the same order of magnitude, although the dynamics after 5 h for the FFA simulation show a slight increase followed by a decrease in FFA concentration compared to the experimental values that show a steady decrease. As postu-
lated earlier, given the decrease in viscosity of the reaction system as time progresses, it is believed that the rate of FFA consumption increases during the reaction and hence the steady decrease in the FFA concentration seen in the experimental data. Given the rate constants are average values in the simulation, it is not possible to capture the behaviour seen. The cumulative frequency distribution plots (Figure 11) are then used to characterize the uncertainty in the model outputs. In this case, at the end of the reaction, the FAME model output has a mean value of 90.83 mass % with a standard deviation of 0.55 mass %. Limitations and Use of the Model. The developed model is only applicable for relatively low methanol concentrations
12
Biotechnol. Prog., 2014, Vol. 00, No. 00
Figure 10. Uncertainty analysis for Exp 13. The experimental values (*) are overlaid on the 500 Monte–Carlo simulations (-).
Figure 11. Cumulative distribution function for the 500 Monte–Carlo simulations (Exp 13) at time 24 h.
(
In this article, a kinetic model for the enzymatic transesterification of rapeseed oil with methanol using CalleraTM Trans L (a liquid formulation of a modified Thermomyces lanuginosus lipase) was developed from first principles. We base the model formulation on a PingPong Bi-Bi mechanism. Methanol inhibition, along with the interfacial and bulk concentrations of the enzyme was also modeled. The model was developed to describe the effect of different oil compositions, as well as different water, enzyme, and methanol concentrations, which are relevant conditions needed for process evaluation, with respect to the industrial production of biodiesel. The developed kinetic model, coupled with a mass balance of the system, was fitted to and validated on experimental results for the fed-batch transesterification of rapeseed oil. The confidence intervals of the parameter estimates, along with the identifiability of the model parameters were presented. The predictive capability of the model was tested for a case using 0.5% (wt. Enzyme/wt. Oil), 0.5% (wt. Water /wt. Oil) and feeding 1.5 times the stoichiometric amount of methanol in total over 24 h. For this case, an optimized methanol feeding profile that constrains the amount of methanol in the reactor was computed and the predictions experimentally validated. Monte-Carlo simulations were then used to characterize the effect of the parameter uncertainty on the model outputs, giving a C 2014 American Institute biodiesel yield, based on the mass of oil, of 90.8 6 0.55 mass %. V of Chemical Engineers Biotechnol. Prog., 000:000–000, 2014 Keywords: enzymatic biodiesel, kinetic modeling, transesterification, uncertainty analysis
Introduction For enzymatic biodiesel production, the conventional biocatalyst is usually immobilized in order to improve enzyme recovery, which combined with increased stability, allows for reuse.1 However, the immobilization carrier, as well as the immobilization process significantly contributes to the price of the biocatalyst, which further necessitate the reuse of the enzymes for the process to be competitive.1–3 Within the last few years, the use of liquid lipase formulations for enzymatic biodiesel production has resulted in a significant reduction in the biocatalyst cost.4–6 Furthermore, compared to the conventional alkali catalysts used to produce biodiesel on an industrial scale, the use of an enzymatic catalyst has the advantage that low-quality feedstock and waste oils that have a high free fatty acid (FFA) content can be treated. This is due to the fact that lipases are able to esterify the FFA contained in waste oils to esters as well as transesterify the acyl-glycerides in the oil.7 This further reduces operating costs of the enzyme-catalysed biodiesel process. Nevertheless, when developing an industrial enzymatic biodiesel process a few issues related to the biocatalyst need to be addressed; methanol inhibition, deactivation at high methanol concentrations, the limited lifespan of the lipase and the relatively high cost of the enzyme per kilogram of oil treated
Correspondence concerning this article should be addressed to J. K. Huusom at [email protected]. C 2014 American Institute of Chemical Engineers V
compared to traditional chemical catalysts (e.g., sodium hydroxide).8–10 All of the issues outlined can be mitigated by operating the enzyme-catalysed biodiesel process in an optimal manner. In terms of developing and optimizing the enzymecatalyzed biodiesel process efficiently, process modeling is a valuable tool to help focus the experimental work needed for process understanding and to support further process development.11 For example, a mechanistic process model of the system can be used to simulate the process performance over a wide range of conditions and to optimize how the process is operated. Hence, a process model can further be used to develop control strategies to mitigate enzyme inactivation and improve enzyme stability. Likewise, modelling can help in the formulation of innovative process designs and configurations.12 Integral to the development of a mechanistic model of any process, is the availability of reliable kinetic models. Over the years, various kinetic models for the enzymatic transesterification of vegetable oils have been proposed.2,13–18 However, only the model by Lv et al. deal with model development for a liquid lipase formulation.2 Likewise, in terms of engineering design purposes such as reactor configuration, optimization, and control, the gap we noticed in the literature for kinetic models for enzymatic biodiesel production is that there was a lack of emphasis placed on: 1. Use of the proposed kinetic models in the combined act of process optimization and experimental validation of the process optimization. This would enable the modeler to 1
2
Biotechnol. Prog., 2014, Vol. 00, No. 00
ascertain how well the model could be extended outside the model fitting. 2. Statistical analysis of the working bounds of the model. This remains a weak point in the credibility of these kinetic models and hence their applicability for engineering design purposes. In this contribution we address the two aforementioned points so that that we can perform model-based engineering design. In using the model for engineering design, it is desired that the kinetic model describing enzymatic biodiesel production, can predict the concentration of all the major species in the reaction. It is also essential to be able to characterize how the process responds to changes in the process conditions over the entire course of the reaction for changes in: 1. Alcohol/oil molar ratios: This is important given the need to balance the amount of methanol needed to shift the final equilibrium conversion of fatty acid methyl esters (FAME, biodiesel) while minimizing the effects of methanol inhibition and deactivation.1,15,19 2. Concentration of reactants in the reactor: The oil composition from batch to batch can vary. The oil composition then needs to be characterized so as to ascertain when the reaction has reached within specification. 3. Enzyme loading: The amount of biocatalyst added will affect reaction time, enzyme efficiency, and potential reuse of the enzyme. 4. Area of the oil–water interface: The liquid lipase CalleraTM Trans L (a liquid formulation of a modified Thermomyces lanuginosus lipase) being used is interfacially activated and hence this phenomena needs to be taken into account when modeling the system.20,21 None of the aforementioned models take into consideration all of the process conditions outlined for a liquid lipase formulation (As mentioned previously, the model by Lv et al. deal with model development for a liquid lipase formulation.2 However, their model formulation does not enable the evaluation of changes in enzyme concentration and the area of the oil-water interface). Hence in this work we: Develop and validate a dynamic mechanistic model from first principles for the transesterification of rapeseed oil with methanol using CalleraTM Trans L, which takes into consideration the effects of the process conditions outlined. Test and experimentally validate the predictive capabilities of the model by optimizing the methanol feed profile. Evaluate the model structure for identifiability of the kinetic parameters (Identifiability Analysis-Correlation matrix and Collinearity index) and characterize the uncertainty in the model outputs due to the uncertainty in the parameter estimates (Uncertainty Analysis-Monte-Carlo Simulations). The article is organized as follows: The proposed reaction mechanism is presented, followed by the experimental and numerical methods used. Subsequently, the results from parameter estimation and identifiability analysis are discussed. Finally, the model is used to predict an optimal methanol feeding profile and the uncertainty analysis is used to characterize the uncertainty in the model outputs.
Reaction Mechanism Model formulation The mathematical model describing the transesterification reaction in the biphasic oil–water system with a liquid lipase,
Figure 1.
Diagrammatic representation of the enzyme at the oil– water interface. The polar phase contains water, methanol, glycerol, and the free enzyme (Ebulk). The nonpolar phase contains the oil components along with the biodiesel formed. At the interface is the penetrated enzyme (E) and the acyl enzyme complex (EX).
CalleraTM Trans L was formulated on the basis of the following assumptions: 1. The reaction proceeds via a Ping-Pong Bi-Bi mechanism.13,15,18,22 2. Alcohol inhibition is competitive.15,18 3. Deactivation due to the alcohol could be ignored at low methanol concentrations.15,18 4. The interfacial and bulk concentrations of the substrate and products are the same (mass transfer from the bulk to the interface is instantaneous).15,18 5. Acyl migration is neglected.18 6. All reaction steps are reversible.18,23 7. The reactor mixture is homogenous (ideal mixing) and the density of the mixture is constant.18 Of the various kinetic models reported,2,13–18 the mathematical formulation of the kinetic model presented by Fedosov and coworkers fulfil the first three of the process conditions outlined in the introduction section. Since their model was developed for an immobilized enzyme, it does not describe the behavior of the liquid lipase at the oil-water interface. Hence, we extend their work by modeling the oil– water interfacial area and present a fully mechanistic formulation for the transesterification reaction using a liquid lipase. In their formulation, various pseudo-components were introduced to imitate the phase boundaries of the system, for the transesterification of rapeseed oil using Novozym 435 (immobilized Candida Antarctica lipase B).18 Key to how we describe the reaction mechanistically for the liquid lipase, is the interaction of the enzyme at the oil–water interface. A schematic illustration of the oil water interface along with the enzymes and its complexes is presented in Figure 1. The triglycerides (T), diglycerides (D), monoglycerides (M), biodiesel (BD), and FFA occupy the nonpolar phase while the bulk enzyme (Ebulk), Water (W), Methanol (CH), and glycerol (G) occupy the polar phase. The lipase used in this study exhibits a pronounced interfacial activation and the reaction is assumed to proceed exclusively at the interface.24 By including the interfacial enzyme concentration (E), the reaction scheme proceeds as shown in Table 1 and Figure 2. The specific interfacial area per unit volume of the oil water mixture (aT [m2/m3]) can be represented as:
Biotechnol. Prog., 2014, Vol. 00, No. 00
3
Table 1. Kinetic Mechanism for the Enzymatic Transesterification i
Reactions
1 2 3 4 5 6 7
Ebulk 1 Af $ E T 1 E $ E.T E.T $ EX 1 D D 1 E $ E.D E.D $ EX 1 M M 1 E $ E.M E.M $ EX 1 G
8 9 10
EX 1 W $ FA 1 E EX 1 CH $ BD 1 E CH 1 E $ E.CH
Figure 2.
Rate of Reaction (ri)
Enzyme in bulk absorbed at the interface In reactions 2, 4, and 6 the penetrated enzyme can react with the substrate to form an enzyme substrate complex E.T, E.D, or E.M (Ping) In reactions 3, 5, and 7 the enzyme substrate complex forms the Acyl enzyme complex and releases the first product D, M or G (Pong) Intermediate steps for the reactions were grouped together given interest is in the overall rate. The acyl enzyme complex can then react with water or methanol (Pong) and then release the second product FA or BD (Ping). Reversible competitive methanol inhibition.
k1 ½Ebulk ½Af 2k21 ½E k2 ½T ½E2k22 ½ET k3 ½ET2k23 ½EX ½D k4 ½D ½E2k24 ½ED k5 ½ED2k25 ½EX ½M k6 ½M ½E2k26 ½EM k7 ½EM2k27 ½EX ½G k8 ½EX ½W2k28 ½FA ½E k9 ½EX ½CH2k29 ½BD ½E k10 ½CH ½E2k210 ½ECH
Conceptual scheme for the overall reaction mechanism describing the enzymatic transesterification.
aT 5
6 Vp : ds V
(1)
where ds is the Sauter mean diameter of the droplets in the system [m], Vp is the polar volume [m3] and V is the bulk volume [m3].25,26 An enzyme coverage Ae, of 2.98 3 107 m2/mol, was used in this study, which denotes the interfacial area that is required for the adsorption of 1 mol enzyme to the interface.20 It is assumed that all forms of the adsorbed enzyme molecules, meaning both the free enzyme as well as all
forms of the enzyme complexes and the inhibition complex, occupy the same area at the interface.20 Given Ae, it is possible to calculate the free specific interfacial area, af [m2/m3] as shown in Eq. 2. af 5aT 2Ae :ðE1EX1ET1ED1EM1ECHÞ
(2)
The free specific interfacial area can then be expressed as a volumetric concentration (Af [mol/m3]), by using the enzyme coverage to estimate a theoretical upper limit of the moles of enzyme molecules that can occupy the free interfacial area.
4
Biotechnol. Prog., 2014, Vol. 00, No. 00 Af 5af =Ae
(3)
In the model formulation, all concentrations used are lumped concentrations (concentration of the component in the entire reaction volume). The mass balance for the system can then be combined with the kinetics to give the system of ordinary differential equations presented below. dð½TVÞ
. dt
dð½DVÞ
52Vðr2 Þ
.
dt
dð½MVÞ
dð½FAVÞ
dð½GVÞ
.
dt
.
dt
.
dt
dð½WVÞ
.
5Vðr1 1r8 1r9 2r2 2r4 2r6 2r10 Þ
dt
dð½EXVÞ
52Vðr9 1r10 Þ
dt
. .
.
5Vðr2 2r3 Þ
dt
dð½E:DVÞ
.
dt
dð½E:MVÞ
5Vðr5 2r6 Þ
.
dt
dð½E:CHVÞ
dð½Ebulk VÞ
dðVp Þ
dðVÞ
.
dt
dt
dt
dt
.
5Vðr6 2r7 Þ
.
.
(4)
5Vðr3 1r5 1r7 2r8 2r9 Þ
dt
dð½E:TVÞ
Fed-batch experiments
5Vðr8 Þ
52Vðr8 Þ
dt
dð½EVÞ
5Vðr9 Þ
5Vðr7 Þ
.
dð½CHVÞ
CalleraTM Trans L with a hydrolytic activity of 1 3 105 LU/g was kindly donated by Novozymes A/S (Bagsværd, Denmark). One LU is defined as the activity required to produce 1 lmol butyric acid in the hydrolysis of tributyrin under standard conditions (pH 7.5, 0.2M substrate).24
5Vðr5 2r6 Þ
dt
dð½BDVÞ
Biocatalyst
5Vðr3 2r4 Þ
.
5Vðr10 Þ
52Vðr1 Þ
For the 13 experiments (Table 3), the water and enzyme content were varied from 3 to 7 and 0.1 to 0.5 wt % oil, respectively. In all the experiments 1.5 equivalents (Eq.) of methanol was reacted with the rapeseed oil. One equivalents corresponds to the stoichiometric amount of alcohol needed to convert all fatty acid residues in the oil to biodiesel (i.e., 1 mol oil:3 mol alcohol). The reaction was carried out in a 0.25 L glass reactor with a tank diameter of 55 mm (T) and 2 baffles, each 0.18 3 T wide. The reactor was immersed in a water bath with temperature control (Julabo Labor-technik GmbH, Seelbach, Germany) maintained at 35 C. A rushton turbine (impeller diameter 0.44 T), spinning at 1400 rpm provided the mixing. Initially 0.2 Eq methanol was charged with the oil in the reactor. When the reaction mixture reached the reaction temperature, the amount of water and enzyme to be used in the experiment (Table 2), was then added to the reactor and methanol feeding started. Methanol feeding was provided by a KNF STEPDOS .03 pump (KNF Neuberger AB, Stockholm, Sweden), calibrated prior to each experiment. Sample preparation
5RG 1RW
5ðFa Þ
The measurement vector ym,i is then shown in Eq. 5 where ym is the measurement matrix [mass %], x are the state variables [mol/L], V is the bulk volume and rmm is the relative molecular mass of component i.
xi V rmmi ym;i 5 X5 3100; x V rmmi i51 i
(Brïndby, Denmark). Absolute methanol (99.8%, technical grade) was purchased from VWR Bie & Berntsen A/S (Herlev, Denmark). n-Heptane (99%), acetic acid (99%), isopropanol (99%), and tert-butyl methyl ether (99.8%) for HPLCAnalysis were obtained from Sigma-Aldrich A/S (Brïndby, Denmark).
0 TAG DAG ym 5 MAG FAME FFA
0 T D and x5 M BD FA (5)
Experimental Materials and Methods Chemicals Rapeseed oil was obtained from Emmelev A/S (Otterup, Denmark) and oleic acid was purchased from Sigma-Aldrich
Fifty-microliter samples were taken from the reactor and mixed with 500 lL solvent A (acetic acid and n-heptane 4:1000 v/v – mobile phase). Samples were then centrifuged at 14,500 rpm for 5 min and 10 lL of the supernatant was mixed with 990 lL of solvent A prior to the HPLC analysis. HPLC analysis Forty-microliters of the prepared sample was injected in the HPLC (Ultimate 3000, Dionex A/S, Hvidovre, Denmark) for analysis of triglycerides (TAGs), diglycerides (DAGs), monoglycerides (MAGs), FFAs and fatty acid methyl esters (FAME). The separation of the different compounds was carried out with a cyanopropyl column (0.25 3 0.004 m) (DisR , Cyano, Sigma Aldrich A/S, Brïndby, Denmark), covery V U3000 auto-sampler, TCC-3000SD column oven, U3400A R Charged Aerosol quaternary pump modules and a CoronaV Detector (Thermo Scientific Dionex, Chelmsford, MA). A binary gradient program was employed for the separation of the different compounds using Solvent A, Solvent B (99.6% v/v tert-butyl methyl ether and 0.4% v/v acetic acid) and iso-propanol as Solvent C.27,28 The detection of the different compounds after separation with the column was carried out R Charged Aerosol Detector from Thermo Sciby a CoronaV entific Dionex (Chelmsford, MA) with nitrogen gas at a
Biotechnol. Prog., 2014, Vol. 00, No. 00
5
Table 2. Experiments for the Data Fitting, Validation, and Optimization
Parameter Fitting
Validation
Optimization
Exp.
Methanol Feed Rate [Eq./h]
Initial Dose Methanol [Eq]
Water [wt.% oil ]
Enzyme [wt.% oil ]
1 2 3 4 5 6 7 8 9 10 11 12 13
0.06 0.06 0.06 0.06 0.06 0.1 2hrs. 0.06 2hrs. 0.06 2hrs. 0.06 0 0 0 3hrs. 0.02
0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.4 0.4 0.4 0.525
3 3 3 5 5 5 5 5 5 5 7 7 5
0.1 0.2 0.3 0.2 0.5 0.3 0.2 0.5 0.3 0.5 0.2 0.5 0.5
0.185 first 0.185 first 0.185 first
0.152 first
thereafter thereafter thereafter
thereafter
pressure of 241 kPa. The composition of the reaction samples was reported on a mass percentage basis, relative to the sum of quantified mass of the five analyzed components (TAG, DAG, MAG, FFA, and FAME). From previous experiments the standard deviation of the measurements, for TAG, DAG, MAG, FFA, and FAME were found to be 0.40, 0.75, 0.18, 0.28, and 0.26 mass %, respectively. Given that the oil composition can vary from different batches of oil, the initial concentration of the oil at the start of each experiment is measured.
Numerical Methods Simulation environment R The model was implemented and simulated in MatlabV (The Mathworks, Natick, MA). All the methods for performing the identifiability, sensitivity, and uncertainty analysis were built on the toolbox based on the work by Sin et al.29 The following sections give further details of the methods used.
Parameter estimation, confidence intervals, and identifiability analysis Parameter Estimation and Confidence Intervals. The 20 unknown kinetic constants (k12k10, k212k210), were estimated by fitting the model equations with full time course data using experiments 1–8, which covered the span of our operating conditions (Table 3). To judge the quality of the fitting, the validation data sets were chosen so that the operating conditions fall within the parameter fitting dataset (Experiments 9) and also to evaluate how the model predicts the initial rates of the components (Experiments 10–12). The differential equations were solved using a stiff variable order solver based on numerical differentiation formulas (ode15s). For the parameter fitting, the squared-sum of the relative errors between the simulated and experimental values for TAG, DAG, MAG, FAME, and FFA were minimized using fminsearch (based on a simplex search algorithm).30 To quickly assess the quality of the data fitting, the histogram of residuals were used to examine the underlying statistical assumptions of the residuals having zero mean and being normally distribution. Scott’s method, was used to determine the number of bins and is based upon the sample standard deviation and number of data points.31 A robust method to calculate realistic confidence intervals for the parameter estimates h^ was made by using a bootstrap method.32,33 To perform the bootstrap analysis, the residuals from the parameter estimates were used to generate synthetic data, which was subsequently used in a Monte–Carlo method
10,000 times to generate simulated data (Increasing the sample size from 5,000 to 10,000 did not change parameter distribution significantly). The simulated data was then used to generate a set of re-estimated parameters (matrix size [10,000 3 20]), where the mean of the distribution was used as the mean parameter estimate. The 95th and 5th percentiles of the re-estimated parameters were then used as the upper and lower bounds of the confidence intervals for the parameters estimates, respectively.34 Identifiability Analysis. The identifiability method based on the work by Brun et al. was used to ascertain, which parameters could actually be identified from the available experimental data given the model structure.35 The method follows three main steps. (1) Calculation of the sensitivity matrix, (2) scaling of the sensitivity matrix, and (3) calculation of the collinearity index for the subset of parameters. Step 1: Calculate sensitivity matrix S: The sensitivities of the model outputs (T, D, M, BD, and FA) to the parameters were calculated by the direct differential method.36 S has dimensions N 3 m (N is the number of experimental data points and m is the number of parameters). The sensitivities of the model outputs were placed in the sensitivity matrix, S: @xi S5 @hj h^
(6)
where i is the total number of observations for the five model outputs for the 8 experiments, and j 5 1:m for each of the parameters. Step 2: Scaling of the sensitivity matrix, S: The nondimensional sensitivity matrix si,j and the normalised sensitivity matrix S^i;j were then computed: Si;j 5Si;j
Dhj si;j and S^i;j 5 sci ksj k
(7)
Where the mean estimate of hj is used for Dhj, the mean of the experimental observations for each model output (i 5 1:5) is used for sci and kSi;j k is the Euclidean norm of the jth column of si,j. Step 3: Calculation of the collinearity index for the different subset of parameters: The collinearity index, ck for the subset of parameters k (k 5 2:M) is then: 1 ck 5 pffiffiffiffiffi kk
(8)
where kk is the smallest eigenvalue of ^s kT^s k with ^s k being a N 3 k submatrix of S^i;j whose columns correspond to the
FAME
0.00 0.90 3.09 6.19 8.32 14.61 18.79 22.66 68.00 70.00 71.00 71.50
0.00 10.75 16.97 22.12 26.00 30.50 33.77 37.36 73.44 75.50 76.92
0.00 7.73 16.85 24.45 29.85 36.58 39.35 44.54 49.00 52.20 77.11 79.30 81.10
Time(min)
0 30 60 90 120 180 240 300 1260 1380 1440 1500
0 30 60 120 180 240 300 360 1320 1380 1440
0 30 60 90 120 150 180 240 300 360 1200 1320 1440
Exp. 1
Exp. 4
Exp. 7
5.75 10.84 12.48 12.16 11.11 9.79 8.79 7.60 6.31 5.48 2.91 2.99 2.96
8.10 14.97 16.80 16.50 15.00 14.00 12.50 11.00 4.51 4.40 4.32
5.73 9.86 11.74 14.44 15.69 17.89 18.91 18.50 5.90 5.64 5.76 5.71
FFA
84.52 50.69 29.50 18.71 12.36 8.68 6.46 4.44 3.49 2.85 1.09 1.19 1.23
85.49 48.58 31.71 18.51 13.00 10.78 9.33 7.92 2.12 2.18 2.23
92.07 77.76 62.90 51.76 42.38 30.55 23.53 19.12 1.17 1.05 0.98 0.88
TAG
Mass %
9.58 26.66 31.24 31.34 30.27 26.33 24.70 22.18 19.82 18.03 7.86 7.32 7.08
5.97 22.08 28.02 32.32 33.00 31.35 29.88 29.11 11.66 11.00 10.58
2.18 11.16 21.28 25.76 30.77 33.63 34.46 34.32 14.00 12.80 11.75 11.88
DAG
0.16 4.08 9.94 13.33 16.41 18.62 20.71 21.24 21.39 21.43 11.02 9.19 7.63
0.44 3.62 6.50 10.95 14.56 14.90 15.10 15.03 8.00 7.00 5.95
0.02 0.32 0.99 1.86 2.84 3.33 4.31 5.40 12.67 12.13 11.50 11.14
MAG
Exp.8
Exp.5
Exp. 2
0 15 30 45 60 120 180 240 300 360 420 480 540
0 30 60 120 180 240 300 360 1320 1380 1440
0 30 60 90 120 180 240 300 1260 1320 1380 1440 1500
Time(min)
Table 3. Experimental Data Used in the Parameter Fitting and Model Validation
0.00 7.10 9.68 11.86 14.71 26.75 32.00 40.33 46.71 51.97 55.10 60.33 65.75
0.00 9.91 16.99 22.53 29.59 34.00 38.59 41.00 85.46 88.51 89.50
0.00 2.94 6.33 10.04 12.80 20.11 24.91 31.74 76.50 78.32 80.20 81.00 81.20
FAME
0.74 9.00 10.84 11.77 11.72 11.17 9.69 8.38 7.07 5.72 5.06 4.22 3.81
7.55 16.02 16.91 16.90 15.50 14.00 13.00 12.46 5.50 5.00 4.60
5.73 11.47 15.25 17.60 18.69 19.59 17.73 16.39 5.06 4.80 4.65 4.65 4.64
FFA
97.80 48.72 38.75 34.23 27.87 14.91 9.93 7.46 6.13 5.44 4.68 4.24 3.99
84.36 39.92 24.11 14.90 11.34 11.00 10.12 9.42 1.48 1.26 1.22
92.07 67.70 51.62 39.99 32.54 23.46 18.01 14.79 1.30 1.13 1.01 1.05 1.03
TAG
Mass %
1.40 29.07 32.71 32.01 34.09 29.70 31.06 26.00 23.31 21.43 20.89 19.03 16.11
6.26 28.80 32.65 34.62 33.00 31.00 29.50 28.22 4.87 3.63 3.09
2.18 17.21 25.23 29.75 32.56 33.30 34.82 32.71 11.42 10.67 10.20 9.01 9.00
DAG
0.06 6.11 8.01 10.14 11.61 17.47 17.31 17.83 16.78 15.44 14.28 12.18 10.35
1.83 5.35 9.34 12.36 12.32 11.50 10.82 10.07 1.98 1.60 1.55
0.02 0.67 1.57 2.62 3.41 3.54 4.54 4.37 5.50 5.08 4.28 4.20 4.10
MAG
Exp. 13
Exp. 6
Exp. 3
0 15 30 45 60 120 180 240 300 360 420 480 540
0 30 60 90 120 180 240 300 360 1320 1380 1440
0 30 60 90 120 180 240 300 1260 1320 1380 1440 1500
Time(min)
0.00 16.30 20.34 23.54 25.68 34.44 44.36 51.97 55.04 60.50 64.29 67.35 71.08
0.00 10.91 20.22 23.06 26.00 31.00 37.00 42.87 48.00 86.85 88.00 89.00
0.00 2.87 6.54 11.79 16.12 20.10 29.22 35.01 86.00 86.22 87.26 87.91 88.62
FAME
0.22 8.31 9.48 9.79 10.13 9.97 8.59 7.29 6.07 5.69 5.53 5.38 5.08
8.03 15.70 17.05 16.49 16.18 13.56 10.81 9.11 7.72 3.22 3.10 3.00
5.73 11.78 16.18 19.15 20.04 21.00 18.26 15.96 4.40 4.35 4.30 4.29 4.20
FFA
98.10 34.23 27.67 21.94 18.09 9.71 5.84 3.85 2.99 2.92 3.20 3.19 3.21
86.59 44.50 27.58 21.29 16.08 11.71 8.23 6.71 5.10 1.16 1.27 1.26
92.07 63.61 48.14 37.10 30.06 22.25 18.08 15.14 0.99 1.03 0.97 0.90 1.00
TAG
Mass %
1.33 30.79 30.09 31.20 30.35 26.18 21.46 17.61 16.49 14.55 13.62 12.92 11.62
5.24 24.73 28.45 30.80 32.08 31.25 29.66 27.00 24.34 5.50 4.96 4.29
2.18 20.90 27.07 29.27 30.96 32.00 31.10 30.19 7.00 5.89 5.43 4.97 4.53
DAG
0.35 10.37 12.43 13.54 15.75 19.69 19.75 19.28 19.41 16.34 13.36 11.16 9.02
0.14 4.16 6.70 8.36 10.51 14.18 15.69 17.28 15.63 4.00 3.84 3.27
0.02 0.85 2.08 2.69 3.00 4.00 3.80 3.69 1.90 1.90 1.89 1.93 1.38
MAG
6 Biotechnol. Prog., 2014, Vol. 00, No. 00
0.44 1.26 2.70 4.83 0.37 1.03 2.73 3.78 5.72 13.33 17.60 24.39 87.26 76.44 64.61 49.86
9.43 8.01 6.71 5.82 4.72 3.81 3.23 2.66 2.20 1.70 1.31 1.10 0.89 0.89 0.67 15.24 12.03 11.85 10.54 8.46 7.34 6.65 5.10 3.99 3.36 2.68 2.58 1.77 1.77 1.31
DAG TAG
3.76 3.33 2.30 2.60 2.36 1.76 0.98 1.14 1.14 0.73 0.88 0.63 0.71 0.71 0.31 600 660 720 780 840 900 960 1020 1080 1140 1200 1260 1320 1380 1440
MAG
Exp. 12
0 5 10 20
0.00 5.44 9.37 18.00
5.85 8.47 8.83 11.31
87.44 69.67 57.30 39.80 6.65 7.37 8.97 10.57 0.00 1.82 6.10 11.39 0 5 10 20
6.27 15.17 21.80 26.06
7.95 7.04 5.67 4.84 4.12 3.51 3.25 2.77 2.36 2.13 1.92 1.66 1.54 1.38 1.27 4.91 4.95 4.65 4.64 4.20 4.07 3.40 4.08 3.71 3.57 3.15 2.81 2.72 3.02 2.71
FFA FAME
72.80 74.40 77.46 79.27 81.50 83.42 85.37 85.96 87.96 89.07 90.50 91.76 92.41 92.54 93.57 600 660 720 780 840 900 960 1020 1080 1140 1200 1260 1320 1380 1440
Time(min)
3.02 3.02 2.94 2.74 2.48 2.35 1.70 1.90 1.66 1.50 1.06 0.84 0.73 0.75 0.53 3.69 3.32 2.91 2.96 3.16 3.05 2.28 2.46 2.13 1.72 2.72 2.11 2.67 2.67 2.14
Results and Discussion
TAG FFA
Uncertainty analysis: Monte–Carlo simulations In this work, the Monte–Carlo technique was used for the evaluation of uncertainty in the kinetic parameters, on the model outputs. This method offers global results due to the large number of model evaluations performed using the randomly sampled parameters, to obtain the distribution of the model outputs. The method samples from the input parameter space, h and generates model outputs, y. The Monte– Carlo analysis of uncertainty involved three steps: Step 1: Specifying parameter input uncertainty: The confidence intervals of the parameters estimates, h^ were used as the upper and lower bounds of a given kinetic parameter. Step 2: Sampling input uncertainty: Latin hypercube sampling with correlation control was used to sample within the input parameter space given the dense stratification over the range of each sampled variable.37 Samples were selected from the input parameter space, where each sample, hi contained one value for each input parameter creating a [m 3 Nlh] matrix. Where m stands for the total number of model parameters, and Nlh is the total number of Latin-Hypercube samples. Step 3: Simulating the model using the sampling matrix: Nlh dynamic simulations were then performed using the [m 3 Nlh] sampled input matrix. Each simulation result was stored in a [Nt 3 Nu 3 Nlh] size array where, Nt is the length of the discrete time series and Nu is the number of model outputs. The complete Monte–Carlo results provide a cumulative distribution function for each output variable at each time instant. The uncertainty of the model outputs were then represented using mean and percentile calculations.
Parameter estimates and confidence intervals
0.89 0.74 1.11 2.21 4.95 10.92 14.87 21.51 84.38 76.76 70.71 59.93 0.00 1.47 2.58 4.09 0 5 10 20 Exp. 10
9.78 10.12 10.73 12.26
0.19 5.77 11.21 14.20 16.90 18.22 18.52 19.55 17.66 16.45 3.12 2.30 1.72 8.26 27.09 31.07 30.66 27.57 24.62 23.94 20.87 19.18 18.10 4.67 3.50 2.89
DAG TAG
84.61 44.78 25.19 15.48 10.79 8.11 6.26 4.42 4.35 3.92 1.05 0.86 0.71 0.00 8.14 17.38 25.58 32.17 37.60 41.33 47.10 51.00 54.19 87.01 89.02 90.36 0 30 60 90 120 150 180 240 300 360 1200 1320 1440
FFA FAME Time(min)
Exp. 9
6.93 14.23 15.16 14.07 12.58 11.44 9.95 8.05 7.81 7.34 4.15 4.32 4.31
MAG
Exp. 11
Time(min)
Mass % Mass % TABLE 3. Continued
parameters in k. Brun et al. determined an empirical threshold of kk being below 15, which is used in this study.35
67.88 73.31 76.22 78.07 81.30 84.03 86.86 88.65 90.55 92.49 92.41 93.57 93.96 93.96 95.56
DAG
7
FAME
Mass %
11.32 10.59 9.27 8.51 7.69 6.65 6.28 5.29 4.31 3.74 3.38 2.93 2.60 2.32 1.93
MAG
Biotechnol. Prog., 2014, Vol. 00, No. 00
The histogram of the residuals for the fitting of Exp. 1–8 were used to assess the quality of the model fitting. For the proposed model, it can be seen in Figure 3 that the histogram is not skewed which gives an indication that the complexity and choice of model is appropriate. Also the residuals have a mean of 20.05 mass % (approximately zero mean) and standard deviation of 2.64 mass %. This signifies that 95% of the residuals are within 20.05 6 5.28 mass %. This result is reasonable, given that a mass balance on the acyl groups for the experimental data close within 3 mass %. The parameter estimates for the proposed model are shown in Table 4 along with the confidence intervals and correlation matrix. Generally, the narrower the confidence interval, the higher the quality of the parameter estimate. The confidence intervals for the kinetic parameters k22, k3, and k26, (the reverse kinetic constants for formation of the TAG enzyme substrate complex (E.T), the rate of DAG production from E.T and the reverse kinetic constants for formation of the MAG enzyme substrate complex (E.M), respectively) deviate more than 40% from the mean estimates signifying a low sensitivity of the model outputs to those parameters. This may be due to the data set not having sufficient information given that the intermediate enzyme substrate complexes were not measured.
1.00 20.26 1.00 20.53 0.61 1.00 0.41 0.01 20.51 1.00 20.33 0.17 0.69 20.07 1.00 20.64 0.08 0.28 20.45 0.04 1.00 0.69 20.03 20.43 0.50 20.25 20.81 1.00 1.00 0.19 20.31 20.23 20.06 20.16 0.07 0.00 1.00 20.01 20.34 0.11 0.12 20.26 20.03 0.36 20.24 1.00 0.18 20.19 20.61 0.18 0.72 20.36 0.78 0.46 20.65 1.00 20.17 0.05 20.16 20.01 20.25 20.35 0.10 20.33 0.08 0.04 1.00 0.02 0.29 0.18 0.17 20.27 20.02 0.30 20.44 0.14 0.36 20.35 1.00 20.04 20.10 20.55 20.23 0.19 0.51 20.13 20.57 0.49 20.48 20.33 0.38 1.00 20.03 20.20 0.50 0.26 20.11 20.26 20.10 20.21 0.09 20.06 0.18 0.01 20.13 1.00 20.03 0.46 20.46 0.14 20.52 20.30 0.06 0.44 20.18 20.61 0.70 20.29 20.48 0.54 1.00 20.20 20.06 20.21 0.25 20.17 0.20 0.12 0.09 20.18 0.07 0.31 20.45 0.11 0.27 20.26 1.00 0.51 20.77 0.18 20.32 0.40 0.01 0.47 0.20 0.02 20.30 0.00 0.42 20.58 0.18 0.40 20.49 1.00 0.47 0.16 20.37 0.46 20.52 0.24 20.08 0.80 20.04 20.08 20.40 20.09 0.63 20.40 0.72 0.15 20.44 1.00 20.32 20.07 20.15 20.07 20.32 0.01 20.01 20.02 20.15 0.19 0.01 20.17 20.03 20.22 0.02 20.25 0.18 0.00
k10 k29 k9 k28 k8 k27 k7 k26 k6 k25 k5 k24 k4 k23 k3 k22 k2 k21 k1
1.00 0.60 20.81 20.38 20.13 20.21 20.40 0.30 20.12 0.08 20.70 0.14 0.10 0.21 20.12 20.54 0.04 20.68 0.03 0.26 – 冑 – 冑 – 2 冑 – – 冑 – – – 冑 冑 冑 冑 冑 冑 104 106 106 104 107 107 106 107 107 105 105 105 106 100 104 106 104 105 102 1024 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
Upper
5.88 7.78 1.96 1.63 3.32 2.51 4.15 1.59 1.79 2.21 1.01 8.41 9.41 5.86 2.84 4.84 2.85 2.44 4.08 5.08 104 100 106 103 107 107 106 107 107 105 104 105 106 100 104 106 104 105 102 1024
Lower
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
Correlation Matrix
Identifiable parameter Subset Confidence Intervals
3.30 4.26 1.45 8.78 1.53 1.75 2.68 1.06 1.36 1.35 8.31 4.14 4.84 4.45 2.00 2.66 2.08 1.66 2.69 3.05 104 100 106 104 107 106 106 107 107 105 104 107 106 100 104 106 104 105 102 1024 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
Having generated parameter estimates, key to the analysis is to find out how specific the parameter estimates actually are, given the experimental data used. One simple method is to look at the correlation between the parameters (Table 4). In this study 60.75 was used to signify highly correlated parameters. This value is chosen based on the correlation value of the inhibition constant (20.81) given the inhibition
Parameter
Correlated and identifiable parameters
Mean Estimate
For Exp. 9–12 the histogram of the residuals (Figure 4) for the validation data set is slightly skewed to the left indicating the model tends to underestimate the concentrations. The performance of the parameter estimates over the entire time course of the reaction is illustrated in Figure 5, using validation data set Exp. 9. The proposed model captures the dynamics for the five components over the entire course of the reaction, although the prediction for FFA and MAG shows some deviation from the experimental data. We also investigated the performance of the model in the prediction of the concentrations of TAG, DAG, MAG, FAME, FFA, and CH for the first 20 min of the reaction (Figure 6), for various enzyme and water concentrations (Note: the water concentrations for Exp. 11 and 12 are outside the range used for the model fitting but are included to investigate the extrapolability of the mechanistic model). The model follows the expected trends of the experimental data. For example, the FAME production increases, with increasing water concentration. The available interfacial area is larger, which means there is an increased chance for substrates to react and hence an increase in the rate of FAME production.38 Likewise, as the enzyme concentration increases from 0.2 to 0.5 wt % enzyme, the production of FAME increases as expected. The same reasoning extends for the other components plotted. However, at the higher water concentrations, the model tends to over predict the amount of FFA produced, reducing the FAME production. The model mismatch observed may be due to process phenomena not taken into account. For example, the viscosity of the reaction media changes one order of magnitude over the 24 h.39 Hence the parameter estimates are average values of the rate constants over the entire course of the reaction. Likewise, the uncertainty in the parameter estimates plays a part given some of the parameter estimates are strongly correlated.
Table 4. Parameter Estimates, Related Confidence Intervals and Correlation of the Parameters for 10000 Bootstrap Samples. Highlighted Values Represent Highly Correlated Values
Figure 3. Histogram of residuals for the fitting of the proposed model for Exp. 1–8. The distribution has a mean of 20.05 mass % and a standard deviation of 2.64 mass %.
4.95 6.60 1.69 1.11 2.07 2.20 3.41 1.33 1.55 1.81 9.13 5.43 7.06 4.93 2.36 3.51 2.54 2.05 3.23 4.39
k210
Biotechnol. Prog., 2014, Vol. 00, No. 00
k1 [L/mol.s] k21 [1/s] k2 [L/mol.s] k22 [1/s] k3 [1/s] k23 [L/mol s] k4 [L/mol s] k24 [1/s] k5 [1/s] k25 [L/mol s] k6 [L/mol s] k26 [1/s] k7 [L/mol s] k27 [L/mol s] k8 [L/mol s] k28 [L/mol s] k9 [L/mol s] k29 [L/mol s] k10 [L/mol s] k210 [1/s]
8
Biotechnol. Prog., 2014, Vol. 00, No. 00
constants are usually strongly correlated.40 For two highly correlated parameters, the change in model output due to changing one of the correlated parameters, can be compensated for by an appropriate change of the other parameter value, preventing a unique estimate of the parameter value. This can be due to the model structure or the similarity of the parameters of the underlying biological system.41 What this may signify is that enzyme activities can be modified by changing one of the correlated parameters. For example k10 is highly correlated with k210. The inhibition constant has a negative coefficient for the correlation value. If one parameter value increases, the other decreases. Hence a modification the enzyme structure that affects inhibition can potentially have a twofold effect. If the forward rate where the dead end complex E.CH is reduced, then the rate of disassociation of E.CH will increase. The correlation seen amongst some of
9
the parameters should be expected given the complex parallel and sequential reactions occurring simultaneously. The plot of the Collinearity index in Figure 7 shows that the Collinearity index increases with the number of parameters and that only a maximum of 10 of the 20 parameters can be identified with the available experimental data. In Table 4 one potential subset was identified, that takes into account the parameters that are correlated. The ticked (冑) parameters were the ones estimated and the others were fixed. Hence, the procedure is iterative, although in this case it only gives a reduction in the squared-sum of the relative errors between the simulated and experimental of 0.01%. It should be noted that fixing parameters, while estimating others, results in reasonable parameter values rather than “true parameter values.” 35 Given the uncertainty in the parameter estimates, we then look at how the model can be used for engineering purposes. Uncertainty analysis: Monte–Carlo simulations
Figure 4.
Histogram of residuals for the fitting of the proposed model for Exp. 9–12. The distribution has a mean of 20.005 mass % and a standard deviation of 2.33 mass %.
The uncertainty in the model outputs for the typically measured variables (TAG, DAG, MAG, FAME, and FFA) during the transesterification reaction are seen in the spaghetti plots of the 500 Latin hypercube samples in Figure 8. When the experimental values are overlaid on the Monte– Carlo simulations, we can see that over the complete course of the reaction, most of the experimental values fall within the bounds of the spaghetti plots. The narrow prediction bands for FAME and TAG reflect the robustness of the predictions for those model outputs over the entire course of the reaction, while the wide bands observed for FFA and MAG show the need for a more accurate estimate of the parameters in order to obtain more certain model predictions. Using the cumulative frequency distribution plots (Figure 9) it is possible to put bounds on the model predictions, which can give the modeler some insight into the reliability of the model to make predictions. Take for example the FAME predictions. At the end of the reaction the model output has a mean value of 87.9 mass % with a standard deviation of 0.64 mass %.
C O L O R Figure 5. Comparison of the model fitting (2) to the Validation data set Exp 9 (*).
10
Figure 6.
Biotechnol. Prog., 2014, Vol. 00, No. 00
Comparison of the model fitting for the first 20 min of the reaction for experiments 10 (* experimental, - simulated), 11 ( experimental,- simulated) and 12 (䉮 experimental, . . . simulated).
Given the uncertainty in the model parameters, the model gives excellent predictions of the FAME and TAG values and shows the deficiencies in the FFA and DAG predictions during the course of the reaction. We hence see this method as a valuable tool to gauge the robustness of a model to parameter uncertainty. Engineering application of the model given the parameter uncertainty For process development, a reliable kinetic model can be used, for example, to simulate and evaluate variations in feed composition, alternative reactor configurations and feeding strategies to mitigate methanol inhibition to name a few. The results from the simulations can then be tested experimentally. Below we investigate a methanol feeding strategy to mitigate inhibition and the uncertainty in the model outputs due to the parameter estimates. Control of Methanol Feeding. One hurdle to industrial implementation of enzymatic biodiesel production is inhibition and deactivation of the biocatalyst by the alcohol substrate. Simulations can be used to devise an optimal feeding policy. We followed a methanol feeding strategy similar to the one proposed by Samukawa et al.42 They found that they could increase the reuse of the immobilized enzyme (a clear indication of a reduction in enzyme deactivation), by keeping the methanol content in the reactor below the concentration that gave the highest initial rate of FAME production (CHcriti42 In our previous work we found CHcritical to be 0.525 cal). 43 Eq. We then combined the model of the system (Eq. 4) with the objective function in Eq. 9. By minimizing the objective function in Eq. 9 we ensure that the methanol concentration in the reactor never goes above the critical value CHcritical at each time step ti by manipulating the methanol feed Fa. To simplify the experimental procedure only two step changes in the methanol feed rate were used. min JEq 5ðfCHEq gti 2fCHcritical gti Þ2
(9)
Fa
The optimization results are now compared to the experimental results as illustrated in Figure 10. The optimization
Figure 7. Plot of how the collinearity index varies with all parameters (Top plot) and how the collinearity index varies with the identifiable parameters for a threshold value of 15 (bottom plot).
objective to constrain the amount of methanol in the reactor was validated experimentally (Exp. 13) and compared to the Exp. 9 where the methanol concentration is over 0.6 Eq. for the last 7 h of the reaction (Figure 8). The optimised case respects the constraint, with the methanol concentration never going above the critical value of 0.525 Eq. for the entire reaction. However, for the optimized case, there was a 2% reduction in the final biodiesel concentration compared to Exp. 8, which had the highest FAME conversion. The FAME equilibrium concentrations at the end of the reaction can be increased by increasing the methanol concentration at the end of the reaction. However, we are then exposing the enzyme to higher concentrations of methanol, which potentially reduces the number of times the enzyme could be reused. What is interesting is the trade-off between downstream processing to bring the final biodiesel concentration within specifications and potential increase in enzyme reuse. However, this analysis is beyond the scope of this article.
Biotechnol. Prog., 2014, Vol. 00, No. 00
11
Figure 8. Uncertainty analysis for the validation experiment (Exp 9). The experimental values (*) are overlaid on the 500 Monte– Carlo simulations (2).
Figure 9. Cumulative distribution function for the 500 Monte–Carlo simulations (Exp 9) at time 24 h.
Characterization of the Model Uncertainty for the Optimization. The same parameter input uncertainty used for the validation data set was used for the uncertainty analysis given we are still operating within the calibrated range of the model. In Figure 10 the narrow prediction bands for TAG, DAG, and FAME reflects the robustness of the predictions. When the experimental values are overlaid on the 500 Monte–Carlo simulations, we can see over the course of the reaction that most of the model outputs fall within the bounds of the spaghetti plots except most notably for the FFA prediction. The concentrations for the FFA predictions are on the same order of magnitude, although the dynamics after 5 h for the FFA simulation show a slight increase followed by a decrease in FFA concentration compared to the experimental values that show a steady decrease. As postu-
lated earlier, given the decrease in viscosity of the reaction system as time progresses, it is believed that the rate of FFA consumption increases during the reaction and hence the steady decrease in the FFA concentration seen in the experimental data. Given the rate constants are average values in the simulation, it is not possible to capture the behaviour seen. The cumulative frequency distribution plots (Figure 11) are then used to characterize the uncertainty in the model outputs. In this case, at the end of the reaction, the FAME model output has a mean value of 90.83 mass % with a standard deviation of 0.55 mass %. Limitations and Use of the Model. The developed model is only applicable for relatively low methanol concentrations
12
Biotechnol. Prog., 2014, Vol. 00, No. 00
Figure 10. Uncertainty analysis for Exp 13. The experimental values (*) are overlaid on the 500 Monte–Carlo simulations (-).
Figure 11. Cumulative distribution function for the 500 Monte–Carlo simulations (Exp 13) at time 24 h.
(
