Kinetic study of lipase-catalyzed glycerolysis of palm olein using Lipozyme TLIM in solvent-free system

Diacylglycerol (DAG) and monoacylglycerol (MAG) are two natural occurring minor components found in most edible fats and oils. These compounds have gained increasing market demand owing to their unique physicochemical properties. Enzymatic glycerolysis in solvent-free system might be a promising approach in producing DAG and MAG-enriched oil. Understanding on glycerolysis mechanism is therefore of great importance for process simulation and optimization. In this study, a commercial immobilized lipase (Lipozyme TL IM) was used to catalyze the glycerolysis reaction. The kinetics of enzymatic glycerolysis reaction between triacylglycerol (TAG) and glycerol (G) were modeled using rate equation with unsteady-state assumption. Ternary complex, ping-pong bi-bi and complex ping-pong bi-bi models were proposed and compared in this study. The reaction rate constants were determined using non-linear regression and sum of square errors (SSE) were minimized. Present work revealed satisfactory agreement between experimental data and the result generated by complex ping-pong bi-bi model as compared to other models. The proposed kinetic model would facilitate understanding on enzymatic glycerolysis for DAG and MAG production and design optimization of a pilot-scale reactor.


Introduction
The prevalence of overweight and obesity owing to the imbalance between energy intake and energy expenditure has triggered attention to the exploitation of an effective approach such as the development of structured-lipid, starch or protein-based fat replacer as well as food supplements and medication, to reverse the fasting growing obesity trend. Recent disclosure of diacylglycerol (DAG)-enriched oil has therefore become a subject of growing interest worldwide a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 due to its anti-obesity properties. Previous research findings revealed that DAG possesses the ability to reduce serum triacylglycerol (TAG), enhance β-oxidation activity and prevent abnormal accumulation of fat [1][2][3]. In addition, DAG in synergistic with monoacylglycerol (MAG) can act as excellent emulsifiers and stabilizers in food and pharmaceutical industry [4]. MAG alone has also been widely employed in textile processing, production of plastics and formulation of oil for machinery on account of its excellent lubricant and plasticizing properties. Therefore, the demand for both DAG and MAG has skyrocketed in recent years [5][6][7].
Mixture of DAG and MAG can be produced chemically. However, the major barrier of the chemical process is extreme operating condition required to accelerate the overall reaction rate with the use of an inorganic catalyst. This eventually results in the formation of undesired by-products (oxidized compounds and soap) and purification of products is therefore needed [8]. The breakthrough in enzyme technology in fats and oils industry for the modification of edible oil has drawn attentions worldwide. Enzymatic approach provides an alternative option due to its mild processing conditions, high regioselectivity of lipases and less environmental impact [9,10]. Various enzyme-mediated methods such as hydrolysis of TAG [11][12][13], glycerolysis of TAG [14][15][16][17], esterification of fatty acids and glycerol [11,15] have been reported for DAG and MAG manufacturing. Among these enzymatic methods, glycerolysis of TAG had been proven to be efficient for the preparation of DAG and MAG due to the huge amounts of glycerol surplus from biodiesel industries as cheap reactants and high production yield [18][19][20].
Earlier we have analyzed enzymatic glycerolysis for DAG production, employing different commercial enzymes namely Lipozyme RMIM, Lipozyme TLIM and Novozyme 435. Lipozyme TLIM was found to surpass others on cost performance basis as it is cheaper and it showed greater affinity towards DAG synthesis which translates to high feasibility of commercialization [17]. Kinetic evaluation and simulation of biochemical process is important in the design of scale-up experiment in order to gain an insight into reaction mechanism of enzymatic glycerolysis. Ping-pong bi-bi and ordered-sequential bi-bi mechanisms are the common mathematical models in the literatures used to represent the lipase-catalyzed reaction regardless of the type of bioconversion [13,[20][21][22][23][24][25][26][27][28]. Nevertheless, to the best of our knowledge, kinetic evaluation of lipase-catalyzed glycerolysis of palm olein in solvent-free system remains vague Therefore, in this paper, three kinetic models namely simple ternary complex, simple ping-pong bi-bi and complex ping-pong bi-bi were proposed based on unsteady-state assumption to simulate the overall glycerolysis process and comparison among the kinetic mechanisms was emphasized.

Glycerolysis reaction
Enzymatic glycerolysis reaction was conducted in solvent-free system as described in earlier work [17]. The reactions were performed by reacting 100 g palm olein in the presence of enzyme at glycerol to enzyme (G/E) mass ratio of 1 in a 250 ml conical flask. Three enzyme loads (3, 5 and 8 wt % oil) were studied. The reaction mixture was incubated at 55˚C and shaken at 200 rpm in a water bath shaker. Samples were taken at specific time intervals where each sample was stored at -20˚C prior to analysis. All reactions were carried out in duplicate.

Reversed-Phase HPLC analysis
Acylglycerol compositions were determined by high-performance liquid chromatography (HPLC) (Gilson, France) equipped with a refractive index detector model 2410 (Waters, USA), using two commercially packed LiChrospher1100 RP-18 column end capped (250mm) with 5 μm particle size in series. The glyceride compositions were eluted with acetone/acetonitrile (70:30) for TAG determination and acetone/acetonitrile (50:50) for free fatty acid (FFA), MAG and DAG determination, both at a flow rate of 1 ml/min. The TAG and DAG were identified according to Swe [29] and Ghazali [30]. Calibration curves were constructed, and the results were given as the weight ratio of total glycerides.

Kinetic modeling
In present work, a time course analysis was conducted to determine the rate constants for each reaction. The rate constants of forward and backward reaction were determined individually so that the effect of each step on lipase-catalyzed glycerolysis process could be better understood.
Three models were considered in this study: a. Ternary complex model; b. Simple ping-pong bi-bi model; and c. Complex ping-pong bi-bi model A material balance was performed for each component involved in the model. Two assumptions were made while performing the material balance: a. Mass transfer limitation in the reaction system was negligible (both internal and external mass transfer resistance can be ignored when considering Lipozyme IM with porous support and sufficient stirring or shaking speed and; b. All the reactions were elementary.
c. No loss of enzyme activity throughout the reaction d. No significant amount of free fatty acid is formed as the glycerolysis was carried out in micro-aqueous environment.
(A) Simple ternary complex model. Referring to Eqs 1-3, simple ternary model proposed in this study suggested that enzymatic glycerolysis reaction follows ordered-sequential bi-bi mechanism in which substrate TAG was initially bound to the enzyme (E) to form a binary complex (E.TAG). Glycerol (G) was then attached to the binary complex before any reaction takes place to form intermediate product of ternary complex (E.G.TAG). Instantaneously, products DAG and MAG were formed and released with restoration of enzyme particles.
where k 1 is the forward reaction rate for TAG and E to form binary complex (E.TAG); k 2 is the reverse reaction of E.TAG complex to cleave into TAG and E; k 3 indicates the forward reaction rate for E.TAG complex to form ternary complex (E.G.TAG); k 4 represents the reverse process rate for ternary complex to reform binary complex; k 5 is forward reaction for DAG and MAG formation.
Considering the steps described in Eqs 1-3 as elementary reactions, the rate equation for the substrates, products and complexes were given as: The sum of all the intermediate enzyme complexes is equal to the enzyme loading as shown in the following equation: (B) Simple ping-pong bi-bi model. Simple ping-pong bi-bi model was also adopted to simulate the kinetic mechanism of enzymatic glycerolysis (Eq 12-15). Binding of the substrate TAG to enzyme (E) to form a binary complex (E.TAG) is the initial step of the reaction, followed by the formation of DAG and a modified enzyme complex (E.FFA). Glycerol (G) is then bound to the modified enzyme complex (E.FFA) forming another binary complex (E.FFA).G.
Finally, MAG is released and the enzyme is regenerated.
where k 1 and k 2 relate to forward reaction coefficient for binary complex (E.TAG) formation and reverse rate constant for TAG and E reformation, respectively; k 3 is the forward reaction of E.TAG to form main product DAG coinciding with the formation of modified enzyme complex (E.FFA); k 4 indicates the reaction between E.FFA complex and G to form secondary binary complex (E.FFA).G; k 5 shows the reverse process which will form E.FFA and G; k 6 is forward reaction coefficient for MAG formation. The entire rate equation for the substrates, products, and complexes were represented as: (C) Complex ping-pong bi-bi model. The complex ping-pong bi-bi model is an extension of the simple ping-pong bi-bi model where this mechanism suggests that some of the diacylglycerol produced become a substrate to the enzyme (E) forming a binary complex (E. DAG), followed by the release of monoacylglycerol (MAG) and a modified enzyme complex (E.FFA) as shown.
where k 1 to k 6 are constants for the reactions as described in simple ping-pong bi-bi model; k 7 and k 8 are both the forward and reverse reactions of enzyme complex E.DAG; k 9 indicates the formation of MAG and E.FFA complex which will react with G for the production of MAG. The rate equation for the substrates, products, and complexes were given as:

Parameter estimation
The rate constant parameters were determined by fitting the model equations into experimental data using software Matlab 7.1. Non-linear regression was used to determine the rate constant parameters by minimizing the sum of square errors (SSE) (Eq 41). Root-mean-square deviation (RMSD) and Chi-squared test are also crucial tools used to evaluate the goodness-offit between both experimental results and model-generated data (Eqs 42 and 43). High-order Runge-Kutta 4,5 was used to solve the ordinary differential equation (ODE). The rate constants obtained from the non-linear regression may be sensitive to initial guesses. At least three different initial guesses were used to ensure the consistency of the rate constants.

SSE
where X calc,j,i is the calculated weight fraction of species i at j number of data; X exp,j,i is the experimental weight fraction of species i at j number of data, and N and Z are the number of experimental data and number of constants used in the proposed kinetic models, respectively.

Results and discussion
The rate constants were initially estimated using experimental data of 3 wt-% enzyme load with assumption being made on the independence of the reaction coefficient on enzyme load. The values of the rate constants and statistical evaluation of kinetic models with respect to experimental data are tabulated in Tables 1 and 2, respectively. Statistical analysis in kinetic modeling is a powerful tool for model discrimination as it enables elucidation of the reliability and accuracy of the estimating equations for distinct kinetic models. As depicted in Table 2, SSE, RMSD and chi-square χ 2 for both 3 wt % and 5 wt % enzyme loads as tested with the three models proposed were found to approach zero, indicating tight fit of the models proposed to the data. Nevertheless, the three models showed relatively high SSE and RMSD as enzyme concentration increased to 8 wt % because excessive enzyme may result in enzyme aggregation, rendering the exposure of the enzyme active site to substrate molecules [20] and therefore improvement and re-estimation of parameters is required. Figs 1-3 illustrate the comparisons between the simulated results and experimental data for enzyme load of 3 wt %, 5 wt %, and 8 wt %, respectively. Even though all models accurately predict the TAG composition for both 3 wt % and 5 wt % enzyme load, complex ping-pong bibi model was found to be able to provide a better estimation on the experimental data with regard to DAG and MAG for reactions treated with 3 wt % and 5 wt % enzyme load as opposed to simple ping-pong bi-bi and ternary complex model (Figs 1 and 2). Based on the statistical data, complex ping-pong bi-bi model was considered as the preferable model as it gave the lowest SSE and RMSD for the two enzyme loads tested (Table 2). However, all models (based on the rate constants obtained using 3 wt-% enzyme load) obviously predict fairly well the acylglycerol compositions at 8 wt % enzyme load, as illustrated in Fig 3. The results agreed with Al-Zuhair's [24] hypothesis that low enzyme concentration model was not suitable for simulating process with high enzyme concentration as the model tends to diverge at high enzyme load.
A new set of rate constants were then revised for the 8 wt-% enzyme load as shown in Table 3. Since the active area of the enzyme was less utilized for TAG at high enzyme concentration as compared to low enzyme concentration, the binding rate between TAG and E, k 1 , was significantly reduced from 0.012 to 0.008 hour -1 g -1 enzyme.  Table 2. Statistical evaluation of models proposed for different enzyme loads. The simulation results with the revised rate constants for 8 wt-% enzyme load for the three models were shown in Fig 4. The revised rate constants were now able to predict the results better for the 8 wt-% enzyme load. Again, the simulation results showed that the complex ping-pong bi-bi model performs better, looking at the fit of the kinetic model to the experimental results (Fig 4). Additional evidence that supports the validity of the complex ping-pong bi-bi model is the statistical analysis of the models (Table 3). Both SSE and RMSD were determined to be 0.0090 and 0.1638, respectively for complex ping-pong bi-bi mechanism. Chi- squared test was evaluated to be 3.8 x 10 −3 , indicating good agreement with observed experimental results. Overall, these simulation results suggested that the glycerolysis reaction followed the complex ping-pong bi-bi mechanism where the DAG was first produced prior to  formation of MAG as the products. The DAG formed earlier then become a substrate for the formation of binary complex (E.DAG) and subsequently led to the formation of MAG. Reaction constant or rate coefficient is crucial in chemical kinetics as it reveals the reaction velocity of each mechanism, giving an insight of a particular biochemical process. Based on the estimated rate constants for ping-pong bi-bi model, it can be concluded the formation of primary binary complex (E.TAG) is slow (low k 1 ). It could be hypothesized that high energy barrier decreases the formation progress. Once the unstable primary complex was formed, it dissociates immediately into first product (DAG) and transforms into modified enzyme complex (E.FFA) instead of reforming the substrates (TAG and E) as k 3 is several times higher than k 2 . Instantly, glycerol molecules show high affinity towards the modified enzyme complex and bind to it to become secondary binary complex (E.FFA.G) (high k 4 ). This unsteady complex decomposes into stable MAG product without delay and enzyme particle is restored instead of following the reverse pathway for modified enzyme complex reformation (k 6 >>> k 5 ). The DAG formed earlier was then bound with enzyme to form enzyme complex (E.DAG) at slow rate due to low value of reaction coefficient k 7 . The unstable enzyme complex (E.DAG) would tend to reverse its reaction pathway in which both enzyme and DAG are regenerated instead of proceeding with the forward direction, producing more MAG and modified enzyme complex (E.FFA) (k 8 > k 9 ).

Conclusions
The enzymatic glycerolysis reaction between TAG and G were successfully modeled based on simple ternary complex, simple ping-pong bi-bi and complex ping-pong bi-bi mechanisms. Based on the statistical analysis, complex ping-pong bi-bi model was proven to give a better and satisfactory agreement between experimental data and model results. Two sets of estimated rate constants were developed for the kinetic models, depending on the enzyme concentration. Accessibility of reactant mixtures to the active surface of the immobilized enzyme should be taken into consideration when enzyme loading is high (! 8 wt%) which was indicated by low k 1 value. The present study provides valuable information and better understanding of the kinetic mechanism for lipase-catalyzed glycerolysis, allowing the use of the model for process optimization in near future.