 Research
 Open Access
 Published:
Development of modified HCH1 kinetic model for longterm enzymatic cellulose hydrolysis and comparison with literature models
Biotechnology for Biofuels volume 12, Article number: 34 (2019)
Abstract
Background
Enzymatic hydrolysis is a major step for cellulosic ethanol production. A thorough understanding of enzymatic hydrolysis is necessary to help design optimal conditions and economical systems. The original HCH1 (Holtzapple–Caram–Humphrey–1) model is a generalized mechanistic model for enzymatic cellulose hydrolysis, but was previously applied only to the initial rates. In this study, the original HCH1 model was modified to describe integrated enzymatic cellulose hydrolysis. The relationships between parameters in the HCH1 model and substrate conversion were investigated. Literature models for longterm (> 48 h) enzymatic hydrolysis were summarized and compared to the modified HCH1 model.
Results
A modified HCH1 model was developed for longterm (> 48 h) enzymatic cellulose hydrolysis. This modified HCH1 model includes the following additional considerations: (1) relationships between coefficients and substrate conversion, and (2) enzyme stability. Parameter estimation was performed with 10day experimental data using αcellulose as substrate. The developed model satisfactorily describes integrated cellulose hydrolysis data taken with various reaction conditions (initial substrate concentration, initial product concentration, enzyme loading, time). Mechanistic (and semimechanistic) literature models for longterm enzymatic hydrolysis were compared with the modified HCH1 model and evaluated by the corrected version of the Akaike information criterion. Comparison results show that the modified HCH1 model provides the best fit for enzymatic cellulose hydrolysis.
Conclusions
The HCH1 model was modified to extend its application to integrated enzymatic hydrolysis; it performed well when predicting 10day cellulose hydrolysis at various experimental conditions. Comparison with the literature models showed that the modified HCH1 model provided the best fit.
Background
Biomass is the only renewable energy resource that can be directly converted to liquid fuels and chemicals. The largest biomass feedstock is lignocellulose, which is found globally in many forms. Converting lignocellulose into biofuels could relieve shortages of liquid fuels and reduce dependence on fossil energy.
In the United States, ethanol is the dominant biofuel, which is usually produced from corn, an important food for animals and humans. To prevent food shortages, cellulosic ethanol is an attractive alternative. In general, there are four major steps for cellulosic ethanol production: pretreatment, hydrolysis, fermentation, and separation. Among these processes, hydrolysis accounts for a large portion (~ 30%) of the total costs [1]. To compete with corn ethanol and petroleumderived gasoline, enzymatic hydrolysis needs optimization and cost reduction [2]; therefore, a thorough understanding of enzymatic hydrolysis is necessary to help design optimal conditions and economical systems.
During the last several decades, various theoretical and empirical models have been developed to simulate enzymatic hydrolysis of cellulose [3,4,5,6]. Because they lack a theoretical foundation, empirical models cannot be applied beyond the range of the data; therefore, this paper only focuses on mechanistic (and semimechanistic) models.
In 1984, Holtzapple et al. [3] proposed a generalized mechanistic model for cellulose hydrolysis termed the HCH1 (Holtzapple–Caram–Humphrey1) model. Figure 1 shows the reaction mechanism of the HCH1 model [3]. As shown in the figure, free enzyme (E^{f}) adsorbs onto a free cellulose site (\(G_{x}^{f}\)) to become adsorbed enzyme (E^{a}), and then complexes with the cellulose to become enzyme–substrate complex (EG_{x}). This complex catalyzes the hydrolysis of the cellulose site to obtain soluble product (G_{s}) with reaction rate k. All enzyme species can complex with product to become inhibited enzyme (\(G_{s} E^{f} ,\;\; G_{s} E^{a}\), and G_{s}EG_{x}). For simplicity, the productbinding constant (β) is assumed to be the same for all the enzyme species. In addition, the adsorption constant (δ) and the complexing constant (η) are assumed not to be affected by the binding of product to the enzyme [3].
The ratelimiting step is the hydrolysis; therefore, the reaction velocity (V) is proportional to the concentration of uninhibited enzyme–substrate complex (EG_{x}). To express the reaction velocity in terms of known variables, substitutions can be made for EG_{x} using material balances of substrate and enzyme species, thus, yielding the HCH1 model (Eq. 1). The detailed model development is described in [3].
where G_{x} is the cellulose concentration (g/L, equivalent to glucose), G_{1} is the glucose concentration (g/L), G_{2} is the cellobiose concentration (g/L, equivalent to glucose), E is the enzyme concentration (g/L), α is the lumped adsorption constant (\(\alpha = \frac{\eta \delta }{\eta + 1}\), g/L), κ is the lumped kinetic constant (\(\kappa = \frac{k}{\eta + 1}\), h^{−1}), β_{1} is the glucose binding constant (L/g), β_{2} is the cellobiose binding constant (L/g), ɛ is the number of cellulose sites covered by adsorbed or complexed enzyme (dimensionless), i is the fraction of total enzyme that is active (dimensionless), and \(\varphi\) is the fraction of total cellulose sites which are free (dimensionless).
Unlike the classic Michaelis–Menten model, the HCH1 model includes a parameter ɛ that describes the number of reactive sites covered by the enzymes [3, 7]. Furthermore, the HCH1 model includes noncompetitive inhibition by sugar products.
Unlike empirical models that apply only in the range where the data were taken, the HCH1 model is mechanistic (Fig. 1) and, therefore, has broader applicability. As a mechanistic model, it applies to individual enzyme species; however, it has also been applied successfully to an enzyme cocktail in which the mixture is treated as a single “lumped” enzyme. Using the initialrate data for pretreated biomass hydrolyzed by an enzyme cocktail, Brown et al. [7] compared mechanistic models and showed that the HCH1 model provided the best fit to experimental data.
Previous studies show that, at high degrees of conversion, the hydrolysis rate drops by 2–3 orders of magnitude [8, 9]. The following factors contribute to the decreasing hydrolysis rates: (1) enzyme deactivation, (2) product inhibition, (3) decreased substrate reactivity, (4) decreased substrate accessibility, and (5) decreased synergism between cellulases [10]. In shortterm enzymatic hydrolysis, these factors are not important and, therefore, are usually not incorporated into models that predict initial rates. However, in longterm batch saccharification, the reaction time is usually 3 to 5 days. As the reaction proceeds, the coefficients in shortterm enzymatic hydrolysis models, such as HCH1, may change because of the enumerated factors above. To describe longterm integrated enzymatic hydrolysis, the initialrate models must be modified.
In this study, the original HCH1 model was modified to describe 10day enzymatic cellulose hydrolysis with commercial enzyme cocktail CTec2. The HCH1 mechanism (Fig. 1) applies to individual enzymes in the cocktail; however, modeling each enzyme component is exceedingly complex. Understanding the kinetics of each enzyme component would be useful when optimizing the cocktail; however, this study uses a cocktail with defined components. Our approach is to treat the enzyme cocktail as a single “lumped” enzyme; hence, the resulting “lumped” parameters reflect the collective kinetics of the cocktail, not the individual components. To describe longterm enzymatic hydrolysis, this study investigates the relationships between the “lumped” parameters in the HCH1 model and substrate conversion. The sensitivities of parameters in the modified model were analyzed. In addition, literature models for longterm (> 48 h) enzymatic hydrolysis were summarized and compared to the modified HCH1 model.
Materials and methods
Materials
Substrate
The substrate used for all experiments was αcellulose (SigmaAldrich, C8002). Compositional analysis showed that the substrate contained glucan 78.5% and xylan 14.4% [11].
Enzyme
The enzyme used in this study was Novozymes Cellic^{®} CTec2 (lot# VCPI 0007), a blend of aggressive cellulases with high levels of βglucosidases and hemicellulases that degrade lignocellulose into sugars [12]. The protein concentration was determined to be 294 mg protein/mL with Pierce BCA assay [11]. Before use, the enzyme solution was diluted ten times with deionized (DI) water.
Citrate buffer
To maintain relatively high enzyme activity, citrate buffer (0.1 M) with a pH of 4.8 was used in enzymatic hydrolysis experiments. To prepare the buffer, citric acid monohydrate and trisodium citrate dihydrate were added to DI water.
Antibiotics
To prevent the growth of contaminating microorganisms that could consume produced sugars, an antibiotic cocktail was added. To prepare the solutions, tetracycline powder was dissolved in an aqueous solution of 70% ethanol at 10 g/L and cycloheximide powder was dissolved in DI water at 10 g/L.
Enzymatic hydrolysis
In the enzymatic hydrolysis experiments, desired amounts of αcellulose, glucose, and DI water together with 125 mL citrate buffer, 2 mL tetracycline solution, and 1.5 mL cycloheximide solution were added to a 1L centrifuge bottle in sequence and then preheated. When the mixture reached 50 °C, enzymes were added. Then, the centrifuge bottle (total working volume of 250 mL) was placed in the incubator at 50 °C for 10 days with an axial rotation of 2 rpm. Liquid samples of 0.5 mL were taken periodically and submerged in boiling water for 20 min to deactivate the enzymes (note: the volume of liquid sample is small relative to the total slurry volume, so it is assumed to have a negligible impact on substrate concentration). Then, to determine the glucose concentration, the samples were filtered and analyzed by a highperformance liquid chromatography (HPLC), which was equipped with a pair of deashing guard columns (BioRad MicroGuard deashing cartridges, 30 mm × 4.6 mm) and an HPLC carbohydrate analysis column (BioRad Aminex HPX87P, 300 mm × 7.8 mm).
Selection of hydrolysis conditions
Experiments for model fitness
Based on our previous study [11], 16 enzymatic hydrolysis conditions were tested including four different substrate concentrations (40, 60, 80, and 100 g/L), two different enzyme loadings (2 and 5 mg protein/g of dry biomass (mg/g)), and two different initial glucose concentrations (0 and 33 g/L).
Experiments for model predictions
Three enzymatic hydrolysis conditions—which were different from the conditions used for model fitness—were tested for model predictions: (1) substrate concentration: 40 g/L, enzyme loading: 1 mg/g, initial glucose concentration: 0 g/L; (2) substrate concentration: 70 g/L, enzyme loading: 1 mg/g, initial glucose concentration: 28 g/L; (3) substrate concentration: 100 g/L, enzyme loading: 5 mg/g, initial glucose concentration: 28 g/L.
Enzyme stability
Enzyme stability was measured by quantifying the soluble protein concentration over the course of 20 days. In this experiment, the desired amount of CTec2 was added to the preheated mixture of citrate buffer, DI water, and antibiotic cocktail. The resulting solution was placed in the incubator at 50 °C. Samples of 0.5 mL were taken periodically and then centrifuged at 13,000 rpm for 10 min. The protein concentration of the supernatant was measured by the Pierce BCA method.
Modification of HCH1 model
Simulation of enzyme stability
Wallace et al. [13] reported that unproductive binding with lignin and thermal deactivation may play a significant role in enzyme deactivation. Considering the substrate used in this study is ligninfree, we assume that enzyme deactivation is solely due to thermal deactivation. RosalesCalderon et al. [14] observed that the protein concentration of a mixture of glucanase and βglucosidase dropped 34% after incubating at 50 °C for 4 days. It was hypothesized that the enzyme proteins suffered a structural change at 50 °C, which led to protein aggregation and precipitation. Additives, whose concentration was assumed constant and proportional to the initial enzyme protein concentration, were supposed to be present in the cocktail to stabilize the enzyme protein against aggregation. Equation 2 incorporates the presence of additives and is proposed to model protein stability [14]. The development of Eq. 2 is described in detail by RosalesCalderon et al. [14].
where E is the native enzyme protein concentration (g/L), E_{0} is the initial enzyme protein concentration (g/L), and k_{1} and k_{2} are the rate constants (h^{−1}).
The stability of CTec2 with three different initial concentrations was tested. Equation 2 was used to fit the experimental data.
Modification of HCH1 model
The HCH1 model was modified by the following steps:
 Step 1::

Use an empirical equation (Eq. 3) to fit the experimental data of the 16 enzymatic hydrolysis conditions (“Experiments for model fitness” section) with high accuracy. This smoothed version of the data provides the reaction rates needed to fit the parameters in the modified HCH1 model of enzymatic hydrolysis.
$$\begin{aligned} { \frac{{\left. {{\text{d}}[G_{1} } \right]}}{{{\text{d}}t}} = \frac{{3.7798\left( {\left[ {G_{x}^{0} } \right]  \left[ {G_{1} } \right]} \right)^{0.6410} \left( {\frac{{\left[ {E_{0} } \right]\left( {0.0574\left[ {E_{0} } \right] + 0.4370\exp \left( {  t\left( {0.4370 + 0.0574\left[ {E_{0} } \right]} \right)} \right)} \right)}}{{0.4370 + 0.0574\left[ {E_{0} } \right]}}} \right)^{0.8500} }}{{1 + 0.0247\left[ {G_{1} } \right]^{1.1579} }}} \\, \end{aligned} $$(3)where \(\begin{array}{*{20}c} { G_{x}^{0} } \\ \end{array}\) is the initial cellulose concentration (g/L, equivalent to glucose).
Equation 3 was developed based on the integrated version of Eq. 2 and an empirical model for batch fermentation [15]. Detailed development of this equation is described in Additional file 1. To fit the parameters, Eq. 3 was solved with the ode45 routine in MATLAB and nonlinear optimization was achieved by the fmincon routine. The objective function was the sum of square errors (SSE), which is the sum of the squared difference between experimental data and the predicted value [16]. The optimal set of parameters corresponds to the smallest SSE value. This empirical correlation of the data provided a coefficient of determination R^{2} = 0.994.
 Step 2::

Divide substrate conversion (from 0 to 1) into 50 equal segments. Using Eq. 3, calculate the reaction rate at each conversion and get a 16 × 50 data set.
 Step 3::

Determine product inhibition.
The inhibition parameter i in the HCH1 model was calculated by determining the initial velocities with and without inhibitor (Eq. 4) [17]. To estimate this value, the same enzyme and cellulose concentrations should be used.
$$\begin{array}{*{20}c} { i = \frac{{V_{{{\text{with}}\;{\text{inhibitor}}}} }}{{V_{{{\text{no}}\;{\text{inhibitor}}}} }}} \\ \end{array} = \frac{{\left[ {\frac{{\kappa \left[ {G_{x}^{0} } \right]\left[ E \right]}}{{\alpha + \left[ {G_{x}^{0} } \right] + \varepsilon \left[ E \right]}}} \right]i}}{{\frac{{\kappa \left[ {G_{x}^{0} } \right]\left[ E \right]}}{{\alpha + \left[ {G_{x}^{0} } \right] + \varepsilon \left[ E \right]}}}}$$(4)The inhibition of enzymatic hydrolysis by cellobiose was not considered in this study, because CTec2 contains a high level of βglucosidase that rapidly converts produced cellobiose into glucose; the cellobiose peak was not found in the HPLC results.
For a single inhibitor, the inhibition parameter i is expressed in Eq. 5 and the glucose binding constant β_{1} is calculated with Eq. 6.
$$\begin{array}{*{20}c} { i = \frac{1}{{\left. {1 + \beta_{1} [G_{1} } \right]}}} \\ \end{array}$$(5)$$\begin{array}{*{20}c} { \beta_{1} = \frac{{\left( {1  i} \right)}}{{i\left[ {G_{1} } \right]}}} \\ \end{array}$$(6)  Step 4::

Use the HCH1 model to fit the 16 reaction conditions at each conversion, and determine the bestfit coefficients κ, α, and ɛ.
 Step 5::

Determine the relationship between parameter κ and conversion x.
Figure 2a presents the relationship between parameter κ in the HCH1 model and substrate conversion x. The data were obtained from Steps 1–4. As shown in the figure, κ drops very fast in the beginning and then stabilizes after conversion reaches 0.38. Nidetzky and Steiner [18] assumed that cellulose consists of an easily hydrolysable part and a recalcitrant part. Based on their twosubstrate hypothesis, the authors derived a mathematical model to describe the kinetics of cellulose hydrolysis. According to the simulation results, the obtained rate constant for easily hydrolysable cellulose was much higher than that of recalcitrant cellulose. Using αcellulose as substrate, they determined that the fraction of easily hydrolysable cellulose was 0.3 [18]. Figure 2a can be explained by this hypothesis, but the rate constant for the easily hydrolysable cellulose decreases as conversion increases instead of being constant. In our experiments, the fraction of easily hydrolysable cellulose (0.38) is close to the result in [18].
Equation 7 was developed to describe the relationship between parameter κ and conversion x.
$$\begin{array}{*{20}c} { \kappa = \frac{{k_{3} }}{{\left( {1 + x^{{k_{4} }} } \right)^{{k_{5} }} }} + k_{6} } \\ \end{array} ,$$(7)where x is the substrate conversion, k_{3}, k_{4}, k_{5}, and k_{6} are the parameters.
The conversion x in the denominator is used to describe the negative effect of conversion on the rate constant. The parameter \(k_{6}\) is considered as the rate constant for recalcitrant cellulose. The parameter k_{3} is used to describe the difference in rate constants between the easily hydrolysable cellulose (initial) and recalcitrant cellulose (height of the curve). The parameters \(k_{4} \;{\text{and}}\;k_{5}\) are used to describe the decrease rate of the rate constant (steepness of the curve) for the easily hydrolysable part. To fit the data, the MATLAB curve fitting tool was used and a coefficient of determination R^{2} = 0.998 was acquired. The values of parameters k_{3}, k_{4}, k_{5}, and k_{6} were determined in this step.
 Step 6::

Determine the relationship between parameter ε and conversion x.
Figure 2b shows the relationship between parameter ε in the HCH1 model and conversion x. As shown in the figure, parameter ε has a very narrow range (0–0.5) over the entire conversion and remains almost unchanged (nearly zero) at conversion 0.1–0.95. Therefore, in this study, parameter ε is assumed not to be affected by conversion and its optimal value should be close to zero. Brown and Holtzapple [19] reported that if [\(\begin{array}{*{20}c} { G_{x}^{0} } \\ \end{array}\)]/[E_{0}] > 35, assuming ε = 0 would not introduce considerable error (< 1%) (note: in our study, [\(\begin{array}{*{20}c} { G_{x}^{0} } \\ \end{array}\)]/[E_{0}] ≥ 200). The parameter ε is needed only at high enzyme loadings. In industrialscale saccharification, considering the high cost of enzymes, the enzyme dosage must be very low; therefore, if modeling a commercial process, the value of ε can be set as zero.
 Step 7::

Determine the relationship between parameter α and conversion x.
The parameter α in the original HCH1 model may be expressed by Eq. 8, which is related to enzyme adsorption:
$$\begin{array}{*{20}c} {\alpha = \frac{{\left[ {E^{f} } \right]\left[ {G_{x}^{f} } \right]}}{{\left[ {E^{a} } \right] + \left[ {EG_{x} } \right]}}} \\ \end{array} .$$(8)Kumar and Wyman [20] showed that glucose addition and enzyme dosage can affect the percentage of cellulase adsorption. Therefore, besides the impact of conversion, the effects of glucose and enzyme concentration on the value of α were tested. Using the bestfit coefficients κ and ɛ obtained from Step 4, two optimal α values corresponding to two initial glucose concentrations were determined by fitting the data (eight data at each initial glucose concentration) from Step 2 at each conversion with the HCH1 model (Fig. 3a). Another two optimal α values corresponding to two enzyme concentrations were determined by repeating this procedure at each conversion (Fig. 3b). As shown in Fig. 3, as the reaction proceeds, the value of α increases and then is unchanged when the conversion reaches a certain point. It is obvious that high initial glucose concentration and low enzyme dosage improve the value of α significantly over the entire conversion range. Equation 9 was proposed to describe the relationship among α, conversion x, enzyme concentration E, and glucose concentration G_{1}. As shown in Fig. 3, all four curves show an “S” shape; therefore, the sigmoid function—which has an Sshaped curve—was used. The core structure of Eq. 9 is a sigmoidal function that describes the relationship between parameter α and conversion x (note: a_{2} and a_{3} are the parameters of the sigmoid function). In addition, because of the significant effect of glucose and enzyme concentrations on the value of α, the terms [G_{1}] and [E] were included in the numerator and denominator of the sigmoid function, respectively. The parameter a_{1} was added to describe the weight of terms [G_{1}] and [E].
$$\begin{array}{*{20}c} { \alpha = \frac{{a_{1} \left[ {G_{1} } \right]}}{{\left[ E \right]\left( {1 + { \exp }\left( {  a_{2} x + a_{3} } \right)} \right)}}} \\ \end{array} ,$$(9)where a_{1}, a_{2}, and a_{3} are the parameters.
 Step 8::

Modify HCH1 model.
Summarizing the proposed equations, Eq. 10 is the modified HCH1 model, where k_{1}, k_{2}, k_{3}, k_{4}, k_{5}, k_{6}, a_{1}, a_{2}, a_{3}, ε, and β_{1} are parameters. Estimates for k_{1}, k_{2}, k_{3}, k_{4}, k_{5}, k_{6}, and β_{1} were determined in previous steps. In this step, the optimal values of a_{1}, a_{2}, a_{3}, and ε were determined by simultaneously fitting the experimental data of the 16 enzymatic hydrolysis conditions (Section: Experiments for model fitness) with Eq. 10.
$$\frac{{\left. {d[G_{1} } \right]}}{dt} = \frac{{\kappa \left[ {G_{x} } \right]\left[ E \right]i}}{{\alpha + \varphi \left[ {G_{x} } \right] + \varepsilon \left[ E \right]}},$$where,
$$i = \frac{1}{{\left. {1 + \beta_{1} [G_{1} } \right]}}$$$$\begin{array}{*{20}c} {\varphi = \frac{{\left[ {G_{x} } \right]  \alpha  \varepsilon \left[ E \right] + \sqrt {\left( {\left[ {G_{x} } \right]  \alpha  \varepsilon \left[ E \right]} \right)^{2} + 4\alpha \left[ {G_{x} } \right]} }}{{2\left[ {G_{x} } \right]}}} \\ \end{array}$$$$ \frac{{\left. {d[E} \right]}}{dt} = k_{1} \left[ E \right]  k_{2} \left( {\left[ {E_{0} } \right]  \left[ E \right]} \right)\left[ {E_{0} } \right]$$$$\kappa = \frac{{k_{3} }}{{\left( {1 + x^{{k_{4} }} } \right)^{{k_{5} }} }} + k_{6}$$$$\begin{array}{*{20}c} { \alpha = \frac{{a_{1} \left[ {G_{1} } \right]}}{{\left[ E \right]\left( {1 + { \exp }\left( {  a_{2} x + a_{3} } \right)} \right)}}} \\ \end{array} .$$(10)Integration of the differential equations described in Eq. 10 was performed using the same numerical methods described in Step 1.
Sensitivity analysis
Local sensitivity analysis
Local sensitivity analysis assesses the local impact of variation in input factors on model outputs. To do this analysis, the direct differential method [21] was used by calculating the sensitivity indices (Eq. 11). The sensitivities of parameters k_{1}, k_{2}, and β_{1} were not analyzed, because their values were obtained from independent experiment or calculation.
where \(S_{{p_{j} }}\) is the nondimensional sensitivity index of the jth parameter, y is the glucose concentration (g/L), and p_{j} is the jth parameter in the parameter vector p.
Global sensitivity analysis
Local sensitivity only analyzes the sensitivity of a solution close to the optimal parameter set. In contrast, global sensitivity analysis assesses the sensitivity of the model for the full range of possible parameter values [16]. In addition, global sensitivity indices can evaluate the combined impact of multiple parameters on model output.
To calculate a global sensitivity index, a normally distributed search of parameter space using the Monte Carlo method was performed and subsequent analysis of the variance in the model outputs was used. In this study, two global sensitivity indices were calculated: firstorder index and totaleffect index [16, 22]. The firstorder index measures the effect of the parameter of interest alone on the output variance. The totaleffect index accounts for not only the effect of the parameter of interest, but also interactions between the other parameters and the parameter of interest at any order.
Model evaluation
The modified HCH1 model was compared with the literature models for longterm enzymatic hydrolysis. As described in Step 1, the same differential equation integration method, nonlinear optimization constraint algorithm, and objective function were used. The Akaike information criterion (AIC) was used to evaluate model quality for the experimental data. The corrected version of AIC (Eq. 12) was used, because the number of observations was not large enough:
where N is the number of observations, P is the number of model parameters, and SSE is the sum square error.
Results and discussion
Enzyme deactivation
Figure 4 shows that, after incubating at 50 °C for 20 days, soluble protein concentrations of CTec2 dropped to 74, 77, and 83% of their initial values of 0.15, 0.26, and 0.61 g/L, respectively. This result is consistent with a previous study [14] that shows higher initial protein concentrations favor lower deactivation rates and supports the additive hypothesis. Equation 2 successfully describes the time profiles of CTec2 protein concentration with a coefficient of determination R^{2} = 0.999. The rate constants in Eq. 2 were determined to be k_{1} = 0.0225 h^{−1} and k_{2} = 0.1740 L/(g h). It should be noted that the modified HCH1 model is a “lumped” model, the performance of each enzyme was not modeled individually; therefore, the stability of each component in the enzyme cocktail was not investigated. Equation 2 describes the overall deactivation of the cocktail.
Production inhibition
Table 1 lists the values of glucosebinding constant calculated from various reaction conditions. The eight β_{1} values are very close to each other and have a standard deviation of 6 × 10^{−6} L/g. The mean value 0.0429 L/g is considered to be the “true” β_{1} value and is used for later calculations.
Model validation
Figure 5a shows the experimental data and modified HCH1 model fitting results for enzymatic hydrolysis with 16 reaction conditions (“Experiments for model fitness” section). Table 2 shows the values of the parameters obtained from the previous section. The model simulation provided the coefficient of determination R^{2} = 0.992, which indicates that the modified HCH1 model describes enzymatic hydrolysis of αcellulose very well.
As a comparison, Fig. 5b shows the original HCH1 model fit to the experimental data with 16 reaction conditions (“Experiments for model fitness” section). The value of β_{1} (0.0429 L/g) was obtained from the previous section (product inhibition). The optimal values (α = 2.0776 × 10^{6} g/L, κ = 9.2889 × 10^{5} h^{−1}, and ɛ = 0.9996) were determined (note: because the original HCH1 model was not developed for integrated cellulose hydrolysis, these parameter values are not be meaningful). The model simulation provided the coefficient of determination R^{2} = 0.947. The calculated SSE and AICc are listed in Table 3.
Model predictions
The modified HCH1 model (Eq. 10) was used to predict the experimental results of the three conditions described in “Experiments for model prediction.” The parameter values were obtained from the fitness of the 16 conditions (Table 2). The experimental and predicted results are shown in Fig. 5c. The simulation provided the coefficient of determination R^{2} = 0.991, which indicates that the modified HCH1 model predicts enzymatic hydrolysis of αcellulose with high accuracy.
Sensitivity analysis
To explore the controlling factors in the proposed model at different hydrolysis stages, local and global sensitivity analyses were performed. Figure 6a shows the parameter sensitivity indices from local sensitivity analysis of the modified HCH1 model over the course of 10 days. As shown in the figure, the sensitivity of k_{3} drops to nearly 10% of its initial value at day 10. The sensitivity of k_{4} increases first and reaches up to 0.4 at around day 1, and then slightly decreases from day 2 to day 10. For the parameters about α, the sensitivity of a_{1} (absolute value) increases as the hydrolysis time increases. The sensitivities of a_{2} and a_{3} only change within the first several reaction days, and then are close to zero after day 3. The sensitivity of ε is close to zero during the entire reaction time.
Figure 6b, c shows the global sensitivity analysis results of the modified HCH1 model. According to the figures, the firstorder indices and totaleffect indices of all variables are almost identical at any time, which means that the variance in this model is not related to any interaction between parameters. At the initial stage of hydrolysis, the variance in the model output only depends on k_{3} and k_{6}. Then, the sensitivity index of k_{3} decreases very fast during the first 2 days, whereas k_{4} increases up to 0.6. From day 2 to day 10, the effects of k_{6} and a_{1} on the model increase. The variables a_{2}, a_{3}, and ε do not show significant effects on the variance in model predictions.
According to Fig. 6, the local and global sensitivity analyses of the modified HCH1 model show a similar trend during the entire reaction time. Figure 7 shows the sensitivity indices calculated from both analyses at day 10. The rankings of the eight sensitivity indices from both analyses are almost the same (k_{6} > a_{1} > k_{4} > k_{5} > k_{3} > a_{2} ≈ a_{3} ≈ ε).
The sensitivity analyses not only determine which parameters have the most influence on model results, but also verify the assumption in Step 6 that the parameter ε is not needed at low enzyme loadings. These analyses provide direction for further modification of the HCH1 model to apply it to realworld lignocellulose that contains lignin.
Model comparison
Based on the methodology used, the published mechanistic and semimechanistic models for cellulose and lignocellulose can be broadly divided into two classes: Michaelis–Menten and enzymeadsorption models [10]. The models following Michaelis–Menten kinetics can also be divided into two subclasses: full Michaelis–Menten models (all rate equations follow Michaelis–Menten kinetics, including the steps of cellulose to cellobiose, cellulose to glucose, and cellobiose to glucose) and partial Michaelis–Menten models (only the step of cellobiose to glucose follows these kinetics). Models employing enzyme adsorption typically use Langmuir adsorption isotherms or the help of kinetic equations [10]. Some literature models incorporate both enzyme adsorption and Michaelis–Menten kinetics.
In this study, the published models for longterm enzymatic hydrolysis of cellulose and lignocellulose were fit to the experimental data using the numerical methods described in Step 1. Some models do not consider product inhibition. To make a fair comparison, these models were only fit to experimental conditions with no initial sugar added (0 g/L initial glucose; four substrate concentrations × two enzyme loadings). Some models teased out fine details in the elementary reaction steps and included some variables that were not determined in this study, such as exocellulase concentration and associated enzyme concentration [23,24,25]. These models are not included in this section. Table 3 summarizes the number of observations and parameters, calculated SSE and AICc values, and the methodology used for the published models. According to the table, the modified HCH1 model has the least SSE and AICc values, which indicates that this model provides the best fit for longterm enzymatic hydrolysis of αcellulose.
Conclusion
The original HCH1 model was modified to extend its application to integrated enzymatic hydrolysis; it performed well when fitting 10day cellulose hydrolysis at various experimental conditions. Local and global sensitivity analyses were performed to determine the controlling parameters at different hydrolysis stages. Mechanistic (and semimechanistic) literature models for longterm enzymatic hydrolysis were compared with the modified HCH1 model and evaluated by AICc. Comparison results show that the modified HCH1 model provides the best description of enzymatic cellulose hydrolysis. The “lumped” modified HCH1 model developed in this study has a simpler form and fewer parameters than mechanistic models of each enzyme component. When each enzyme is modeled separately, the kinetics is extremely complex with the potential to overparameterize. For the specific commercial enzyme cocktail used in this study, excellent fits to the data were obtained without the need to model each enzyme component individually.
Abbreviations
 mg/g:

mg protein/g of dry biomass
 AIC:

Akaike information criterion
 SSE:

sum of square errors
 DI water:

deionized water
 M–M:

Michaelis–Menten kinetics
 Ads:

adsorptionbased approach
References
 1.
Wallace R, Ibsen K, McAloon A, Yee W. Feasibility study for colocating and integrating ethanol production plants from corn starch and lignocellulosic feedstocks (revised). Golden: National Renewable Energy Lab; 2005.
 2.
Wooley R, Ruth M, Glassner D, Sheehan J. Process design and costing of bioethanol technology: a tool for determining the status and direction of research and development. Biotechnol Prog. 1999;15:794–803.
 3.
Holtzapple MT, Caram HS, Humphrey AE. The HCH1 model of enzymatic cellulose hydrolysis. Biotechnol Bioeng. 1984;26:775–80.
 4.
Kadam KL, Rydholm EC, McMillan JD. Development and validation of a kinetic model for enzymatic saccharification of lignocellulosic biomass. Biotechnol Prog. 2004;20:698–705.
 5.
Fenila F, Shastri Y. Optimal control of enzymatic hydrolysis of lignocellulosic biomass. Resour Effic Technol. 2016;2:S96–104.
 6.
Zhou J, Wang YH, Chu J, Luo LZ, Zhuang YP, Zhang SL. Optimization of cellulase mixture for efficient hydrolysis of steamexploded corn stover by statistically designed experiments. Bioresour Technol. 2009;100:819–25.
 7.
Brown RF, Agbogbo FK, Holtzapple MT. Comparison of mechanistic models in the initial rate enzymatic hydrolysis of AFEXtreated wheat straw. Biotechnol Biofuels. 2010;3:6.
 8.
Bommarius AS, Katona A, Cheben SE, Patel AS, Ragauskas AJ, Knudson K, Pu Y. Cellulase kinetics as a function of cellulose pretreatment. Metab Eng. 2008;10:370–81.
 9.
Hong J, Ye X, Zhang YH. Quantitative determination of cellulose accessibility to cellulase based on adsorption of a nonhydrolytic fusion protein containing CBM and GFP with its applications. Langmuir. 2007;23:12535–40.
 10.
Bansal P, Hall M, Realff MJ, Lee JH, Bommarius AS. Modeling cellulase kinetics on lignocellulosic substrates. Biotechnol Adv. 2009;27:833–48.
 11.
Zentay AN, Liang C, Lonkar S, Holtzapple MT. Countercurrent enzymatic saccharification of cellulosic biomass. Biomass Bioenergy. 2016;90:122–30.
 12.
Novozymes: cellulosic ethanol—novozymes Cellic^{®} CTec2 and HTec2—enzymes for hydrolysis of lignocellulosic. 2010. http://www.shinshuu.ac.jp/faculty/engineering/chair/chem010/manual/Ctec2.pdf. Accessed 29 Sep 2018.
 13.
Wallace J, Brienzo M, GarcíaAparicio MP, Görgens JF. Lignin enrichment and enzyme deactivation as the root cause of enzymatic hydrolysis slowdown of steam pretreated sugarcane bagasse. New Biotechnol. 2016;33:361–71.
 14.
RosalesCalderon O, Trajano HL, Duff SJ. Stability of commercial glucanase and βglucosidase preparations under hydrolysis conditions. PeerJ. 2014;2:e402.
 15.
Fu Z, Holtzapple MT. Anaerobic mixedculture fermentation of aqueous ammoniatreated sugarcane bagasse in consolidated bioprocessing. Biotechnol Bioeng. 2010;106:216–27.
 16.
Ordoñez MC, Raftery JP, Jaladi T, Chen X, Kao K, Karim MN. Modelling of batch kinetics of aerobic carotenoid production using Saccharomyces cerevisiae. Biochem Eng J. 2016;114:226–36.
 17.
Holtzapple M, Cognata M, Shu Y, Hendrickson C. Inhibition of Trichoderma reesei cellulase by sugars and solvents. Biotechnol Bioeng. 1990;36:275–87.
 18.
Nidetzky B, Steiner W. A new approach for modeling cellulase–cellulose adsorption and the kinetics of the enzymatic hydrolysis of microcrystalline cellulose. Biotechnol Bioeng. 1993;42:469–79.
 19.
Brown RF, Holtzapple MT. Parametric analysis of the errors associated with the Michaelis–Menten equation. Biotechnol Bioeng. 1990;36:1141–50.
 20.
Kumar R, Wyman CE. An improved method to directly estimate cellulase adsorption on biomass solids. Enzyme Microb Technol. 2008;42:426–33.
 21.
Peri S, Karra S, Lee YY, Karim MN. Modeling intrinsic kinetics of enzymatic cellulose hydrolysis. Biotechnol Prog. 2007;23:626–37.
 22.
Sobol IM. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math Comput Simul. 2001;55:271–80.
 23.
Shang BZ, Chang R, Chu JW. Systemslevel modeling with molecular resolution elucidates the ratelimiting mechanisms of cellulose decomposition by cellobiohydrolases. J Biol Chem. 2013;288:29081–9.
 24.
CruysBagger N, Alasepp K, Andersen M, Ottesen J, Borch K, Westh P. Rate of threading a cellulose chain into the binding tunnel of a cellulase. J Phys Chem B. 2016;120(25):5591–600.
 25.
Jeoh T, Cardona MJ, Karuna N, Mudinoor AR, Nill J. Mechanistic kinetic models of enzymatic cellulose hydrolysis—a review. Biotechnol Bioeng. 2017;114:1369–85.
 26.
Drissen RE, Maas RH, Van Der Maarel MJ, Kabel MA, Schols HA, Tramper J, Beeftink HH. A generic model for glucose production from various cellulose sources by a commercial cellulase complex. Biocatal Biotransformation. 2007;25:419–29.
 27.
Fan LT, Lee YH. Kinetic studies of enzymatic hydrolysis of insoluble cellulose: derivation of a mechanistic kinetic model. Biotechnol Bioeng. 1983;25:2707–33.
 28.
Liao W, Liu Y, Wen Z, Frear C, Chen S. Kinetic modeling of enzymatic hydrolysis of cellulose in differently pretreated fibers from dairy manure. Biotechnol Bioeng. 2008;101:441–51.
 29.
Gusakov AV, Sinitsyn AP, Klyosov AA. Kinetics of the enzymatic hydrolysis of cellulose: 1. A mathematical model for a batch reactor process. Enzyme Microb Technol. 1985;7:346–52.
 30.
Philippidis GP, Smith TK, Wyman CE. Study of the enzymatic hydrolysis of cellulose for production of fuel ethanol by the simultaneous saccharification and fermentation process. Biotechnol Bioeng. 1993;41:846–53.
 31.
Shen J, Agblevor FA. Kinetics of enzymatic hydrolysis of steamexploded cotton gin waste. Chem Eng Commun. 2008;195:1107–21.
 32.
Zhang Y, Xu JL, Xu HJ, Yuan ZH, Guo Y. Cellulase deactivation based kinetic modeling of enzymatic hydrolysis of steamexploded wheat straw. Bioresour Technol. 2010;101:8261–6.
 33.
RosalesCalderon O, Trajano HL, Posarac D, Duff SJ. Modeling of oxygen delignified wheat straw enzymatic hydrolysis as a function of hydrolysis time, enzyme concentration, and lignin content. Ind Biotechnol. 2016;12:176–86.
Authors’ contributions
CL, MNK, and MH conceived the study and participated in its design. CL performed all the experiments and drafted the manuscript. All authors were involved in the development of the model, sensitivity analysis, and model evaluation. CL, MNK, and MH revised the manuscript. All authors read and approved the final manuscript.
Acknowledgements
The authors would like to acknowledge Mr. Chiranjivi Botre, Mr. Haoran Wu, Mr. Opeyemi Olokede, Mr. ShenChun Hsu, and Ms. Pallavi Kumari for their advice on this project.
Competing interests
The authors declare that they have no competing interests.
Availability of supporting data
All data generated or analyzed in the present study are included in this article.
Consent for publication
All authors obtained consent for publication.
Ethics approval and consent to participate
Not applicable.
Funding
This work was supported by the Michael O’Connor Chair II Endowment. The open access publishing fees for this article have been covered by the Texas A&M University Open Access to Knowledge Fund (OAKFund), supported by the University Libraries and the Office of the Vice President for Research.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Author information
Additional file
13068_2019_1371_MOESM1_ESM.doc
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Liang, C., Gu, C., Raftery, J. et al. Development of modified HCH1 kinetic model for longterm enzymatic cellulose hydrolysis and comparison with literature models. Biotechnol Biofuels 12, 34 (2019) doi:10.1186/s1306801913715
Received
Accepted
Published
DOI
Keywords
 Enzymatic cellulose hydrolysis
 Simulation
 HCH1
 Model comparison
 Sensitivity analysis