Expanded Isotopomer Model

An isotopomer model was constructed in two phases. First, a central metabolic isotopomer model that accounts for 85 reactions including glycolysis, the TCA cycle, the pentose phosphate pathway, oxidative phosphorylation, pyruvate metabolism, and anaplerotic reactions was derived from the iJR904

*E. coli *reconstruction [

16]. This initial model was equivalent in reaction content to commonly used isotopomer models for

*E. coli *[

25,

26].

An expanded model was then constructed that includes both central and biosynthetic pathways. The iMC1010 metabolic network [

19] was evaluated to determine which reactions can sustain non-zero fluxes during growth on glucose, acetate, or lactate when only certain by-products are allowed to be secreted (acetate, formate, D-lactate, pyruvate, succinate, glycerol, CO

_{2}, and ethanol). Blocked reactions, which must have zero net flux at steady state, were subsequently omitted from consideration. Groups of reactions that could be merged together without affecting model results (e.g. linear pathways) were combined in order to reduce the number of variables. Large sets of biosynthesis reactions that produce phospholipids, nucleotides, co-factors were also combined, since there no experimental measurements existed for these high-carbon metabolites. However, by-products resulting from high-carbon metabolite production (e.g. CO

_{2}, formate, succinate, fumarate, and pyruvate) that could enter back into the metabolic network were tracked. Of the original 932 reactions in the complete metabolic iMC1010 network, nearly a third were represented in the biosynthetic isotopomer model, either individually or as grouped reactions.

The final isotopomer model accounts for a total of 313 irreversible reactions, including 278 which track carbon. Inclusion of these additional pathways is likely important for accurate assessment of the flux-resolving power of

^{13}C experiments both within and beyond central metabolism [

7]. A complete listing of the reactions and metabolites in the biosynthetic network can be found in the Additional File

1.

Monte Carlo Sampling Approach

To compute possible flux distributions of the

*E. coli *model, the network was sampled using a Markov Chain, Monte Carlo (MCMC) sampling algorithm (see Methods). The steady-state mass balance and uptake rate constraints for the metabolic network create a convex hyperspace that contains all biochemically feasible steady-state flux distributions [

27]. Monte Carlo sampling generates a set of flux distributions that are spread uniformly throughout the feasible space. The inclusion of

^{13}C experimental data reduces the feasible space in which the true flux state must lie by requiring that the IDV calculated from the putative flux distribution must match the experimental data within error. While the space of feasible flux distributions depends only on reaction stoichiometry, the space of resulting simulated IDVs differs depending on the input substrate labeling pattern. Hence, different labeling patterns can have differing abilities to resolve each reaction flux.

Here, we used Monte Carlo sampling of flux distributions to analyze the degree to which reaction fluxes can be determined by steady-state ^{13}C labeling experiments in terms of several possible experimental objectives. For example, one possible experimental objective is to determine whether a particular reaction has a flux above or below a specified value. For this objective, a well-designed labeling pattern would be one in which flux distributions that have an objective reaction flux greater than the specified value can be easily distinguished from flux distributions with an objective reaction flux less than the specified value. As seen in Figure , a hypothetical experiment 1 produces measurement distributions which overlap whereas experiment 2 shows greater separation. If one were interested in differentiating between the two partitions, experiment 2 would be much preferable. This method allows for the scoring of any label for any given experimental objective without first knowing the true cellular flux distribution *v*.

**Generating and Evaluating **^{13}**C Experimental Hypotheses**

An experimental hypothesis is defined as a partition of the sampled flux distribution set. While many possible hypotheses could be considered, two rational hypotheses were studied. The first case attempts to elucidate whether a reaction has high or low flux (hi-lo). The solution space is partitioned into all points with *v*_{j }>threshold versus *v*_{j }<threshold. A different hypothesis is generated for each reaction *j*. The threshold was chosen to be the median of all *v*_{j }so that half of all points would be in each of the two partitions. The second set of hypotheses tested consisted of biologically relevant flux ratios. For each point the ratio of two reactions, *v*_{i}/*v*_{j}, was determined to be above or below some threshold that formed a partition.

Intuitively, a hypothesis score should be high if the isotopomer distributions coming from one partition are distinguishable from distributions in the other partition. While there are several ways of doing this, we chose a heuristic metric based on a Z-score, which is commonly used to determine the difference between two samples. A Z-score was calculated for each fragment (element) of the calculated MDV for each simulated flux distribution:

where

and

are the average fragment enrichments for the upper and lower partitions, respectively,

and

are the variances of fragment enrichments for the upper and lower partitions, respectively, and the

*α *is a constant equal to 0.014.

*α *is on the order of magnitude of the uncertainty in measurements. This slight modification to the standard Z-score puts a lower bound on the expected experimental variation. The Z-score of each fragment is added together to give the Z-score of the experiment.

Using this approach, candidate flux states were sampled uniformly and experimental hypotheses tested. Z-scores were calculated for the hi-lo hypothesis corresponding to 1) individual reactions 2) reaction ratios and 3) two 'random' reactions or ratios. Random hypotheses were tested to estimate the level of noise associated with the set of flux distributions. Raw and normalized Z-scores are given in the Additional File

1. Z-scores varied from the level of noise to a maximum of

*>*20-fold the level of noise.

To illustrate the differences in label-dependent reaction resolving capacity, two sets of Z-scores corresponding to [1-^{13}C] glucose and [6-^{13}C] glucose are plotted in Figure . Lighter colors indicate higher Z-scores and ease of measurement. In this case, [6-^{13}C] glucose scores higher at measuring the pentose phosphate pathway and most of lower glycolysis, whereas [1-^{13}C] glucose glucose scores much higher at measuring the glyoxylate shunt. The results suggest that there is no single label that yields a high score for all experimental objectives. For example, the exchange of formate (EX_for) could be easiest measured with a [1,2-^{13}C] glucose; however, this labeling pattern is bested by [1-^{13}C] glucose for the measurement of reaction formyltetrahydrofolate deformylase (FTHFD) (Figure ). This non-universality of labels is in line with expectations, as it has been previously shown that the choice of labels can affect the flux resolution. For many reactions, the best experiment that could be performed involves hypothetical (non-commercially available) labels. One example is the ratio of phosphofructokinase (PFK) flux to fructose bisphosphate aldolase (FBA) flux. The best label for determining this ratio is [1,2,3-^{13}C] glucose (Z = 28.0), which gives a much higher Z-score than the best commercially available label, [1,2-^{13}C] glucose (Z = 18.5). Thus, there may be motivation to synthesize compounds with more complex labeling patterns than commonly used. Additionally, there are certain reactions which are predicted to be difficult to measure with any labeling pattern. For example, the Z-scores for each possible labeled glucose substrate for the reaction pyruvate oxidase (POX) all lie within the level of noise, as determined by the comparison with random hi-lo experiment Z-scores.

In addition to label-specific reaction flux elucidation properties, ^{13}C experiments show a clear pathway bias regardless of labeling pattern. The maximum Z-score of all labeling patterns was found for each reaction, giving a metric for the maximum potential for reaction flux determination using ^{13}C-labeled glucose (Figure ). Then, the fraction of reactions that had a maximum potential at least twice the noise level was found for each subsystem (Figure ). Reactions that were stoichiometrically fixed by the measured constraints on acetate, glucose, D-lactate, oxygen, growth rate, and ATP maintenance were also categorized by subsystem. Stoichiometrically fixed (i.e. constraint-determined) reactions have a confidence interval of zero, and thus are label-independent and receive no additional knowledge from ^{13}C experiments. It was found that histidine, valine, leucine, and isoleucine metabolism fluxes are completely identified solely based on the flux constraints. On the other hand, prior constraints fix none of the fluxes in central carbon metabolic systems such as glycolysis, citric acid cycle, pentose phosphate pathway, and anaplerotic reactions; however, fluxes in these pathways are all predicted to be identifiable with a ^{13}C experiment using optimal labeling patterns for each reaction. This result is expected as these identifiable pathways are the typical pathways being studied using ^{13}C analysis. Other pathways, such as cysteine, threonine, and lysine metabolism, are completely identifiable through a combination of prior stoichiometric constraints combined with well-chosen ^{13}C experiments. However, many of the remaining subsystems have a fraction of reactions that cannot be determined using any ^{13}C-labeling pattern of glucose. In particular, no additional information can be obtained from ^{13}C-labeled glucose experiments about certain biosynthetic pathways, nucleotide salvage pathways, reductive citric acid cycle reactions, and certain alternate pathways peripheral to glycolysis, such as an alternate pathways from DHAP to D-lactate. Measuring metabolites other than amino acids may give more information on these pathways. Note that in this discussion of identifiability, the Z-score metric indicates that an experiment can significantly reduce the confidence interval of a particular reaction but does not specifically predict the value of the confidence interval. Confidence intervals are directly calculated for experimental data sets in a later section and compared to the Z-scores for the same labeling patterns.

Dimensionality of Isotopomer Data

The Monte Carlo sampling approach enables the determination of the dimensionality of simulated ^{13}C experiments for a particular substrate labeling pattern. The dimensionality gives an indication of the degree to which a particular substrate labeling pattern can specify the free dimensions inherent in a network structure, given a set experimental error. In an extreme case, if all the data falls on one point (zero dimensions), no additional information is given from the data. Similarly, a dimensionality of one indicates that the data can specify one degree of freedom. Singular value decomposition (SVD) is a data reduction technique that allows the estimation of data dimensionality (Figure ). A data matrix *M *of size (*n*_{fragments }x *n*_{points}), consisting of all sample points generated from Monte Carlo Sampling, is decomposed into *M *= *U *· Σ · *V *^{T }where *U *and *V *are orthonormal bases and Σ is a diagonal matrix containing singular values in descending order. The singular values are effectively weightings that describe the information content of the corresponding vectors in *U *and *V *towards reconstructing the full matrix *M*. A partial reconstruction of *M *is possible by taking only a subset of the singular values greater than some threshold. These thresholds have a direct interpretation as the uncertainty with which a data point can be measured. For example, a threshold cutoff of 0.01 indicates that the remaining uncertainty of the data falls within 0.01 or 1% error in the measurement of isotope enrichment.

To determine the dimensionality of the isotopomer data, SVD was performed on ^{13}C fragments derived from uniformly sampled flux distribution sets for several glucose labels. The results are summarized in Figure . Globally, the choice of glucose labels affects the dimensionality of the resulting isotopomer data set. At the 1% (0.01) threshold, the label with the highest dimensionality was hypothetical [1,2,5-^{13}C] glucose with 73 dimensions. The three labels for which experimental data was measured in the subsequent section, [1-^{13}C] glucose, [6-^{13}C] glucose, and 20% [U-^{13}C] glucose, had dimensions 53, 53 and 35, respectively, at this cutoff. These values are all significantly lower than the best label, and, in particular, the uniform labeled experiment only produces half of the dimensionality as the optimal experiment. This result is significant. While 139 dimensions (the number of undetermined dimensions for the model used) are required to specify a unique flux vector, the dimensionality of the ^{13}C data for each label is significantly lower. The best labeling experiment specifies just over half (73/139 = 0.52) the degrees of freedom required, and 20% [U-^{13}C] glucose only specifies about one fourth of the possible degrees of freedom (0.26). It is worth noting that SVD is a linear operation used to approximate properties of a non-linear system and the true degrees of freedom may be even lower than reported. SVD serves as a useful upper bound on the dimensionality of data for non-linear systems, but the difference between SVD dimensionality and true dimensionality may grow to be unacceptable for large systems. For the system studied here, SVD was found to be of practical use.

Experimental Validation

In order to assess the agreement of computationally predicted flux elucidation capacity with experimental data, we took fluxomic measurements for three labeling patterns in *E. coli*. Flux distributions that best explain each set of ^{13}C data were calculated using a non-linear optimization problem:

The function Error(v) is a score of how well a given flux distribution fits the experimental data. It is defined as:

where measured _{i }is the measured fractional enrichment of fragment *i*, fragment _{i}(*v*) is the computed fractional enrichment of fragment *i *as a function of the flux distribution *v*, and *α *= 0*:*014 is the standard deviation of the fragments as calculated from experimental replicates.

Reaction flux confidence ranges were then computed for all reactions using all three sets of ^{13}C data and all combinations thereof. Confidence intervals for reaction rates were computed by maximizing and minimizing the value of each reaction in turn subject to a slightly relaxed score.

where *c*_{i }= (0, 0, ...0, 1, 0...0) is a vector of all zeros with a 1 in position i, and *Error*_{max }was set based on the confidence value. Because different data sets provide different levels of consistency, Error_{max }was chosen to be 30 units greater than the minimum error found.

These intervals were compared with Z-scores calculated through Monte Carlo methods to assess the ability of the Z-scores to predict the size of experimental reaction ranges in a label-specific manner (Figure ). The Z-scores were found to be correlated with the relative flux ranges in a statistically significant manner (Student's T-test, p *<*8.6 × 10^{-34}). A receiver operating characteristic (ROC) curve suggests that the Z-scores can identify with both sensitivity and specificity the reactions that can be elucidated in a label-specific manner, with better performance predicting ranges that are restricted more by data (Figure ). These findings indicate that the Z-score is indeed a useful predictor of the degree a flux range will be constrained by a particular ^{13}C experiment and provide experimental support for the computational approach taken.

The number of reactions elucidated at particular confidence intervals was then found (Figure ). Using different labels provides different levels of reaction confidence (Additional File

1). Including no

^{13}C data generates the largest flux ranges (lower black line), while adding

^{13}C data reduces the ranges and shifts the curve left. With almost no exception, including one experiment yields larger confidence intervals than any combination of two carbon sources which in turn is a larger range than including all three sets. Of the single experiment curves, the 20% [U-

^{13}C] glucose curve provides notably worse ranges than the other two experiments, consistent with the finding that 20% [U-

^{13}C] substrate provides data with the smallest number of dimensions.

At a reaction confidence of 1 mmol · gDW^{-1}·h^{-1 }(a relatively non-stringent cutoff), 85 reactions are specified simply from uptake rate data without any ^{13}C data. Performing the least informative ^{13}C experiment, using 20% [U-^{13}C] substrate, yields 105 reactions that meet the confidence criterion, whereas the combination of all three ^{13}C experiments yields 125 reactions that meet the criterion. In other words, performing all three experiments will increase the number of elucidated reactions by 40 reactions or about 50%. As the model used contains 278 carbon-tracking reactions and reaction groups, the increase in knowledge at 1 mmol · gDW^{-1}·h^{-1 }confidence from 85 to 125 reactions from using ^{13}C data indicates that a large gap in the knowledge remains. It seems apparent that other methods must be developed to obtain flux information at the genome scale from single experiments, as would ultimately be desirable. However, as noted in the above section, the majority of reactions that are elucidated by ^{13}C-labeled glucose experiments lie in central metabolic pathways, which tend to be both of high interest and not-specifiable by constraints alone.