Current limitations in quantitatively predicting biological behavior hinder our efforts to engineer biological systems to produce biofuels and other desired chemicals. Here, we present a new method for calculating metabolic fluxes, key targets in metabolic engineering, that incorporates data from 13C labeling experiments and genome-scale models. The data from 13C labeling experiments provide strong flux constraints that eliminate the need to assume an evolutionary optimization principle such as the growth rate optimization assumption used in Flux Balance Analysis (FBA). This effective constraining is achieved by making the simple but biologically relevant assumption that flux flows from core to peripheral metabolism and does not flow back. The new method is significantly more robust than FBA with respect to errors in genome-scale model reconstruction. Furthermore, it can provide a comprehensive picture of metabolite balancing and predictions for unmeasured extracellular fluxes as constrained by 13C labeling data. A comparison shows that the results of this new method are similar to those found through 13C Metabolic Flux Analysis (13C MFA) for central carbon metabolism but, additionally, it provides flux estimates for peripheral metabolism. The extra validation gained by matching 48 relative labeling measurements is used to identify where and why several existing COnstraint Based Reconstruction and Analysis (COBRA) flux prediction algorithms fail. We demonstrate how to use this knowledge to refine these methods and improve their predictive capabilities. This method provides a reliable base upon which to improve the design of biological systems.
Systems biology aims to understand and predict how a cell’s behavior emerges from the interaction of its molecular parts [1–3]. Determination of metabolic fluxes (i. e. , the number of metabolites traversing each biochemical reaction per unit time [4,5]) is crucial to this effort because they map how carbon and electrons flow through metabolism to enable cell function [2,5]. Metabolic fluxes typically cannot be measured directly, but must be inferred from experimental data through computational algorithms [4]. Among the most popular methods for studying metabolic fluxes are Metabolic Flux Analysis (MFA, [6,7]), Flux Balance Analysis (FBA, [8]) and 13C Metabolic Flux Analysis (13C MFA, [4,5, 9]). MFA calculates fluxes by using a stoichiometric model for the major intracellular reactions and assuming no metabolite accumulation [6]. The inputs are extracellular fluxes obtained through measurements of external concentrations of metabolites such as glucose or lactate as a function of time. If more flux measurements are available than degrees of freedom, the system is said to be overdetermined and a unique solution can be obtained. In the opposite case, the system is underdetermined and several flux profiles are compatible with the experimental data. MFA has been used to study fluxes in (e. g.) chinese hamster ovary [10], S. cerevisiae [11] and hybridoma cells [12]. Present-day FBA enhances MFA by expanding the network to include all reactions in metabolism, or at least as many as can be inferred from the genome through a metabolic reconstruction that yields a genome-scale stoichiometric model [13]. Since the degrees of freedom for such a model are usually over a hundred and measured fluxes are usually an order of magnitude less, the system is grossly underdetermined. Hence, fluxes are determined through linear programming (LP) by assuming that metabolism is tuned, due to evolutionary pressure, to maximize growth rate (typically; but see Schuetz et al. [14] for other suggested alternatives). The use of this objective function to interpret stoichiometric models is a key feature of FBA, even when stoichiometric models were not genome-scale. FBA forms the basis of a family of flux analysis methods named COnstraint-Based Reconstruction and Analysis (COBRA) methods [15], some of which can be used to produce flux predictions which are often used in bioengineering. These predictions can be full (when fluxes are determined without data from the actual experiment [16]) or partial (when some data from the experiment, like glucose consumption, is used for the prediction [17]), qualitative (e. g. , prediction of whether an organism will grow or not under given conditions [18]) or quantitative (e. g. , the value of the flux is predicted [19]). Full predictions are particularly useful for bioengineering purposes since they enable quick testing of the consequences of engineering approaches [20]. These methods have been used to facilitate the large-scale industrial production of 1,4-butanediol, a commodity chemical used to manufacture over 2. 5 million tons annually of high-value polymers [21]. Recently, this rationally developed strain was used for a 5 million pound commercial production [22], and BASF has licensed this strain for future production of renewable 1,4-butanediol [23]. FBA has applications beyond bioengineering, including the prediction of ratios of microbial species in a simple microbial community [24] and providing valuable insights into tumor cell metabolism [25,26]. 13C MFA improves on MFA by using data obtained from 13C substrate labeling experiments together with a limited reaction stoichiometry and measured extracellular fluxes to measure intracellular fluxes. In these experiments, the organism under study is grown on 13C labeled substrate and the labeling pattern (i. e. , the fraction of molecules with 0,1, 2, … 13C atoms incorporated or Mass Distribution Vector, MDV [27]) is measured for a set of metabolites. Since the labeling pattern is highly dependent on the flux profile, it is possible to back-calculate the fluxes that best explain the measured labeling pattern if we know the fate of each carbon atom (carbon transitions, see [9]) for all reactions in the model. This involves a nonlinear fitting problem where the fluxes are the parameters. The usual approach is to consider only a small subset of all metabolism assumed to be comprehensive enough to explain the labeling of the measured metabolites, typically central carbon metabolism [28,29]. Hence the model exhibits relatively few degrees of freedom that can be fully constrained by the labeling data. 13C MFA has found applications in metabolic engineering [30], biotechnology [31] and biomedicine [32]. Each of these methods involves its own advantages and disadvantages. Unlike 13C MFA, FBA uses the comprehensive description of metabolism contained in a genome-scale model and explicitly takes into account the system-wide balances of metabolites that can be crucial for host engineering [33]. Because of the exhaustive description of metabolism incorporated in genome-scale models, they often point towards completely unexpected regions of metabolism involved in the studied processes. For example, these models have been used to show that biosynthesis and degradation of heme compensates for the lack of a functional TCA cycles in cancer cells [26]. Furthermore, FBA can be used in combination with COBRA methods to make full predictions. 13C MFA, on the other hand, is a descriptive method for determining the metabolic fluxes compatible with the accrued experimental data but does not postulate general principles that can be used to make predictions for experiments that have not been performed. However, 13C MFA does not rely on maximum growth assumptions, the general applicability of which has been questioned [14,34,35] and shown to be inaccurate for engineered strains that are not under long-term evolutionary pressure [17]. Moreover, the comparison of measured and fit labeling patterns provides a degree of validation and falsifiability that FBA does not possess: an inadequate fit to the experimental data indicates that the underlying model assumptions are wrong. In contrast, FBA produces a solution for almost any input. Several attempts to combine the complementary virtues of 13C MFA and FBA have been reported. For example, they have often been combined to test new FBA-based methods, since 13C MFA is considered to be the most authoritative determination of fluxes. Both Segrè et al [16] and Yizhak et al [19] used 13C MFA to validate MOMA and IOMA, respectively. More recently, Schuetz et al [14], used 13C MFA-derived fluxes to compare predictions from an MFA model containing ∼ 100 reactions and using different objectives, and also to demonstrate the applicability of pareto optimality to predict fluxes [35]. Choi et al [36] combined both methods by using flux ratios obtained from 13C MFA to constrain FBA for genome-scale models through the use of artificial metabolites. Although genome-scale models were not used, Suthers et al [27] presented 13C MFA for a large scale model containing 350 reactions. The OPENflux open source software by Quek et al [37], allows for certain reactions in 13C MFA to be used only for stoichiometric modeling purposes. More recently, Chen et al [38] for the first time modeled the same E. coli strain under the same conditions using FBA and 13C MFA, using some of the information of the latter to constrain the former. In a similar vein, Kuepfer et al [39] constrained a S. cerevisiae genome-scale model with fluxes obtained from 13C MFA and determined the flux distribution by using the minimization of the overall intracellular flux as the objective function. However, to date there has been no attempt to use the data from 13C labeling experiments to constrain fluxes for a genome-scale model without assuming that metabolism is evolutionarily tuned to optimize an objective function (such as the growth rate optimization typically used in FBA). The traditional underlying assumption in the field of 13C MFA is that the degrees of freedom of the model must carefully match the amount of information obtained from the labeling patterns. Indeed, if the mathematical formulation of 13C MFA were that of a linear programming problem (as for FBA) it would be pointless to try to constrain over a hundred degrees of freedom with the approximately 50 measurements of labeling used in this study, for example. However, 13C MFA is a nonlinear fitting problem (with fluxes being the parameters) and these problems behave very differently for the underdetermined case, a special case of what is known as “sloppy” models in statistical mechanics [40–43]. These underdetermined nonlinear fits exhibit some degrees of freedom which are highly constrained (the parameters cannot be changed without noticeable measurable effects) and other degrees of freedom barely constrained at all (parameters can be changed freely, see section 3. 3 in Brown et al [40] for a concrete example). Even if all degrees of freedom are not fully determined, the model can still be used to test hypotheses effectively [44]. These characteristics have been shown to be of general nature and apply to a variety of nonlinear problems appearing in systems biology [44], insect flight and variational quantum wave functions [45], interatomic potentials [46] and a model of the next-generation international linear collider [42]. In this paper, we present a systematic and rigorous framework to take full advantage of the nonlinear nature of the flux fitting problem and find all fluxes compatible with the 13C labeling data for a genome-scale model. In order to constrain fluxes effectively, we make the biologically relevant assumption that metabolic flux flows from core metabolism (defined below) to peripheral metabolism and does not flow back. The simultaneous use of 13C labeling data and genome-scale models produces a new computational approach combining the advantages of both FBA and 13C MFA: two-scale 13C Metabolic Flux Analysis (2S-13C MFA, see Fig 1). 2S-13C MFA determines fluxes for a full genome-scale model, taking into account the system-wide balances of metabolites. However, 2S-13C MFA does not rely on maximum growth assumptions: instead it uses the data obtained from 13C labeling experiments to constrain feasible fluxes. The use of this data is shown to constrain glycolytic and pentose phosphate pathway fluxes 8–50 fold more effectively than using only measured extracellular fluxes. We use this new method to compare different predictive methods and show where and why they fail. Based on that information, we develop a new predictive method that is able to produce a full quantitative prediction of 48 labeling measurements, going beyond the usual qualitative (e. g. grow/no grow) predictions.
In 13C MFA, most reactions in the cell’s enormous metabolic network have a very limited contribution to the 13C labeling of the observed metabolites (regardless of whether these are free metabolites [47] or proteogenic amino acids [4,9]). Thus, since the complexity of the problem of flux measurement from labeling information scales nonlinearly with the number of reactions, the metabolic network in 13C MFA consists of a minimal set of reactions (core set) that most influence labeling patterns (typically central carbon metabolism but may include other extra reactions, such as those describing protein turnover [48]), stemming from the literature or the researcher’s experience. This approach is able to convincingly explain labeling patterns for amino acids and intracellular metabolites for model organisms (e. g. E. coli [49,50] and S. cerevisiae [51,52]) under well-studied conditions (e. g. glucose feed). The good fits to experimental data support, to a good approximation, the underlying assumptions that carbon precursors flow from a central core metabolism to peripheral metabolism and do not flow back. However, it would be desirable to generate the minimal core network systematically through a computational method, so the technique can be generally applied to non-standard cases: e. g. non-model organisms, bioengineered strains, alternative non-standard carbon sources, human cells or microbial communities. Most crucially, such a method would explicitly test the assumption that reactions not included in the minimal network do not significantly influence labeling. As an added benefit, it would take advantage of the significant community efforts that have been put into developing and improving genome-scale metabolic models [53]. Another important benefit of integrating 13C MFA with genome-scale metabolic models is to check for consistency of the inferred core fluxes with peripheral metabolism. The central metabolic network needs to produce not only carbon precursors for peripheral metabolism, but also ATP and reducing equivalents. Genome scale networks offer the possibility to track every reaction consuming and producing energy or reducing equivalents, and therefore can be used effectively to study the interplay between peripheral and core metabolism. On one hand, since core metabolism involves the reactions with largest fluxes in metabolism, once these are set by the 13C labeling data, peripheral metabolism is expected to be highly constrained. On the other hand, changes in the usage of core intermediates in peripheral metabolism will also have significant effects on core fluxes since they need to provide these carbon intermediates, on top of ATP and reducing equivalents. While this has been already predicted in terms of small changes in biomass composition having significant effects on central carbon fluxes [54], the effect is an order of magnitude more important in bioengineered strains [55]. In such cells, peripheral metabolism can be vastly altered, and it is crucial to have a method that can take the diverse changes into account in a systematic manner. 2S-13C MFA addresses both of these issues and provides a complete rigorous framework for estimating fluxes using 13C labeling data in the context of a genome scale metabolic network. The algorithm is comprised of a set of optimization problems to be applied sequentially, as shown in the outline in Fig 2. Our approach divides the metabolites and reactions in a genome-scale model into two groups to be modeled at different scales of resolution (Fig 3A), in the spirit of the multi-scale approach in engineering and physics [56]. For “core” metabolites and reactions both stoichiometry and carbon labeling are tracked. For the remaining “non-core” metabolites and reactions, only stoichiometry is tracked and their contribution to the core set labeling is ignored. Crucially, the algorithm recursively adjusts the division into core and non-core reactions according to the error in the labeling fitting. This procedure is illustrated with a small illustrative network (see Fig 3B and 3C) throughout the following subsections. Notice that, because of this core readjustment and the inclusion of other reactions needed to fully explain labeling [48], the core set needs not be the same as is usually referred to as central carbon metabolism. 2S-13C MFA is supported by two independent data sets: 13C labeling experimental data and genome-scale stoichiometry. We will now show that the method is robust to errors in the 13C labeling experimental data and more robust than FBA to stoichiometric errors. Fluxes calculated through 2S-13C MFA are robust with respect to the experimental accuracy of the 13C labeling data. We show this by generating new sets of 13C data where the new labeling is randomly chosen within the experimental error. The calculated profiles do not change for this new labeling significantly from the initial profiles, as can be observed in Fig 7. Cofactor ambiguities in the reconstruction of genome-scale models introduce errors in determining fluxes which are much smaller than for FBA. Reconstruction errors are still possible in genome-scale models, in spite of the care used to develop them and the continous improvement in the reconstruction for each new release [13]. We studied the resulting flux profiles after a reconstruction error was simulated (by changing NADPH to NADH dependence for a high flux reaction, G6PDH2r) and found that the change in fluxes was much less severe than for FBA (Fig 8). We found similar results for other reactions (S16 Fig). We attribute this difference to FBA relying heavily on stoichiometric information whereas 2S-13C MFA can additionally count on 13C labeling data to constrain fluxes. Furthermore, changes in the biomass composition reaction do not affect calculated fluxes significantly (S17 Fig). 2S-13C MFA can produce cofactor balance information which cannot be produced using either FBA (does not provide confidence intervals bounded by 13C data) or 13C MFA (does not consider all possible cofactor-producing reactions in the cell metabolism). Discerning cofactor production and consumption can be crucial to understand cell behavior. Not only does cofactor balancing have an essential bearing on tumor cell death [65], but it can also have a significant impact on final production for bioengineered cells [33,55,66,67]. Core metabolism and biosynthetic pathways are linked through cofactor balancing, so an expansion of 13C MFA to include genome-scale stoichiometry such as that provided by 2S-13C MFA is of great interest for bioengineering. While some previous 13C MFA efforts [27,28,37] did take cofactor balances into account, this has not been done for all metabolites in a genome-scale network. This capability is illustrated (Fig 6) with the main electron carriers in metabolism, NADPH and NADH, responsible for redox balances and, ultimately, for respiration. A majority of NADPH production (47–55%, obtained by normalizing data from Fig 6) in wild type E. coli is produced by the pentose phosphate pathway (reactions G6PDH2r and GND: 25–27%) and the TCA cycle (reaction ICDHyr: 22–29%). Notice that, unlike FBA, 2S-13C MFA produces confidence intervals bounded by 13C data: for example, if the flux through reaction GND were to be more than 3. 2 mMol/gdw/hr (30% of total NADPH creation according to upper left part of Fig 6), the computationally derived labeling for (e. g.) m = 0,4 for malate (mal-L), for m = 1 for phosphoenolpyruvate (pep) and for m = 1,2 for ribose 5-phosphate (r5p) would differ from the measured value by more than the CEMS measuring error (see Eq 23 and S1–S3 Figs). The transhydrogenase reaction THD2, transferring electrons from NADH to NADP, also contributes significantly (∼ 44%) but, unfortunately, the flux through this reaction is only very loosely constrained by the current data set. The main NADPH sink involves glutamate production (reaction GLUDy: 29–65%), with fatty acid biosynthesis (C181SN) in a distant second place (7%) and a collection of other reactions at equally low levels (< 12%). NADPH consumption and production are radically altered when the gene encoding the initial enzyme in glycolysis (PGI) is knocked out (Fig 6, right upper panel). This genetic manipulation forces the rerouting of all glycolytic flux into the PPP and a lower growth rate (0. 17–0. 23 vs 0. 83–0. 9 h-1). 2S-13C MFA allows us to study the impact of this flux rerouting on the rest of the cell metabolism. For example, the overall NADPH demand falls three fold (vs a ∼ four fold decrease in growth rate). Furthermore, the demand is not equally met by the pentose phosphate pathway and the TCA cycle anymore: in the Δpgi, the PPP enzymes G6PDH2r and GND produce most of the NADPH (81%) with a much less important role played by the TCA cycle (13–19%). Interestingly, the absolute flux through the PPP for the pgi KO and the wild type is very similar. NADH is generated mostly by the glyceraldehyde-3-phosphate (36–37%, GAPD), malate (10–27%, MDH), pyruvate (15–25%, PDH) and 2-Oxoglutarate dehydrogenases (7–11%, AKGDH), which provide 79–98% of all NADH. Consumption is dominated by the ubiquinone-8 NADH dehydrogenase (NADH6,43–83%). After pgi is knocked out, NADH demand decreased ∼ 4 fold, with relative production by MDH undergoing the biggest drop (from 10–27% to 5–7%). These results are consistent with previous 13C MFA calculations for NADPH production [68], which validates this approach. Furthermore, they also provide information on secondary metabolism and NADH balances in a systematic manner. A similar analysis can be done for acetyl-CoA, a key metabolite and common bottleneck in metabolic engineering [69,70] or any other relevant metabolite such as ATP, AMP, Coenzyme A or FADH, which would provide detailed information on the organism’s underlying physiology. The 2S-13C MFA results pinpoint which (non measured) metabolites are expected to be detected in the extracellular medium based on the values of exchange fluxes as constrained by the 13C labeling experiments (S13 Fig). Knowledge and quantification of the full range of excreted metabolites (also known as exometabolome or metabolic footprint [71]) is desirable for a full understanding of the biochemical impact of the cell on its environment. For metabolic engineering purposes, this knowledge provides important clues as to how to close the carbon balance and whether toxic compounds are being produced in the fermentation. For human metabolism, better prediction of extracellular fluxes can yield improved metabolic predictions when integrated with physiologically-based pharmokinetic models [72]. The exchange fluxes predictions typically have large confidence intervals, but these intervals are much smaller than for FBA (see S13 Fig). Metabolites expected to be detected in the medium are those whose exchange fluxes have net positive maximum and minimum values for a period of time sufficiently long so as to reach detection limits. For the E. coli strains considered here, urea, glycolate (glyc, S14 Fig), fumarate (fum, S15 Fig) and acetaldehyde (acald) are the non-typical metabolites expected to be present in the medium (acetate is already measured). This prediction of atypical metabolites is of particular interest in light of the recent discovery of extended overflow metabolism [73]. A full prioritized list can be obtained from the exchange flux information and used to direct mass spectrometry, NMR or vibrational spectroscopy efforts to find the missing metabolites until carbon balance is met. For future improvement of intracellular metabolic flux predictions, constraints introduced by extracellular metabolite measurements are very effective and usually easy to measure, but it is necessary to know which metabolites to look for and 2S-13C MFA provides precisely that type of insights. 2S-13C MFA produces nearly the same results as 13C MFA for central carbon metabolism (see e. g. S4 Fig and S18 Fig). This similarity is not surprising since 2S-13C MFA is designed to mimic 13C MFA for this part of metabolism (see “Limiting flux to core reactions” section). The only difference for the current data set can be found in the TCA cycle flux, as described below. These differences arise because genome-scale models account for fluxes to biomass in a more detailed and realistic manner and because they do not rule out unexpected metabolic routes compatible with the available data. Flux through the TCA cycle is lower in the 13C MFA solution (S4 Fig vs S18 Fig) because of an inaccurate account of fluxes to biomass in the 13C MFA model. Specifically, the large fluxes draining acetyl-CoA into biomass and to the exterior of the cell as acetate (each a third of the pyruvate dehydrogenase flux, PDH) imposed in the original publication [47] have been overestimated in this particular case. While the typical biomass function used in 13C MFA involves a specific stoichiometry of central carbon intermediates converted directly to biomass, this is only an approximation since some of the metabolites required for biomass growth are not represented in the minimal network model and need to be substituted by their requirements in terms of intermediates present in the minimal network (acetyl-CoA, in this case). These effects require significant effort to be accurately incorporated into a small-scale model, but they are elegantly handled by the genome-scale model. In this case a large flux of acetyl-CoA to biomass is assumed in the 13C MFA solution, while the 2S-13C MFA solution determines the requirements on core metabolism by setting the biomass flux to the measured growth rate, and automatically determining the fluxes needed to provide the metabolites present in the biomass equation through the comprehensive stoichiometry network reflected in Eq 2. Hence, genome-scale models provide an automatic and detailed accounting of biomass requirements. The difference in TCA cycle flux extends to the glyoxylate shunt in some cases in which the isocitrate lyase (ICL) is active for the 2S-13C MFA (S5 and S6 Figs) solution, but not in the 13C MFA solution due to a limited 13C MFA model. In the 2S-13C MFA solution, the objective function (Eq 1) becomes slightly lower by diverting flux into the ICL and then shuttling it out of the core set from glyoxylate (glx) into glycolate. This route is not included in the 13C MFA model and activating the ICL for this model would produce a large amount of flux through the malate synthase (MLS). This MLS flux would significantly deteriorate the fit to malate labeling, so the glyoxylate shunt remains inactive for 13C MFA case. The shuttling of glyoxylate into glycolate is unexpected and, in fact, could be the result of a numerical artifact since the only labeling data available in the full TCA cycle is that of malate (mal-L) and small errors in the labeling measurements (or their confidence intervals) for this metabolite can lead to erroneous solutions. The glyoxylate shunt is known to be inactive under the given conditions and one could use this information and constrain its flux to zero. However, we aim to produce a general method; of use under conditions where this information may not be available (e. g. exotic feeds or bioengineered cells). Hence, we decided not to constrain the glyoxylate flux in order to show how the genome-scale model constrained by measured data can produce testable and falsifiable consequences to detect this type of errors. In this case, the shuttling of glyoxylate through glycolate results in glycolate being exported out to the medium (S14 and S15 Figs) and detectable through (e. g.) MS methods. If no glycolate is found, that extracellular flux can be set to zero and the glyoxylate shunt flux will decrease to zero, so as not to deteriorate the fit of the malate labeling pattern. Alternatively, the availability of labeling patterns for additional metabolites (fumarate, fum and glyoxylate, glx) would confirm or deny the ICL activity. In this way, 2S-13C MFA can fruitfully use available data to suggest and test unexpected metabolic activity; such as the surprising heme degradation in cancer cells with a non-functional TCA cycle [26]. Our goal in developing 2S-13C MFA is to make a clear distinction between highly reliable constraints such as those induced by 13C labeling data, measured extracellular fluxes, carbon transition information and reaction stoichiometry, and reasonable hypotheses (such as growth rate optimization) that may not be universally applicable. 2S-13C MFA has been developed to provide a self-consistent and unbiased determination of the range of metabolic fluxes compatible with available experimental data. Unlike other COBRA methods, the fluxes obtained by 2S-13C MFA are backed up by the extra validation provided by the correct fit of 48 relative labeling measurements (see S1–S3 Figs). Hence, we use it here as a reference to compare various COBRA predictive methods based on different hypotheses. This procedure permits us to determine which method predicts fluxes most accurately. We compared flux profile predictions for pyk and pgi gene knock outs calculated at three different time points (5,6 and 7 hrs for wild type and pyk KO, and 16,21 and 23 hrs for pgi KO) using six different methods: FBA maximizing either growth and ATP production [74], Minimization of Metabolic Adjustment (MOMA [16]), Regulatory On/Off Minimization (ROOM [75]) and two new methods developed in this paper (Fig 9 and S19 and S20 Figs). These new methods are 13C MOMA and 13C ROOM, which leverage 2S-13C MFA flux profiles obtained for the wild type strain to improve flux predictions (see S3 Text), similar in spirit to the approach by Kuepfer et al [39]. Growth maximizing FBA provides predictions that are not quantitatively accurate, but seem approximately right for most cases (S19 and S20 Figs). This fact probably explains its success in metabolic engineering applications. Since FBA is typically used to make partial flux predictions by using some measured flux information for the predicted experiment, we tested this variation as well. Once fluxes for growth rate and excreted metabolites were constrained to the measured values, pgi predictions became extremely accurate (S20 Fig). However, when a full prediction (i. e. , not using any data from the predicted experiment) was sought, FBA noticeably failed for this KO due to its inability to predict the drop in growth rate. ATP maximizing FBA fails most noticeably for the pyk KO, probably because the PYK reaction is involved in ATP production and its elimination significantly changes the ATP balance when ATP is to be maximized. MOMA, ROOM, 13C MOMA and 13C ROOM flux predictions fail most blatantly for the pgi knockout strain because growth rates change radically and the method tries to maintain the previous flux levels (S19 Fig). Hence, we can expect these methods to offer good results only when changes expected in glucose intake and growth rates are relatively small. An improvement of these methods could be obtained if only the relative flux profiles are used for the prediction in the algorithm and the growth rate is separately obtained. Since 13C MOMA and 13C ROOM use 2S-13C MFA flux profiles from the wild type strain, we expected that they would more accurately predict fluxes than would MOMA and ROOM. We find that is the case for fluxes in glycolysis and PPP, but the TCA cycle flux values are less accurate (see S19 Fig). This phenomena is most likely due to our initial flux profile not being very accurate for the TCA cycle (large confidence intervals, see S4–S6 Figs) since many fewer labeling measurements are available for TCA metabolites (only malate vs eight metabolites available for glycolysis and PPP). Labeling data for more metabolites next to branch points (e. g. fumarate, glyoxylate, isocitrate, alphaketoglutarate) would help reduce the flux confidence intervals for this area of metabolism. These comparisons illustrate how 2S-13C MFA can be used to test COBRA methods, see where they fail and why, and use this information to improve them. The combination of 13C labeling data and genome-scale models with COBRA prediction methods allows us to make some predictions we cannot do using either 13C MFA or FBA. The ability to make accurate predictions of metabolism is of fundamental importance to make metabolic engineering a more predictive discipline. The 13C MOMA predictions of glycolysis and PPP fluxes for the pyk KO at the 5-hr time point surpass those of all other methods (Fig 9 and S19 Fig). In fact, the flux predictions are precise enough that we can even predict with reasonable accuracy the metabolite labeling to be expected from that strain at that time point (Fig 10). This is a prediction of directly measured data instead of a derived measurement such as flux. No information from the pyk strain at 5 hrs was used. We used the 2S-13C MFA flux calculations from the wild type strain at 5 hrs to obtain fluxes compatible with those measurements and then, through 13C MOMA, obtain the predicted fluxes if pyk were to be knocked out. The corresponding labeling patterns were then derived. Notice that this is very different from the partial predictions of labeling that 13C MFA routinely produces, where the labeling for the present experiment is predicted. Hence, this constitutes a full quantitative prediction of metabolite labeling from a 13C labeling experiment, which has not been reported before, to our knowledge. The prediction of metabolite labeling is particularly relevant as validation because they refer to a directly measured quantity (i. e. the MDV obtained from the CE-TOFMS). Fluxes are derived quantities relying on a variety of implicit assumptions (i. e. the two-scale approximation, metabolic pseudo-steady state, no accumulation of intermediate metabolites, genome-scale model completeness and accuracy, cell homogeneity…). The real test that these assumptions are not severely violated and that the method provides reliable flux profiles is to use them to predict directly measured quantities for other experiments (in this case, labeling patterns). Moreover, this example shows that coupling 13C labeling data with COBRA methods opens the possibility to go beyond qualitative predictions (e. g. grow/no grow). In this manuscript, we have shown how to maximize the information obtained from 13C data to constrain genome-scale models, and that once core metabolism is set by 13C labeling data information, the rest of metabolism is generally highly constricted. As is a usual behavior in “sloppy” nonlinear fitting problems [40], some fluxes are very effectively constrained and some others only loosely. Confidence intervals obtained through 2S-13C MFA immediately identify these two types of fluxes and show that the use of 13C labeling experimental data produce much narrower confidence intervals than those produced by FBA. The method is generally applicable to any genome-scale model or feed, and can be used to expand the use of 13C-based flux analysis beyond the customary cases to tackle non-standard feeds, exotic organisms and systems described by large stochiometric matrices such as the human metabolic network [76], those derived from adding macromolecular synthesis [77] or microbial communities [24]. 2S-13C MFA produces similar results as 13C MFA for the region where the latter is valid: central carbon metabolism. 2S-13C MFA, however, extrapolates the constraints induced by the 13C labeling data to a genome-scale model, providing fluxes not only central carbon metabolism but also for peripheral metabolism. This was illustrated firstly by a detailed description of NADPH and NADH production and consumption and, secondly, by predicting unmeasured metabolites expected in the extracellular medium. 2S-13C MFA does not use an evolutionary optimization principle (such as growth rate optimization) but, rather characterizes all flux profiles compatible with the experimental data. The extra validation gained by matching the measured labeling values is used to test the validity of the maximization hypothesis by comparing these results with predictions obtained through FBA and other COBRA methods. The comparison of flux profiles predicted with COBRA methods and those obtained through 2S-13C MFA provided not only a ranking of accuracy for predictions but also insight as to how to improve predictive methods. An improved version of MOMA using 2S-13C MFA profiles as a starting point is able to predict the outcome of 48 direct measurements of metabolite labeling for a pyk KO experiment using only data from a different experiment involving the wild type. This capability shows that using 13C labeling experimental data enables accurate predictions beyond qualitative cases (e. g. grow/no grow). This method represents another step in the effort to make bioengineering a more predictable endeavour. The new method hinges on a simple assumption: flux flows from the core set to peripheral metabolism and does not flow back. This assumption is supported by the good fits obtained in general by 13C MFA methods thus far. There might be situations where this two-scale assumption is not applicable and these will be pinpointed by unacceptable fits to the labeling data. The core set of reactions, however, is flexible and can be enlarged as needed to provide acceptable fits to the labeling data. Hence, phenomena such as protein turnover [48,78] or cell scavenging in stationary phase can be included by adding the appropiate reactions to the core. The availability of carbon transitions for genome-scale models [79] facilitates a systematic core enlargement. In spite of its simplicity, FBA-based modelling has already exhibited significant success: only three elements are used in this modelling scheme (genome-scale stochiometry, measured extracellular fluxes and an optimization principle) but they have been successfully used to rationally engineer strains used for large-scale industrial production [21–23]. However, certainly not every single flux in a flux profile obtained through FBA can be trusted. 2S-13C MFA unites the informative constraints of 13C labeling experiments with genome-scale stoichiometry to improve the determination of internal metabolic fluxes and set confidence intervals based on experimental data. 2S-13C MFA completes and improves 13C MFA by enforcing a global balance of metabolites instead of balancing only a few chosen metabolites. We believe it will be a tool of extreme utility in bioengineering, at a time when a variety of different frameworks for flux prediction for genome-scale models are becoming available [80,81]. Furthermore, we think that its widespread use to determine metabolic fluxes will affect our understanding of fundamental biological problems [82] beyond bioengineering.
The two-scale approximation assumes that non-core reactions do not contribute directly to the labeling of core metabolites, since carbon precursors flow from core metabolism into peripheral metabolism and do not flow back. This is represented in terms of a genome-scale model by limiting to zero the flux of reactions flowing into core metabolism (see Algorithm 1). The first step in 2S-13C MFA (Fig 2) hence consists in taking each reaction that has a product in core metabolism and setting the upper bound to zero. However, it may be the case that this extra constraint makes it impossible to meet the measured growth rate (we check this by solving the corresponding FBA problem). In that case, setting the upper bound to a fraction of the glucose uptake rate is tested (first 0. 05 and then 0. 2 for this case). Since the labeling of core metabolism can be impacted by reversible reactions with reactants included in the core set as well, we cover this case by limiting the lower bound of the reaction to zero or the lowest value that permits growth. The impact of the reactions that could not be set to zero will be checked later through External Labeling Variability Analysis (ELVA, Figs 2 and 4). The input for the first step (Fig 2) is the genome-scale model with the carbon transitions for the core set integrated in it. The output consists of the genome-scale model with lower and upper bounds modified by this “Limiting flux to core” procedure. A detailed description of this first step in the diagram shown in Fig 2 can be found in the pseudo code in Algorithm 1. Algorithm 1. “Limiting flux to core” pseudo code for each reaction j flowing into core: limits = [0,0. 05,0. 2]*glucose_uptake limit = limits. next () goOn = True while goOn: if reaction j has forward flux: ub[j] = min (ub[j], limit) else if reaction j has backward flux: lb[j] = sign (lb[j]) *min (abs (lb[j]), limit) solve FBA problem goOn = (FBA problem has no solution) and (limit is not the last value in limits) limit = limits. next () Where limits. next () obtains next value in list limits (the first one if uninitiated), and glucose_uptake is the value of the glucose uptake rate. Solve FBA problem refers to finding the solution to the problem given by equations 1–3 in S1 Text. has forward flux refers to the reaction having a possible positive flux (i. e. positive upper bound, ub) flowing into the core set and has backward flux refers to the reaction having a possible negative flux (i. e. negative lower bound, lb) flowing into the core set for reversible reactions. sign (lb[j]) is the sign of the lower bound lb[j]. ub[j] and lb[j] denote upper and lower bounds for reaction j, respectively. The second step in 2S-13C MFA (Fig 2) involves fitting the measured metabolite labeling by solving the optimization problem in Eqs 1–7, where the upper (ubj) and lower bounds (lbj) have been limited by the previous step (“Limiting flux to core”). 2S-13C MFA is a hybrid of FBA and 13C MFA (see Fig 1 and S1 Text) where the stoichiometry constraint is applied to the full genome-scale network, as in the case of FBA, and the labeling constraints [58] are applied only to the core set of reactions and metabolites, as is the case for 13C MFA. Unlike previous efforts [36,39], these constraints are enforced simultaneously, instead of sequentially using the results of 13C MFA to constrain the FBA problem. This simultaneous approach is more rigorous than using slack coefficients (δ in Kuepfer et al [39]), does not need to invoke an optimization principle to calculate fluxes and allows for the global metabolite balance to affect core metabolism fluxes. Furthermore, it is also self-consistent whereas in the sequential approach one might find that fluxes that flow into core metabolism are active, even though they were not taken into account to do the initial 13C MFA fit. In the notation of Suthers et al [58] (GAMS files available in S1 Code): Minimize OF= (∑e∈Emeasm∈Me (femexp−femΔem) 2/| Me |) /| Emeas | (1) Subject to: ∑ j S i j v j = 0 ∀ i ∈ I N, j ∈ J (2) l b j ≤ v j ≤ u b j ∀ j ∈ J (3) ∑ m ∈ M e f e m = 1 ∀ e ∈ E c o (4) Σe′∈Eco ( (Σl|EMMe′−>el>0EMMe′−>elVl) fe′m) + (Σl|Sil*<0Sil*Vl) fem=0∀m∈Me, e∈Ei, i∈IcoN (5) f e m = ∑ w ∈ W e m ∏ n = 1 | E e | f e n m n ∀ m ∈ M e, e ∈ E c o c (6) v j = ∑ l ∈ J c o B m a p j l V l ∀ j ∈ J (7) where: S e t s _ I ≡ { i }: Set of all metabolites. I c o ⊂ I: Set of core metabolites. I c o N ⊂ I c o: Set of non-exchange core metabolites. J: Set of fluxes. J c o ⊂ J: Set of core fluxes. J B: Set of fluxes with backward and forward fluxes differentiated, e. g. PGI_f, PGI_b, PGL…etc. J c o B ⊂ J B: Set of core fluxes for J B. E = { e }: Elementary Metabolite Units (EMUs). E c ⊂ E: Combined EMUs. E i ⊂ E: EMUs from metabolite i ∈ I. E c o c ⊂ E c: Core combined EMUs. E e ⊂ E: EMUs that produce combined EMU e. E c o ⊂ E: EMUs corresponding to core metabolites. E m e a s ⊂ E: EMUs corresponding to measured EMUs. W e m: Set of every possible mass isotopomer multiplet of E e that produce the mass isotopomer m of e. M e: m values for MDV of emu e: 0,1, ⋯, # of carbons ine. P a r a m e t e r s _ E M M e ′ - > e l = 1 k ife ′ produces e through reactionl ∈ J c o B, 0 otherwise. See Suthersetal[58]. S i j: Stoichiometry matrix. S i l *: Stoichiometry matrix with backward and forward fluxes differentiated. u b j, l b j: Upper and lower bounds for reaction j. f e m e x p ∈ [ 0,1 ]: Experimentally measured MDV for emu e from metabolite m. Δ e m: Measurement error for f e m e x p. m a p j l = 1 * g l u c u p t ifl corresponds to forward flux of j. = - 1 * g l u c u p t ifl corresponds to backward flux of j. V a r i a b l e s _ v j: Flux value of reaction j ∈ J, in mmol/gdw/h. V l: Flux value of reaction l ∈ J c o B, normalized to glucose input rate. f e m ∈ [ 0,1 ]: Mass isotopomer fraction (MDV) for emu e from metabolite m. Notice that S i j* is not the same as Sij, since J and JB are slightly different sets of fluxes. In fact: S i l * = S i j if l is the forward version of j. S i l * = - S i j if l is the backward version of j. (8) Notice that Eqs 2 and 3 are the traditional FBA contraints and that Eqs 4–6 are the 13C labeling constraints, but they have been limited to core metabolites and reactions. The mapping (Eq 7) converts the fluxes from the 13C MFA description (higher resolution) where fluxes are normalized to the glucose uptake rate (glucupt) into the FBA description (lower resolution). We describe the 13C MFA description as of higher resolution because the 13C labeling data can pick up differences in forward and backward fluxes (JB set), whereas the purely stochiometric approach of FBA can only constrain net fluxes (J set). For all 2S-13C MFA calculations all input flux was supposed to be routed through GLCpts (glucose transport via PEP: pyr pts [57]). Optimization problems were run N = 30 times and the one with lowest objective function picked. The value of N was chosen so as to avoid relative minima of OF. A plot of how the OF saturated as N increased can be found in S21 Fig in the supplementary material. In addition to the OF, the more typical sum-of-squares residuals (SSR) for the fits have been included in the legends of S1–S3 Figs. However, standard good-of-fitness metrics, such as those proposed by Antoniewicz et al [85], are not applicable to 2S-13C MFA. By using genome-scale models, the number of estimated free fluxes (p) is higher than the number of independent measurements (n = 48 − 9 = 39 in this case) and the χ2 (n − p) distribution of a negative number of degrees of freedom is then not properly defined. This apparent paradox arises because of the implicit assumptions in the χ2 statistics approach to these types of fits. The null model assumes that each of the terms in the SSR are independent random variables, hence the degrees of freedom are the number of terms in the SSR [86]. Nonetheless, we know that, for 13C MFA, the terms in the SSR are typically not independent: the labeling pattern (MDV) of related compounds (i. e. amino acids arising from the same precursors) are very similar. Hence, Antoniewicz et al [85], decided to use as degrees of freedom the number of individual measurements minus the number of estimated free fluxes. This proposal seems reasonable for standard 13C MFA, but breaks down for genome-scale models with a much larger number of degrees of freedom. One could choose a different null model for the χ2 statistics (and previous approaches have shown that the standard χ2 goodness-of-fit approach is probably too conservative, see Fig 3A in [27]) but the crux of the matter is that we believe that using p-values < 0. 05 as an absolute arbitrary threshold for significant vs insignificant results is too simplistic, as do other biological researchers [87] or R. A. Fisher himself [88]. Hence, what we report here (in S1–S3 Figs) is what we believe is a better (and more intuitive) way to estimate how good a fit is: the average objective function normalized to the measurement error (Eq 1). This is the answer to the question: how different are my fits from the experimentally measured values, measured in units of the experimental error? This is in accordance in spirit with the suggestion of presenting measures of significance without arbitrary thresholds [87]. Notice that none of the objective functions (OF) in S1–S3 Figs is smaller than one, indicating that the difference between experiment and theory cannot be explained through experimental error. This may be because the experimental error for the labeling pattern was underestimated or because the model fails to explain the full labeling pattern. These results are in line with the general trend that fits to intracellular metabolites tend to be worse [89] than fits to proteogenic amino acids [49]. The discarding of the 0. 05 p-value criterium does not imply that the considerable effort employed in obtaining excellent fits to data [49,51] goes unrewarded. Under the 2S-13C MFA method, a bad fit to the experimental data results in a larger value for δe m in Eq 23 and wider confidence intervals. Hence, worse fits beget less flux resolution. Finally, we think that the best validation of the of a flux fit is using the flux distribution to predict the results of another experiment, as we did in Figs 9 and 10. Once a core set is chosen, ELVA establishes the maximum impact of non-core metabolism in the labeling of the measured metabolites. In order to do so, only core metabolism is considered and the impact of non-core metabolism is represented through “inflow” metabolites and reactions (see Fig 3C for an example). Inflow reactions agglomerate all non-core reactions flowing into (or out of) a particular core metabolite (see S22 Fig in supplementary material) and are assigned trivial carbon transitions (e. g. abc --> abc for a three carbon metabolite). Inflow metabolites are dummy metabolites with the same number of carbons as the involved core metabolite. ELVA constrains fluxes to the solution obtained in the “Fit data” step (Eq 10) and maximizes and minimizes the computational MDV for each m in each of the labeled metabolites. Since fluxes are fixed for all reactions and labeling is fixed for all metabolites except inflow metabolites, the optimization problem in Eqs 9–15 quantifies the maximum and minimum effect that this unknown labeling (since it comes from non-core metabolism) could have on the labeling pattern for the measured metabolites (fem∀m, e ∈ Emeas): Min/max f e m ∀ m ∈ M e, e ∈ E m e a s (9) Subject to: V j = V ¯ j j ∈ J c o e x t B (10) ∑ j S i j * V j = 0 ∀ i ∈ I c o e x t N, j ∈ J c o e x t B (11) l b j ≤ V j ≤ u b j ∀ j ∈ J c o e x t B (12) ∑ m ∈ M e f e m = 1 ∀ e ∈ E (13) ∑e′∈E ( (∑j|EMMe′−>ej>0EMMe′−>ejVj) fe′m) + (∑j|Sij*<0Sij*Vj) fem=0∀m∈Me, e∈Ei, i∈IN (14) f e m = ∑ w ∈ W e m ∏ n = 1 | E e | f e n m n ∀ m ∈ M e, e ∈ E c o e x t c (15) where symbols are as explained before, with the addition of: S e t s _ I c o e x t B: Set of extended metabolites. J c o e x t B: Set of extended fluxes with backward and forward fluxes differentiated. P a r a m e t e r s _ V ¯ l: solutions to the problem given by equations 1 - 7 J c o e x tB and Icoext are the set of reactions and metabolites obtained after expanding the core set to meet stoichiometry requirements (Eq 11) as is explained in supplementary S22 Fig. The input for this step is the flux profile obtained from the data fit and the output is an ELVA plot (Fig 4) used to decide whether the solution is self-consistent or not. Once a self-consistent core set has been determined through the recursive procedure in Fig 2, flux ranges compatible with the experimental data are obtained through the following optimization problem: Min/max V j ∀ j ∈ J (16) Subject to: ∑ j S i j v j = 0 ∀ i ∈ I N, j ∈ J (17) l b j ≤ v j ≤ u b j ∀ j ∈ J (18) ∑ m ∈ M e f e m = 1 ∀ e ∈ E c o (19) Σe′∈Eco ( (Σl|EMMe′−>el>0EMMe′−>elVl) fe′m) + (Σl|Sil*<0Sil*Vl) fem=0∀m∈Me, e∈Ei, i∈IcoN (20) f e m = ∑ w ∈ W e m ∏ n = 1 | E e | f e n m n ∀ m ∈ M e, e ∈ E c (21) v j = ∑ l ∈ J c o B m a p j l V l ∀ j ∈ J (22) (f e m - f e m e x p) 2 ≤ δ e m 2 ∀ e ∈ E m e a s, m ∈ M e (23) where symbols are as explained before, with the addition of: P a r a m e t e r s _ δ e m: Maximum error allowed for f e m. δ e m = Δ e mifΔ e m > ϵ e m and δ e m = 1. 1 * ϵ e m if Δ e m < = ϵ e m where ϵ e m = f e m f i t - f e m e x p from the solution to equations1-7. The fluxes vj and Vj are initialized to the values obtained from solving Eqs 1–7. The results of the minimization and maximization give the flux upper and lower bound compatible with the experimental data from the 13C labeling experiments. This procedure is similar to the FluxRange procedure in Suthers et al [27], with the exception that we use set constraints for each measured data point (Eq 23) instead of only for the objective function. This approach guarantees that fluxes produce labeling patterns within the experimental error in a much more efficient way than a monte carlo approach. For example, consider these reactions which conform a futile cycle: NDPK 1: g d p c + a t p c < = > g t p c + a d p c ADK 1: a m p c + a t p c < = > 2. 0 a d p c ADK2: g t p c + a m p c < = > g d p c + a d p c and for which the 13C labeling data cannot constrain the flux values in any fashion. This is shown correctly in the confidence intervals obtained with the method expressed through Eqs 16–23: NDPK1: [ - 497. 4 - 500 ] ADK1: [ - 497. 4 - 500 ] ADK2: [ - 498. 0 - 499. 3 ] whereas if we perform a Monte Carlo, where labeling data is randomly chosen within the experimental error (100 instances), we do not obtain the proper range: NDPK 1: [ - 490. 8 - - 489. 7 ] ADK 1: [ 500. 0 - 500. 0 ] ADK2: [ - 498. 0 - - 498. 0 ] The input for this step is the flux profile obtained in the “Fit data” step and the outputs are the flux profile with corresponding confidence intervals, as shown in S4–S12 Figs.
While metabolic fluxes constitute the most direct window into a cell’s metabolism, their accurate measurement is non trivial. The gold standard for flux measurement involves providing a labeled feed where some of the carbon atoms have been substituted by isotopes with higher atomic mass (13C instead of 12C). The ensuing labeling found in intracellular metabolites is then used to computationally infer the metabolic fluxes that produced the observed pattern. However, this procedure is typically performed with small metabolic models encompassing only central carbon metabolism. The genomic revolution has afforded us easily available genomes and, with them, comprehensive genome-scale models of cellular metabolism. It would be desirable to use the 13C labeling experimental data to constrain genome-scale models: these data constrain fluxes very effectively and provide in the labeling data fit an obvious proof that the underlying model correctly explains measured quantities. Here, we introduce a rigorous, self-consistent method that uses the full amount of information contained in 13C labeling data to constrain fluxes for a genome-scale model where underlying assumptions are explicitly stated.
A Method to Constrain Genome-Scale Models with 13C Labeling Data