One of the most obvious phenotypes of a cell is its metabolic activity, which is defined by the fluxes in the metabolic network. Although experimental methods to determine intracellular fluxes are well established, only a limited number of fluxes can be resolved. Especially in eukaryotes such as yeast, compartmentalization and the existence of many parallel routes render exact flux analysis impossible using current methods. To gain more insight into the metabolic operation of S. cerevisiae we developed a new computational approach where we characterize the flux solution space by determining elementary flux modes (EFMs) that are subsequently classified as thermodynamically feasible or infeasible on the basis of experimental metabolome data. This allows us to provably rule out the contribution of certain EFMs to the in vivo flux distribution. From the 71 million EFMs in a medium size metabolic network of S. cerevisiae, we classified 54% as thermodynamically feasible. By comparing the thermodynamically feasible and infeasible EFMs, we could identify reaction combinations that span the cytosol and mitochondrion and, as a system, cannot operate under the investigated glucose batch conditions. Besides conclusions on single reactions, we found that thermodynamic constraints prevent the import of redox cofactor equivalents into the mitochondrion due to limits on compartmental cofactor concentrations. Our novel approach of incorporating quantitative metabolite concentrations into the analysis of the space of all stoichiometrically feasible flux distributions allows generating new insights into the system-level operation of the intracellular fluxes without making assumptions on metabolic objectives of the cell.
Metabolic fluxes give immediate insights into the metabolism of a cell [1], [2]. Metabolic flux analysis has proven to be useful, for example for the determination of enzyme functions [3], for the identification of regulatory mechanisms in response to environmental perturbations [4], or as a tool in metabolic engineering [5]. The most common method to quantify metabolic fluxes uses labeled substrates, and the measured label distribution in intracellular metabolites is interpreted together with measured uptake and production rates by means of a metabolic network model [6]. Despite successful quantification of fluxes with flux analysis in different conditions, the method has several limitations. For example, it is limited to the main branches in central carbon metabolism and fluxes cannot be resolved per compartment [7], despite compartmentation being a highly relevant aspect of eukaryotes [8]. Moreover, today' s flux analysis rests on a number of a priori assumptions, e. g. on reaction reversibilities or on relevant parts of the network [7], [9]. To improve flux quantification, we need additional constraints on the possible flux distributions in a metabolic network. In stoichiometric network analysis, the metabolic network is modeled as a collection of biochemical reactions where all internal metabolite concentrations are assumed constant [10]. Next to the typical constraints, such as uptake and excretion rates, reaction reversibilities and maximum flux capacities, the field recently began to incorporate thermodynamic information, whereby statements on feasibility of reaction fluxes or flux distributions can be made based on calculation of changes in Gibbs energy using metabolite concentrations [11]–[13]. For example, using flux balance analysis (FBA) and related approaches, metabolite concentrations were used as additional constraints to predict fluxes in the non-compartmentalized organism E. coli [14], [15] or in a model of liver metabolism [16]. Here, we develop a novel approach to integrate metabolite data into metabolic network flux analysis, to get additional insight into the compartmentalized flux physiology of Saccharomyces cerevisiae. The method combines network embedded thermodynamic (NET) analysis [12], elementary flux mode (EFM) analysis [17]–[19], and experimentally determined metabolome data. We employ EFM analysis instead of flux balance analysis because the collection of the generated flux modes can yield insight into all feasible flux distributions, as compared to the single thermodynamically feasible flux solution that is obtained with thermodynamically constrained FBA. Additionally, assumptions on a metabolic objective function of the cell are not required. Our new approach to analyze the compartmentalized central metabolic network of S. cerevisiae using quantitative metabolite data acquired under glucose batch growth conditions allowed us to generate novel insight into the system-level causalities underlying the intercompartmental redox metabolism. Specifically we show that the and NADH concentrations in the cytosol and the mitochondrion do not allow for the ethanol-acetaldehyde redox shuttle to be active under the investigated condition. Further, we identified a number of maximal reaction activities that could be used as constraints for flux analysis or FBA. We envision that our method becomes a useful tool to unravel system-level insights about a complex metabolic system from metabolome data.
The mass balanced flux solution space of a stoichiometric metabolic network can be described with a non-negative linear combination of its EFMs: (1) where any flux distribution is a sum of EFMs with coefficients. As we show in the proof provided in Text S1, a thermodynamically feasible flux distribution only consists of thermodynamically feasible EFMs: (2) where the thermodynamically feasible flux distribution is only composed of EFMs from the feasible set, . The mathematical proof demonstrates that by eliminating infeasible EFMs we do not loose feasible flux distributions, because any feasible flux distribution can be composed of only feasible EFMs. Specifically, in the hypothetical case that an infeasible EFM is part of a feasible flux distribution it must involve a cancellation or directionality change of a specific reversible reaction. In this case, the flux distribution can be decomposed into one or more feasible EFMs, and a feasible combination of an infeasible EFM with another EFM. The combination of the infeasible EFM and another EFM must then be either a feasible EFM by itself, or it must be possible to achieve the feasible combination by other feasible EFMs. Hence, Eq. (2) allows us to exclude thermodynamically infeasible EFMs from the complete set of EFMs without excluding thermodynamically feasible flux distributions. It is important to note that the resulting flux solution space defined by the feasible EFMs can still contain infeasible flux solutions because it is possible that the combination of multiple feasible EFMs leads to an infeasible flux distribution. Exploiting that EFMs allow us to exclude thermodynamically infeasible EFMs, we aimed at developing an approach to generate novel insights into the complex flux physiology of the central metabolism of the yeast S. cerevisiae. Therefore, we assembled a 230 reaction stoichiometric network of its central carbon metabolism and amino acid synthesis pathways (cf. Materials and Methods) encompassing the cytosolic and mitochondrial compartments and many parallel pathways. With our approach we aim to obtain additional insight into the metabolic network operation, therefore we build upon current knowledge by defining the reversibilities in our model as they are defined in the original model [21]. First, we needed to calculate all EFMs. As the EFM calculation is computationally demanding, we initially applied steps to constrain the mass-balanced solution space as much as possible upfront, before we started enumerating EFMs (Fig. 1). Thus, in a first step, we performed flux variability analysis (FVA) [22] on the basis of the measured uptake and production rates of external compounds and biomass, to determine the reversibility of each reaction under the investigated physiological conditions, that is, for growth on glucose. From FVA, we obtain a minimum and maximum achievable flux for each reaction. A reaction is reversible if both a negative and positive flux can be achieved, else it is unidirectional. Reactions that have a minimum and maximum flux that are either positive or negative, and cannot be inactive, are reactions that are always active in the respective direction. Using this approach, we could classify 67 initially reversible reactions as unidirectional (Fig. 2, Dataset S1). Next, we employed measured metabolite concentrations from glucose batch cultures, and NET analysis to identify additional reaction irreversibilities [12], [23]. For the metabolite data, we assembled published and unpublished data from glucose batch experiments, and generated a consensus data set to define lower and upper concentration limits for 55 metabolites (see Materials and Methods and Dataset S2). For NET analysis, we used the reaction activities inferred from FVA. With the consensus metabolite data set, we obtained constraints on the reversibility of three additional reactions (see Fig. 2). Another iteration of FVA and NET analysis with the obtained constraints as input did not yield any further constraints. For the obtained condition-specific constrained metabolic network, we computed the EFMs and obtained 71. 266. 960 EFMs. Using NET analysis and the consensus metabolite data set, we classified 38. 420. 207 (54%) EFMs as feasible and 32. 846. 753 (46%) EFMs as infeasible. Assuming that the flux solution space in the metabolic network of an organism can be approximated by the number of EFMs, our result shows that roughly at most half of the solution space is thermodynamically feasible. With the finding that 54% of all the EFMs are thermodynamically feasible we reduced the number of EFMs that can constitute a thermodynamically feasible flux distribution considerably. In a first analysis step towards generating insights into the flux distribution, we searched the feasible EFM set for reactions that only use a subset of the possible reaction directions compared to the complete set of EFMs. For each reaction in each EFM we determined whether a backward, inactive or forward reaction activity was used. Then, the possible reaction activities for each reaction were compared between the complete set of EFMs and the feasible set of EFMs. Here, we found that the oxaloacetate transport from the mitochondrion to the cytosol is never used in an EFM of the feasible set, meaning that it has to be inactive under the investigated growth condition. Indeed, we find no contradicting evidence for the prediction when comparing with the experimental observation that a knock-out of the corresponding gene OAC1, whose translated protein is responsible for the respective oxaloacetate transport reaction, does not have an effect on growth rate under glucose batch conditions [24], [25]. Further, as we found that all EFMs with acetaldehyde transport out of the mitochondrion are infeasible, we conclude that during growth on glucose, acetaldehyde can only be transported into the mitochondrion. It is important to note that the EFMs with active oxaloacetate transport, or acetaldehyde transport out of the mitochondrion, are not infeasible because of the metabolite concentration constraints on the respective single reaction only, since single reaction infeasibilities are removed in the first NET analysis step before EFM generation (see Fig. 1). Instead, as we will show later, the infeasibility is the result of a system of coupled reaction activities, where all individual reactions need to be thermodynamically feasible simultaneously. Next, we wanted to test whether the feasible set of EFMs differs from the infeasible set in terms of reaction rates. Such a comparison is possible by normalizing the reaction rates in each EFM to the glucose uptake rate of the EFM (all EFMs have glucose uptake, EFMs without glucose uptake are internal cycles and they were removed because they are physiologically meaningless [26]). Principal component analysis (PCA) of the normalized EFMs shows a clear difference between the feasible and infeasible EFMs in principal component 2 (PC 2 in Fig. 3). The reactions with the highest loadings in this component are alcohol dehydrogenase in both the cytosol and the mitochondrion (ALCD2x, ALCD2m) and acetaldehyde and ethanol transport to the mitochondrion (ACALDtm, ETOHtm), and these reactions are likely involved in causing thermodynamic infeasibilities. Note, that although there is a separation between feasible and infeasible EFMs because there are no feasible EFMs in the center of the graph, by combination of feasible EFMs it could still be possible to obtain a feasible flux distribution that would be projected in this area of the PCA. Therefore, the loadings of PC1 are not considered. Next, we searched for the highest and lowest rate of each reaction in the complete set of EFMs and in the feasible set of EFMs to define the flux ranges that can be achieved in terms of flux per unit of glucose uptake. Any flux value in this range can in principle be achieved through a combination of EFMs. When comparing the flux ranges that can be realized by the feasible EFMs with the flux ranges of the complete set of EFMs, we find that eight reactions cannot assume the full range for thermodynamic reasons (see Fig. 4), with four of these reactions already having shown high loadings in the second principal component (see Fig. 3). These quantitative flux constraints result from metabolite concentrations and thermodynamics, and can be applied as constraints in flux balance analysis.
In this work, we demonstrated that an in vivo thermodynamically feasible metabolic flux distribution is only composed of thermodynamically feasible elementary flux modes. This EFM property allowed us to integrate EFMs and NET analysis into a novel approach to study the system-level properties of complex metabolic networks on the basis of quantitative metabolome data. As exemplified with a compartmentalized model of central metabolism in S. cerevisiae and cell-averaged metabolome data generated under glucose batch conditions, 46% of the 71. 3 million EFMs were found thermodynamically infeasible, leading to direct insights into reaction directionalities, to constraints in several metabolic rates, and to the identification of reaction patterns that must be inactive due to a thermodynamic infeasibility. This work builds on earlier work that integrated FBA, thermodynamics and quantitative metabolite data [14], [15], and extends it by using EFMs, allowing us to identify the reasons underlying the infeasibilities without making a priori assumptions on the metabolic objectives of the cell, such as maximization of biomass production, as is the case with FBA. Additionally, in our study we considered a compartmented metabolic network of S. cerevisiae to analyze cell-averaged metabolite data. Notably, the results we obtain are directly related to compartmentation, as can be seen from the identified infeasibility patterns that involve both compartments. The identified infeasibility of the ethanol-acetaldehyde redox shuttle has been previously identified using NET analysis [12] and manual consideration of the system. In this work we demonstrate that by using the flux patterns that are obtained from the EFMs we systematically identify such infeasibility patterns. Although the system-level constraints and their underlying causes can be rationalized without using EFMs, we need the generated EFMs to determine the infeasible patterns that are part of a stoichiometrically balanced flux distribution. In addition, because the number of EFMs can be considered approximately proportional to the flux solution space, we find that roughly half of the flux solution space is thermodynamically infeasible due to systems of reaction activities. With the recently developed new group contribution method to estimate thermodynamic properties on a genome-scale [35], the recently increased availability in thermodynamic properties through experimental methods [36], the advances in quantitative metabolomics [37] and the now available methods to calculate EFMs also for large stoichiometric network [38], [39], we envision that the here presented approach will be helpful to shed light on metabolic flux physiology in more complex metabolic system such as higher cells simultaneously growing on multiple carbon substrates, where the applicability of classical flux analysis methods are still rather limited.
We used experimental data on metabolite concentrations for Saccharomyces cerevisiae obtained from four independent experiments with at least two replicates [40]–[42] with equal growth medium but under different cultivation conditions (bioreactor, shake flask, 96-well). Based on the data from the independent experiments, we constructed a consensus data set, where for each metabolite a minimum and maximum concentration was defined. The minimum and maximum concentrations were determined from all the replicates of measurements for each metabolite. To reduce the effect of outliers on the ranges, when more than 3 replicates were available, we removed the values higher than the third quartile +1. 5 IQR (inter quartile range), and values lower than the first quartile −1. 5 IQR. Physiological data was obtained for S. cerevisiae on glucose as carbon source from one of the four experiments. In Dataset S2 we describe the details of the experimental conditions of the data sets, the obtained concentration ranges and physiological data. The stoichiometric metabolic network model describes the core central carbon metabolism of S. cerevisiae in the cytosol and mitochondrion with 230 reactions and 218 metabolites (see Dataset S3), and was developed on the basis of the genome-scale metabolic model iND750 that contains 1149 reactions and 646 metabolites [21]. For our model, we selected the cytosolic and mitochondrial reactions belonging to glycolysis/gluconeogenesis, pentose-phosphate pathway (PPP), TCA cycle, anaplerosis, pyruvate metabolism, and oxidative phosphorylation. The reversibility of each reaction was taken from Duarte et al. [21]. A cytosolic malate synthase was added to complement the glyoxylate shunt in the cytosol [43]. A citrate synthase was added to the cytosol since this is supported by localization studies [44]. To allow the model to synthesize all amino acids that are required for biomass, we added the following pathways: For L-alanine, two biosynthetic routes from pyruvate were included: cytosolic and mitochondrial alanine transaminase reactions, which were assumed to solely produce but not degrade L-alanine [45], [46]. Furthermore, L-glutamate could be produced via three alternative pathways: cytosolic or mitochondrial NADP-dependent glutamate dehydrogenase from alpha-ketoglutarate or mitochondrial NAD-dependent glutamate synthase from alpha-ketoglutarate and glutamine [47]. For glycine synthesis, we implemented three pathways such that it could be synthesized in the mitochondria via (i) alanine-glyoxylate transaminase [48], or in the cytosol by (ii) glycine hydroxymethyltransferase from serine [49], or (iii) from L-threonine via threonine aldolase [50]. As the latter reaction was assumed to be reversible, it could also be used to produce L-threonine, and such it constitutes a second possibility to produce L-threonine next to the linear pathway from L-aspartate. For all other amino acids, the model contains only one linear cytosolic pathway consisting of consecutive enzymatic reaction steps. Here, no alternative paths exist or they are excluded based on biochemical literature as it was done also to construct models for -based flux analysis [45], [46]. The model further includes transport reactions across the mitochondrial membrane for metabolites that participate in reactions in both the cytosol and the mitochondria. Additional transport reactions that were not contained in iND750 (i. e. for L-glutamate, alpha-ketoglutarate, homocitrate, glyoxylate, and 2-oxobutanoate) were added to properly connect additionally included alternative pathways for amino acid synthesis to the metabolic network. The biomass composition was adopted from iND750 besides that trehalose and glycogen were discarded since carbohydrate storage was not considered in our model. Lumped reactions for synthesis of the remaining biomass constituents, i. e. lipids, nucleotides, and cell wall components from the corresponding precursors were determined based on the biomass composition as provided in iND750. In the model, carbon molecules that can be exchanged with the environment are glucose, glycerol, pyruvate, acetate, ethanol, succinate, and. The model is not proton balanced. The reason for this is that it is close to impossible to do the proton balancing correctly (e. g. , for transport reactions). Thus, we did not want to add any potentially wrong constraints on the model and therefore did not account for proton balancing, with one exception. We only balanced the protons around the respiratory chain by replacing the cytosolic protons produced and consumed in the reactions CYOR_u6m, CYOOm and ATPS3m by a unique species “hcyt”. As a result, ATP generated in ATPS3m can only occur through the respiratory chain. With NET analysis we can determine the feasibility of a flux distribution based on ranges for the concentrations of the involved metabolites. A flux distribution is thermodynamically infeasible when one or more reaction activities conflict with the calculated Gibbs energy of reaction range (s). Conversely, a flux distribution is feasible when no conflicts are found. It is important to note that a metabolite concentration is constrained in NET analysis by any reaction that has the metabolite as a reactant. Therefore, a flux distribution can be infeasible due to propagated constraints in a pathway. The NET analysis implementation constrains metabolite concentrations of metabolites that occur in multiple compartments as a sum of the compartment specific concentrations, corrected for their volume. Therefore, the compartmental distribution of such metabolites is left free. The compartmental volume fractions of the cytosol and mitochondrion are set to 0. 35 and 0. 1, respectively. The NET analysis approach is similar to other thermodynamic analysis approaches [14], [15]. A main difference from other approaches is that with NET analysis we aim at checking flux distributions for thermodynamic feasibility, and at estimating ranges of Gibbs reaction energies and metabolite concentrations. For NET analysis we used the concentration ranges defined from the experimental data. For all other metabolites in the network we assumed a default range with a minimum concentration of 0. 0001 mM and a maximum of 120 mM, except for carbon dioxide (“co2tot”), phosphate (“pi”) and diphosphate (“ppi”) that were constrained to a range of 1 mM to 100 mM [51], [52], and oxygen (“o2”) that was constrained to 0. 001 mM to 0. 1 mM [53]. By using such large metabolite concentration ranges we account for the noise in the metabolite concentration data and uncertainties in Gibbs energies of formation. Typical uncertainties in formation energies are in the order of 0. 02–2 kJ/mol [35], which are overshadowed by variations in metabolite concentration data. The compartmental pH values were set to 5,6. 5 and 7 for the external, cytosolic and mitochondrial environment, respectively [54], [55]. The ionic strength was set to 0. 15 M for all compartments. For the correct consideration of transport thermodynamics in NET analysis, we defined the specific transported species for each transport reaction where possible, and calculated transport reaction values according to Jol et al. [56]. Computations for FVA, NET analysis and EFM generation were done using MATLAB (The Mathworks). For optimization of FVA problems we used the LINDO API library (LINDO Systems Inc.) and for NET analysis we used anNET [23] in combination with the LINDO global solver. To generate EFMs we used the Java implementation from Terzer and Stelling [19] on a quad-core system (3 GHz) with 128 GB memory. To test the thermodynamic feasibility of each EFM we used anNET, which was modified to run in an automated way on a cluster of computers encompassing on average 60 CPUs (3 GHz). Testing EFMs for thermodynamic feasibility was computationally intensive and took approximately 14 days. To find the reaction activity patterns that cause infeasibility, we considered each infeasible EFM separately and performed an iterative analysis. In NET analysis of the EFM, we removed consecutively each reaction' s activity constraint from the NET analysis optimization and determined the feasibility. If the activity pattern became thermodynamically feasible, the reaction activity that was removed was identified as part of the pattern. Then we continued with removing the next reaction activity, while keeping the activities that were identified as part of the pattern. We continued this process for all the reaction activities. This process led to a set of reaction activities, which is a subset of the activities in the analyzed EFM, of which the removal of one activity leads to a feasible system. The set of reaction activities is only infeasible as a whole system, and no single reaction can be marked infeasible by itself. All possible infeasible reaction patterns may not be found when multiple patterns are present in an EFM, because the order of reaction activity removal determines which pattern is found. The obtained infeasible patterns cover all the infeasible EFMs. Flux balance analysis with maximization of biomass production, maximization of ATP production and maximization of the ratio of ATP production over the sum of squared fluxes was performed according to Schütz et al. [28] with the physiological data used for FVA as constraints. The computations were done using MATLAB (The Mathworks) using the LINDO API library (LINDO Systems Inc.) for optimization.
Fluxes in metabolic pathways are a highly informative aspect of an organism' s phenotype. The experimental determination of such fluxes is well established and has proven very useful. To address some of the limitations of experimental flux analysis, such as when the cell is divided in multiple compartments, stoichiometric modeling provides a valuable addition. The approach that we take is based on stoichiometric modeling where we consider the thermodynamic feasibility of many different possible routes through the metabolic network of Saccharomyces cerevisiae using experimentally determined metabolite concentrations. We show that next to conclusions on single biochemical reactions in the metabolic network, we obtain system-level insights on thermodynamically infeasible flux patterns. We found that the compartmental concentrations of and NADH are the causes for the system-level infeasibilities. With the current advances in quantitative metabolomics and biochemical thermodynamics, we envision that the presented method will help gaining more insight into complex metabolic systems.
systems biology
model organisms
yeast and fungal models
biology
computational biology
microbiology
genetics and genomics
System-Level Insights into Yeast Metabolism by Thermodynamic Analysis of Elementary Flux Modes