Curation and Analysis of a Saccharomyces cerevisiae Genome-Scale Metabolic Model for Predicting Production of Sensory Impact Molecules under Enological Conditions

One approach for elucidating strain-to-strain metabolic differences is the use of genome-scale metabolic models (GSMMs). To date GSMMs have not focused on the industrially important area of flavor production and, as such; do not cover all the pathways relevant to flavor formation in yeast. Moreover, current models for Saccharomyces cerevisiae generally focus on carbon-limited and/or aerobic systems, which is not pertinent to enological conditions. Here, we curate a GSMM (iWS902) to expand on the existing Ehrlich pathway and ester formation pathways central to aroma formation in industrial winemaking, in addition to the existing sulfur metabolism and medium-chain fatty acid (MCFA) pathways that also contribute to production of sensory impact molecules. After validating the model using experimental data, we predict key differences in metabolism for a strain (EC 1118) in two distinct growth conditions, including differences for aroma impact molecules such as acetic acid, tryptophol, and hydrogen sulfide. Additionally, we propose novel targets for metabolic engineering for aroma profile modifications employing flux variability analysis with the expanded GSMM. The model provides mechanistic insights into the key metabolic pathways underlying aroma formation during alcoholic fermentation and provides a potential framework to contribute to new strategies to optimize the aroma of wines.


Introduction
Saccharomyces cerevisiae is employed to produce many specialized metabolites for the beverage, food, cosmetics, and pharmaceutical industries [1]. When producing alcoholic beverages, in particular wines, sensorial properties are of major importance to quality. Therefore, it is beneficial for the wine industry to be able to systematically optimize the production of key aroma metabolites. However, in order to control flavor compound formation and optimize the quality of wine, a better understanding of the aroma formation during fermentation is essential. This is not a simple task, as the metabolic pathways leading to flavor and aroma impact molecules are complex and not fully understood. Despite the large amount of work using systems biology approaches to understand aroma production during fermentation, we are still only beginning to understand how metabolic changes impact aroma production [2,3]. Achieving a better understanding will allow more control over alcoholic beverage production.
(GSMMs) have been developed for yeast that show the potential of these simulations to predict the impact of genetic changes on cell phenotype [33][34][35][36][37]. In addition, Vargas et al. [38] applied a GSMM to study wine fermentation over time. While earlier metabolic modeling studies that focused on simulating yeast metabolism successfully predicted flux distributions, they were not genome scale representations of yeast metabolism (fewer than 100 reactions) [30,32,39], making it difficult to study individual fluxes that are likely related to overall aroma formation. Despite the increasingly expanded coverage of yeast metabolic network models over the last few years, all published GSMMs of yeast still contain a partial or incomplete representation of the formation of aroma impact molecules. More specifically, GSMMs up to this point have not included MCFA ester formation, as well as production of important sulfur-containing compounds such as dimethyl sulfide. Larger GSMMs have been used to simulate yeast fermentations, but studies applying them have mostly lacked exploration of enological fermentations where nitrogen was not the growth limiting nutrient and/or the system was not anaerobic [36,37,39]. Vargas et al. [38] did apply an earlier yeast GSMM model to dynamically model a wine fermentation, although their study did not cover aroma formation. While there have been at least 12 GSMMs for yeast released since the first one was published in 2003 by Förster et al. [33], the available GSMMs still lack a comprehensive coverage of essential metabolic pathways needed to investigate aroma formation during alcohol fermentation.
In this study, we present an expanded genome-scale metabolic model of S. cerevisiae to incorporate aroma forming pathways including parts of the Ehrlich pathway, pathways related to ester synthesis, and parts of the sulfur reduction pathways. A widely used published reconstruction of yeast metabolism referred to as iMM904 [40,41] was taken as a framework to incorporate known aroma forming reactions, as well as relevant chemical, non-enzymatic conversions and transport reactions. Our expanded GSMM was validated using published data from a study that employed a commercial wine yeast strain EC1118 under enological conditions (synthetic grape juice/must medium using an industrial wine strain and anaerobic fermentation) [30]. Furthermore, we applied the model to understand the underlying metabolic differences regarding aroma production under two different growth conditions. Here, we illustrate the advances and remaining limitations of this expanded GSMM to simulate metabolism and production of aroma impact molecules under enological conditions.

Model Curation and Development
A current model was evaluated as a base model for expansion, a widely adopted stoichiometric model of the S. cerevisiae metabolic network iMM904 [40,41]. The most recent GSMM, Yeast 8 [42], was also evaluated as a base model for this work. While predicted fluxes of primary metabolites with an expanded Yeast 8 model were similar to those using the iMM904 base model, prediction of secondary metabolites was far superior with the iMM904 base model (data not shown). Furthermore, Heavner and Price [43] concluded that yeast GSMM performance and predictions varies according to simulation environments and evaluation metrics. Thus, all subsequent work was performed using the expanded iMM904 model. iMM904 was expanded by incorporating more fusel alcohol production routes into the Ehrlich pathway, adding pathways for acetate and medium-chain fatty acid (MCFA) ester formation reactions, and expanding sulfur reduction pathways. An overview illustrating the key additions and modifications to the base model to construct the expanded model (iWS902) is represented in Figure 1 and Table 1, and comprehensive details are provided in the Supplementary Information Tables S1-S3. Uptake exchange reactions for the sugar, essential amino acids, and other nutrients were set according to measured flux values from a chemically defined medium during anaerobic nitrogen-limited fermentations found in the literature (Supplementary Table S1 from Quiros et al. [30]). No strictly anaerobic growth was predicted with flux balance analysis with this particular medium composition. Therefore, an extremely low O 2 exchange flux of −0.01 mmol g DW −1 h −1 was assumed for our simulations. This is a reasonable assumption because Aceituno et al. demonstrated Figure 1. Representation of the the components of the yeast genome-scale metabolic model (GSMM) related to aroma formation. Each colored section shows metabolism related to a respective class of aroma compounds, including the expanded pathways related to amino acid and fatty acid degradation and sulfur reduction pathways, as well as formation of fusel alcohols, esters, and sulfur compounds in yeast model metabolism. The entire metabolic network utilized is shown in the Supplementary Material ( Figure S2).
Sulfur-containing compounds are also formed during alcoholic fermentation. These compounds are diverse in character from simple thiols (e.g., methanethiol) to sulfides (e.g., hydrogen sulfide) to thioesters. Given that there is an extensive variety of sulfur-containing compounds that can arise during alcoholic fermentation and that many volatile sulfur pathways in yeast are still poorly understood [8], we mainly focused on the assimilation of sulfur and its links to amino acid metabolism. Sulfur compounds are produced by the sulfate-reduction sequence pathway where, by means of sulfate permease (SUL1 or SUL2), sulfate is taken up from the medium and reduced to sulfite, followed by further reduction to sulfide through MET5 and MET10, respectively [9,58]. Next, excess sulfide is either used to form hydrogen sulfide, which is excreted from the cell or becomes incorporated into the amino acid synthesis pathway. Sulfide combines with O-acetyl-homoserine to form homocysteine [5]. From here, methionine-γ-lyase can convert the reaction from methionine to a methanethiol, which can subsequently serve as a precursor to some thioesters and sulfides ( Figure 1). We extended the coverage of known sulfur producing pathways in our GSMM by adding pathways for methionol, methionyl acetate, and dimethyl disulfide. The other parts of sulfur metabolism were already included in the GSMM.
In order to properly model metabolic impacts of aroma formation, we incorporated and added to the reconstruction known aroma forming reactions, as well as relevant chemical, non-enzymatic, conversions and transport reactions. In total, we integrated 58 reactions, 43 metabolites, and three enzyme-catalyzed reaction-associated genes involved in flavor formation into the base GSMM to form the newly expanded model iWS902. iWS902 contains eight compartments: cytosol, endoplasmic reticulum, extra organism, Golgi apparatus, mitochondria, nucleus, peroxisome, and vacuole. Also, iWS902 is annotated according to BiGG database convention (http://bigg.ucsd.edu/). related to aroma formation. Each colored section shows metabolism related to a respective class of aroma compounds, including the expanded pathways related to amino acid and fatty acid degradation and sulfur reduction pathways, as well as formation of fusel alcohols, esters, and sulfur compounds in yeast model metabolism. The entire metabolic network utilized is shown in the Supplementary Material ( Figure S2).

Flux Balance Analysis
FBA is a widely used approach for studying biochemical networks. This approach allows prediction of the flux of metabolites through the network [46,47]. It is a technique based on mass balances. For m metabolites and n reactions, a m × n stoichiometric matrix S can be formulated and, if we neglect the accumulation of metabolites (pseudo-steady-state assumption), a mass balance for all metabolites can be written as where v is the flux vector [mmol g DW −1 h −1 ]. Furthermore, upper and lower bounds for each flux can be included to define the reversibility of each reaction, as well as to account for experimental observations of extracellular metabolite fluxes (for example, nutrient utilization or product formation).
Here, constraints on the exchange reactions of the amino acids and nutrients taken up from the medium are taken from those measured in Quiros et al. [30]. The uptake and secretion fluxes are expressed in mmol g DW −1 h −1 (note that a negative flux for an exchange reaction means uptake of the corresponding metabolite). In our study, we simulated two different growth conditions described by Quiros et al. [30] (both in chemostats at steady state using a defined medium with a glucose feed concentration of 280 g/L): (i) temperature of 16 • C with a growth rate, µ, of 0.1 h −1 (Case I) and (ii) temperature of 28 • C with a µ of 0.25 h −1 (Case II). In addition, we simulated Quiros et al. [30] (in chemostats at steady state using a defined medium with a glucose feed concentration of 240 g/L): (i) temperature of 16 • C with a growth rate, µ, of 0.1 h −1 (Case III) and (ii) temperature of 28 • C with a µ of 0.25 h −1 (Case IV) (results shown in Figure S1 of the Supplementary Material). In these studies, the growth rate was fixed by the dilution rate for each case, as the dilution rate is equal to the growth rate for a chemostat. Though not typical of grape juice that contains relatively equal concentrations of glucose and fructose, the defined medium used in this study used only glucose. Both cases had a yeast assimilable nitrogen level (YAN) of 220 mg/L in the chemostat feed. To assess the predictive power of the model, we performed FBA using only measured experimental consumption rates [30] to set constraints. Production rates for key metabolites were not constrained; therefore, predictions could be used to evaluate the predictive performance of our expanded GSMM. However, when assessing differences in predicted intracellular metabolic fluxes for the two cases, we also constrained secretion fluxes, in addition to the consumption rates, using data from Quiros et al. [30]. A more detailed overview of model constraints and conditions are represented in the Supplementary Material (Table S4). Since metabolic models are usually highly underdetermined (there are many more reactions than metabolites), a minimized or maximized objective function is applied to the model to solve for all fluxes. Here, an optimization problem was solved to maximize growth rate (µ), a biologically relevant objective function, with linear programming as shown below where LB and UB are vectors containing the lower and upper bounds, respectively, for all of the fluxes in units of mmol g DW −1 h −1 .

Flux Variability Analysis
Flux variability analysis (FVA) is a technique which identifies the allowable range of flux values through a given reaction by finding its maximum and minimum flux values for a given optimal (i.e., maximum) objective value [48]. This analysis method begins by finding the optimal value of the objective function, as specified in Equation (2), for a given set of constraints. The objective function is set for each reaction in the network and minimum and maximum flux values for each reaction are computed while using the constraints related to the optimal solution [49]. The FVA problem can be constructed as follows: where v i is the flux value for each reaction, both maximum and minimum flux values of a reaction were calculated to determine the full range. Z 0 is the optimal (maximum) value calculated by FBA. λ, whose value can be between 0 and 1, is a parameter to govern whether the analysis is carried out at a suboptimal network state (λ < 1) or at an optimal state (λ = 1). The range of a flux i (δ i ) is determined in Equation (4) as the difference between the maximum and the minimum value: For this analysis, we constrained the model using the same set of experimental uptake and secretion fluxes as in the FBA (see Section 2.2) (constrained case). Additionally, we performed another FVA where all substrate and production rates were unconstrained (unconstrained case). Since it is of interest in this study to maximize aroma production and the complexity of the aroma profile of wines, it is necessary to increase the flux of certain aroma compounds. It must be noted that FBA allows for optimization of a single flux. Therefore, this problem was resolved by adding an extra reaction to the network whereby the 1-propanol, 3-methyl-1-butanol, 2-methyl-1-propanol, 2-methyl-1-butanol 2-phenylethanol, indole-3-ethanol, and 4-(2-hydroxyethyl) phenol were lumped into a reaction using functions in Cobra Toolbox. Note: these lumped reactions were only added for the FVA. Subsequently, a transport reaction from the cytosolic lumped product to the same extracellular product was added. This allows for optimization of all seven fusel alcohols at the same time.

Computing Environment
Modeling was performed in MATLAB ® 2018b (The MathWorks, Inc., Cambridge, MA, USA) using Cobra Toolbox version 3.0 [50][51][52] and implemented on a Windows 10 (Microsoft Corporation, Redmond, WA, USA) Intel ® (Intel Corporation, Santa Clara, CA, USA) Core™ i7-7500 CPU @ 2.70 GHz-2.90 GHz processor. Git version 2.19. and curl version 7.61.1 were installed before cloning COBRA with Github and initializing COBRA in Matlab. The solver used for solving the linear and mixed integer problems was Gurobi version 9.0 (Gurobi Optimization, LLC, Beaverton, OR, USA). The GSMM was imported into MATLAB, as an SBML file, and curated using Cobra Toolbox. For network visualization, the model was exported into Omix TM version 1.9.34 (Omix Visualization GmbH & Co. KG, Lennestadt, Germany) [53] and Escher (The Regents of the University of California, Oakland, CA, USA) [54].

Coverage of Aroma Pathways in the Model (GSMM)
When examining current GSMMs, it becomes apparent that some essential aroma impact molecules are missing. Therefore, we addressed this deficiency through literature curation and pathway incorporation into an existing GSMM (see Figure 1 and Table 1 for details). The existing model we used contained five relevant (higher) alcohols and their respective esters, which are shown in Table 1. Higher alcohols, which form the highest concentrations of aroma impact molecules, are produced predominately via the Ehrlich pathway [10,55]. This three-step pathway ( Figure 1) features (i) a transamination reaction between an amino acid and α-ketoglutarate, (ii) a decarboxylation reaction facilitated by several pyruvate decarboxylases where an α-keto acid is converted to a branched-chain aldehyde, and (iii) a reduction step catalyzed by alcohol dehydrogenases reducing the aldehyde to a higher alcohol [4,56]. Higher alcohols and their respective pathways that were added to the GSMM include 1-propanol, tyrosol, 1-hexanol, and methionol, which have aroma attributes ranging from sweet, floral, and fruity (tyrosol) to meaty and oniony (methionol).
Two different classes of esters, acetate esters and fatty acid ethyl esters, are formed depending on their use of acetyl-CoA or acyl-CoA, respectively ( Figure 1). Acetate esters are compounds such as ethyl acetate, isoamyl acetate, isobutyl acetate, 2-methylbutyl acetate, hexyl acetate and 2-phenylethyl acetate. Many of these acetate esters are contained in existing GSMMs with the exception of hexyl acetate, which we added to our model. Members of the second class of esters, i.e., fatty acid ethyl esters, are formed from a reaction between a bonded CoA-MCFA and ethanol as illustrated (see Figure 1). Fatty acid ethyl esters include compounds such as ethyl butanoate, ethyl hexanoate, ethyl octanoate, and ethyl decanoate [57].
Sulfur-containing compounds are also formed during alcoholic fermentation. These compounds are diverse in character from simple thiols (e.g., methanethiol) to sulfides (e.g., hydrogen sulfide) to thioesters. Given that there is an extensive variety of sulfur-containing compounds that can arise during alcoholic fermentation and that many volatile sulfur pathways in yeast are still poorly understood [8], we mainly focused on the assimilation of sulfur and its links to amino acid metabolism. Sulfur compounds are produced by the sulfate-reduction sequence pathway where, by means of sulfate permease (SUL1 or SUL2), sulfate is taken up from the medium and reduced to sulfite, followed by further reduction to sulfide through MET5 and MET10, respectively [9,58]. Next, excess sulfide is either used to form hydrogen sulfide, which is excreted from the cell or becomes incorporated into the amino acid synthesis pathway. Sulfide combines with O-acetyl-homoserine to form homocysteine [5]. From here, methionine-γ-lyase can convert the reaction from methionine to a methanethiol, which can subsequently serve as a precursor to some thioesters and sulfides ( Figure 1). We extended the coverage of known sulfur producing pathways in our GSMM by adding pathways for methionol, methionyl acetate, and dimethyl disulfide. The other parts of sulfur metabolism were already included in the GSMM.
In order to properly model metabolic impacts of aroma formation, we incorporated and added to the reconstruction known aroma forming reactions, as well as relevant chemical, non-enzymatic, conversions and transport reactions. In total, we integrated 58 reactions, 43 metabolites, and three enzyme-catalyzed reaction-associated genes involved in flavor formation into the base GSMM to form the newly expanded model iWS902. iWS902 contains eight compartments: cytosol, endoplasmic reticulum, extra organism, Golgi apparatus, mitochondria, nucleus, peroxisome, and vacuole. Also, iWS902 is annotated according to BiGG database convention (http://bigg.ucsd.edu/). A comprehensive table of reactions (enzyme commission number, reaction equation, etc.) (Table S3), their associated genes (Table S1), and metabolites (Minimum Information Required in the Annotation of Models (MIRIAM), composition, charge number, compartment, etc.) ( Table S2) that were added to the model are shown in the Supplementary Material (Supplementary Tables S1-S3).

Validation of Model Predictions
FBA (see Section 2) was used to assess the predictive capability of our modified GSMM model to simulate yeast fermentation behavior under two relevant growth conditions for the industrially relevant strain EC1118 in a chemostat, Cases I and II as defined in the Section 2 (for more detail see Supplementary Table S4). These growth conditions likely mimic a stage of a wine fermentation corresponding to the end of exponential growth phase for each case. The use of a chemostat facilitated direct measurement of the growth and consumption/production rates of key metabolites as they are constant over time. Since the maximum rate of production of key volatile higher alcohols [59] and other aroma precursors occurs just at the onset of stationary phase as cells are leaving exponential phase, these data were chosen for evaluating model predictions of aroma and other compounds produced by the yeast.
GSMM predictions are shown with the experimental data for the two enological fermentation cases ( Figure 2). Here it is illustrated that our model, using the set of constraints developed in this study, succeeded in predicting the formation of many key metabolites. Furthermore, our simulations predicted soundly for major metabolites as well as biomass growth in Case I while predicting somewhat less accurately for biomass growth in Case II and acetate formation for both cases. More specifically, the biomass growth rate was predicted to be 24% less than the experimentally determined biomass growth rate in Case II and 10% less than the experimentally determined biomass growth rate in Case I. Despite underestimating biomass growth, we were able to successfully predict metabolites formation for both cases. More specifically, ethanol, carbon dioxide, and glycerol were predicted well, while carbon dioxide production was predicted to be 29% higher than the experimental flux for Case II. The model predicted relatively successfully the formation of some metabolites including the production of volatile aroma metabolites. While for both cases the model predicted production of organic acids (succinate and lactate) close to experimental production rates, acetate production in both cases was significantly under-predicted. Applying FBA to our model, we effectively predicted the production rates of higher alcohols, particularly; 1-propanol, which was added to the GSMM, is predicted rather closely to experimental production rates for both conditions. Furthermore, the simulated production rates of other higher alcohols (3-methylbutanol (isoamyl alcohol), isobutanol, 2-phenyl ethanol, and 2-methyl-1-butanol) show considerable agreement with their experimental production rates for Case I, although there is a slight over prediction for many of these higher alcohols for Case II. GSMM predictions for Case III and Case IV are shown in Figure S1 of the Supplementary Material. as they are constant over time. Since the maximum rate of production of key volatile higher alcohols [59] and other aroma precursors occurs just at the onset of stationary phase as cells are leaving exponential phase, these data were chosen for evaluating model predictions of aroma and other compounds produced by the yeast.
GSMM predictions are shown with the experimental data for the two enological fermentation cases (Figure 2). Here it is illustrated that our model, using the set of constraints developed in this study, succeeded in predicting the formation of many key metabolites. Furthermore, our simulations predicted soundly for major metabolites as well as biomass growth in Case I while predicting somewhat less accurately for biomass growth in Case II and acetate formation for both cases. More specifically, the biomass growth rate was predicted to be 24% less than the experimentally determined biomass growth rate in Case II and 10% less than the experimentally determined biomass growth rate in Case I. Despite underestimating biomass growth, we were able to successfully predict metabolites formation for both cases. More specifically, ethanol, carbon dioxide, and glycerol were predicted well, while carbon dioxide production was predicted to be 29% higher than the experimental flux for Case II. The model predicted relatively successfully the formation of some metabolites including the production of volatile aroma metabolites. While for both cases the model predicted production of organic acids (succinate and lactate) close to experimental production rates, acetate production in both cases was significantly under-predicted. Applying FBA to our model, we effectively predicted the production rates of higher alcohols, particularly; 1-propanol, which was added to the GSMM, is predicted rather closely to experimental production rates for both conditions. Furthermore, the simulated production rates of other higher alcohols (3-methylbutanol (isoamyl alcohol), isobutanol, 2-phenyl ethanol, and 2-methyl-1-butanol) show considerable agreement with their experimental production rates for Case I, although there is a slight over prediction for many of these higher alcohols for Case II. GSMM predictions for Case III and Case IV are shown in Figure S1 of the Supplementary Material.

Analysis of Active Aroma Reactions in Metabolic Network
After validating the model with existing data, it was then possible to use the model to understand the differences in metabolic states of strain EC1118 between Cases I and II. Metabolic differences are visually depicted in network maps of our two phenotypes from FBA simulations for core carbon metabolism (Figure 3), as well as sulfur and amino acid degradation metabolism ( Figure 4). Additionally, more complete and integrated networks featuring metabolic flux distributions within central carbon metabolism and amino acid metabolism are shown in the Supplementary Material (see Figures S3 and S4). For central carbon metabolism the most apparent differences are concerning the pyruvate-oxaloacetate node which is a major branch that functions as a linking point between glycolysis, gluconeogenesis, and the TCA cycle. Our model predicted that there is a difference in the distribution of fluxes regarding the oxaloacetate shuttle where, for Case I,

Analysis of Active Aroma Reactions in Metabolic Network
After validating the model with existing data, it was then possible to use the model to understand the differences in metabolic states of strain EC1118 between Cases I and II. Metabolic differences are visually depicted in network maps of our two phenotypes from FBA simulations for core carbon metabolism (Figure 3), as well as sulfur and amino acid degradation metabolism ( Figure 4). Additionally, more complete and integrated networks featuring metabolic flux distributions within central carbon metabolism and amino acid metabolism are shown in the Supplementary Material (see Figures S3 and S4). For central carbon metabolism the most apparent differences are concerning the pyruvate-oxaloacetate node which is a major branch that functions as a linking point between glycolysis, gluconeogenesis, and the TCA cycle. Our model predicted that there is a difference in the distribution of fluxes regarding the oxaloacetate shuttle where, for Case I, there is a preference for αketoglutaric acid formation, whereas for Case II, there is a preference for the formation of citric acid. The model also predicts differences in the manner in which acetic acid forms (Figure 3). For Case I, the simulated metabolism favors using both aldehyde dehydrogenase and acetyl CoA synthase to produce acetate. However, for Case II, the simulated metabolic flux distribution prefers solely the aldehyde dehydrogenase route.    For amino acid metabolism, there are apparent differences in simulated flux distributions when observing utilization pathways of leucine, valine, and tryptophan ( Figure 4). There is a slight difference in flux distribution between our two cases where simulated routes of leucine and valine in Case I chose a different path using IPPSm (mitochondrial isopropylmalate synthase) instead of using IPPS (cytosolic 2 isopropylmalate synthase) to form α-ketoisovaleric acid. This, in turn, impacts the predicted formation of fusel alcohols such as isobutanol. For tryptophan metabolism in Case II, the model predicts that there is flux to the production of indole-3-aldehyde, which is a precursor to tryptophol, in addition to the flux to tryptophan. The model predicts no flux being diverted to tryptophol for Case I. However, overall, there was not much difference in flux distribution of many of amino acid pathways responsible the fusel alcohols alcohol between the two cases, although the fluxes are generally higher in Case II by a factor of two.
When examining sulfur metabolism, there is a stark contrast in predicted flux distribution between the two cases, especially when examining the flow and utilization of sulfate and the sulfur-containing amino acids cysteine and methionine, which are highly involved in sulfur metabolism. In Case I, it was predicted that sulfate is further metabolized and then coupled with serine metabolism to form cysteine, whereas in Case II, flux to serine is higher and sulfate is not utilized. In addition, the model predicts fluxes directed toward formation of undesirable sulfur compounds such as hydrogen sulfide and methanethiol for Case I. However, the model predicts formation of methanethiol only for Case II.

Exploring Network Flexibility via FVA
In order to be able to modulate aroma production in yeast fermentations, it is important to know which enzymes and reactions are essential and contribute most to their production. Therefore, we applied flux variability analysis (FVA) to identify these enzymes and reactions. In brief, the method computes the minimal and maximal flux values that a metabolic reaction can take without affecting the cellular objective(s). Total fusel alcohol production was chosen to be maximized here, as an example, since these aroma molecules are formed in the highest quantity (excluding ethanol) during yeast fermentation and serve as precursors to other aroma molecules such as acetate esters. To determine the flux range, we first computed the production of the maximal total fusel alcohol production using FBA. Next, the total fusel alcohol production was fixed at the value obtained from FBA and used to run FVA. FVA was performed for two scenarios: (1) the model was constrained based on experimental data on uptake and secretion rates under Case I conditions from literature [30] and (2) the model was unconstrained mimicking a "super-rich" medium (see Section 2.3). The resulting flux ranges from FVA for every reaction can be divided into three groups: (I) essential reactions, where the minimal and maximal flux is non-zero; (II) non-essential reactions, where the minimal flux is zero and maximum value is non-zero; and (III) inactive reactions, where both minimum and maximum fluxes are zero.
There are some discernible contrasts in flux variability illustrated between the two scenarios regarding Group I aroma-associated reactions ( Figure 5). The total amount of Group I aroma-associated variable fluxes is different where there are 19 of these type of reactions in the Constrained Case ( Figure 5a) and 10 reactions in the Unconstrained Case (Figure 5b). Since the Unconstrained Case represents an extreme rich medium it is expected that in this case the number of essential reactions is lower. However, 10 of the essential reactions overlap for the two cases. For the Constrained Case, all of these Group I essential aroma reactions contained a zero to small flux range (<0.2687 mmol gDW −1 h −1 ) ( Figure 5). In fact, 14 out of 19 reactions have a zero-flux range representing essential reactions that when slightly changed, the objective, i.e., aroma production, will immediately change as well.

Discussion
Studies have pressed the need for understanding how altering the fermentation environment induces changes in secondary metabolism, especially when impacting organoleptic properties of products [10]. This illustrates the demand for suitable tools to aid in understanding the role of yeast metabolism in aroma formation especially with regard to higher alcohol and ester synthesis. We improved the GSMM for S. cerevisiae for studying wine fermentation by including aroma forming pathways previously not contained in existing S. cerevisiae GSMMs. This expanded GSMM permitted the comparison of yeast metabolism, specifically aroma formation under varying

Discussion
Studies have pressed the need for understanding how altering the fermentation environment induces changes in secondary metabolism, especially when impacting organoleptic properties of products [10]. This illustrates the demand for suitable tools to aid in understanding the role of yeast metabolism in aroma formation especially with regard to higher alcohol and ester synthesis. We improved the GSMM for S. cerevisiae for studying wine fermentation by including aroma forming pathways previously not contained in existing S. cerevisiae GSMMs. This expanded GSMM permitted the comparison of yeast metabolism, specifically aroma formation under varying fermentation conditions. While this is an important advance over previous models, further curation and adaptation may be necessary to exhaustively predict yeast-derived aroma compounds, such as those derived from phenolics, which can sometimes be important to flavor. Despite this, one can apply approaches found in this study to expand GSMMs in order to investigate aroma in wine or other alcoholic fermentations. It was also noted here that this model still requires small fluxes of oxygen into the cells in order to predict cell growth and fermentation. Since other GSMMs [42,44,60] have been developed where this is not necessary, it would be useful to continue development on this model to incorporate this behavior, along with the ability to accurately predict yeast-derived aroma production.
Using FBA, our expanded GSMM generated relatively accurate predictions for the secondary metabolite formation including the production of volatile aroma metabolites. The model did, however, underestimate the growth rate in both Case I and II (10% in Case I and 24% in Case II). This could be due to several factors. First of all, the model biomass equation may not reflect the most precise medium conditions used in the experiment [61]. Second, there is genetic variation between wine yeast strains that affect the utilization of nitrogen which is needed for biomass growth [62]. This highlights the need for more environmental, as well as genetically specific biomass equations, including one for wine applications. Co-factor balance is confirmed based on our observation that the model with the given constraints predicts growth of biomass. Without co-factor balance, growth cannot be predicted. While the model predicted production of organic acids (succinate and lactate) close to experimental production rates for both simulated cases, acetate production in both cases was under-predicted. Acetate production is known to be coupled to additional ATP generation within yeast cells. This only happens if the redox balance is restored somewhere else in the metabolic network. In addition, the under-prediction of the acetate flux could be related to an under estimation of the maintenance coefficient, reaction-ATPM (see GSMM). This could point to competition with pyruvate dehydrogenase linked reactions within the network and respiratory pyruvate dehydrogenase (PDC5) linked reactions where PDCs facilitate more pyruvate from the pathway than necessary diverting it toward ethanol [63]. In the study by Quiros et al., the flux toward oxaloacetic acid production (OAA) was not significantly impacted by the variables investigated in the study. However, our results indicate the routes of OAA production shifted from mitochondrial production to cytosolic when going from Case I to Case II. This is interesting as it points to environmental effects (e.g., temperature) influencing the mechanism of acetate production as has been suggested in other studies [64]. High osmolarity resulting from high sugar concentrations in the feed induce glycerol production [65,66]. In addition, it is clear from these results that the model is predicting more carbon flux directed towards ethanol and glycerol production than to cell growth compared to experimental observations, which may indicate a need for further model curation.
Our GSMM describes production of volatile aroma compounds by linking aroma formation pathways to central carbon metabolism where precursors are synthesized. The production of fusel alcohols (under anaerobic conditions) can and will be used as extra capacity to regenerate NAD+. This results in adding more flexibility in the network in terms of redox co-factor regeneration. Moreover, it has been proposed that ester production enables control of intracellular redox balance and some esters are involved in the preservation of plasma membrane fluidity [67,68]. Intracellular esterification of MCFAs could also help yeast cells to get rid of toxic MCFAs through diffusion out of the plasma membrane [69]. In any case, extracellular ester production is believed to play a significant evolutionary role in leading to spread of yeast throughout the environment as esters attract insects by indicating the presence of rotting fruits [70]. Our expanded GSMM was able to predict the formation of some these desirable esters (isoamyl acetate, isobutyl acetate, and 2-methylbutyl acetate) (data not shown), where, under the Case I, there was more predicted ester production than Case II. This result is supported by studies that have indicated higher alcohol and consequently higher acetate production under higher temperature conditions [71,72]. Moreover, Saerens et al. [22] showed that alcohol O-acetyltransferase 1 (ATF1) and alcohol O-acetyltransferase 1 (ATF2) expression is positively correlated with temperature and thus would boost acetate ester production under certain temperature conditions. Production of sulfur compounds from sulfur reduction pathways is also represented in our model simulations ( Supplementary Figures S1 and S2). There are some contrasts in flux distribution throughout sulfur metabolism when comparing the flux distributions found for our two cases. This can be due to the Case I having methionine or cysteine-limited conditions, thus leading to O-acetylserine and O-acetylhomoserine limitations, which results in excess sulfide that is diverted to hydrogen sulfide [73]. On the other hand, Case II could contain non-limiting cysteine or methionine concentrations which would lead to sulfide reacting with O-acetylserine or O-acetylhomoserine to form homocysteine [5].
In order to investigate the metabolic capability and the potential for suggesting metabolic modification, the solution space of our expanded GSMM was evaluated using FVA. The flexibility of the GSMM network was explored for a constrained case and an unconstrained case. When examining the unconstrained case, it is observed that there are three Group I essential aroma reactions; ALCD5xi, 3MOBDC, and PYRDC4 which have zero flux difference. These reactions are also essential with zero flux range in the constrained case, suggesting unconditional essentiality with zero flux spans. Moreover, in the genes associated with the three most essential reactions are PDC1, PDC5, PDC6, ADH1, ADH2, ADH4, ADH5, and SFA1. Many of the Group I aroma-associated reactions are closely connected to the Ehrlich pathway and several are reactions we added to the GSMM e.g., aldehyde dehydrogenases (via propanol and tyrosol). Based on these results, we propose that Group I reactions with zero flux span are of particular interest for metabolic engineering of aroma production, since these reactions will change aroma production once their fluxes are perturbed even slightly. This applies especially to the unconditionally essential ones with zero flux spans as these are environmentally independent (i.e., both present in the constrained case and the unconstrained case). In addition, our results are supported by studies that have increased higher alcohol production by modulating alcohol dehydrogenases (ADH1, ADH2) to increase higher alcohol production [74,75]. Taken together, we can conclude that, under the given experimental conditions, Group I aroma-associated reactions can serve as potential candidates for steering aroma profiles.
Even with some of these limitations, we have shown that this GSMM incorporating aroma production pathways is able to predict aroma compound production under different growth conditions and helps to explain the metabolic fluxes underlying these differences. These predictions were made using constraints relevant to winemaking or enological conditions, including mostly anaerobic and nitrogen-limited conditions. With this model, it should be possible now to predict and understand aroma metabolism in yeast over the course of an entire wine fermentation and as a function of different commercial yeast strains, thus creating a tool to modify fermentation behavior through culture conditions and strain selection or development. Although this study was specific to wine aroma formation, it is expected that the work presented here can be applied to understand and study aroma formation in other types of alcoholic fermentations as well.

Conclusions
In this study, a GSMM for S. cerevisiae was expanded by adding metabolic pathways that lead to the formation of various aroma compounds. Many aroma compounds producing pathways particularly related to ester and sulfur formation were incorporated into a widely used GSMM for S. cerevisiae. This study is distinctive from previous aforementioned studies in that it applies a GSMM to comprehensively study yeast metabolic impacts on aroma formation under enological conditions. This expanded model was used for examining aroma formation under enological conditions with input of published experimental data for two growth conditions. By applying FBA, fermentation performance such as biomass growth rates and product secretion rates were predicted and compared favorably to experimental values for a commercial yeast strain. Furthermore, we applied the extended model to compare key differences in metabolic flux distributions between the two growth condition cases. Finally, we predicted key essential reactions that are of particular interest for metabolic engineering purposes aiming to modulate aroma production. Therefore, we have illustrated how this newly expanded GSMM can be a helpful tool for winemakers and researchers alike that are interested in understanding aroma formation during wine formation.
Supplementary Materials: The following are available online at http://www.mdpi.com/2227-9717/8/9/1195/s1, Figure S1. Comparison of GSMM simulations (Red bars) with experimental data [30] (Blue bars). Model predicted and experimental values for biomass growth rates, primary metabolite production fluxes, and secondary metabolite production fluxes represent conditions for Case III and Case IV. Experimental data are shown with standard deviation error bars; Figure S2. Comprehensive metabolic network map of yeast (S. cerevisiae); Figure S3. Network map of central carbon metabolism (left) and amino acid metabolism (right). Network fluxes were determined using data from Case I where the network illustrates the optimal solution from FBA. The fluxes were scaled between −1 and +1 for ease of visualization; Figure S4. Network map of central carbon metabolism (left) and amino acid metabolism (right). Network fluxes were determined using data from Case II where the network illustrates the optimal solution from FBA. The fluxes were scaled between −1 and +1 for ease of visualization; Table S1. Genes added to current GSMM; Table S2. Metabolites added to current GSMM; Table S3. Reactions added to current GSMM; Table S4. Experimentally measured consumption/production rates of the different metabolites analyzed from Quiros et al. 2013 and constraints used in the study.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.