Modeling Tool for Studying the Influence of Operating Conditions on the Enzymatic Hydrolysis of Milk Proteins

Systematic modeling of the enzymatic hydrolysis of milk proteins is needed to assist the study and production of partially hydrolyzed milk. The enzymatic hydrolysis of milk proteins was characterized and evaluated as a function of the temperature and protease concentration using Alcalase, Neutrase and Protamex. Modeling was based on the combination of two empirical models formed by a logarithmic and a polynomial equation to correlate the kinetic constants and the operating conditions. The logarithmic equation fitted with high accuracy to the experimental hydrolysis curves with the three proteases (R2 > 0.99). The kinetic constants were correlated with the operating conditions (R2 > 0.97) using polynomial equations. The temperature and protease concentration significantly affected the initial rate of hydrolysis, i.e., the kinetic constant a, while the kinetic constant b was not significantly affected. The values for the kinetic constant a were predicted according to the operating conditions and they were strongly correlated with the experimental data (R2 = 0.95). The model allowed for a high-quality prediction of the hydrolysis curves of milk proteins. This modeling tool can be used in future research to test the correlation between the degree of hydrolysis and the functional properties of milk hydrolysates.


Introduction
The use of hydrolyzed milk proteins is a very attractive practice directed to children's food not only for the hypoallergenicity of this product, but also for the immunomodulatory action of enzyme-generated peptides [1]. The hydrolyzed proteins are obtained by the proteolytic action of proteases, which reduces the molecular weight of proteins hydrolyzing the peptide bonds, thus decreasing the allergenicity caused by epitopes [2,3]. An increase in the degree of hydrolysis and the levels of low molecular weight peptides are generally correlated with the reduced antigenicity of milk proteins [4]. Wroblewska, et al. [5] found that an increase in the degree of hydrolysis (DH) from 15.3% to 15.9% after a two-step process with Alcalase and papain showed a more effective decrease in immunoreactivity despite the nonsignificant increase in the DH. This is a key control parameter during the production of partially hydrolyzed formulas, which could be preferred over extensively hydrolyzed formulas because of their better taste and nutritional value [2]. Partially hydrolyzed formulas are more orally tolerated and less expensive than extensively hydrolyzed formulas [6]. The DH is a useful parameter to characterize the progression of the proteolytic reaction and it is influenced by the operating conditions, traditionally pH, temperature, protease dose and protease-to-substrate ratio [7]. As mentioned above, control of the DH is a major concern during the formulation of hypoallergenic milk. The correlation of DH with hypoallergenicity is a common strategy to determine the adequate DH to be achieved. In addition, the DH achieved depends on the operating conditions and the protease used. The characterization of the hydrolysis process regarding the operating conditions must be evaluated before allergenicity tests. Theoretical modeling by the Michaelis-Menten equation has been proposed by various authors with different protein substrates and proteases [8][9][10][11][12][13][14][15]. The logarithmic equation obtained by Márquez-Moreno and Fernández-Cuadrado [16] has been used by different authors [17][18][19][20][21][22][23][24] due to its simplicity and successful fitting to the experimental data. Empirical modeling based on response surface methodology (RSM) has been applied to protein hydrolysis [25][26][27][28][29]. In a recent publication, a two-level structured model was used as a new methodology to characterize the enzymatic hydrolysis of proteins [30]. The advantages of this methodology are (i) a good description of the hydrolysis process, (ii) an operating definition of kinetic constants, (iii) excellent fitting to experimental data, (iv) the simultaneous exploration of multiple variables (operating conditions) and (v) very good predictability.
The objective of this work was to characterize the enzymatic hydrolysis of milk proteins regarding the operating conditions with different proteases using this recently published modeling methodology. The proposed methodology can be used to correlate the DH with immunoreactivity reduction during the hydrolysis of milk proteins. The information obtained will be useful to produce hypoallergenic milk by designing and setting the operating conditions, thus achieving the proper DH during the process.

Materials
Milk powder from Hormel Foods USA was used. The milk contained 33.7% (w/w) protein and 5.3% (w/w) moisture. Commercial proteases Alcalase 2.4 L (serine endoprotease from Bacillus licheniformis), Neutrase (metallo-endoprotease from Bacillus amyloliquefaciens) and Protamex (a protease with endoprotease and exopeptidase activity from Bacillus licheniformis and Bacillus amyloliquefaciens) supplied by Novozymes (Bagsvaerd, Denmark) were used. Analytical grade quality reagents were used in all experiments.

Hydrolysis Curves
The milk was prepared by mixing the low-heat milk powder with deionized water at a 1:10 proportion, resulting in a 3.1% (w/w) protein solution. The mixture was subjected to magnetic agitation and the homogenized mixture is hereafter delineated milk. The milk was loaded in a glass vessel and thermoregulated by a water bath at a constant temperature. After thermal equilibrium, the pH was fixed at 6.5. The protease previously diluted with deionized water was added to start the hydrolysis reaction. Different protease amounts were added to 40 g of milk. The release of free α-amino groups was followed by the pH-stat technique (G20 Compact Titrator, Mettler-Toledo) to obtain the hydrolysis progress. The pH was constantly monitored and 0.5 N NaOH was added to maintain the pH of 6.5 during the 60 min of reaction. Calculations were performed according to Equation (1): In the above equation, V is the added NaOH volume, N is the NaOH concentration, V T is the total reaction volume and α is the average degree of dissociation of the α-NH groups (an average value of 0.2 was used for all conditions). The total α-amino groups in milk proteins were determined by quantification with the o-phthaldialdehyde method (OPA) reported by Nielsen [31] after total hydrolysis in 6 N HCl at 110 • C for 24 h. The DH in each experiment was calculated according to Equation (2): where h T is the content of α-amino groups per protein mass in milk (10.7 meq/g prot ) and M P is the mass of protein in the reaction mixture (1.24 g). The release of α-amino groups against time was plotted to build the hydrolysis curves.

Experimental Design and Statistics
A central composite circumscribed (CCC) design was used to set the number and conditions of each hydrolysis experiment. Each experiment corresponded to a hydrolysis curve with fixed reaction conditions of temperature (T) and protease concentration (E). The T and E values resulting from the design were between 48-62 • C and 19-231 mAU, respectively (Table 1). A total of 10 experiments resulted from the CCC design shown in Table 2. Each experiment was randomly run as one replicate. The central point was replicated in Experiments 5 and 10 ( Table 2). Each experiment consisted of a hydrolysis curve where the product concentration (α-NH) was registered against the reaction time. The autotitrator registered one data point per second, resulting in 3600 points of product concentration (α-NH) as a function of time for a 1 h experiment. The first modeling level corresponds to the description of each hydrolysis curve performed by the estimation of the kinetics constants a and b from Equation (3). The logarithmic equation proposed by Márquez-Moreno and Fernández-Cuadrado [16] in terms of DH can be rewritten in terms of product concentration: where P is the concentration of released α-amino groups during milk protein hydrolysis, t is the reaction time and a and b are the kinetic constants of the model. These constants were estimated by the nonlinear least-squares method by fitting Equation (3) to the experimental data. The determination coefficient (R 2 ) and the standard error (se) were used to characterize the fitting quality according to Equation (4): whereσ 2 is the estimator of σ 2 calculated from the mean sum of squares and C jj is the diagonal element of the dispersion matrix (X T X) −1 . The second modeling level correlated each of the kinetic constants a and b with the reaction conditions. Equation (5) resulted from considering the variables T (x 1 ) and E (x 2 ) (operating conditions) and the combination of linear, square and interaction components.
A multivariable regression was performed to estimate the regression coefficients (β i ). A t test was used to calculate the statistical significance of each individual regression coefficient. This information was used to eliminate the nonsignificant regression coefficients and reduce the polynomial model. An analysis of variance (ANOVA) was used to test the significance of the regression (Tables 3-5). Additional hydrolysis experiments were performed to test the predictability of the model ( Table 6). The prediction range was calculated from the confidence intervals according to Equation (6): where x 0 corresponds to the vector of the experimental points considered in the validation experiment and x 0 T (X T X) −1 x 0 is the prediction variance used to calculate the prediction error of the kinetic constants a and b.

Results and Discussion
The description of the enzymatic hydrolysis of milk proteins was assessed through a two-level mathematical model using a previously published methodology [30]. The first-level model corresponds to the logarithmic equation used to model a wide spectrum of hydrolysis curves with different protein sources and proteases. The effect of E and T on the reaction performance can be observed in Figure 1. The hydrolysis curves generated similar results for Alcalase and Protamex. On the other hand, Neutrase exhibited high efficiency at the initial phase of the hydrolysis curve (high a values) and a sharp decrease in the hydrolysis rate after 5 or 10 min of reaction. Consequently, lower degrees of hydrolysis were obtained with Neutrase compared to Alcalase and Protamex.

Results and Discussion
The description of the enzymatic hydrolysis of milk proteins was assessed through a two-level mathematical model using a previously published methodology [30]. The firstlevel model corresponds to the logarithmic equation used to model a wide spectrum of hydrolysis curves with different protein sources and proteases. The effect of E and T on the reaction performance can be observed in Figure 1. The hydrolysis curves generated similar results for Alcalase and Protamex. On the other hand, Neutrase exhibited high efficiency at the initial phase of the hydrolysis curve (high a values) and a sharp decrease in the hydrolysis rate after 5 or 10 min of reaction. Consequently, lower degrees of hydrolysis were obtained with Neutrase compared to Alcalase and Protamex. These results agree with previous observations concerning Neutrase thermal stability. Neutrase activity has an optimal temperature range of 45 °C to 55 °C and its activity sharply decreases at temperatures over 55 °C [32]. The half-life was estimated to be approximately 10 min at 60 °C and a pH of 7. This important decrease in Neutrase activity explains the observed shape of the hydrolysis curve. The logarithmic equation was fit to determine the kinetic constants a and b in each experiment. These results are shown in Table 2 for Alcalase, Neutrase and Protamex. A high fitting quality was obtained for all the hydrolysis curves (R 2 > 0.99). This result agrees with previous fittings of hydrolysis curves [22,30]. As stated previously, the kinetic constant a corresponds to the initial slope of the hydrolysis curve, i.e., ⁄ at = 0 [30,33]. Thus, it can be predicted that the kinetic constant a increases when T and E increase. It can be observed in Figure 1 that faster hydrolysis rates and larger reaction extensions were obtained at higher temperatures and protease concentrations. The highest values for the kinetic constant a and the lowest values for the kinetic constant b were obtained with Neutrase. This is reflected in a high initial rate and an abrupt decrease after 5 or 10 min of reaction. Despite this high initial rate, Neutrase produced the lowest DHs. On the other hand, Alcalase and Protamex produced lower initial rates but higher DHs than Neutrase. This result can be explained by the These results agree with previous observations concerning Neutrase thermal stability. Neutrase activity has an optimal temperature range of 45 • C to 55 • C and its activity sharply decreases at temperatures over 55 • C [32]. The half-life was estimated to be approximately 10 min at 60 • C and a pH of 7. This important decrease in Neutrase activity explains the observed shape of the hydrolysis curve. The logarithmic equation was fit to determine the kinetic constants a and b in each experiment. These results are shown in Table 2 for Alcalase, Neutrase and Protamex. A high fitting quality was obtained for all the hydrolysis curves (R 2 > 0.99). This result agrees with previous fittings of hydrolysis curves [22,30]. As stated previously, the kinetic constant a corresponds to the initial slope of the hydrolysis curve, i.e., dP/dt at t = 0 [30,33]. Thus, it can be predicted that the kinetic constant a increases when T and E increase. It can be observed in Figure 1 that faster hydrolysis rates and larger reaction extensions were obtained at higher temperatures and protease concentrations. The highest values for the kinetic constant a and the lowest values for the kinetic constant b were obtained with Neutrase. This is reflected in a high initial rate and an abrupt decrease after 5 or 10 min of reaction. Despite this high initial rate, Neutrase produced the lowest DHs. On the other hand, Alcalase and Protamex produced lower initial rates but higher DHs than Neutrase. This result can be explained by the higher thermal stability of Alcalase and Protamex. Valencia et al. [15] reported an Alcalase activity loss of 20% after 3 h of reaction during the hydrolysis of 1.7% (w/w) salmon muscle proteins.
According to the obtained results, the kinetic constants a and b can characterize the different catalytic efficiencies that proteases exhibit during milk protein hydrolysis. The second modeling level was based on the correlation between the kinetic constants (a, b) and the reaction conditions T and E. The polynomial equation contained T and E as the independent variables. The multivariable regression results are shown in Tables 3-5 for Alcalase, Neutrase and Protamex, respectively. The nonsignificant coefficients (β i ) were eliminated from the polynomial equation to obtain a reduced model for each protease (Tables 3-5). The reaction conditions significantly affected the kinetic constant a for all proteases. A lower correlation between the kinetic constant a and operating conditions was obtained with Neutrase (R 2 = 0.8746) compared to Alcalase and Protamex. This can be explained by the lower thermal stability of Neutrase. A better correlation can be observed at lower temperatures for Neutrase (40-50 • C). The kinetic constant b was poorly correlated with the operating conditions when Alcalase and Protamex were used (R 2 = 0.7972 and R 2 = 0.8693, respectively). For this reason, the kinetic constant b was fixed at 0.0410 and 0.0423 for Alcalase and Protamex, respectively, corresponding to the central value. The reaction conditions impacted the kinetic constant b to a lesser extent than the kinetic constant a. Average values of b for different operating conditions have been used in previous publications [16,19,21,22,24]. This result is consistent with previous findings, where the kinetic constant b depended exclusively on the substrate concentration [30]. In the present study, the protein concentration was constant because the same milk formulation was used in all experiments. This explains the poor correlation of the kinetic constant b with the T and E. Regression analysis and the analysis of variance, shown in Tables 3-5. The response surfaces from the reduced models for the kinetic constant a were plotted against T and E in Figure 2.     The ANOVA results shown in Tables 3-5 indicated that these models were significant in describing the kinetic constant a as a function of operating conditions, with the exception of the kinetic constant a for Neutrase (p > 0.01). This result was explained by the thermal inactivation of Neutrase. A better fitting can be achieved at lower temperatures for Neutrase. Considering the results obtained with Alcalase and Protamex, the polynomial models can be used to calculate the values of the kinetic constant a from a set of operating conditions (T, E) and plot the hydrolysis curves considering the fixed values for the kinetic constant b. The predictability of these models was tested with additional experiments, which are presented in Table 6. The models for kinetic constant a were used to calculate their predicted values, while the experimental values were calculated from the fitting of the hydrolysis curves from validation experiments. The predicted vales of the kinetic constant a were correlated against their experimental values in Figure 3. Good agreement was obtained for both Alcalase (R 2 = 0.9830) and Protamex (R 2 = 0.9295) independently and with both proteases together (R 2 = 0.9528). For Neutrase, high prediction errors were obtained, thus, this protease was excluded from validation experiments. This error is explained by significant thermal inactivation of Neutrase. The predicted hydrolysis curves are plotted in Figures 4 and 5 for Alcalase and Protamex, respectively.  The predicted hydrolysis curves were in good agreement with the experimental curves, despite the errors observed in Table 6 for the kinetic constant a, reaching approximately 20% in some cases. Indeed, the predicted DH after 60 min of reaction resulted in less than 5% error for most of the cases (Table 7). This result is an important source of validation for the present model because it can not only be used to predict the DH, but also to study the reactor performance and design to achieve a certain DH value. As noted in a previous article [30], the prediction quality of this model can still be improved by using a different experimental design and/or expanding the studied range of the operating conditions. Furthermore, this methodology allowed us to compare the performance of three different commercial proteases. It was observed that in the range of temperatures studied, Neutrase was rapidly inactivated, achieving a lower DH after 60 min of reaction compared to Alcalase and Protamex. Protamex showed an initial rate twice that observed for Alcalase. This performance is reflected in the values of the kinetic constant a (Table 2). Nevertheless, Alcalase reached almost the same values of DH as Protamex after 60 min of reaction. This performance is reflected in the lower values of the kinetic constant b obtained with Alcalase ( Table 2). The performance of Protamex can be explained by a possible higher thermal inactivation and/or product inhibition than Alcalase. In resolution, this methodology can be used to evaluate the performance of any protease under different operating conditions. Different milk sources or formats can be evaluated using this The predicted hydrolysis curves were in good agreement with the experimental curves, despite the errors observed in Table 6 for the kinetic constant a, reaching approximately 20% in some cases. Indeed, the predicted DH after 60 min of reaction resulted in less than 5% error for most of the cases (Table 7). This result is an important source of validation for the present model because it can not only be used to predict the DH, but also to study the reactor performance and design to achieve a certain DH value. As noted in a previous article [30], the prediction quality of this model can still be improved by using a different experimental design and/or expanding the studied range of the operating conditions. Furthermore, this methodology allowed us to compare the performance of three different commercial proteases. It was observed that in the range of temperatures studied, Neutrase was rapidly inactivated, achieving a lower DH after 60 min of reaction compared to Alcalase and Protamex. Protamex showed an initial rate twice that observed for Alcalase. This performance is reflected in the values of the kinetic constant a (Table 2). Nevertheless, Alcalase reached almost the same values of DH as Protamex after 60 min of reaction. This performance is reflected in the lower values of the kinetic constant b obtained with Alcalase ( Table 2). The performance of Protamex can be explained by a possible higher thermal inactivation and/or product inhibition than Alcalase. In resolution, this methodology can be used to evaluate the performance of any protease under different operating conditions. Different milk sources or formats can be evaluated using this methodology. Considerations about the nitrogen content must be prevented for the proper calculation of the DH. In the case that autotitration is not available, a different analytical technique can be used for the quantification of released α-NH, such as the TNBS [34] or OPA [31] methods. Sample withdrawal is required in these methods in order for the sample to be mixed with the proper reagent. If parameter conditions other than protease dose needs to be tested, such as milk concentration, pH or E/S, a new factor can be added to the polynomial equation after including this parameter in the experimental design. Moreover, different experimental designs can be tried if better statistical properties are achieved in the experimental configuration. Thus, the actual methodology offers not only reliability, but also versatility.
are achieved in the experimental configuration. Thus, the actual methodology offers not only reliability, but also versatility.   Table 6.  Table 6. The logarithmic equation has largely been demonstrated to be a reliable model because of its good description, explanation, and prediction. The proposed modeling tool will improve the study of the correlation between DH and allergenicity during the enzymatic hydrolysis of milk proteins. An intermediate correlation between DH and allergenicity is further needed. In this way, DH can be used as an objective and quality control parameter during the hydrolysis process, thus allowing the reaction to be stopped at the desired allergenicity. This improvement will allow us to establish the operating conditions to produce a specific DH in the milk protein hydrolysate and, consequently, a specific functional property, for example, a specific allergenicity reduction.
The logarithmic equation has largely been demonstrated to be a reliable model because of its good description, explanation, and prediction. The proposed modeling tool will improve the study of the correlation between DH and allergenicity during the enzymatic hydrolysis of milk proteins. An intermediate correlation between DH and allergenicity is further needed. In this way, DH can be used as an objective and quality control parameter during the hydrolysis process, thus allowing the reaction to be stopped at the desired allergenicity. This improvement will allow us to establish the operating conditions to produce a specific DH in the milk protein hydrolysate and, consequently, a specific functional property, for example, a specific allergenicity reduction.  Table 6.

Conclusions
The performance of the enzymatic hydrolysis of milk proteins was characterized and mathematically modeled against the temperature and protease concentration. The enzymatic hydrolysis performance was readily and accurately predicted using this modeling methodology to compare the protease efficiency and temperature effect. The temperature and protease concentration significantly affected the initial rate of hydrolysis, while the reaction extension was not affected. The model allows the possibility of studying the effect of additional variables, such as operating conditions, by simply adding variables to the polynomial equation and using an adequate experimental design. This methodology can be used in future research to test the effect of hydrolysis conditions on the reduction in antigenicity of hydrolysates. In this way, antigenicity can be correlated with operating conditions and the process can be designed to obtain hypoallergenic milk products.  Table 6.

Conclusions
The performance of the enzymatic hydrolysis of milk proteins was characterized and mathematically modeled against the temperature and protease concentration. The enzymatic hydrolysis performance was readily and accurately predicted using this modeling methodology to compare the protease efficiency and temperature effect. The temperature and protease concentration significantly affected the initial rate of hydrolysis, while the reaction extension was not affected. The model allows the possibility of studying the effect of additional variables, such as operating conditions, by simply adding variables to the polynomial equation and using an adequate experimental design. This methodology can be used in future research to test the effect of hydrolysis conditions on the reduction in antigenicity of hydrolysates. In this way, antigenicity can be correlated with operating conditions and the process can be designed to obtain hypoallergenic milk products.