Reduction of Off-Flavor Generation in Soybean Homogenates: A Mathematical Model

The generation of off-flavors in soybean homogenates such as n-hexanal via the lipoxygenase (LOX) pathway can be a problem in the processed food industry. Previous studies have examined the effect of using soybean varieties missing one or more of the 3 LOX isozymes on n-hexanal generation. A dynamic mathematical model of the soybean LOX pathway using ordinary differential equations was constructed using parameters estimated from existing data with the aim of predicting how n-hexanal generation could be reduced. Time-course simulations of LOX-null beans were run and compared with experimental results. Model L2, L3, and L12 beans were within the range relative to the wild type found experimentally, with L13 and L23 beans close to the experimental range. Model L1 beans produced much more n-hexanal relative to the wild type than those in experiments. Sensitivity analysis indicates that reducing the estimated Km parameter for LOX isozyme 3 (L-3) would improve the fit between model predictions and experimental results found in the literature. The model also predicts that increasing L-3 or reducing L-2 levels within beans may reduce n-hexanal generation. Practical Application: This work describes the use of mathematics to attempt to quantify the enzyme-catalyzed conversions of compounds in soybean homogenates into undesirable flavors, primarily from the compound n-hexanal. The effect of different soybean genotypes and enzyme kinetic constants was also studied, leading to recommendations on which combinations might minimize off-flavor levels and what further work might be carried out to substantiate these conclusions.


Introduction
Soybean (Glycine max) is an important ingredient in many processed foods. However, their use may be limited by the generation, following homogenization, of so-called off-flavors (Lei and Boatright 2001;Lozano and others 2007). The presence of these flavors or aromas may cause variability in product quality, particularly beverages made using crushed soybeans. It is believed that in many plants the off-flavors are generated by a range of aldehydes known as C6-aldehydes or green leaf volatiles (GLVs) (Matsui 2006;Pulvera and others 2006). These GLVs are generated by the action of enzymes in the lipoxygenase (LOX) pathway (Feussner and Wasternack 2002). Of particular interest for food manufacturers is the means of controlling one such GLV, n-hexanal.
LOX is ubiquitous in plants and catalyses the dioxygenation of polyunsaturated fatty acids (PUFAs) containing (1Z,4Z)pentadiene systems (Brash 1999;Liavonchanka and Feussner 2006). Soybean seeds are rich in 2 such suitable substrates, linoleic acid (LA) and α-linolenic acid (Liu 1997). Before beans are crushed, LOX and LA are separated within the cell, but following homogenization are mixed and begin to react to form the products of the LOX pathway, notably n-hexanal. The products of the LOX reaction are the 13-and 9-hydroperoxides of the PUFA substrate (Brash 1999). Soybean seeds have 3 LOX isozymes (L-1, L-2, and L-3), which vary in product positional specificity and kinetic behavior (Axelrod and others 1981). If LA is the substrate, then certain specific 13-hydroperoxide isomers formed go on to be cleaved by hydroperoxide lyase, a CYP74, to form n-hexanal and 12-oxo-(Z)-9-dodecenoic acid (Matoba and others 1985b;Pulvera and others 2006). The 2 LOX product isomers that go on to form n-hexanal, 13HOD-S(Z,E), and 13HOD-R(Z,E) (see nomenclature), are formed in different proportions, as well as at different rates, by the 3 LOX isozymes (Andre and Funk 1986;Fukushige and others 2005). At higher pH, n-hexanal may itself be converted to n-hexanol by alcohol dehydrogenase (ADH), which has less impact on flavor quality (Matoba and others 1989).
There have been several studies using LOX-null strains of soybean, which aim to elucidate the mechanism by which n-hexanal and other GLVs are formed (Matoba and others 1985a(Matoba and others , 1985b(Matoba and others , 1989Pulvera and others 2006). In particular, studies using various existing soybean strains that lack one or more LOX isozymes show how different soybean genotypes produce different concentrations of n-hexanal following homogenization (Takamura and others 1991;Zhuang and others 1991;Nishiba and others 1995). From some of these studies, kinetic parameters can be derived that can be used to model the biological system mathematically.
While many laboratory studies for soybean n-hexanal generation have been published, no mathematical model of the LOX pathway was found that uses the existing data to try and describe observed behavior. This work has created a mathematical model of the soybean LOX pathway at near-neutral pH that produces time-course predictions of n-hexanal concentration comparable with existing laboratory data. Also the model has been used to predict and compare the effect of removing one or more LOX isozymes from the model. Further analysis may then suggest possible methods to reduce the generation of n-hexanal following soybean homogenization, or experiments that will help determine the source of off-flavor generation.

Materials and Methods
Kinetic parameters Bild and others (1977) gives values for K m and V max for purified L-1 at both pH 6.8 and pH 9.0. Unfortunately, explicit values and K m or V max under similar conditions for L-2 or L-3 were not found in the literature. Takamura and others (1991) gave measurements at pH 6.5 of relative LOX activity for wild-type beans containing all 3 isozymes, and for 3 lines of mutant beans, each only possessing just one of the isozymes. By assuming that the ratios of these results correspond to the relative V max of the 3 LOX isozymes, estimates of the different V max parameters were made. The same K m parameters were used for L-2 and L-3 as for L-1. The significance of the selection of K m values was later tested by sensitivity analysis (Figure 6). The parameters used in the model for the 3 LOX reactions are given in Table 1.
The relative activity of Hydroperoxide Lyase (HPL) against the different LA hydroperoxide isomers is given by Matoba and others (1985b). This is given relative to the value of 100% for 13HOD-S(Z,E), with the only other isomer found to have significant activity as substrate for HPL being 13HOD-R(Z,E) with a mean relative activity of 13.5%. The V max parameter used of 19 μmol·mg/min protein for 13HOD-S(Z,E) was the mean value given for wild-type beans. V max for 13HOD-R(Z,E) was then estimated as 13.5% of this value. A K m value of 40 to 60 μmol is given, and the mid value of 50 μmol has been used in the model for both active substrates.

Initial concentrations
For simplicity, the initial concentration of all metabolites in the model, except the active enzymes and the initial LA substrate, is set to zero. An estimate of initial LA concentration of 67 nmol/mL has been used. This is based on data for the fatty acid composition of typical soybeans given by Liu (1997) and scaled relative to the total protein concentration given by Matoba and others (1985a) of 10 to 20 mg protein/mL. To account for the difference in V max for purified enzyme and for whole bean homogenate, the assumption was made that in wild-type beans L-1, L-2, L-3 are present in similar quantities (Takamura and others 1991) and that each constitutes approximately 1% of total protein (Brash 1999). The final set of V max and K m estimates is given in Table 1.

Product ratios
Data for the ratio of 13-and 9-HOD isomers produced by each of the LOX isozymes is taken from Andre and Funk (1986) for the L-1 isozymes, and from Fukushige and others (2005) for L-2 and L-3. For L-1, product ratios are given for pH 7 and are assumed to approximate values for pH 6.8, while values for L-2 and L-3 are given at pH 6.8. The product ratios for each LOX isozyme used in the model are shown in Table 2.

Modeling methods
The availability of K m and V max kinetic parameters for L-1 and HPL (Bild andothers 1977, Matoba andothers 1985b), and the relative activity of the 3 LOX isozymes (Takamura and others 1991), leads to the use of a mathematical model with Michaelis-Menten kinetics to simulate the biochemical network of the LOX pathway. Copasi (http://www.copasi.org) is a desktop software application, free for noncommercial use, designed for the modeling of such biochemical networks, and was used to create a deterministic model using a set of ordinary differential equations. Time-course data of the concentrations of the different metabolites in the model could then be generated for further analysis, for different parameter sets, and null-bean simulations. The biochemical network modeled is shown in Figure 1.

Assumptions
For the model, the 1st assumption is that the activity of all 3 LOX isozymes is governed by Michaelis-Menten kinetics, and that once the product is disassociated from the LOX enzyme, the reaction is irreversible. Second, there is an assumption that each of the 8 possible 13-or 9-hydroperoxide isomers is produced at a different, fixed ratio for each of the 3 LOX isozymes. Also assumed is that oxygen is freely available for the LOX reaction, and is not rate limiting. For the HPL reaction, Michaelis-Menten kinetics were also used. It is assumed that HPL is present in similar quantities, has similar substrate specificity, and kinetic parameters, for all soybeans, regardless of which LOX enzymes are present. Only one substrate, free LA, is present in the model. Therefore, Table 2-Ratios of product isomers in the model for each LOX isozyme at pH 6.8 and pH 9.0. Data from Andre and Funk (1986) and Fukushige and others (2005).

L-1
L-2 L-3 (% of total (% of total (% of total Product product) product) product) 0.4 1.6 9.7 there is an assumption that n-hexanal is only produced via the LOX pathway, and that only free LA is a suitable substrate. The final initial assumption is that there are no feedback or inhibitory effects from any of the products or by-products generated by the model over the time course.
To model the predicted behavior of wild-type beans that contain all 3 LOX isozymes, the combined action of the 3 separate reactions is assumed to have an additive effect. Each reaction L-1, L-2, and L-3 competes for the LA substrate, and each produces the same range of 8 possible products in different proportions. The transient concentrations of these 8 products are the sums of those produced by all 3 reactions. In addition, by setting V max for one or more LOX isozymes in the model to zero, the effect of using real beans lacking one of more LOX isozymes can be simulated.
The 2nd step in the model pathway is the cleaving of 13HOD-S(Z,E) or 13HOD-R(Z,E) by HPL to form n-hexanal and 12oxo-(Z)-9-dodecenoic acid using the parameters given by Matoba and others (1985b). Model simulations were run to give time- course data for the generation of n-hexanal for all the different combinations of LOX-null beans.
The model is (will be) available for download in SBML format at http://www.biomodels.net.

Results and Discussion
Primary LOX products Figure 2 shows the concentration of products after 100 min predicted by the combined LOX model without the HPL reactions included. Also shown is the predicted concentration of the primary LOX products after 100 min for each of the possible combinations of LOX isozymes. In each case, with the exception of simulated L-3 only beans (L 3 , see nomenclature), the majority of the product formed is 13HOD-S(Z,E). The products formed in greatest concentrations by the simulated L 3 beans are 9HOD-S(Z,E) and 9HOD-R(Z,E). The preferred substrate of HPL, 13HOD-S(Z,E), was modeled in the next step of the LOX pathway, and so the relative concentrations in which it is formed in different bean genotypes are significant. In the model, the L 2 beans formed the most 13HOD-S(Z,E), closely followed by the L 12 , L 23 , and L 123 (wild type) beans, then at around one-third lower concentrations the L 1 beans, the L 13 beans, and finally the L 3 beans producing by far the lowest concentration.

Null-bean simulation compared with experimental data
Several sets of experimental results are available in the literature that gives n-hexanal generation for different LOX-null beans under various different experimental conditions. Table 3 summarizes the results by showing the n-hexanal generated by each LOX-null bean relative to the wild type for the different experimental results, alongside the simulated values.
Model L 2 , L 3 , and L 12 values all fell within the range of values obtained experimentally. L 13 and L 23 beans were also fairly close to the range of experimental values with 0.70 and 1.00 model compared to experimentally derived ranges of 0.38 to 0.63 and 0.73 to 0.97, respectively. The L 1 beans show the greatest difference between the model and the experimental data, with a model value of 0.80 times the wild-type n-hexanal concentration compared with a much lower range of 0.06 to 0.50 times the wild type in the lab experiments. Where values for L 0 beans were given in experimental data, there was a background level of nhexanal generation observed of 0.06 to 0.22 times the wild-type value.
Qualitatively at least, the differences between simulated null beans were mostly fairly similar to those observed in the null-bean experiments. In the literature, the soybean genotype that consistently produces the most n-hexanal is that which only contains the L-2 isozyme (L 2 ) (Takamura and others 1991;Nishiba and others 1995;Pulvera and others 2006). This is also true of the model. This is to be expected as n-hexanal is derived mainly from 13HOD-S(Z,E) that is produced in highest proportions by L-2. Since there is a finite amount of substrate, when other LOX isozymes are present more of it is converted into other products that are not converted into n-hexanal. Similarly, it is to be expected that L 3 beans produce relatively low amounts of n-hexanal, as L-3 generates the least 13HOD-S(Z,E) of the 3 isozymes. In the model, the L 3 beans do produce the least n-hexanal, and within the range concentrations relative to the wild type found in the literature. Experimentally, however, while L 3 beans produce very low amounts of n-hexanal, the L 1 beans also produce similarly low or lower amounts (Matoba and others 1985b;Takamura and others 1991;Nishiba and others 1995;Pulvera and others 2006). This is not true of the model results, in which L 1 beans produce more than 7 times the amount of n-hexanal as the L 3 beans. The simulated beans in which just one isozyme is missing produce n-hexanal in a similar order to that found in most experiments surveyed, with L 12 producing the most, followed by L 23 , wild-type beans, and finally L 13 producing the least (Matoba and others 1985a; Zhuang and others 1991).
It is unclear why beans containing just the L-1 isozyme should produce similarly low amounts of n-hexanal as the beans containing just L-3. Takamura and others (1991) measures LOX activity for L 1 beans as more than 3 times that for L 3 beans at pH 6.5, and the proportion of 13HOD-S(Z,E) is given as 57.4% for L-1 (Andre and Funk 1986) compared to just 6.8% for L-3 (Fukushige and others 2005). Other possibilities for low n-hexanal generation by L 1 beans is that other enzymes affecting generation are absent or present in different quantities, or that there is less free fatty acid substrate present. Matoba and others (1985) does show that there is some variation in fatty acid composition between bean genotypes, but there is no data for L-1 only beans. It may also be that the kinetics of the reactions does not proceed as assumed in the model and that some feedback or inhibitory effects of the by-products generated may occur.

Time-course results compared with experimental data
Using the estimate of initial free LA of 67 nmol/mL, estimated using the value of 10 to 20 μmol/mL protein content given by Matoba and others (1985a) results in a model prediction of around 50 nmol/mL n-hexanal produced by wild-type beans. The value given in the paper is around 0.3 nmol/mL for wildtype Suzuyakata beans, more than a 10-fold decrease. This may be due to the degradation of free LA by other pathways, or an overestimation of initial LA or LOX concentrations. In order to compare the relative n-hexanal generation for different bean types between model and experimental data, the data were therefore plotted relative to the maximum wild-type value for both model predictions and experimental data, respectively. Figure 3A shows the comparison between different simulated LOX-null beans, and Figure 3B shows the experimental results (taken directly from the paper itself, which gives no standard errors and a pH range of 6.5 to 7.0). Comparing the overall pattern of the results, some correlation is observed. The L 13 (L-2 null) beans show clearly the lowest generation of n-hexanal in both model and experiment, with 0.4 and 0.6 times the peak wild-type value, respectively. The experimental data then show (in order of increasing peak n-hexanal) the L 23 , the L 123 (wild type), and the L 12 beans steadily increasing from an initial value to a peak after 60 min of around 0.8 (L 23 ) and 1.2 (L 12 ) times the peak wild-type value. The model predicts a similar order, with L 12 beans producing slightly more n-hexanal than the wild type, and L 23 slightly less, although the difference between the 3 curves appears smaller. There are 2 main qualitative differences between the experimental and model datasets. The 1st is the curves representing the L 2 beans. Experimentally, these beans consistently produced the most n-hexanal throughout the time course, reaching a maximum value of around 1.8 times the wild-type value after 60 min. In contrast, in the model, while the peak value for L 2 is still the highest predicted, it is much closer to the values predicted for the L 12 , wild type, and L 13 beans, and for much of the time course, the predicted n-hexanal is lower than that of other beans. The experimental data for all beans show an initial base level of n-hexanal that is not accounted for in the model. For simplicity, initial n-hexanal concentration is set to zero in the model, but the experimental data give a mean initial concentration of around 0.14 times peak wild-type n-hexanal or 0.4 nmol/mL. In addition, Table 3 shows that in experiments where beans with none of the LOX isozymes present (L 0 ), n-hexanal is still generated in concentrations comparable to the null beans L 1 and L 3 . It is possible that this is a steady state value of n-hexanal concentration prior to homogenization, when LA and LOX are isolated. After homogenization, LA and LOX are mixed, with the end result that n-hexanal concentration increases. Alternatively, there may be another metabolic pathway by which n-hexanal is generated. 0.07 n/a n/a n/a 0.40 0.10 0.29 L 12 1.04 0.89 1.17 0.59 n/a n/a n/a L 13 0.70 0.63 0.33 0.38 n/a n/a n/a L 23 1.00 0.73 0.90 0.97 n/a n/a n/a L 0 0.00 n/a n/a n/a 0.22 0.10 0.06 Figure 4 shows the same data but with a direct comparison between model and experimental values for each type of bean. Again there appears to be a rough correlation, with the exception of the L 2 beans, and the presence of initial concentrations of nhexanal observed experimentally in all bean types. In the present model, 13HOD-S(Z,E) and 13HOD-R(Z,E) concentrations are only depleted by the activity of HPL. In reality, these molecules may by degraded or metabolized by other pathways, and the rates of these pathways will have an effect on n-hexanal concentration. An experiment in which 13HOD-S(Z,E) and 13HOD-R(Z,E) levels are recorded while HPL activity is blocked or otherwise inhibited would enable an estimation of these rates, and enable their inclusion in the model. If we assume that these parameters and those relating to HPL remain the same regardless of the presence or absence of the 3 LOX isozymes, then while varying them will have an effect on the absolute concentrations of n-hexanal, they should have little or no effect on the relative concentrations obtained for the different null bean simulations.
The model is based on data at or near-neutral pH. Some data are however available at higher pH. In particular, it is shown that purified L-1 enzyme has peak activity at around pH 9.0 to 9.5 (Axelrod and others 1981), and as may be expected that LOX activity for L 1 beans is also much higher at pH 9.5 (Takamura and others 1991). Andre and Funk (1986) do give data for the relative concentrations of product isomers produced by L-1 at pH 9, but the precise values for each isomer is not given explicitly, and no statistical errors are given. While Fukushige and others (2005) give excellent data for isomer ratios for L-2 and L-3 at pH 6.8, no data are given on how the ratios differ at higher pH. A consistent set of experiments in which product ratios for all 3 isozymes are determined at a range of pHs would be useful in extending the model to simulate behavior at different pH. The behavior of HPL and estimates for the K m of the LOX isozymes at increasing pH would also be required to extend the model in this way. Matoba and others (1989) suggests a means by which n-hexanal is converted to n-hexanol at higher pH by the activity of ADH. In this experiment, n-hexanal peaked after around 15 min after homogenization before returning to near its initial value. Since n-hexanol has been reported as having a less negative effect on soybean flavor, increasing this conversion of n-hexanal to nhexanol may be desirable.

Sensitivity analysis
Changing the kinetic parameters of each of the 3 LOX isozymes in the model affects the peak n-hexanal concentration for simulated null beans to varying degrees ( Figure 5). Most of the variation was seen by reducing L-3 K m (increasing V max ), with peak n-hexanal of around one-fifth that predicted using the default value following a 1000-fold decrease. Increasing L-3 K m 1000fold only results in a very slightly higher peak n-hexanal concentration. Increasing L-2 K m (reducing V max ) 1000-fold reduces peak n-hexanal concentration to around 0.6 times the value for the default parameter, while decreasing 1000-fold has relatively little effect. Increasing or decreasing L-1 K m (decreasing or increasing V max ) has the least effect relatively of the 3 isozymes, with peak n-hexanal remaining in a similar range for all the parameter values tested. Reducing or decreasing V max values of the 3 LOX isozymes in the model should have a similar effect as changing individual LOX protein concentration and so more precise data for the concentration of the different LOX isozymes in the different bean types would be useful in testing the model further. The analysis would suggest that increasing levels of L-3 or reducing L-2 would have the greatest effect on reducing n-hexanal generation.
Clearly, altering the parameter estimates used for the LOX reactions will affect model predictions for n-hexanal generation. Figure 6 shows how different parameter values will affect the comparison between "wild type" and each of the simulated null Figure 5-Simulated peak n-hexanal concentration against varying K m parameter values for L-1, L-2, and L-3 isoenzymes in "wild type" beans. Peak n-hexanal is given relative to the value at default K m (0.49 mM). Figure 6-Peak n-hexanal values relative to "wild type" for each type of simulated null bean for increasing K m parameter values for L-1, L-2, and L-3 isoenzymes. The final set of values marked "Exp" represent the mean experimental values found in the literature (Table 3).
Vol. 75, Nr. 7, 2010 r Journal of Food Science R137 beans, and how these predictions compare with experimentally obtained values. The set of experimental values (marked "Exp") represents the mean of the values in Table 3 sourced from the literature. The pattern of experimental results appears to be high generation of n-hexanal for L 2 beans (around 1.5 time wild type), followed by L 23 and L 12 beans (both around 0.9 times wild type), then L 13 (around 0.5 times wild type), with the lowest n-hexanal generation in the L 1 and L 3 beans (both less than 0.3 times wild type). When comparing to the default model values (K m = 0.49 mM), the major differences are that values are much closer to the wild-type beans for both L 2 and L 1 beans. L 23 , L 13 , and L 12 simulated beans all generate slightly more n-hexanal than observed experimentally, while simulated L 3 generates less than the mean experimental value.
The best fit between model and experiment in this set of results is given by decreasing L-2 K m by a factor of 10. This increases the L 2 n-hexanal near to the experimental value of around 1.5 times wild type, while increasing the L 3 n-hexanal nearer to the experimental value. However, this also has the effect of increasing the value for L 1 beans nearer to the L 23 beans, and the L 12 bean n-hexanal increases along with the L 2 beans, so the fit is still not perfect. Increasing L-1 K m slightly may improve the model by reducing n-hexanal for the L 1 and L 13 beans toward the experimental values, while leaving values for the other beans unaffected. Reducing L-1 K m from the default value of 0.49 mM does not improve model performance as it increases the value for both the L 1 and L 13 beans closer to the wild-type value, away from the experimental values. Increasing L-2 K m reduces the fit between the model and experimental values, as it simultaneously increases L 1 and reduces L 2 n-hexanal. Reducing L-2 K m has little effect on relative levels between the different beans. More sophisticated parameter fitting and optimization may be required in order to find the best fit between the available experimental results and model predictions.

Conclusion
While the model does partially reflect the behavior observed experimentally, there is some room for improvement. In particular, the model predicts that increasing L-3 or reducing L-2 levels within beans may reduce n-hexanal generation, as may be expected. On the other hand, the difference between the behavior of simulated and laboratory L 1 beans is not accounted for. Sensitivity analysis for LOX K m values suggests that parameter sets that compare more favorably with experimental results could be found by parameter fitting or optimization techniques. The results of Matoba and others (1989) suggest that including the effects of increasing pH and the inclusion of the reaction converting nhexanal to n-hexanol by ADH at higher pH will be vital in using the model to predict optimum conditions for reducing n-hexanal generation during soybean processing.