Abstract
Background
Metabolic fluxes provide invaluable insight on the integrated response of a cell to environmental stimuli or genetic modifications. Current computational methods for estimating the metabolic fluxes from ^{13}C isotopomer measurement data rely either on manual derivation of analytic equations constraining the fluxes or on the numerical solution of a highly nonlinear system of isotopomer balance equations. In the first approach, analytic equations have to be tediously derived for each organism, substrate or labelling pattern, while in the second approach, the global nature of an optimum solution is difficult to prove and comprehensive measurements of external fluxes to augment the ^{13}C isotopomer data are typically needed.
Results
We present a novel analytic framework for estimating metabolic flux ratios in the cell from ^{13}C isotopomer measurement data. In the presented framework, equation systems constraining the fluxes are derived automatically from the model of the metabolism of an organism. The framework is designed to be applicable with all metabolic network topologies, ^{13}C isotopomer measurement techniques, substrates and substrate labelling patterns.
By analyzing nuclear magnetic resonance (NMR) and mass spectrometry (MS) measurement data obtained from the experiments on glucose with the model microorganisms Bacillus subtilis and Saccharomyces cerevisiae we show that our framework is able to automatically produce the flux ratios discovered so far by the domain experts with tedious manual analysis. Furthermore, we show by in silico calculability analysis that our framework can rapidly produce flux ratio equations – as well as predict when the flux ratios are unobtainable by linear means – also for substrates not related to glucose.
Conclusion
The core of ^{13}C metabolic flux analysis framework introduced in this article constitutes of flow and independence analysis of metabolic fragments and techniques for manipulating isotopomer measurements with vector space techniques. These methods facilitate efficient, analytic computation of the ratios between the fluxes of pathways that converge to a common junction metabolite. The framework can been seen as a generalization and formalization of existing tradition for computing metabolic flux ratios where equations constraining flux ratios are manually derived, usually without explicitly showing the formal proofs of the validity of the equations.
Background
From microorganisms to animals and plants, cells adjust their metabolic operations to fulfill the demand of energy and biosynthetic precursors. In nature this is a challenging task because substrate availability is typically limited and often changing in its composition. To ensure viability on a broad palette of chemically heterogeneous substrates, cells have developed intertwined enzymatic networks that also confer robustness against genetic mutations. Optimum redistribution of molecular fluxes in metabolism is achieved by multilevel regulation circuits. In recent years, experimental measurement of in vivo metabolic fluxes has attracted much attention. For example, in biotechnology metabolic fluxes are utilized to lead rational strain engineering, whereas systems biologists assess fluxes to unravel targets and mechanisms of metabolic regulation.
Metabolic fluxes are often estimated using flux balance analysis (FBA) [1,2]. In FBA, fluxes are solved by fixing some objective for the metabolism of an organism, such as maximal growth. Then, a corresponding optimization problem is solved by using the stoichiometry of the metabolic network as a constraint to the optimization. FBA is a viable method for studying the metabolic capabilities of an organism, but as a method for estimating metabolic fluxes it has some weaknesses. First, selecting the correct objective for the metabolism is far from trivial [3], especially when mutant strains or behaviour in exceptional environmental conditions is analyzed. Second, there can be many biologically interesting flux distributions that give an optimal solution to the optimization problem of FBA.
A more direct method for experimental determination of the metabolic fluxes is to feed an organism with ^{13}C labelled substrate, observe the fate of ^{13}C atoms in the cell at isotopomeric steady state with mass spectrometry (MS) or nuclear magnetic resonance (NMR) measurements, and then infer the metabolic fluxes from the measurements. The rationale behind these ^{13}C tracer experiment is that, often alternative pathways between metabolites in the network manipulate the carbon backbones of the metabolites differently, thus inducing different ^{13}C labelling patterns to metabolites. Then, constraints to fluxes complementary to the basic stoichiometric constraints can be derived by measuring the relative abundances of different labelling patterns in the metabolites.
The main difficulty in applying the procedure in practice is that current measurement techniques only can produce incomplete information about relative abundances of different ^{13}C labelling patterns, the isotopomer distributions, of some metabolites, usually protein bound amino acids in the network, and no isotopomer information at all for many intermediate metabolites of interest [46]. This imposes a highly nonlinear dependency between the measured isotopomer distributions of the metabolites and the metabolic fluxes, which is very challenging to solve both computationally and statistically.
Currently, two main approaches for ^{13}C metabolic flux analysis exist. In the global isotopomer balancing approach, the problem of estimating metabolic fluxes from the isotopomer measurements is formulated as a nonlinear optimization problem, where candidate flux distributions are iteratively generated until they fit well enough to the experimental ^{13}C labelling patterns [711]. Global isotopomer balancing is a versatile approach that can be applied with all network topologies, substrate labelling distributions and with all measurement techniques – also in short time scales where isotopomeric steady state is not reached [1214]. However, due to the nonlinearity of the problem, it is hard to make sure that the optimization has converged to a global optimum and that this optimum is unique [15]. Also, to apply the global isotopomer balancing approach successfully, one usually needs comprehensive information on the uptake and production rates of external metabolites, as well as about biomass composition of the cell. This information can be hard to obtain, especially in largescale experiments where dozens to hundreds of mutants or different organisms are comparatively analyzed [16,17].
A metabolic flux ratio approach (METAFoR) [4,18,19] for ^{13}C metabolic flux analysis has traditionally relied more on the expertise of a domain specialist than advanced computational techniques. In metabolic flux ratio analysis, the aim is to write linear equations constraining the ratios of fluxes producing the same metabolite. The equations are manually derived by domain experts, by careful (and tedious) inspection of metabolic networks. The motivation for the approach is that, in many cases, the knowledge about the flux ratios already offers enough information about the response of an organism to its environment.
The ratio of competing fluxes or pathways producing the same metabolite is easy to understand, and in many cases easier to estimate reliably than all the fluxes in the network – some interesting flux ratios might be obtainable from scarce measurement data or from the incomplete model of metabolic network that would not allow a reliable estimation of a complete flux distribution using global isotopomer balancing. Flux ratios can also be obtained without knowing the uptake and production rates of external metabolites of the biomass composition. And, if enough nonredundant flux ratios are identified, it is possible to use this information to construct and solve an equation system for the full flux distribution of the metabolic network [2022].
As a downside, manually derived flux ratio equations depend heavily on the topology of a metabolic network, measurement capabilities and substrate labelling distributions. Thus, each time a new organism or new mixture of substrates are introduced, flux ratio equations have to be verified and new ones possibly derived. To date, flux ratio equations are manually derived for central carbon metabolisms of three model organisms on glucose, S. cerevisiae [17,23], B. subtilis [24] and for Escherichia coli [18,19]. Recently, flux ratio equations of S. cerevisiae were modified for Pichia pastoris grown on glycerol and on glycerol/methanol mixtures [25,26]. Facilitating the process of deriving flux ratio equations for other organisms and other substrates clearly calls for automatic tools. Also, many times the (simplifying) assumptions made by the expert in the derivation and solution of flux ratio equations, are not reported in detail. Thus, it is often nontrivial to verify the correctness of given flux ratio equations.
In this article we present a novel automatic framework for deriving flux ratios when the measurement data and the model of metabolic network are given as input. The framework is based on the graph algorithmic flow analysis of metabolite fragments in the metabolic network [27] and on the interpretation and manipulation of both NMR and MS data with vector space techniques [21]. The goal of our work is to combine the good aspects of global isotopomer balancing and manual flux ratio analysis: like global isotopomer balancing, our framework is systematic and can be applied with all network topologies, substrates and substrate labelling distributions and with all current isotopomer measurement techniques. Thus, laborious and errorprone manual inspection of metabolic network models and the tailoring of the equation systems constraining the fluxes separately for each experimental setting required in manual flux ratio techniques can be avoided. On the other hand, during the automated construction of flux ratios we resort to linear optimization techniques only, combined with graph algorithms of polynomial worst case time complexity. Thus, our framework is computationally efficient and avoids problems with local and multiple optima frequently met in global isotopomer balancing. The tradeoff of this philosophy is, however, the requirement of measuring isotopomer distributions of metabolites more rigorously to obtain full flux distribution. Given insufficient measurements, our framework can solve the flux ratios only for some, but not necessary for all metabolites in the network. We expect that, especially as measuring isotopomers of intermediate metabolites becomes more routine, our framework will be an attractive method for ^{13}C flux analysis.
Results and Discussion
In this section we demonstrate the applicability of the presented framework by empirical results. We show that our automatic and systematic framework is able to reproduce flux ratios previously determined by a manual analysis from NMR and GCMS isotopomer measurements of protein bound amino acids of S. cerevisiae and B. subtilis on glucose. Thus, we can conclude that the presented framework is powerful enough to provide interesting flux ratio information in the well studied experimental settings. Furthermore, we show that the framework can be applied to study less known experimental conditions without any further effort by discovering the flux ratios that can be estimated when B. subtilis is grown on malate instead of glucose. The results of this analysis show that our framework can detect profound effects the change of external substrate can have to the flux ratio computations. The results indicate that our framework is a good tool to study flux ratios of microbes in different experimental conditions – a claim that will try to validate with more experiments in our further work.
We obtained NMR and GCMS labelling data, where isotopomer distributions of protein bound amino acids of S. cerevisiae and B. subtilis grown on different conditions were measured. Then, available flux ratios were computed with the presented framework. Models of metabolic networks applied in the analysis can be found from the supplementary material of this article: additional files 1 and 2 contain the SBML model file and a visualization of the model of S. cerevisiae, while additional files 3 and 4 contain the same information for B. subtilis. In the models, some simplifications common to ^{13}C metabolic flux analysis were applied by pooling metabolites whose isotopomer pools can be assumed to be fully mixed (cf. [28]). Pooling of metabolites was carried for (i) the three pentosephosphates in PPP, (ii) phopshotrioses between GA3P and PEP in glycolysis, and (iii) oxaloacetate and malate in the TCA cycle. In these cases, pooling is justified by the existence of fast equilibrating, bidirectional reactions between the listed intermediates and the empirical evidence that their isotopic labeling is not distinguishable with the current analytical tools. Cofactor metabolites were excluded from the model as cofactor specificities and activities are not accurately known for many reactions.
Additional file 1. Model of central carbon metabolism of S. cerevisiae in SBML format. Additional file 1 contains the stoichiometric model of central carbon metabolism of S. cerevisiae used in the experiments of the article. The model is provided as a text file containing SBML compliant descriptions of metabolites and reactions in the model. Metabolites are identified by their KEGG LIGAND codes.
Format: TXT Size: 43KB Download file
Additional file 2. Illustration of the model of central carbon metabolism of S. cerevisiae. Bidirectional reactions in PPP pathways are depicted as unidirectional for better readability.
Format: PDF Size: 22KB Download file
This file can be viewed with: Adobe Acrobat Reader
Additional file 3. Model of central carbon metabolism of B. subtilis in SBML format. Additional file 2 contains the stoichiometric model of central carbon metabolism of B. subtilis used in the experiments of the article. The model is provided as a text file containing SBML compliant descriptions of metabolites and reactions in the model. Metabolites are identified by their KEGG LIGAND codes.
Format: TXT Size: 36KB Download file
Additional file 4. Illustration of the model of central carbon metabolism of B. subtilis. Bidirectional reactions in PPP pathways are depicted as unidirectional for better readability.
Format: PDF Size: 19KB Download file
This file can be viewed with: Adobe Acrobat Reader
The bulk of the carbon mappings of reactions in the metabolic network were provided by ARM project [29]. Carbon mappings from amino acids to their precursors were conform to [4] and [23]. Before the analysis of the real measurement data, the correctness of the implementation of the framework was empirically verified by estimating flux ratios for junction metabolites in the metabolic network of S. cerevisiae from the artificial data generated by the 13CFLUX software [8].
NMR measurements from S. cerevisae on glucose
In the first experiment we analyzed NMR isotopomer measurement data from protein bound amino acids of S. cerevisiae that was grown on uniformly labelled glucose (see Section Experimental NMR and GCMS methods for more details on experimental settings).
From the 15 measured amino acids we were able to estimate flux ratios for seven junction metabolites: oxaloacetate, PEP, glycine and serine on cytosol and for oxaloacetate, acetylCoA and pyruvate in mitochondria. Furthermore, an upper bound for a ratio of GA3P molecules that have visited the transketolase reaction was obtained by manually simplifying the model to imitate the previously reported ways to manually compute the corresponding upper bound (cf. [4]). (The structural analysis of the metabolic network model described in Section Structural analysis of isotopomer systems can help in discovering such simplifications, but they also need some expert insight. As the simplifications are currently not done automatically, the systematical framework is unable to validate them.)
The computed flux ratios were compared with the manually derived ratios [23], when the assumptions made in the manual derivation of flux ratios were consistent with the model used. In all cases, automatically computed flux ratios agreed well with the manually derived ratios (Table 1). Differences between the estimations can be explained by numerical instabilities and by differences in computational procedures: in manually derived ratios the estimations are based on the breakage of a single bond in different routes leading to a metabolite while in our framework more isotopomer information is possible utilized in the estimation.
Table 1. Estimated flux ratios from NMR measurements of S. cerevisiae.
GCMS measurements from B. subtilis on glucose
In the second experiment we analyzed GCMS isotopomer measurement data from protein bound amino acids of Bacillus subtilis that was grown on uniformly labelled glucose (see Section Experimental NMR and GCMS methods for more details on experimental settings).
In comparison to eukaryotic S. cerevisiae, the metabolic network of prokaryotic B. subtilis lacks cellular compartments. Thus, from the point of view of ^{13}C metabolic flux analysis, there are fewer interesting junction metabolites in the central carbon metabolism of B. subtilis where the flux ratios can be estimated. From the GCMS measurements of 14 amino acids we were able to compute flux ratios for four junction metabolites – oxaloacetate, pyruvate, PEP and glycine – when [U13C]glucose was used as a carbon source. Furthermore, an upper bound for a ratio of GA3P molecules that have visited transketolase reaction was obtained by manually simplifying the model of the metabolic network. Excluding pyruvate, we were able to compute the same ratios with [113C]glucose as a carbon source.
We compared the computed flux ratios to ones obtained with the software FiatFlux [30] that is based on the manually derived analytic equations for computing flux ratios. Currently, manually derived flux ratio equations for [113C]glucose as a carbon source exist only for PEP and for the upper bound to the flux through oxidative pentose phosphate pathway. In general, the flux ratios computed with different methods from the same data and with the same assumptions about the topology of metabolic network were in good agreement (Table 2). (As a data cleaning procedure, we removed from [113C]glucose data the mass distributions of fragments whose fractional enrichment deviated more than 5% from the assumed fractional enrichment of 20% in [U13C]data. This was done because differences in fractional enrichments can be tracked in uniformly labelled data where the fractional enrichment of each fragment is know a priori, but not in positionally labelled data, where the fractional enrichment of a fragment depends on the network topology and the fluxes. This kind of irregularities are in general caused by noise in fragments with low intensity or by coeluting analytes with overlapping fragment masses.) Again, differences between the flux ratios estimated by different methods can be explained by numerical instabilities and by differences in isotopomer information applied during the estimation. Variation in the estimated flux ratios between repeated experiments (six repetitions for [113C]glucose experiment, four repetitions for uniformly labelled glucose experiment) was negligible.
Table 2. Estimated flux ratios from GCMS measurements of B. subtilis.
In silico calculability analysis of B. subtilis on malate
One of the strengths of the presented framework is that it is able to automatically produce metabolic flux ratios also when other external labelled substrates than commonly used glucose are fed to organisms. To demonstrate this ability, we applied our framework to predict what flux information would be available, if we feed B. subtilis with malate.
We applied the in silico calculability analysis (see Section Calculability analysis) to examine which flux ratios are calculable in the best case from GCMS measurements of amino acids, when B. subtilis is grown on [U13C]labelled malate. Interestingly, our fragment flow analysis revealed that – with the reaction reversibilites in the applied model – the isotopomer distributions of GA3P, PEP and pyruvate depend only on the isotopomer distribution of the fragment containing the first three carbons of oxaloacetate, but not on the relative fluxes producing these metabolites. Thus, isotopomer balances for GA3P, PEP and pyruvate reduce to mass balances and the ratios of fluxes producing these metabolites cannot be estimated. This somewhat surprising phenomenon is due the fact that the rearrangements of carbon chains occurring in PPP pathway will affect only to the carbon fragments that will be recycled in the upper metabolism but not the carbon fragments that can flow back GA3P, PEP and pyruvate from PPP (we modelled a reaction from GA3P to F6P as unidirectional one).
Preliminary experiments with GCMS data from B. subtilis grown on [U13C]labelled malate agreed with the results of fragment flow analysis: constraints to the isotopomer distributions of fluxes entering to these metabolites were identical within the limits of measurement accuracy. On the other hand, our framework was able to estimate for example the TCAcycle activity also when B. subtilis is grown on malate, just as predicted by the calculability analysis.
Conclusion
In this article we introduce a systematic and analytic framework for ^{13}C metabolic flux ratio analysis. At the heart of the framework lie the techniques for flow analysis of a metabolic network and for manipulating isotopomer measurements as linear subspaces. These techniques facilitate the efficient and analytic computation of the ratios between the fluxes producing the same junction metabolite. The framework can be seen as a generalization and formalization of existing analytic methods for computing metabolic flux ratios [23,30,31] where equations constraining flux ratios are manually derived. Like the recent methods to improve the speed of the simulation of isotopomer distributions in the global isotopomer balancing framework [10,32], our framework relies on graph algorithms. However, both our goals and applied techniques are quite different from these approaches. In [10] and [32] connected components of isotopomer graphs are discovered to divide the simulation of isotopomer distributions to smaller subtasks. In our framework, flow analysis techniques are applied to discovered metabolite fragments with equivalent isotopomer distributions in every isotopomeric steady state.
Our experiments with NMR and MS data show that the framework is able to provide relevant information about metabolic fluxes, even when only constraints to the isotopomer distributions of proteinbound amino acids are measured.
Thanks to recent advancements in measurement technology improving the feasibility of mass isotopomer measurements of intermediate metabolites [13,33], we expect that the full power of the framework will be harnessed in near future. Measurements from intermediates will make it possible to use larger models of metabolic networks and estimate flux ratios more accurately, without simplifying assumptions about the topology of the metabolic network or directionality of the fluxes. However, these measurements will not be easy to conduct, because of the low abundances of intermediates in the cell. Thus, systematic methods for experimental planning and data quality control are required. The presented framework provides powerful tools for these tasks. First, the framework facilitates time saving in silico calculability analysis.
Second, the manipulation of isotopomer measurements as linear subspaces offers a natural way for comparing measurements from different metabolites that contain overlapping information to detect inconsistencies in the measurements: it is enough to compare propagated isotopomer information in the fragments that belong to the same equivalence class. Third, as MS isotopomer measurement techniques have to be developed separately for different intermediate metabolites or metabolite classes, it will be very useful to select a small subset of intermediates that gives enough information about the interesting metabolic fluxes with least experimental effort. In future research, we want to tackle this problem by generalizing our earlier experimental planning method [34], to all measurement data and to realistic measurement error models.
As the presented framework for ^{13}C metabolic flux analysis only resorts to linear optimization techniques, it is not always able to provide as much information about the metabolic fluxes as the global isotopomer balancing frameworks utilizing more powerful, nonlinear optimization techniques [8,35], that do not necessarily converge to the global optimum. On the other hand, some flux ratios might be computable from the scarce data or incomplete model of metabolic network that does not allow global isotopomer balancing. The differences in the practical performance of different approaches require further research. We see these alternative approaches as complementary ones. A very nice goal would be an integration of our work with global isotopomer balancing: our analytic flux ratios could speed up and direct the optimization process of global isotopomer balancing, that would then fill in the flux ratios possibly missed by our framework.
Methods
In this section we describe a systematic and analytic computational framework for ^{13}C metabolic flux analysis. At the end of the section we also shortly describe the experimental method that were applied to produce the isotopomer measurement data that was analyzed in Section Results and Discussion.
The overall goal of our computational framework is to automatically infer from the available isotopomer measurement data produced by MS, MSMS or NMR techniques, an equation system constraining the fluxes. The crucial question is to derive as many nonredundant equations as possible, ideally constraining the flux distribution to a point solution, or in general, as lowdimensional convex set as possible.
In short, the framework consists of the following steps:
1. The model of the metabolic network of an organism is constructed by selecting a set of biochemical reactions operating in the organism and by designating them to correct cellular compartments;
2. Structural analysis of the isotopomer system is conducted, consisting of the following steps:
(a) Flow analysis of the metabolic network is conducted in order to discover equivalent fragments, fragments of carbon backbones of metabolites that will have the same theoretical isotopomer distribution, regardless of the fluxes.
(b) Independence analysis of fragments is conducted to find statistically independent carbon subsets from metabolites, that is, subsets that have been at some point separated along every pathway able to producing them and that have flux invariant isotopomer distributions. This guarantees that the isotopomer distribution of their union assumes the form of a product distribution.
(c) In silico calculability analysis is performed to test if the available measurement techniques and substrate labellings make it in principle possible to obtain the required flux information.
3. Wetlab isotopomer tracer experiments are conducted and constraints to isotopomer distributions are measured;
4. The fluxes of the network are estimated. The process consists of the following steps:
(a) Isotopomer measurement data is propagated in the metabolic network model from the measured metabolites to unmeasured ones according to the equivalences discovered in step 2.
(b) An equation system tying the isotopomer data and the fluxes together is constructed and solved, either to obtain a flux distribution for the metabolic network as a whole, or for a single junction metabolite to obtain the ratios of fluxes producing it.
(c) The statistical analysis of obtained fluxes or flux ratios is carried out.
In the following we first formalize the ^{13}C flux analysis problem and then detail the computational steps above.
Preliminaries
In ^{13}C metabolic flux analysis the carbon atoms of metabolites are of special interest. We denote with M the set of carbon locations M = {c_{1}, ..., c_{k}} of a kcarbon metabolite. By M = k we denote the number of carbons in M. Fragments of metabolites are subsets F = {f_{1}, ..., f_{h}} ⊆ M of the carbons of the metabolite. A fragment F of M is denoted as MF. A metabolic network G = (, ) is composed of a set = {M_{1}, ..., M_{m}} of metabolites and a set = {ρ_{1}, ..., ρ_{n}} of reactions that perform the interconversions of metabolites.
With isotopomers we mean molecules with similar element structure but different combinations of ^{13}C labels. Isotopomers of the molecule M = {c_{1}, ..., c_{k}} are represented by binary sequences b = {b_{1}, ..., b_{k}} ∈ {0, 1}^{k }where b_{i }= 0 denotes a ^{12}C and b_{i }= 1 denotes a ^{13}C in location c_{i}. Molecules that belong to the bisotopomer of M are denoted by M(b). Isotopomers of metabolite fragments MF are defined in an analogous manner: a molecule belongs to the F(b)isotopomer of M, denoted MF(b_{1}, ..., b_{h}), if it has a ^{13}C atom in all locations f_{j }that have b_{j }= 1, and ^{12}C in other locations of F. Isotopomers with equal numbers of labels belong to the same mass isotopomer. We denote mass isotopomers of M by M(+p), where p ∈ {0, ..., M} denotes the number of labels in isotopomers belonging to M(+p).
The isotopomer distribution D_{M }of metabolite M gives the relative abundances 0 ≤ P_{M}(b) ≤ 1 of each isotopomer M(b) in the pool of M such that
The isotopomer distribution D_{MF }of fragment MF and the mass isotopomer distribution of mass isotopomers M(+p) are defined analogously: D_{MF }of metabolite M gives the relative abundances 0 ≤ P_{MF}(b) ≤ 1 of each isotopomer MF(b) and gives the relative abundances 0 ≤ P_{M}(+p) ≤ 1 of each mass isotopomer M(+p).
Reactions are pairs ρ_{j }= (α_{j}, λ_{j}) where α_{j }= (α_{1j}, ..., α_{mj}) ∈ ℤ^{m }is a vector of stoichiometric coefficientsdenoting how many molecules of each kind are consumed and produced in a single reaction eventand λ_{j }is a carbon mapping describing the transition of carbon atoms in ρ_{j }(see Figure 1). Metabolites M_{i }with α_{ij }< 0 are called substrates and with α_{ij }> 0 are called products of ρ_{j}. If a metabolite is a product of at least two reactions, it is called a junction. If α_{ij }< 0, a reaction event of ρ_{j }consumes α_{ij} molecules of M_{i}, and if α_{ij }> 0, it produces α_{ij} molecules of M_{i}. Bidirectional reactions are modelled as a pair of reactions.
Figure 1. An example of a metabolic reaction. In the reaction ρ_{j}, a fructose 1,6bisphosphate (C_{6}H_{14}O_{12}P_{2}) molecule is produced from glycerone phosphate (C_{3}H_{7}O_{6}P) and glyceraldehyde 3phosphate (C_{3}H_{7}O_{6}P) molecules. Carbon maps are shown with dashed lines. Glyceraldehyde 3phosphate is equivalent to the gray fragment of fructose 1,6bisphosphate while glycerone phosphate is equivalent to the white fragment.
A pathway p in network G from metabolite fragments {F_{1}, ..., F_{k}} to fragment F' is a sequence of reactions such that a composite carbon mapping , defined by p maps the carbons of {F_{1}, ..., F_{p}} to the carbons of F'.
For the rest of the article, it is important to distinguish between the subpools of a metabolite pool produced by different reactions. Therefore, we denote by M_{ij}, the subpool of the pool of M_{i }produced (α_{ij }> 0) or consumed (α_{ij }< 0) by reaction ρ_{j}. The concept of the subpools of a metabolite pool is illustrated in Figure 2.
Figure 2. An example of subpools of a metabolite pool. Phosphoenolpyruvate (PEP) is produced by two different reactions (ρ_{i }and ρ_{j}), either from Oxaloacetate (OAA) or from glyceraldehyde 3phosphate (GA3P). Thus, PEP has two in flow subpools, PEP from OAA and PEP from GA3P (grey boxes) that are mixed in the common PEP pool (at the bottom).
By M_{i0 }we denote the subpool of M_{i }that is related to the external in flow or external out flow of M_{i}. We call the sources of external in flows external substrates. Subpools of fragments are defined analogously. In ^{13}C metabolic flux analysis, the quantities of interest are the rates or the fluxes v_{j }≥ 0 of the reactions ρ_{j}, giving the number of reaction events of ρ_{j }per time unit. We denote by v the vector [v_{1}, ..., v_{n}] of fluxes, or a flux distribution.
Generalized isotopomer balances
The framework for ^{13}C metabolic flux analysis presented in this article rests on the assumption that the metabolic network is in metabolic and isotopomeric steady state. In the metabolic steady state, metabolite balance, or mass balance
holds for each metabolite M_{i}. Here, β_{i }is the measured external in flow (β_{i }< 0) or external out flow (β_{i }> 0) of metabolite M_{i}. From balance equations (1) defined for every metabolite M_{i }we will obtain a metabolite balancing, or stoichiometric equation system.
In isotopomer steady state, for each isotopomer b of each metabolite M_{i }in the metabolic network the following isotopomer balance holds:
For metabolic flux analysis (1) and (2) bear a fundamental difference: the former cannot be used to estimate fluxes of alternative pathways producing M_{i }while the latter can. However, using (2) is not in general admissible in practice: typically abundances P_{M}(b) of isotopomers are not fully determined from the measurements, and we need to settle for some constraints to the distribution D_{M}. A crucial building block of our framework is the representation of isotopomer measurements as systems of linear equations (c.f. [21])
where s_{bh }is the weight of isotopomer b in the h'th constraint, d_{h }is a value derived from isotopomer measurements, and r is the total number of constraints. We call (3) isotopomer constraints. For a k carbon metabolite, 2^{k }linearly independent isotopomer constraints – one for each isotopomer – are necessary and sufficient to constrain the isotopomer distribution D_{M }to a point solution. A set of isotopomer constraints has a natural matrix representation SD_{M }= d, where S = (s_{bh})_{b,h }is a 2^{k }by r matrix, where 1 ≤ r ≤ 2^{k }is the number of isotopomer constraints (the trivial constraint Σ_{b}P_{M}(b) = 1 by definition always holds).
The use of (3) follows from the simple observation that isotopomer balance (2) implies that each linear combination of isotopomers is balanced. Thus, we can write a new balance equation that constrains the fluxes producing M_{i }as soon as we know the value of the same linear combination of isotopomer abundances for each subpool of M_{i}. We have
where each is a linear combination of the form (3), with coefficients s_{b }that do not depend on j, i.e. they are the same for each subpool M_{ij}. We call (4) a generalized isotopomer balance.
Representing MS and NMR isotopomer measurements as linear constraints
In the following, we will show by examples how MS and NMR data can be represented as isotopomer constraints. Let us first consider mass isotopomer distributions obtained from MS. Here we omit discussion on practicalities such as corrections for natural abundances of ^{13}C isotopes (c.f., [6,36]) and concentrate on the conceptual level. For example, the +2 mass isotopomer of a three carbon metabolite M satisfies
which conforms to (3) by taking s_{011,2 }= s_{110,2 }= s_{101,2 }= 1, and s_{b,2 }= 0 otherwise. Similarly, the coefficients s_{bk }can be derived for all mass isotopomers +k, k = 0, ..., 3.
Isotopomer data originating from Tandem MS, or MSMS, falls into the same representation. Consider, for example, a fragment MF of a threecarbon metabolite, containing the first and second carbon of M. The following equation holds for the mass isotopomer MF(+1):
The equation can be written in terms of the precursor M, but the exact form of the equation depends on the mode of tandem MS. If the full scan mode is used, we obtain
as all precursor molecules M that have exactly one carbon among the first and second location contribute to the fragment mass isotopomer MF(+1). On the other hand, in the daughterion scanning mode a single mass isotopomer, for example M(+2), is selected as the precursor. Then we obtain
as the precursor must always have two ^{13}C atoms in total. We refer the reader to [6,36] for further details. Also NMR ^{13}C isotopomer measurements, where relative intensities of different combinations of ^{13}C and ^{12}C atoms that are coupled to an observed ^{13}C atom are measured, can be expressed as linear combinations of isotopomer abundances. For example, for a threecarbon metabolite M, the following constraints to can be inferred for the labeling pattern 010:
where d(010) is the measured intensity. Rewriting the above as
and denoting s_{b }= d(010), for b ∈ {110, 011, 111}, s_{010 }= d(010)  1 and s_{b }= 0 for b ∈ {000, 100, 001, 101}, the above can be seen to conform to (3). Similar derivation can be used for other isotopomer signals present in the NMR spectrum to obtain the corresponding isotopomer constraints.
Projection of isotopomer measurements to fragments
In our computational framework, it will be necessary to project the measurement data obtained for a metabolite M to its fragments MF and vice versa. In this projection, we want to avoid any unnecessary loss of measurement information, that is, we want to obtain as many linearly independent constraints to the isotopomer distribution of F as possible. For example, if we have measured that a twocarbon metabolite M has the isotopomer distribution
and we need to know the isotopomer distribution of the fragment MF consisting of the first carbon of M, we should obtain
For the general model of isotopomer measurements (4) the projection of measurement information from a metabolite to its fragments can be done by the techniques of linear algebra introduced in [21]. We recapitulate the techniques in the following.
Recall the general form of isotopomer measurement SD_{M }= d, where S denotes a matrix with 2^{k }columns, one column for each isotopomer b of kcarbon metabolite M, and each row h of S corresponds to a measured isotopomer constraint (3). The rows of S span a subspace in a 2^{k }dimensional vector space spanned by all possible isotopomer distributions D_{M}.
Also the metabolite fragments are naturally represented as vector subspaces. Let U_{F }denote a matrix with also a column for each isotopomer M(b) and a row for each isotopomer F(b') of MF, that is,
The rows of U_{F }span another subspace . Any isotopomer distribution D_{MF }lies in this subspace, and hence also any isotopomer constraint S_{F}D_{MF }= d for fragment MF necessarily lies in the same subspace.
In conclusion, the available information about D_{M }is given as its linear projection onto , and anything we can express about D_{MF }in terms of isotopomer constraints is contained within . Thus, any isotopomer constraint for D_{MF }that we can derive from the measurements can be expressed in terms of the vector space intersection .
Thus, to obtain isotopomer constraints for fragment MF from a measurement SD_{M }= d, we need to compute the vector space intersection and project the measurement to . This can be done by standard linear algebra (c.f. [21]). This process gives us as output isotopomer constraints of the required form
Finally, transforming a fragment constraint Y_{F}D_{MF }= d_{F }into an isotopomer constraint SD_{M }= d is easy: we postmultiply the fragment constraint with the matrix U_{MF }: S = Y_{F}U_{MF }and d = d_{F}U_{MF}.
Structural analysis of isotopomer systems
The incomplete nature of ^{13}C measurement data requires us to find ways to use the available data the best way possible. The central concept is to find invariants of isotopomer distributions that remain through the pathways, and allow us to trade or propagate measurement information from one metabolite to another. This allows us to write or augment generalized isotopomer balances for metabolites for which the isotopomer distributions are not completely determined by measurements. Thus, the fluxes are potentially more tightly pinpointed as well.
In particular, we use two techniques: First, flow analysis is used to uncover sets of metabolite fragments that have the same isotopomer distribution regardless of the fluxes. Second, independence analysis of fragments is used to uncover situations where two fragments of the same metabolite induce the product distribution for the isotopomer distribution of their union.
Flow analysis of metabolic networks
The goal of the flow analysis [27] is to partition the fragments of the metabolites in the network to equivalence classes such that fragments in the same equivalence class have identical isotopomer distributions in every steady state. This can be guaranteed if a fragment is produced from another a such a manner that the carbons within the fragment never depart from each other regardless of the pathway that is being used.
Formally, we say that fragment F' dominates fragment F if the following conditions are met
1. F and F' have the equal number of carbons;
2. all carbons of F originate always from the carbons of F';
3. carbon of F' stay connected to each other via all pathways from F' to F;
4. composite carbon mappings are the same in all pathways from F' to F.
Intuitively, a dominated fragment (F) is always produced from its dominator (F') without manipulating the carbon chain of the fragment. Thus, isotopomer distribution of the dominated fragment does not contain any information about the metabolic fluxes. For a fragment F that has no dominators, the transitive closure of the domination relation corresponds to the class of equivalent fragments in the network.
The simplest example of fragment equivalence is the one between a substrate M_{k }and product M_{i }in a single reaction ρ. If the atoms in M_{i}F originate from M_{k}F', then the fragments M_{k}F' in the subpool M_{ij }produced by reaction ρ_{j}, are equivalent with the fragment M_{i}F (Figure 1). Furthermore, if metabolite M_{i }has only one producing reaction ρ_{j}, isotopomer distributions of subpool M_{ij }and M_{i }coincide. Thus, if fragment M_{i}F is produced from a single fragment M_{k}F' of some substrate M_{k }of ρ_{j}, F and F'are equivalent. By transitivity, all fragments in the linear pathway are equivalent.
More complicated case of fragment equivalence is found when a fragment of a junction metabolite is dominated by an upstream fragment (Figure 3). In [27] we show that the the classes of equivalent fragments corresponding the conditions (1–4) can be efficiently computed. Very brie fly, first the metabolic network is transformed to a fragment flow graph that connects substrate metabolite fragments to their product fragments for each reaction in the network. Then, the dominator tree [37,38] of the fragments in the fragment flow graph is constructed. It turns out that the subtrees of this dominator tree correspond to the required fragment equivalence classes (see Figure 4).
Figure 3. An example of fragment equivalence classes in a branched pathway. An example of equivalence classes of fragments in the metabolic network that contains dominated junction fragments ME and MF. Grey and white fragments constitute two equivalence classes. Dashed lines depict carbon mappings.
Figure 4. An example of a fragment flow graph and a dominator tree. A metabolic network (left), the corresponding fragment flow graph (up right) and the subtrees of the dominator tree (down right).
Fragment equivalence classes have many uses [27]. Most importantly, measured isotopomer constraints to fragment F can be directly propagated to another fragment F' in the same equivalence class, by applying the joint carbon mappings between F and F'. This helps in the construction of generalized balance equations (4) where the same isotopomer information is required for each subpool of junction metabolites.
Independence analysis of fragments
A complementary property to fragment equivalence ^{13}C is the statistical independence of fragment isotopomer distributions. Intuitively, if two fragments of the same metabolite are statistically independent, new isotopomer constraints to the union of them can be obtained by taking a product of the isotopomer distributions of the independent fragments.
More formally, the basic question is on what conditions the distribution D_{ME∪F }of union of two fragments will necessarily have the form of a product distribution
where ⊗ denotes the tensor product consisting all terms of the form P_{E∪F}(b) = P_{E}(b')·P_{F}(b"), where b' (resp. b") ranges over all isotopomers of E (resp. F), and b is the isotopomer of ME ∪ F formed by joining E and F.
The utility of fragment independence is in that it gives us constraints to the isotopomer distributions that are complementary to the isotopomer constraints (3) obtained from the measurements.
In general, two criteria need to be satisfied for statistical independence of two fragments ME and MF. First, the fragments need to be structurally independent, meaning that along all pathways producing the metabolite, at some point all carbons of the fragments have resided in different metabolite molecules. This property can be defined in recursive manner. Fragments ME and MF are structurally independent if for all carbon pairs (a, b), a ∈ E and b ∈ F, for all reactions ρ producing M, it holds that
• a and b originate from different reactants of ρ, or
• a and b originate from the same reactant M', and the reactant fragments M'F_{a }and M'F_{b}, where F_{a }= (a), F_{b }= (b), are structurally independent.
The second necessary condition is that the two fragments need to be dominated by some other metabolite fragments in the network. This will make the fragment distributions flux invariant. Together, the two criteria guarantee (6) to hold.
A simple case of statistical independence of fragments is a (subpool) product metabolite M_{i }of a single reaction ρ_{j}, where the fragments M_{i}E and M_{i}F are disjoint and originate from different reactants. The fragments are structurally independent (by originating from different reactants) and are dominated (by reactant fragments of ρ_{j}). Hence (6) holds. The underlying assumption here is that enzymes pick their reactants independently and randomly from the available pools. This case of statistical independence of fragments is depicted Figure 1, where white and grey fragments of Dfructose 1,6biphosphate are statistically independent (in the subpool of reaction ρ_{j}).
Another simple example is a junction metabolite M_{i }that has two or more producers with associated subpools M_{ij}. If M_{ij}E and M_{ij}F are structurally independent in all subpools, M_{i}E and M_{i}F are structurally independent as well. If M_{i}E and M_{i}F are dominated by some fragments in the network, all subpools have the same distribution which takes the form of (6). Without dominance the distribution will in general be a fluxdependent mixture of product distributions . This case of statistical independence is depicted in Figure 5.
Figure 5. An example of statistical independence of fragments. White and grey onecarbonfragments of M_{i }are statistically independent: both fragments are dominated by onecarbonfragments of M, and the fragments are disjoint in every pathway that produce M_{i }from M.
A generalized form of (6), useful for propagation of isotopomer constraints, is derived as follows. Assume independent fragments ME and MF of metabolite M and isotopomer constraints S_{E∪F}D_{ME∪F }= d_{E∪F}, S_{E}D_{ME }= d_{E }and S_{F}D_{MF }= d_{F}, where S = S_{E }⊗ S_{F}. Now, the the h'th constraint for fragment E and g'th constraint for fragment F, written as Σ_{a}s_{ah,E}P_{ME}(a) = d_{h,E }and Σ_{c}s_{cg,F}P_{MF}(c) = d_{g,F}.
Multiplying the constraints, and denoting s_{b }= s_{ah,E}·s_{cg,F }where b is the isotopomer consistent with fragment isotopomers a and c, we get the following equation
for the l'th constraint for E ∪ F. The equations of the above kind can be concisely written in terms of tensors:
From above, if two of the three vectors d, d_{E}, d_{F }are known, the remaining unknown one can be solved.
We note in passing that computing constraints to the metabolite given constraints to the fragments is straightforward application of (7).
Applying fragment independence analysis to flux ratio computation
In our framework, statistical independence of fragments has two uses. We apply it
1. to compute isotopomer constraints for the union of independent fragments, given isotopomer constraints to its independent fragments, and
2. to compute isotopomer constraints for an independent fragment given isotopomer constraints to the other fragments and the metabolite as a whole.
In both cases making use of (6) gives us a larger set of constraints than the vector space and flow analysis approach alone.
Next we describe how (7) generalizes the basic measurement propagation step of the traditional metabolic flux ratio analysis [31]. In the basic case, the flux ratios are solved for two competing pathways p and q, which p cleaves a certain carboncarbon bond b of junction M while the q preserves b intact from the external substrate. (See Figure 6 for an example). This serves also as an example of applying (7) to compute isotopomer constraints for the union of independent fragments.
Figure 6. An example of using fragment independence to obtain new isotopomer constraints under uniform substrate labelling. Constraints to the isotopomer distributions of striped metabolites are assumed to be known, either by direct measurement of measurement propagation. In pathway q = (ρ_{2}, ρ_{4}), the isotopomer distribution of M_{F }molecules will be the same as in M_{E}. In pathway p = (ρ_{1}, ρ_{3}), the isotopomer distribution of M_{F }can be derived by applying fragment independence: the isotopomer distributions of single carbon metabolites produced by ρ_{1 }are known a priori to be equal to the labelling degree of uniformly labelled substrate. As the two carbons of M_{F3 }are produced from two different metabolites, these carbons are statistically independent to each other in the subpool and the isotopomer distribution D() of M_{F }molecules produced by p can be computed by applying Equation 7.
When a uniformly labelled substrate is used, the labelling degree of every carbon in the network is the same (and known a priori) when the system reaches isotopomeric steady state. Thus, the isotopomer distribution of a twocarbon fragment F (metabolite M_{F }in Figure 6) containing bond b can be computed by (7) for pathway p cleaving b, while for pathway q, can be propagated from the external substrate (metabolite M_{E }in Figure 6) using the fragment equivalence classes of the previous section. If we are able to measure (constraints to) the isotopomer distribution of the mixed pool F, we can then automatically derive a generalized isotopomer balance corresponding the manually derived ratio. To use (7) to compute isotopomer constraints for an independent fragment from the known isotopomer constraints to the other independent fragment and to the whole metabolite is complicated by the incompleteness of the measurement data: an arbitrary measurement SD_{ME∪F }= d might not be directly representable via a tensor product S = S_{E }⊗ S_{F}. Instead, we need to first compute isotopomer subspaces for known isotopomer constraints where (7) can be applied.
The detailed description of this technique is rather technical and omitted from this article. Here we give an example of the technique (See Figure 7). We assume that we know the mass isotopomer distributions of metabolite M_{i }(metabolite M_{1 }in Figure 7) and of fragment M_{i}E. We furthermore assume that M_{i}E and M_{i}F are independent. From this information the mass isotopomer distribution of can be solved. To be exact, can be solved from the system containing an equation
Figure 7. An example of using fragment independence to obtain new isotopomer constraints for a reactant. The mass isotopomer distributions of striped metabolites are assumed to be measured. Fragments M_{1}E and M_{2 }belong to the same fragment equivalence class. Thus, D_{m}(M_{1}E) can be derived from D_{m}(M_{2}) by the measurement propagation inside equivalence classes. Furthermore, fragments M_{5}E' and M_{5}F' dominate fragments M_{1}E and M_{1}F, and the bond between M_{1}E and M_{1}F is broken in all pathways producing M_{1 }from M_{5}. Thus, M_{1}E and M_{1}F are statistically independent, and D_{m}(M_{1}F) can be deduced from D_{m}(M_{1}) and D_{m}(M_{1}E) by utilizing Equation 7. Computed D_{m}(M_{1}F) can then be propagated to M_{4}, as M_{1 }and M_{4 }belong to the same fragment equivalence class. Finally, D_{m}(M_{4}) helps to solve the ratios of fluxes entering to M_{3}.
for each mass isotopomer p of M_{i}. To see that (8) conforms to (7), we denote the measured mass isotopomer distribution = S·D(M_{i}) of M_{i }by d_{M }(i.e. rows of coefficient matrix S correspond to different mass isotopomers of M_{i}) and the measured mass isotopomer distributions of E and F by d_{E }and d_{F}. Let E = 2, F = 1 and M_{i} = 3, thus M_{i }= E ∪ F. We have
with the tensor product
As the two matrices are not the same (7) is not directly applicable. However, by summing up the second and the third rows and the fourth and the fifth rows of S_{E }⊗ S_{F }we obtain S. Intuitively, this means that Equation (7) can be applied to compute , when we take into account that the isotopomer constraints corresponding both second and third rows of S_{E }⊗ S_{F }contribute to mass isotopomer M_{i}(+1), while the isotopomer constraints corresponding fourth and fifth rows of S_{E }⊗ S_{F }contribute to mass isotopomer M_{i}(+2). (From the definition of the tensor product we see that, for example, the second row of S_{E }⊗ S_{F }corresponds to isotopomer constraints P_{E}(00)·P_{F}(1) = P_{E}(+0)·P_{F}(+1) and the third row corresponds to the constraints P_{E}(01)·P_{F}(0) + P_{E}(10)·P_{F}(0) = P_{E}(+1)·P_{F}(+0), thus validating our intuitive observation.) When the similiar information for all rows of S_{E }⊗ S_{F }is collected to a linear equation system, we will obtain the following constraints to the mass isotopomer distribution (which in the case of onecarbonfragment MF coincides with the isotopomer distribution D_{MF}):
which is equal to (8).
Calculability analysis
Isotopomer tracer experiments using less common carbon sources can be very costly because of the prices of purposefully labelled substrates. Thus, it is very useful to be able to first conduct in silico calculability analysis to find out, whether it is even in principle possible to obtain the required flux information from the tracer experiment. By analyzing the fragment equivalence classes, it is relatively easy to perform this kind of "structural identifiability analysis" (cf. [39,40] for global isotopomer balancing), that is, to discover the set of junction metabolites for which the flux ratios can be calculated (in the best case) from the given measurement data: it is enough to check what type of isotopomer constraints
can be propagated to each subpool M_{ij }of junction metabolites M_{i }from the measured metabolites (we need to know only coefficients s_{bij}, not the isotopomer abundances d_{ij}). Then, by applying the techniques of computing vector subspace intersection described above, we can compute the maximal number of linearly independent constraints obtainable for the flux ratios of each junction. Thus, it is possible to check before costly and timeconsuming wet lab experiments, whether the experiments even have potential to answer the biological questions at hand. The results of the calculability analysis tell which flux ratios are in principle determinable, given the labelling of external substrates, topology of the metabolic network and the available measurement data. It then depends on the actual flux distribution and the accuracy of the measurements, whether these ratios can be reliably determined from the experimental data.
Estimating the flux distribution of the metabolic network
In the main step of our framework for ^{13}C metabolic flux analysis, the fluxes of the metabolic network are estimated by forming and solving generalized isotopomer balance equations (4). The generalized isotopomer balance equations are based on the isotopomer measurement data that is first propagated in the network to unmeasured metabolites by utilizing the results of the structural analysis presented above.
Measurement propagation
The aim of the propagation of measurement data is to infer from the isotopomer constraints of measured metabolites as many isotopomer constraints as possible to unmeasured metabolites. As a rule of thumb, more constraints the unmeasured metabolites will get more generalized isotopomer balance equations (4) bounding the fluxes can be written.
Fragment equivalence classes can be utilized in the measurement propagation: from isotopomer constraints known for fragment M_{i}F isotopomer constraints for other fragments M_{l}F_{k }in the equivalence class of F can be easily computed. The process is the following:
1. Before measurements are propagated from fragment MF of measured metabolite M to other fragments in the equivalence class of F, isotopomer constraints to F are computed from the constraints measured to the whole metabolite M by using the vector space projection techniques (see Section Projection of isotopomer measurements to fragments).
2. The fragment constraints are propagated to all fragments F' that have been found equivalent to F via the flow analysis technique.
This requires mapping of isotopomers of F to isotopomers of F' by applying the carbon mappings of the reactions along any pathway between F and F'.
3. After the propagation of measurement data inside the fragment equivalence classes, new isotopomer constraints for independent fragments of the same metabolite can be derived, as described in Section Independence analysis of fragments.
Steps 2 and 3 can be iterated until no new isotopomer constraints to the fragments are discovered.
Construction of generalized isotopomer balances
After the propagation step, we have some isotopomer constraints for each subpool j of every junction metabolite M_{i}. (For nonjunction metabolites, isotopomer balance equations do not contain any additional flux information compared to the mass balances.) In the best case we know complete isotopomer distribution , in the worst case we have only trivial constraints stating that the sum of relative abundances of all isotopomers equals one.
Next, a linear equation system containing flux constraints obtained from mass balances (1) and generalized isotopomer balances (4) is constructed.
However, the isotopomer constraints of different subpools do not yet conform to (4) as the matrices S_{ij }are not necessarily the same.
Thus we still need to compute a common subspace ( is spanned by the rows of S_{ij}) of the isotopomer constraints known for each subpool M_{ij }and project subpool constraints to .
This can be done with the same techniques that were previously applied to project measured isotopomer information of a metabolite to its fragments. Let Y_{i }be a matrix with row space . After the projection we obtain isotopomer constraints for each subpool M_{ij }(See Figure 8 for an example).
Figure 8. An example of the computation of the common subspace of isotopomer constraints in different subpools. The mass isotopomer distribution of junction metabolite M_{1 }is assumed to be measured. For the in flow subpools M_{11 }and M_{12 }we obtain isotopomer constraints from the above reactant metabolites by measurement propagation. These propagated constraints must be projected to mass isotopomer to the subspace defined by the mass isotopomer distribution of M_{1 }before generalized isotopomer balances are constructed.
Now the isotopomer constraints of all the subpools lie in the same subspace of and we are ready to write the system of generalized isotopomer balance equations (4) for every junction M_{i}:
that is,
where g_{i }= β_{i}z_{i0}.
Estimating the fluxes
The ratios of (forward) fluxes producing M_{i }can be computed by solving the corresponding Equation (11) augmented with a constraint that fixes the out flow from M_{i }to equal 1. Thus, we obtain flux ratios of junction metabolites without manual derivation of ratio equations, without nonlinear optimization and without knowing intake and outtake rates of external metabolites or biomass composition.
In addition, when the equations (11) of all junction metabolites are combined with the mass balances (1) of nonjunctions, we obtain a linear equation system
constraining the fluxes v of the network that contains a block (junctions) or a row (nonjunctions) A_{k }for each metabolite M_{k}. Measured external fluxes and other known constraints, such as the composition of biomass can also be added to (12) as additional constraints. Additional constraints, like ones derived from gene regulatory information [41] or from thermodynamic analysis of metabolism [4244] can also easily be included to (12).
If (12) is of full rank, the whole flux distribution can be solved with standard linear algebra [45]. Also, more complex, nonlinear methods can be applied to model the effect of experimental errors to the estimated flux distribution [20]. In a common case where the system is of less than full rank, a single flux distribution can not be pinpointed without additional constraints. Instead, (12) defines the space of feasible flux distributions, that are all equally good solutions. In that case we can apply techniques developed for the analysis of stoichiometric matrices to determine as many fluxes as possible [46] from (12). More generally, by linear programming we can obtain maximum (resp. minimum) values for each flux v_{i}:
where and are predetermined minimum and maximum allowable values for v_{i}
Furthermore, it is possible to search for in some sense optimal flux distribution – for example a flux distribution maximizing the production of biomass – from the feasible space defined by (12) by linear programming techniques of flux balance analysis [1,3,47,48]. In that case, isotopomer data constrain the feasible space more than the stoichiometric information would alone do, thus possibly allowing more accurate estimations of the real flux distribution.
Statistical analysis
For an experimentalist, it is important to know how sensitive the obtained estimation of fluxes is to measurement errors. If enough repeated measurements are not available to assess this sensitivity, it has to be estimated by computational techniques. In the global isotopomer balancing framework for ^{13}C metabolic flux analysis, many mathematically or computationally involved methods have been developed to analyze the sensitivity of estimated flux distributions to errors in isotopomer measurements and the sensitivity of the objective function to the changes in the generated candidate flux distributions [4953].
As our direct method for ^{13}C metabolic flux analysis is computationally efficient, we can afford to a simple, yet powerful Monte Carlo procedure to obtain estimates on the variability of individual fluxes due to measurement errors:
1. For each measured metabolite M_{i}: By studying the variability in the repeated measurements, fix the distribution Ω_{i }from which the measurements of M_{i }are sampled.
2. Repeat k times:
(a) For each measured metabolite M_{i}: sample a measurement from Ω_{i}.
(b) Estimate fluxes v_{l }from the sampled measurements.
3. Compute appropriate statistics from the set V = {v_{1}, ..., v_{k}} to describe the sensitivity of fluxes to measurement errors.
Possible statistics that can be applied in the last step of the above algorithm include standard deviation, empirical confidence intervals [53], kurtosis, standard error etc. of each individual flux v_{j }and measures of "compactness" of V, such as (normalized) average distance of items of V from the sample average.
Experimental NMR and GCMS methods
In this section we shortly describe the experimental procedures applied in NMR and GCMS isotopomer measurements that produced the data for Section.
In the first experiment S. cerevisiae was grown in an aerobic glucoselimited chemostat culture at dilution rate 0.1 h^{1}. After reaching a metabolic steady state, as determined by constant physiological parameters 10% of the carbon source in the medium was replaced with fully carbon labelled glucose ([U13C]) for approximately 1.5 residence times, so that about 78% of the biomass was uniformly labelled. 2D [^{13}C, ^{1}H] COSY spectra of harvested and hydrolysed biomass were acquired for both aliphatic and aromatic resonances at 40°C on a Varian Inova 600 MHz NMR spectrometer. The software FCAL v.2.3.0 [19] was used to compute isotopomer constraints for 15 amino acids from the spectra. Detailed description of the cultivation set up can be found in [54] whereas similar ^{13}C labeling set up, NMR experiments and spectral data analysis as were applied here have been described in [55].
In the second experiment B. subtilis was grown on shake flasks containing 50 ml M9 minimal medium. In the experiment, the medium was supplemented with 50 mg/L tryptophan and 3 g/L glucose labelled to the first carbon position ([113C]) (99%; Cambridge Isotope Laboratories) or a mixture of 0.6 g/L fully carbon labelled glucose ([U13C]) (99%; Cambridge Isotope Laboratories) and 2.4 g/L unlabeled glucose as the sole carbon source. Fourteen derivatized amino acids were analyzed for ^{13}C labeling patterns with a series 8000GC combined with an MD800 mass spectrometer (Fisons instruments). More information about the details of the measurement procedure can be found from [20] where identical measurement techniques were applied.
Authors' contributions
AR, JR and EU are the main contributors of the computational methods given in the article. AR implemented the methods and ran the experiments. PJ constructed the models of metabolic networks and contributed to the development of the computational methods. NZ and HM contributed the development of the computational methods. AR, JR, NZ and PJ wrote the manuscript. All authors have read and accepted the contents of the manuscript.
Acknowledgements
This research has been supported by Academy of Finland (SYSBIO programme, grant 207436, SYSFYS project). We thank Uwe Sauer and Merja Penttilä for providing the data and for their comments to the work presented in the article; Markus Heinonen, Esa Pitkänen and Arto Åkerlund for the implementation of software tools supporting the presented framework; Roelco Kleijn, SarahMaria Fendt, Simon Tännler and Martin Rühl for providing the data and for fruitful discussions. We also thank the anonymous reviewers whose comments helped us to improve the manuscript substantially.
References

Varma A, Palsson B: Metabolic flux balancing: basic concepts, scientific and practical use.
Nature Biotechnology 1994, 12(10):994998. Publisher Full Text

Edwards JS, Ibarra RU, Palsson B: In silico predictions of Escherichia coli metabolic capabilities are consistent with experimental data.
Nature Biotechnology 2001, 19(2):125130. PubMed Abstract  Publisher Full Text

Schütz R, Kuepfer L, Sauer U: Systematic evaluation of objective functions for predicting intracellular fluxes in E. coli.

Szyperski T: Biosynthetically directed fractional ^{13}Clabelling of proteinogenic amino acids.
European Journal of Biochemistry 1995, 232(2):433448. PubMed Abstract  Publisher Full Text

Dauner M, Sauer U: GCMS analysis of amino acids rapidly provides rich information for isotopomer balancing.
Biotechnology Progress 2000, 16(4):642649. PubMed Abstract  Publisher Full Text

Rousu J, Rantanen A, Ketola R, Kokkonen J: Isotopomer distribution computation from tandem mass spectrometric data with overlapping fragment spectra.

Schmidt K, Carlsen M, Nielsen J, Viladsen J: Modeling Isotopomer Distributions in Biochemical Networks Using Isotopomer Mapping Matrices.
Biotechnology and Bioengineering 1997, 55(6):831840. Publisher Full Text

Wiechert W, Petersen MöllneyS, de Graaf A: A Universal Framework for ^{13}C Metabolic Flux Analysis.
Metabolic Engineering 2001, 3(3):265283. PubMed Abstract  Publisher Full Text

Wiechert W, de Graaf A: Bidirectional Reaction Steps in Metabolic Networks: I. Modeling and simulation of carbon isotope labeling experiments.
Biotechnology and Bioengineering 1997, 55(1):101117. Publisher Full Text

Antoniewicz M, Kelleher J, Stephanopoulos G: Elementary Metabolite Units (EMU): a novel framework for modeling isotopic distributions.
Metabolic Engineering 2007, 9(1):6886. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

van Winden W, Heijnen J, Verheijen P: Cumulative bondomers: a new concept in flux analysis from 2D [^{13}C, ^{1}H] COSY NMR data.
Biotechnology and Bioengineering 2002, 80(7):731745. Publisher Full Text

Nöh K, Wahl A, Wiechert W: Computational tools for isotopically instationary ^{13}C labeling experiments under metabolic steady state conditions.
Metabolic Engineering 2006, 8(6):554577. PubMed Abstract  Publisher Full Text

Nöh K, Grönke K, Luo B, Takors R, Oldiges M, Wiechert W: Metabolic flux analysis at ultra short time scale: Isotopically nonstationary ^{13}C labeling experiments.
Journal of Biotechnology 2007, 129(2):249267. PubMed Abstract  Publisher Full Text

Young J, Walther J, M A, Yoo H, Stephanopoulos G: An elementary metabolite unit (EMU) based method of isotopically nonstationary flux analysis.
Biotechnology and Bioengineering 2007, 99(3):686699. Publisher Full Text

Ghosh S, Zhu T, Grossmann I, Ataai M, Domach M: Closing the loop between feasible flux scenario identification for construct evaluation and resolution of realized fluxes via NMR.
Computers and Chemical Engineering 2005, 29(3):459466. Publisher Full Text

Blank L, Kuepfer L, Sauer U: Largescale ^{13}Cflux analysis reveals mechanistic principles metabolic network robustness to null mutations in yeast.
Genome Biology 2005, 6(6):R49. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Blank L, Lehmbeck F, Sauer U: Metabolicflux and network analysis in fourteen hemiascomycetous yeasts.
FEMS Yeast Research 2005, 5(6–7):545558. PubMed Abstract  Publisher Full Text

Fischer E, Sauer U: Metabolic flux profiling of Escherichia coli mutants in central carbon metabolism using GCMS.
European Journal of Biochemistry 2003, 270(5):880891. PubMed Abstract  Publisher Full Text

Szyperski T, Glaser R, Hochuli M, Fiaux J, Sauer U, Bailey J, Wütrich K: Bioreaction Network Topology and Metabolic Flux Ratio Analysis by Biosynthetic Fractional ^{13}C Labeling and TwoDimensional NMR Spectrometry.
Metabolic Engineering 1999, 1(3):189197. PubMed Abstract  Publisher Full Text

Fischer E, Zamboni N, Sauer U: Highthroughput metabolic flux analysis based on gas chromatographymass spectrometry derived ^{13}C constraints.
Analytical Biochemistry 2004, 325(2):308316. PubMed Abstract  Publisher Full Text

Rousu J, Rantanen A, Maaheimo H, Pitkänen E, Saarela K, Ukkonen E: A method for estimating metabolic fluxes from incomplete isotopomer information.
Computational Methods in Systems Biology, Proceedings of the First International Workshop, CMSB Volume 2602 of Lecture Notes in Computer Science 2003, 88103.

Sauer U, Hatzimanikatis V, Bailey J, Hochuli M, Szyperski T, Wüthrich K: Metabolic fluxes in riboflavinproducing Bacillus subtilis.
Nature Biotechnology 1997, 15(5):448452. PubMed Abstract  Publisher Full Text

Maaheimo H, Fiaux J, Cakar Z, Bailey J, Sauer U, Szyperski T: Central carbon metabolism of Saccharomyces cerevisiae explored by biosynthetic fractional ^{13}C labeling of common amino acids.
European Journal of Biochemistry 2001, 268(8):24642479. PubMed Abstract  Publisher Full Text

Zamboni N, Sauer U: Knockout of the highcoupling cytochrome aa3 oxidase reduces TCA cycle fluxes in Bacillus subtilis.
FEMS Microbiology Letters 2003, 226(1):121126. PubMed Abstract  Publisher Full Text

Sola A, Maaheimo H, Ylönen K, Ferrer P, Szyperski T: Amino acid biosynthesis and metabolic flux profiling of Pichia pastoris.
European Journal of Biochemistry 2004, 271(12):24622470. PubMed Abstract  Publisher Full Text

Sola A, Jouhten P, Maaheimo H, SanchezFerrando F, Szyperski T, Ferrer P: Metabolic flux profiling of Pichia pastoris grown on glycerol/methanol mixtures in chemostat cultures at low and high dilution rates.
Microbiology 2007, 153(1):281290. PubMed Abstract  Publisher Full Text

Rantanen A, Maaheimo H, Pitkänen E, Rousu J, Ukkonen E: Equivalence of metabolite fragments and flow analysis of isotopomer distributions for flux estimation.
Transactions on Computational Systems Biology VI, Lecture Notes in Bioinformatics 2006, 4220:198220.

Kleijn R, Geertman J, Nfor B, Ras C, Schipper D, Pronk J, Heijnen J, van Maris A, van Winden W: Metabolic flux analysis of a glyceroloverproducing Saccharomyces cerevisiae strain based on GCMS, LCMS and NMR derived ^{13}Clabeling data.
FEMS Yeast Research 2006, 7(2):216231. PubMed Abstract  Publisher Full Text

Arita M: In Silico Atomic Tracing of SubstrateProduct Relationships in Escherichia coli Intermediary Metabolism.
Genome Research 2003, 13(11):24552466. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Zamboni N, Fischer E, Sauer U: FiatFlux – a software for metabolic flux analysis from ^{13}Cglucose experiments.
BMC Bioinformatics 2005., 6(209) PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Szyperski T: ^{13 }CNMR, MS and metabolic flux balancing in biotechnology research.
Quarterly Reviews of Biophysics 1998, 31(1):41106. PubMed Abstract  Publisher Full Text

Weitzel M, Wiechert W, Nöh K: The topology of metabolic isotope labeling networks.
BMC Bioinformatics 2007., 8(315) PubMed Abstract  Publisher Full Text  PubMed Central Full Text

van Winden W, van Dam J, Ras C, Kleijn R, Vinke J, Gulik W, Heijnen J: Metabolicflux analysis of Saccharomyces cerevisiae CEN.PK1137D based on mass isotopomer measurements of ^{13}Clabeled primary metabolites.
FEMS Yeast Research 2005, 5(6–7):559568. PubMed Abstract  Publisher Full Text

Rantanen A, Mielikäinen T, Rousu J, Maaheimo H, Ukkonen E: Planning optimal measurements of isotopomer distributions for estimation of metabolic fluxes.
Bioinformatics 2006, 22(10):11981206. PubMed Abstract  Publisher Full Text

Schmidt K, Nielsen J, Villadsen J: Quantitative analysis of metabolic fluxes in Escherichia coli, using twodimensional NMR spectroscopy and complete isotopomer models.
Journal of Biotechnology 1999, 71(1–3):175189. PubMed Abstract  Publisher Full Text

Rantanen A, Rousu J, Ketola R, Kokkonen J, Tarkiainen V: Computing Positional Isotopomer Distributions from Tandem Mass Spectrometric Data.
Metabolic Engineering 2002, 4:285294. PubMed Abstract  Publisher Full Text

Aho A, Hopcroft J, Ullman J: The Design and Analysis of Computer Algorithms. Addison Wesley; 1974.

Appel A: Modern Compiler Implementation in Java. Cambridge University Press; 1998.

Isermann N, Wiechert W: Metabolic isotopomer labeling systems. Part II: structural identifibiality analysis.
Mathematical Biosciences 2003, 183(2):175214. PubMed Abstract  Publisher Full Text

van Winden W, Heijnen J, Verheijen P, Grievink J: A priori analysis of metabolic flux identifiability from ^{13}Clabeling data.
Biotechnology and Bioengineering 2001, 74(6):505516. Publisher Full Text

Covert M, Schilling C, Palsson B: Regulation of Gene Expression in Flux Balance Models of Metabolism.
Journal of Theoretical Biology 2001, 213(1):7388. PubMed Abstract  Publisher Full Text

Kümmel A, Panke S, Heinemann M: Systematic assignment of thermodynamic constraints in metabolic network models.
BMC Bioinformatics 2006., 7(512) PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Henry C, Broadbelt L, Hatzimanikatis V: ThermodynamicsBased Metabolic Flux Analysis.
Biophysical Journal 2007, 92(5):17921805. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hoppe A, Hoffmann S, Holzhütter H: Including metabolite concentrations into flux balance analysis: thermodynamic realizability as a constraint on flux distributions in metabolic networks.
BMC Systems Biology 2007., 1(23) PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Schwarz H: Numerical Analysis: A Comprehensive Introduction. John Wiley & Sons; 1989.

Klamt S, Schuster S: Calculating as many fluxes as possible in underdetermined metabolic networks.
Molecular Biology Reports 2002, 29(1–2):243248. PubMed Abstract  Publisher Full Text

Bonarius H, Schmidt G, Tramper J: Flux analysis of underdetermined metabolic networks: the quest for the missing constraints.
Trends in Biotechnology 1997, 15(8):308314. Publisher Full Text

Edwards J, Covert M, Palsson B: Metabolic modelling of microbes: the fluxbalance approach.
Environmental Microbiology 2002, 4(3):133140. PubMed Abstract  Publisher Full Text

Möllney M, Wiechert W, Kownatzki D, de Graaf A: Bidirectional Reaction Steps in Metabolic Networks IV: Optimal Design of Isotopomer Labeling Experiments.
Biotechnology and Bioengineering 1999, 66(2):86103. Publisher Full Text

Kleijn R, van Winden W, Ras C, van Gulik W, Schipper D, Heijnen J: 13CLabeled Gluconate Tracing as a Direct and Accurate Method for Determining the Pentose Phosphate Pathway Split Ratio in Penicillium chrysogenum.
Applied and Environmental Microbiology 2006, 72(7):47434754. Publisher Full Text

AraúzoBravo M, Shimizu K: An improved method for statistical analysis of metabolic flux analysis using isotopomer mapping matrices with analytical expressions.
Journal of Biotechnology 2003, 105(1–2):117133. PubMed Abstract  Publisher Full Text

Wiechert W, Siefke C, de Graaf A, Marx A: Bidirectional Reaction Steps in Metabolic Networks: II. Flux Estimation and Statistical Analysis.
Biotechnology and Bioengineering 1997, 55(1):118134. Publisher Full Text

Antoniewicz M, Kelleher J, Stephanopoulos G: Determination of confidence intervals of metabolic fluxes estimated from stable isotopome measurements.
Metabolic Engineering 2006, 8(4):324337. PubMed Abstract  Publisher Full Text

Wiebe M, Rintala E, Tamminen A, Simolin H, Salusjärvi L, Toivari M, Kokkonen J, Kiuru J, Ketola R, Jouhten P, Huuskonen A, Maaheimo H, Ruohonen L, Penttilä M: Central carbon metabolism of Saccharomyces cerevisiae in anaerobic, oxygenlimited and fully aerobic steadystate conditions and following a shift to anaerobic conditions.
FEMS Yeast Research 2008, 8(1):140154. PubMed Abstract  Publisher Full Text

Fiaux J, Cakar P, Sonderegger M, Wüthrich K, Szyperski T, Sauer U: MetabolicFlux Profiling of the Yeasts Saccharomyces cerevisiae and Pichia stipitis.
Eukaryotic Cell 2003, 2(1):170180. PubMed Abstract  Publisher Full Text  PubMed Central Full Text