Subcellular flux analysis of central metabolism in a heterotrophic Arabidopsis cell suspension using steady-state stable isotope labeling.

The presence of cytosolic and plastidic pathways of carbohydrate oxidation is a characteristic feature of plant cell metabolism. Ideally, steady-state metabolic flux analysis, an emerging tool for creating flux maps of heterotrophic plant metabolism, would capture this feature of the metabolic phenotype, but the extent to which this can be achieved is uncertain. To address this question, fluxes through the pathways of central metabolism in a heterotrophic Arabidopsis (Arabidopsis thaliana) cell suspension culture were deduced from the redistribution of label in steady-state (13)C-labeling experiments using [1-(13)C]-, [2-(13)C]-, and [U-(13)C(6)]glucose. Focusing on the pentose phosphate pathway (PPP), multiple data sets were fitted simultaneously to models in which the subcellular compartmentation of the PPP was altered. The observed redistribution of the label could be explained by any one of three models of the subcellular compartmentation of the oxidative PPP, but other biochemical evidence favored the model in which the oxidative steps of the PPP were duplicated in the cytosol and plastids, with flux through these reactions occurring largely in the cytosol. The analysis emphasizes the inherent difficulty of analyzing the PPP without predefining the extent of its compartmentation and the importance of obtaining high-quality data that report directly on specific subcellular processes. The Arabidopsis flux map also shows that the potential ATP yield of respiration in heterotrophic plant cells can greatly exceed the direct metabolic requirements for biosynthesis, highlighting the need for caution when predicting flux through metabolic networks using assumptions based on the energetics of resource utilization.

Extensive subcellular compartmentation, with unique locations for many steps and pathways, as well as the duplication of other steps and pathways in different compartments, adds greatly to the structural complexity of the plant metabolic network (Lunn, 2007;. The need to consider discrete pools of metabolites in specific compartments, and the transporters that link them, complicates the quest for a detailed, predictive understanding of the regulation of plant metabolism; as a result, it remains difficult to manipulate flows of material through the central metabolic network in a predictable way (Carrari et al., 2003;Sweetlove et al., 2008).
The emergence of steady-state metabolic flux analysis (MFA) as a practicable systems biology tool for generating flux maps of the central metabolic pathways in plants offers new opportunities for analyzing plant metabolic phenotypes (Ratcliffe and Shachar-Hill, 2006;Libourel and Shachar-Hill, 2008;Schwender, 2008;Kruger and Ratcliffe, 2009). In this approach, substrates labeled with stable isotopes are introduced into the network, and fluxes are determined by measuring the labeling of the system after it has reached an isotopic and metabolic steady state. Subcellular compartmentation of steps and pathways can be incorporated into the model that describes the redistribution of the label, and flux maps of central carbon metabolism in plant cells typically aim to distinguish between fluxes in the cytosol, mitochondria, and plastids. In principle, compartmented fluxes can be deduced from the labeling patterns of metabolites that are synthesized in a single compartment Shachar-Hill, 2005, 2006;Allen et al., 2007), and this property has been exploited with varying degrees of success in different studies.
Most recent flux maps, such as those for developing embryos of oilseed rape (Brassica napus; Schwender et al., 2006;Junker et al., 2007) and sunflower (Helianthus annuus;Alonso et al., 2007a) or the one for a heterotrophic cell suspension of Arabidopsis (Arabidopsis thaliana; Williams et al., 2008), have well-defined mitochondrial fluxes. This reflects the high sensitivity of the available label measurements to these fluxes, but it also reflects the arbitrary way in which the models assume that some fluxes, for example from citrate to 2-oxoglutarate, are exclusively mitochondrial and others, for example between malate and oxaloacetate, occur in an unspecified compartment. There is less consensus in the description of the pathways of carbohydrate oxidation, with the subcellular compartmentation of glycolysis and the oxidative pentose phosphate pathway (PPP) being described at various levels of detail. For example, a recent model of soybean (Glycine max) cotyledons included complete pathways for glycolysis and the PPP in both cytosol and plastids (Iyer et al., 2008), whereas models for developing sunflower embryos (Alonso et al., 2007a), maize (Zea mays) root tips (Alonso et al., 2007b), and an Arabidopsis cell culture (Williams et al., 2008) confined the oxidative and nonoxidative steps of the PPP to the plastid. The limited compartmentation of the pathways of carbohydrate oxidation in many flux maps can be traced to an analysis of developing oilseed rape embryos (Schwender et al., 2003), which indicated that the labeling data could be satisfactorily explained without duplication of glycolysis and the PPP. Since some of the enzymes of these pathways are known to be present in both compartments, the model implied fast exchange of intermediates across the plastid envelope, giving rise to a situation in which the physically compartmented pathways are functionally uncompartmented (Schwender et al., 2003;Ratcliffe and Shachar-Hill, 2006).
Although the extent to which the steps of the PPP are duplicated in the cytosol and plastids has not been fully established, most of the biochemical evidence suggests that the enzymes for the oxidative steps at least are present in both compartments (Kruger and von Schaewen, 2003). So modeling the PPP as fully duplicated in both cytosol and plastids, or alternatively as exclusively plastidic, ignores the most likely structure of the network. Models based on all three network structures have now been tested in an Arabidopsis cell culture. Building on a steady-state MFA of the same culture (Williams et al., 2008), similar experiments were conducted with multiple substrates and the labeling data sets were simultaneously fitted to models in which the subcellular compartmentation of the PPP was varied. Acceptable fits to the data were found for all three models, emphasizing the importance of building the most appropriate level of subcellular compartmentation into the model at the outset. Moreover, the flux maps showed that the data were consistent with a major contribution from the cytosol to the oxidative steps of the PPP.

Biosynthetic Outputs
The Arabidopsis cell culture was maintained in the dark with Glc as the sole respiratory substrate, and the principal activities of the metabolic network, together with the associated biosynthetic outputs, were assessed by determining the redistribution of radiolabel following metabolism of [U-14 C]Glc (Table I). Overall, the pattern of labeling was comparable to that observed previously for heterotrophic cell cultures of Arabidopsis and is similar to that reported for a wide range of nonphotosynthetic plant material (Kruger et al., 2007a). However, extending the analysis to include components that are frequently ignored during similar fractionation procedures, such as the organic phase obtained during chloroform/methanol fractionation and the medium following removal of cells, allowed a more comprehensive analysis of the fate of metabolized Glc and measurement of the labeling of lipids and ethanol. The uptake of label was linear over 24 h, and the proportion of radioactivity recovered in each of the major chemical fractions was the same after 12 and 24 h of incubation (Supplemental Fig. S1). This implies that the pools of metabolic intermediates had achieved isotopic steady state within this time and that the accumulation of radiolabel in the products reported in Table I reflects the rate of synthesis of these compounds rather than just isotopic equilibration within existing pools. The data show that almost 40% of the metabolized Glc was released as CO 2 , while less than 25% of the carbon was Table I. Metabolism of [U-14 C]Glc by cell suspension cultures of Arabidopsis Samples (5 mL, 0.618 6 0.017 g fresh weight) from 6-d-old cultures grown in Glc at natural abundance 13 C were supplemented with 185 kBq [U-14 C]Glc (yielding a specific activity of 308 6 15 MBq mol 21 ) and incubated for 24 h prior to extraction in methanol and subsequent metabolic fractionation. Total 14 C uptake is the sum of radioactivity in CO 2 , ethanol, and the chloroform-soluble, aqueous methanol-soluble, and methanol-insoluble fractions in each cell culture. Each value is the mean 6 SE of the distribution of radioactivity in four cell cultures. converted into biopolymers (comprising starch, cell wall, protein, and lipid), with the remainder accumulating mainly as soluble sugars, organic acids, and amino acids. While this analysis defines the major output fluxes associated with growth, further information is required to determine the fluxes that generate the precursors needed for biomass production, and this can be provided by steady-state MFA.

Metabolic Network
In MFA, the metabolic network is reduced to a model that can account for the observed redistribution of the label. Typically, this leads to a model in which (1) multiple steps between branch points in the network are represented as a single reversible or irreversible step and (2) subcellular pools of metabolites are treated as single pools in the absence of compartment-specific information about their labeling. The scope of the model is strongly influenced by the choice of label precursor and by the extent and quality of the measurements that define the redistribution of the label after the system has reached an isotopic and metabolic steady state. It has already been shown that statistically robust flux maps of central carbon metabolism can be obtained for heterotrophic Arabidopsis cell suspensions using [1-13 C]Glc as the label source (Williams et al., 2008), and here the aim was to establish whether improved definition of the subcellular compartmention of carbohydrate oxidation could be obtained by combining the results of experiments with [1-13 C]Glc, [2-13 C]Glc, and [U-13 C 6 ]Glc as label precursors.
The existing model, in which the oxidative and nonoxidative reactions of the PPP were confined to the plastids (Fig. 1A), was compared with two alternatives: (1) a model in which the oxidative steps of the PPP were allowed to occur in both the cytosol and plastids, with the nonoxidative steps confined to the plastids (Fig. 1B); and (2) a model in which the complete PPP was present in both compartments (Fig. 1C). Models were constructed in the format used by the steady-state MFA software 13C-FLUX, and full details of the model in which the oxidative steps of the PPP were duplicated in the cytosol and plastids (Figs. 1B and 2) are given in Supplemental File S1. In essence, the model defines the relationship between the carbon atoms in substrates and products for each step in the network so that the underlying fluxes can be deduced from the redistribution of label in the steady-state labeling experiment. The models used here are directly related in most respects to the one used earlier (Williams et al., 2008), and they were constructed in the same way. Several features of the models, and the biochemical networks they represent, should be noted.
First, analysis of the subcellular distribution of flux depends on obtaining information about the steadystate labeling of metabolic intermediates in specific compartments. Here, the labeling of the plastidic pools Figure 1. Alternative metabolic models for the subcellular compartmentation of the PPP. A, Oxidative and nonoxidative reactions of the PPP confined to the plastid. B, Oxidative steps of the PPP in both the cytosol and plastid, with the nonoxidative steps in the plastid. C, The complete pathway in both compartments. Net and exchange fluxes are shown with unidirectional and bidirectional arrows, respectively, and the direction of a unidirectional arrow corresponds to a positive flux in the output of the model. Flux names are given in italics, and the pentose phosphates are grouped into a single indistinguishable pool in all three models. Ery-4-P, Erythrose 4-phosphate; Rbu-5-P, ribulose 5-phosphate; Sed-7-P, sedoheptulose 7-phosphate; TA-C3, the threecarbon substituted enzyme intermediate in the reaction catalyzed by TA; TK-C2, the two-carbon substituted enzyme intermediate in the reaction catalyzed by TK; Xlu-5-P, xylulose 5-phosphate. was deduced from the labeling of starch and from amino acids that are known to be synthesized exclusively in the plastid (Cys, Gly, His, Ile, Leu, Lys, Phe, Ser, Trp, Tyr, and Val). Note that although Ser, following synthesis in the plastid (Ho and Saito, 2001), can be reversibly converted to Gly via Ser hydroxymethyltransferase, the fact that this can potentially occur in several compartments has no effect on the interpretation of the Ser labeling pattern. Similarly, the labeling of cytosolic intermediates was deduced from the labeling of Suc and Ala (Miyashita et al., 2007), and Arg, g-aminobutyrate (GABA), Glu, Gln, and Pro were assumed to be derived from mitochondrial 2-oxoglutarate. In contrast, it was not possible to deduce information about the subcellular labeling of oxaloacetate and malate, and the synthesis of these metabolites was not assigned to a specific compartment.
Second, because the model describes the redistribution of label through the network, expected biochem-ical features of the network may be concealed or omitted. For example, it is possible to synthesize Gly from Ser using Thr aldolase (Joshi et al., 2006b), but using this as the sole route to Gly leads to a substantially worse fit to the data, and adding the Thr aldolase reaction as an extra component of the model in Figure  2 leads to no change in the flux estimates (data not shown). Thus, there is no evidence from the labeling data to justify the inclusion of the Thr aldolase reaction in the model. Theoretical analysis also has a bearing on the construction of the model. For example, flux through the GABA shunt is combined with the flux through the tricarboxylic acid cycle between 2-oxoglutarate and fumarate because these pathways have identical labeling outcomes in the model. Similarly, Suc cycling is omitted because this can only be correctly analyzed in steady-state MFA with information on the labeling of the cytosolic and vacuolar Glc pools (Kruger et al., 2007b). The flux supported by the Output fluxes used to constrain the model are shown with dotted black arrows; net and exchange fluxes deduced from the 13 C label redistribution are shown with solid black unidirectional and bidirectional arrows, respectively. The direction of a unidirectional arrow corresponds to a positive flux in the output of the model, and the dotted gray boxes indicate that the enclosed metabolites were considered as a single pool. Standard abbreviations are used for the amino acids. C3-P, Three-carbon phosphate ester pool; cit, citrate; CO 2in , internal carbon dioxide; fum, fumarate; aKG, 2-oxoglutarate; OAA, oxaloacetate; PEP, phosphoenolpyruvate; TAG, triacylglycerol.
putative Glc-6-phosphatase  is omitted for the same reason, because there is no easy way to measure the labeling of the cytosolic and vacuolar Glc pools, and ignoring the compartmentation of Glc is likely to produce very misleading flux values (Kruger et al., 2007b). In contrast, discrete cytosolic and plastidic pools of Glc-6-P were retained in the model because the glucosyl units of Suc and starch differed appreciably in their patterns of labeling, establishing that they were synthesized from isotopically distinct precursors (Supplemental Fig.  S2, A-C; Supplemental Table S7). Subsequent in silico analysis established that the calculated flux distribution would produce differential labeling in products derived directly from cytosolic and plastidic hexose phosphate pools (Supplemental Fig. S2, D-F).
Third, the steps catalyzed by transketolase (TK) and transaldolase (TA) in the nonoxidative steps of the PPP are represented in terms of half-reactions to allow for the underlying mechanism of the enzymes and their lack of specificity (Kruger and von Schaewen, 2003;Kleijn et al., 2005;Selivanov et al., 2005). Conventional stoichiometric models of the nonoxidative reactions ( Fig. 3A) assume that the label is directed through a subset of the possible reactions and thus misrepresent the underlying fluxes. Adding the extra reactions for an equivalent formulation of the pathway (Fig. 3B) leads to ambiguity in the flux solution (Fig. 3D), but Figure 3. Representation of the oxidative PPP for steady-state MFA. A and B, Conventional (A) and alternative (B) stoichiometries for the nonoxidative reactions of the PPP reflecting the substrate specificity of TK. C, Formulation of the nonoxidative reactions of the PPP in terms of the two-and three-carbon substituted enzyme intermediates formed during the ping-pong mechanisms of TK and TA. D, Steady-state analysis of the fluxes through the reactions catalyzed by TK for a simulated set of 13 C-labeling data. Isotopomer abundances were first deduced for a flux distribution corresponding to the conventional network in A, and then these abundances were used to predict the flux distribution through a metabolic network that permitted flux through both the conventional (A) and alternative (B) pathways. The 100 solutions each correspond to a different proportion of the three exchange reactions catalyzed by TK: triangles, Xlu-5-P + Rib-5-P 4 Tri-3-P + Sed-7-P; circles, Xlu-5-P + Ery-4-P 4 Tri-3-P + Fru-6-P; squares, Sed-7-P + Ery-4-P 4 Rib-5-P + Fru-6-P. E, The result of the same exercise performed with the network (C) defined using half reactions for TK: circles, Xlu-5-P 4 Tri-3-P + Enz-C2; squares, Fru-6-P 4 Ery-4-P + Enz-C2; triangles, Sed-7-P 4 Rib-5-P + Enz-C2. Using half-reactions removes the ambiguity apparent in D and leads to a unique solution for the flux distribution (E). All flux solutions were obtained using the metabolic models defined in Supplemental File S3. this can be avoided by using the half-reactions based on the ping-pong mechanisms of TK and TA (Fig. 3, C and E). Theoretical analysis has shown that this formulation is essential for a correct description of the redistribution of label and that it allows the investigation of parallel pathways supported by different TK isozymes (Kleijn et al., 2005). Accordingly, the halfreaction scheme is used here to maximize the possibility of resolving cytosolic and plastidic contributions to the PPP fluxes.

Data Handling and Mathematical Modeling
Heterotrophic Arabidopsis cell cultures were incubated with [1-13 C]Glc, [2-13 C]Glc, or 10% [U-13 C 6 ]Glc, the label redistribution at isotopic and metabolic steady state was quantified by 13 C-NMR, and flux solutions, constrained by measurements of Glc consumption and biosynthetic outputs (Supplemental Tables S1-S4), were calculated in 13C-FLUX using the protocol summarized in Figures 4 and 5. Label was routinely measured in soluble metabolites, including amino acids, citrate, fumarate, GABA, malate, and Suc, in the amino acids derived from protein hydrolysate and in the Glc derived from starch. The best fit flux solutions presented here were based largely on Suc, starch, and protein hydrolysate labeling data, and full details of the measurements used in the fitting process are given in Supplemental File S2. Measurement errors were generally set at 5% for amino acid data and 15% for other metabolites, but a few measurements were assigned larger deviations as specified in Supplemental File S2.
Label measurements for a particular metabolite were initially grouped together for the fitting process, but subsequently it was established that better fits (lower residua) could be obtained by treating 13 C peaks with multiplet structures as separate subgroups. Simulations in 13C-FLUX showed that the reduction in the residuum could be attributed to the inclusion of additional scaling factors, compensating for discrepancies in the relative intensities of the signals from different carbon atoms within a metabolite, and that the subgrouping procedure did not compromise the accuracy of flux estimation. Subgrouping multiplet signals avoids the need to establish the relative intensities of multiplets from the same metabolite, and it is used routinely when two-dimensional NMR methods, which often mask the true relative intensities, are used to analyze [U-13 C 6 ]Glc labeling experiments (Sriram et al., 2004).
Each of the three labeling experiments was performed twice, and various combinations of experiments were analyzed simultaneously in 13C-FLUX by setting up models in which the reaction network was replicated two, three, or six times. Software limitations, arising from the large size of the network and the quantity of data, placed some restriction on this approach, but a three-substrate model, comprising a subnetwork for each labeled precursor (i.e. [1-13 C]Glc, [2-13 C]Glc, or 10% [U-13 C 6 ]Glc) and two data sets for each subnetwork, worked well, and a single flux solution corresponding to the entire data set could be generated by constraining equivalent fluxes in the three subnetworks to be equal in the model definition file. Ultimately, this procedure allowed a more reliable definition of the fluxes in the network than could be obtained by analyzing different labeling experiments separately (see below).
Preliminary fits of the data to the model were used to refine the model and to identify errors in the labeling data. In this iterative process, 13C-FLUX   was used to generate 50 to 100 solutions for a given model and data set. Changes were then made to either the model or the data set, for example, the addition of a pathway or the correction of an incorrectly annotated data point, and the fitting process was repeated. Comparing the residua for the five best fits provided the basis for accepting or rejecting a change to the model or the data set. For example, the reaction catalyzed by isocitrate dehydrogenase, and the mitochondrial uptake of pyruvate, were both defined as unidirectional because making them reversible had no effect on the fit. Similarly, CO 2 uptake at natural abundance was removed from the model because it had no discernible effect on the fit. In contrast, adding unidirectional plastidic pyruvate uptake improved the fit, as did the inclusion of separate cytosolic and plastidic pools for phosphoenolpyruvate and triose phosphate. However, the latter metabolites were combined to create a single pool of three-carbon phosphate esters in each compartment, since this removed indeterminable fluxes from the network and had no effect on the fit. The iteration process led to the model in Figure 2 and the validated data set in Supplemental File S2, thus providing a secure foundation for testing the extent to which the method could reveal the subcellular compartmentation of carbohydrate oxidation (Fig. 1).

Best Fit Fluxes and Statistical Analysis
The six data sets, comprising over 700 positional isotopomer measurements and 29 biomass fluxes, were initially analyzed in a three-substrate model based on the flux network in Figure 2. Estimates for the free fluxes were obtained by running the optimization routine 1,000 times, generating a correspondingly large set of flux solutions, each with a set of predicted labeling measurements. Some of these solutions failed to satisfy constraints in the model and were identified as infeasible solutions by the software; but the majority were feasible solutions, and these were extracted from the 13C-FLUX output file with a customized software tool (available from the authors on request). The residua for the feasible solutions varied over a wide range (Fig. 6A), and a subset of the solutions corresponding to residua below 3,000 was selected for further analysis. This cutoff was justified both by visual inspection of the distribution of residua ( Fig. 6A) and by principal component analysis of the mean-centered, unit variance-scaled flux solutions (Fig. 6B). The one-dimensional scores plot corresponding to the first principal component confirmed the similarity between the selected feasible solutions and provided further justification for rejecting the solutions with residua above the cutoff.
Monte Carlo simulation with bootstrap sampling of isotopomer measurements provides an effective method for characterizing the size and shape of the flux space (Wiback et al., 2004;Joshi et al., 2006a), and the selected solutions showed considerable variability in the definition of the free fluxes (Supplemental Fig.   S3). Some fluxes were tightly defined, for example, many of those in the tricarboxylic acid cycle, while others were normally distributed over a wider range of values. To reduce the hundreds of selected feasible solutions to a global best fit, the mean values from the Monte Carlo flux space were used as starting free flux values for a final run of the optimizer with the validated data set. This run was performed without bootstrap Monte Carlo sampling, and it resulted in a global best fit solution with the lowest residuum (Table II). These procedures were repeated for the models in which the full PPP was either confined to the plastids Feasible flux solutions with less than 3,000 residua (black circles) form a distinct cluster that is well resolved from those with more than 3,000 residua (white circles). The latter have a more scattered distribution that is similar to that obtained from random combinations of fluxes that produce infeasible flux maps (crosses). Table II. Influence of metabolic network structure on determination of flux estimates Fluxes were determined by fitting isotopomer abundances and biosynthetic output estimates to a family of models differing in the extent of PPP compartmentation ( Figs. 1 and 2). Each model contained a complete plastidic PPP pathway. Free fluxes that varied during the fitting procedure are indicated in boldface. Net fluxes and normalized hyperbolic transformed exchange fluxes are the global best fit estimate 6 SD as determined by EstimateStat and are calculated relative to the rate of Glc uptake (set to unity). A normalized exchange flux for which SD takes the estimated value outside the permissible [0,1] range is considered indeterminate (Ind).  (Fig. 1C), and the final best fit fluxes for the three reaction networks are given in Table II. Statistically reliable flux values were generated for all three models of the PPP (Table II) with only minor differences in residuum. The flux distribution for the steps involving sugar phosphates differed substantially for the three models (Fig. 7), and it was noticeable that the bulk of the flux through the oxidative steps of the PPP switched to the cytosolic pathway when this option was included in the network structure. Thus, while the flux through the oxidative steps was necessarily confined to the plastids in the absence of a cytosolic pathway (Figs. 1A and 7A), the plastidic pathway made only a minor contribution when the cytosolic pathway was included in the model (Figs. 1,B and C,and 7,B and C). The degrees of freedom in the models increased with the number of fluxes, and as expected, this led to a small improvement in the residuum achieved in the fitting procedure and an increase in the SD values for the flux estimates. For example, including the full PPP in both cytosol and plastids (Fig. 1C) produced a marginal improvement in residuum, but it produced a flux map in which 28 of the 33 fluxes shared by all three models had the largest SD values.
The validity of the flux maps was tested by comparing the predicted contribution of the carbon atoms of metabolized Glc with the relative specific activity of 14 CO 2 released from specifically labeled Glc (Fig. 8). Ratios of 14 CO 2 labeling from [1-14 C]-, [2-14 C]-, [3,4-14 C]-, and [6-14 C]Glc are typically used to characterize carbohydrate oxidation (Harrison and Kruger, 2008), and the pattern of 14 CO 2 release predicted by the global best fit flux estimates for each model was not significantly different from the experimentally determined pattern (Fig. 8), implying that each of the flux maps provides an adequate explanation of the observed data for the release of 14 CO 2 from positionally labeled Glc.
Since the models are indistinguishable, it can be concluded that the subcellular compartmentation of the PPP cannot be determined unambiguously from the available labeling data alone. However, based on current biochemical evidence, the most likely model is the one in which the oxidative steps of the PPP are duplicated in the cytosol and the plastids (Figs. 1B and 7B; Kruger and von Schaewen, 2003). The global best fit solution for this network is presented in detail in Supplemental Table S5, and it provides good agreement between the observed labeling data and the isotopomer abundances predicted from the global best fit flux map (Fig. 9A). There is also excellent correspondence between the measured, scaled, and predicted relative fractional abundance of positional isotopomers for individual metabolites, as exemplified by the data on protein-derived Asp/Asn (Fig. 9B). More than 80% of the total residuum was confined to 176 measurements (25% of total), and the flux map obtained after removal of these values was essentially identical to that obtained from the complete data set, suggesting that the poorly fitting measurements did not have a major influence on the flux solution (Supplemental Fig. S4). Furthermore, there was no consistency in the identities of the isotopomers that  Table II. contributed most to the residua between duplicate feeding studies involving identically labeled substrates, implying that the discrepancies between the observed and predicted values were due to poor estimates of measurement error rather than to a failure of the model to account adequately for the pattern of label redistribution in a particular section of the network.
Finally, in silico analysis of the chosen network revealed that the reliability of the flux estimates improved by combining data from labeling studies using [1-13 C]-, [2-13 C]-, and [U-13 C 6 ]Glc. This was particularly noticeable for the oxidative steps of the PPP in the cytosol and the plastids, and it also emphasized the importance of using well-defined labeling measurements (Fig. 10).

Choice of Experimental System and Application of Steady-State MFA
A plant cell suspension was chosen for the investigation to simplify the task of testing different models of the subcellular compartmentation of the PPP. The Arabidopsis cell suspension is experimentally tractable, it has been used successfully for MFA (Williams et al., 2008), and it avoids the difficulties that might arise from multiple cell types in differentiated tissues. Moreover, many of the general metabolic features of the chosen Arabidopsis culture are similar to Figure 8. Validation of Arabidopsis flux map by radiorespirometric analysis. The ratios of the specific activity of 14 CO 2 released from cell cultures supplied with [ 14 C]Glc specifically labeled in the carbon positions indicated were compared with those predicted from the flux maps. The latter were derived from the specific abundance of label in CO 2 determined by simulations in CumoNet using different Glc isotopomers as the input substrate. Experimental data are shown as means 6 SD of ratios from four replicate measurements for each of two independent cell cultures (white bars). The predicted ratios were determined from the global best fit solutions (Table II; Supplemental  Table S5) for the networks in which the PPP is confined to the plastid (model A; light gray bars), the oxidative steps of the PPP are located in both the cytosol and plastids (model B; dark gray bars), or the entire pathway is present in both compartments (model C; black bars). There is no significant difference in the experimentally determined ratios for the two separately grown cultures as assessed by repeated-measures ANOVA, and all the predicted ratios are within the 95% prediction intervals of the ratios determined experimentally. Figure 9. Assessment of goodness of fit of isotopomer data to the global best fit solution of the metabolic network defined in Figure 2. Data from cell cultures labeled separately with [1-13 C]Glc, [2-13 C]Glc, and [U-13 C 6 ]Glc were used in these evaluations. A, Comparison of predicted and measured isotopomer abundances for all metabolites analyzed. B, Comparison of the predicted relative fractional abundance of isotopomers of Asp with measured values, before and after scaling of subgroups of 13 C peaks with multiplet structure. Shading is conserved for individual isotopomers in the pie charts for a given labeled substrate. the heterotrophic metabolism that occurs in roots (Lehmann et al., 2009).
The steady-state MFA protocol used here was broadly similar to existing best practice, but the application of this tool in plant biology is still developing, and several new features were introduced with the aim of further improving the precision and reliability of the flux estimates. First, half-reactions were used to represent the underlying carbon transitions that occur during the steps catalyzed by TA and TK. This approach, which provides a more accurate description of the redistribution of the label by the PPP, has not, to our knowledge, been used previously in steady-state MFA of plants and has only rarely been used in microbial work. Second, Monte Carlo simulations with bootstrap sampling of the isotopomer abundances provided an effective method for exploring the flux space that was consistent with the observed redistribution of label. Stochastic sampling of the starting fluxes is routine in the modeling process (Schwender et al., 2006), and combining this with statistical sampling of the data points improves the fitting procedure (Kelly et al., 1990;Namba, 2004). Surprisingly, the availability of a subroutine to exploit this method in the commonly used 13C-FLUX software appears to have been overlooked. Third, the output from the Monte Carlo simulations was subjected to principal component analysis to identify the set of best fit feasible solutions. This novel approach provided an objective method for analyzing the outputs of the modeling process and an effective route to the global best fit. Finally, redistribution of radioactivity from U-14 C-substrate was used as a convenient method for ensuring that all major metabolic products had been identified and as a sensitive method for determining the output fluxes.

Quantifying the Compartmented Fluxes of Central Metabolism
The subcellular compartmentation of metabolism in eukaryotes complicates steady-state MFA (Roscher et al., 2000), and the problem is acute in plants because of the more extensive duplication of steps and pathways in the metabolic network . Capturing the characteristic compartmentation of the pathways of carbohydrate oxidation in a realistic model of heterotrophic plant metabolism exemplifies the problem, and it would be particularly useful if this could be done by comparing models for networks compartmented to different extents, since this would allow an in vivo assessment of the functional compartmentation of the network. Since most steadystate MFA models are based on overdetermined data sets, reflecting the wealth of stable isotope labeling information that can be obtained from NMR or mass spectrometry (Ratcliffe and Shachar-Hill, 2006), it is easy to add extra steps to the model to allow it to match biochemical reality more closely. However, the relationship between fluxes and measurements is not always obvious, and even if the extra fluxes are potentially resolvable, they may not be determinable with any great confidence.
In practice, many steady-state MFA investigations of the PPP in plants use models in which the compartmentation of the pathway is predefined. For example, Figure 10. Statistical analysis of flux estimates obtained with different steady-state labeling strategies. The SD values of optimum flux estimates were derived from the global best fit solution of the metabolic network defined in Figure 2 using the EstimateStat routine in 13C-FLUX. Isotopomer measurements obtained from analysis of extracts of cell cultures labeled separately with [1-13 C]Glc (1-13 C), [2-13 C]Glc (2-13 C), or [U-13 C 6 ]Glc (U-13 C) were used independently or in combination. Fluxes were expressed relative to the rate of Glc consumption (set to unity), and SD values are coded from black (low; flux well defined) to white (high; flux poorly defined), as indicated by the logarithmic gray scale at the top. The bottom row contains an analysis of data combined from all three feeding strategies (1/2/U-13 C), in which the SD of all isotopomer measurements was reduced to 5%. a single plastidic pathway has been assumed in maize root tips (Dieuaide-Noubhani et al., 1995;Alonso et al., 2007b), tomato (Solanum lycopersicum) cell cultures (Rontein et al., 2002), and Arabidopsis cell cultures (Williams et al., 2008), whereas complete duplication of the pathway in the cytosol and plastid has been assumed in the analysis of cultured soybean embryos (Sriram et al., 2004;Iyer et al., 2008). In some tissues, such as oilseed rape embryos (Schwender et al., 2003) and soybean embryos (Allen et al., 2009), the modeling task is made easier because the observed labeling patterns imply that intermediates in glycolysis and the PPP are in fast exchange between the cytosol and the plastids, ensuring that duplicate pathways in different physical locations function indistinguishably from a single uncompartmented pathway. However, in most of the tissues that have been examined so far, including soybean embryos in some instances (Sriram et al., 2004), the hexose phosphate pools in the cytosol and the plastids are found to be labeled to different extents. Typically, this evidence might be based on differences in the labeling patterns of the glucosyl moieties from Suc and starch (Supplemental Fig. S2), and it is then necessary to consider the subcellular compartmentation of glycolysis and the oxidative PPP explicitly.
Inspection of the labeling patterns in products derived specifically from either the cytosolic or plastidic hexose phosphate pool may provide circumstantial evidence for the operation of the full PPP in both compartments (Krook et al., 1999), justifying the construction of a model with duplicate pathways. However, intuitive interpretation of the data is complicated by the multiple metabolic processes that contribute to the labeling of these metabolites , and a more robust quantitative approach is needed to compare the results of fitting all the labeling data to models in which the pathways are duplicated to different extents. Here, steady-state MFA was used to compare three different models of the oxidative PPP (Fig. 1), and the analysis indicates that several aspects of the protocol contribute to its effective implementation. First, it is necessary to work with overdetermined models so that the statistical significance of the flux solutions can be properly evaluated. Second, for the PPP, it is important to define the carbon transition network in terms of the half-reactions for TA and TK to maximize the chance of being able to discriminate between parallel events in the cytosol and plastid. Third, it is essential to use high-quality data from experiments with several differently labeled substrates in order to deduce fluxes with sufficient confidence to allow a comparison of the models (Fig. 10). Finally, it is necessary to establish a procedure (Figs. 4 and 5) that leads to an optimal best fit based on clearly defined criteria so that the results for different models can be compared.
Analysis of three models for the oxidative PPP ( Fig.  1) emphasized the difficulty of determining the flux distribution with sufficient confidence to be able to discriminate between the models. There was little to distinguish a plastid-only model (Fig. 1A) from the model in which the oxidative steps were duplicated in the cytosol (Fig. 1B), even though the latter is more plausible on the basis of molecular evidence. The increase in the degrees of freedom in the second model (Fig. 1B) produced only a marginal improvement in the residuum obtained for the global best fit, and the main observations that support the model are the nonrandom allocation of flux between the oxidative steps in the two compartments (Supplemental Fig. S3) and the reliability of these flux estimates, as reflected in their confidence intervals (Table II). Adding further degrees of freedom by duplicating the complete oxidative PPP in the two compartments (Fig. 1C) led to a further, modest improvement in the fit and a radically different flux distribution through the pathways of carbohydrate oxidation (Table II), but it also reduced the confidence in the flux values, with most of the fluxes showing an increase in SD. Thus, this model provides a less reliable flux map and can be rejected on the basis that the fit to the available data is less convincing. Such a conclusion is always subject to the proviso that data might become available that would put further constraints on the assessment of the compartmentation of the pathway (Kruger and Ratcliffe, 2009), but currently, the available data imply only a plastidic location for the nonoxidative steps of the PPP in Arabidopsis cells.
In silico analysis of the sensitivity of cumomer abundance to changes in flux through the irreversible oxidative section of the PPP can identify the isotopomer measurements that are best able to define the fluxes through this step in the cytosol and plastids (Supplemental File S4). Although the sensitivities are influenced by the label distribution in the supplied substrate and by the flux distribution through the network, estimates based on the biochemically preferred flux map (Fig. 7B) indicate that isotopomers of cytosolic and plastidic hexose phosphates and compounds derived directly from these pools are among the most diagnostic measures of these PPP fluxes. This discrimination arises in part because, in the absence of plastidic Fru-1,6-bisphosphatase, cytosolic Glc-6-P is the only substrate replenishing the plastidic Glc-6-P pool, and any difference in labeling must therefore stem from plastidic PPP activity. However, intermediates within the reversible nonoxidative section of the PPP, and their derivatives such as His and aromatic amino acids, are also relatively sensitive to perturbations in PPP flux. Thus, measurements of label distribution in these compounds would be expected to further improve the flux estimates and could lead to greater discrimination between the different models.
A similar, but less comprehensive, analysis of the compartmentation of metabolism in sunflower embryos found that models of the PPP that included cytosolic steps were no improvement over a model in which the PPP was confined to the plastids (Alonso et al., 2007a). However, the flux distributions in the alternative models were not discussed, and it is per-haps surprising that the additional degrees of freedom in the models with cytosolic steps led to no improvement in the fit. Irrespective of these points, since all the models gave effectively the same fit to the data, it is unclear on what basis one of them was preferred.
The difficulties encountered in analyzing the compartmentation of the oxidative PPP by steady-state MFA highlight some of the practical limitations of the approach. The objective is to capture a realistic metabolic phenotype, but in practice the flux map is a representation of the movement of the chosen label through the system. Resolving the contribution of chemically distinct pathways between two metabolites in the same compartment, or of duplicated steps and pathways in distinct compartments, is not always possible or even attempted. For example, the steps from citrate to 2-oxoglutarate are routinely assigned to the mitochondrion in steady-state MFA, even though there may be parallel fluxes in the cytosol. Thus, if a function of MFA is to guide metabolic engineering Libourel and Shachar-Hill, 2008), then current flux maps are unlikely to be useful in predicting the result of manipulating a compartment-specific enzyme activity. Similarly, even if MFA provides evidence that a physically compartmented pathway is functionally uncompartmented as a result of fast exchange of intermediates (Schwender et al., 2003;Allen et al., 2009), the resulting metabolic phenotype still conceals the associated compartmental demands for reducing/oxidizing power and energy. It can be concluded that it is essential to push steady-state MFA to its limits to obtain the most informative model for the flux distribution.

Metabolism in a Heterotrophic Arabidopsis Cell Suspension
The network structure used to determine the biochemically preferred flux map (Fig. 7B) differs from the one used in the previous analysis of Arabidopsis cell suspension cultures (Williams et al., 2008) in two main ways. First, the models differ in their treatment of the lower section of glycolysis. In steady-state MFA of plant metabolism, there is usually insufficient labeling information to support a model with separate cytosolic and plastidic pools of both glyceraldehyde 3-phosphate and phosphoenolpyruvate. In the earlier model (Williams et al., 2008), each metabolite was assumed to be in a common uncompartmented pool, implying fast exchange across the plastid envelope, whereas here it was assumed that the two metabolites were at isotopic equilibrium in a three-carbon phosphate ester pool (Fig. 2) within each compartment. This assumption allowed the degree of compartmentation of the pools to be tested by not arbitrarily assuming fast exchange between the cytosol and plastid, and it led to a marked alteration in the flux distribution, since all three flux maps showed a significant net flux from the plastid to the cytosol (Fig. 7) at the level of the three-carbon phosphate ester pool and only a limited exchange flux (Table II). Second, the biochemically preferred model that emerged from the comparison of three different models for the oxidative PPP (Fig. 1) allowed the oxidative steps to operate in both the cytosol and plastids. This compartmentation is supported by estimates of transcript abundance based on publicly available microarray data for genes encoding isozymes of the PPP (Supplemental Fig. S5A) and associated plastid envelope translocators (Supplemental Fig. S5B).
A striking feature of the flux map generated from the biochemically preferred model is that the oxidative steps of the PPP are almost completely confined to the cytosol. Since most biosynthetic processes requiring NADPH are plastidic, this flux distribution implies extensive exchange of reducing power equivalents across the plastid envelope. However, this is unlikely to be problematic, since previous studies on maize root tips suggest that the oxidative PPP in the cytosol and plastids can cooperate in providing NADPH for biosynthesis throughout the cell (Averill et al., 1998). Nevertheless, the dominance of the cytosolic pathway in catalyzing flux through the oxidative section of the PPP is unexpected, and it is unclear whether this is a general characteristic of the organization of carbohydrate metabolism in higher plants, since compartmentation of the PPP is not included in most flux maps. The only other plant tissues for which information is available are rosy periwinkle (Catharanthus roseus) hairy root cultures and developing soybean embryos, and in both of these only 25% to 50% of flux through the oxidative section of the PPP is cytosolic (Sriram et al., 2004(Sriram et al., , 2007Iyer et al., 2008). However, in these systems, total flux through the oxidative PPP is considerably higher than that observed in Arabidopsis cells, and when expressed as a proportion of the rate of sugar consumption, the cytosolic fluxes in all three systems are similar (rosy periwinkle, 42.6%; soybean, 47.4%; Arabidopsis, 37.2%), implying that variation in the compartmentation of PPP flux is due to fluctuations in plastidic activities rather than to a decrease in the significance of the cytosolic flux. Moreover, irrespective of the exact proportion of the oxidative PPP flux occurring in the cytosol, two recent studies demonstrate the importance of the cytosolic activity. First, disrupting the expression of the two genes encoding cytosolic isozymes of Glc-6-P dehydrogenase (G6PD5 and G6PD6) in Arabidopsis through insertional inactivation alters seed lipid content (Wakao et al., 2008). Second, supplementation or replacement of cytosolic Glc-6-P dehydrogenase with an alternative isoform improves the resistance of tobacco (Nicotiana tabacum) leaves to Phytophthora infestans, and this is attributed to changes in the provision of NADPH and stimulation of the hypersensitive defense response (Scharte et al., 2009). Both observations highlight a requirement for the oxidative section of the PPP in the cytosol to support normal metabolic activity.

Energy and Redox Balance in Plant Metabolic Networks
The metabolic activity of the Arabidopsis cell culture appears to be broadly representative of a wide range of heterotrophic plant systems, and one measure of this can be found in its efficiency in converting sugar to biomass components. The proportion of the [U-14 C]Glc metabolized to macromolecules and soluble compounds (Supplemental Table S1) indicates a carbon use efficiency of 60%, in good agreement with the 65% calculated from a comparison of substrate uptake and estimated CO 2 output in the previous study of Arabidopsis cells grown under equivalent conditions (Williams et al., 2008). These estimates of carbon use efficiency are within the range obtained for other heterotrophic plant systems, including a rosy periwinkle hairy root culture (24%; Sriram et al., 2007), maize root tips (42%-47%; Alonso et al., 2007b), developing sunflower embryos (50%; Alonso et al., 2007a), and a tomato cell culture (52%-68%; Rontein et al., 2002). Although considerably higher values have been reported for developing embryos of oilseed rape (85%-95%) and soybean (82%-83%), in these systems there is appreciable reassimilation of respired CO 2 through Rubisco, and light makes significant contributions to the provision of ATP and/or reductant required for biosynthesis (Goffman et al., 2005;Allen et al., 2009). In contrast, the energy and reducing power required for biosynthesis to support growth in heterotrophic Arabidopsis cells must be generated exclusively by oxidation of Glc, the sole respiratory substrate, through the central pathways of carbon metabolism.
Comparison of the oxidation state of the substrates in the culture medium and the products generated by the Arabidopsis cells confirms that the metabolic network is able to operate without an input of reducing power. The primary substrates for cell growth were Glc, nitrate plus ammonium, and sulfate, in which the average oxidation states of carbon, nitrogen, and sulfur are 0, +1, and +6, respectively. The analysis in Supplemental File S5 shows average carbon oxidation states of 0 for carbohydrates, 21.56 for lipid, +0.92 for organic acids, and 21.02 for protein and soluble amino acids in which nitrogen and sulfur have oxidation states of 23 and 22, respectively. The corresponding values for carbon in CO 2 and ethanol, the principal respiratory end products, are +4 and 22, respectively. It follows from the relative partitioning of substrate during growth (Table I) that the average change in oxidation state in metabolic end products is +1.33 per carbon (Supplemental File S5). The analysis also shows that assimilation of nitrate used in amino acid production represents the greatest demand for reducing power in the Arabidopsis cell culture under the growth conditions used in this study.
However, this simple redox analysis ignores the coenzyme specificity of different processes and does not consider the preferential involvement of NADPH in many reductive biosynthetic steps. This limitation can be overcome by using flux maps to evaluate the coenzyme balance in the network. In contrast to the situation in developing soybean embryos (Iyer et al., 2008), inspection of the fluxes determined for the Arabidopsis cells establishes that the oxidative PPP Table III. Coenzyme requirements for the biosynthetic outputs from the central metabolic network in heterotrophic Arabidopsis cells The coenzyme requirements for biosynthesis were determined from the fluxes in Supplemental Table S4 using the reaction stoichiometries in  Supplemental Table S6 and are presented as means 6 SE based on the estimates of label partitioning following metabolism of [U-14 C]Glc in Table I. The coenzyme requirements for lipid assume the synthesis of triacylglycerol containing only oleoyl (C18:0) side chains. Costs for amino acid production are based on the assimilation of equal amounts of NO 3 2 and NH 4 + . The ATP requirement for protein production is calculated from the cost of synthesizing proteins from the component amino acids and is equivalent to 4 ATP per amino acid residue. Accumulating organic acids are withdrawn directly from intermediate pools in the network and incur no additional cost. Generation of coenzymes by the central metabolic network is determined from the net fluxes of the internal interconversions (i.e. excluding output steps) of the model defined in Figure 2 using conventional reaction stoichiometries and ignoring the possible contribution of inorganic pyrophosphate as a phosphoryl donor. NADH and NADPH yields (in italics) assume conversion of isocitrate to 2-oxoglutarate by NADP + -dependent isocitrate dehydrogenase. ATP generation is attributed to substratelevel phosphorylation, and FADH 2 production is assumed to be equivalent to 0.63 NADH. Values for coenzyme generation are means 6 SE of estimates obtained from the best fit flux solutions of 501 Monte Carlo simulations; values in parentheses are determined from the global best fit flux solution of the network in which the oxidative steps of the PPP are located in both the cytosol and plastids, as detailed in Supplemental Table S5 is insufficient to meet the NADPH requirements for biosynthesis predicted by the network (Table III). This observation is the result of two assumptions: (1) the assimilation of equal amounts of nitrate and ammonium (Table III); and (2) the NADPH dependence of Glu synthesis (Supplemental Table S6). In this situation, reductive biosynthesis will require other NADP +dependent processes, and in the network considered here, the demand could be met by the operation of NADP + -isocitrate dehydrogenase (Table III).
In contrast, the production of NADH is more than enough for its utilization in reductive processes, including the demands for ATP synthesis. Indeed, if it is assumed that excess NADH (and FADH 2 ) are reoxidized through the mitochondrial electron transport chain to produce ATP through oxidative phosphorylation (at stoichiometries of 2.3-2.5 and 1.4-1.5 ATP per NADH and FADH 2 , respectively; Hinkle, 2005), then the maximum yield of ATP generated by the metabolic network is over seven times that required to support the flux to biosynthetic products. The actual yield of ATP will depend on the extent of engagement of the alternative oxidase and the activity of mitochondrial uncoupling proteins, but the comparison suggests that the metabolic network has a considerable overcapacity for ATP production. Of course, these calculations ignore the energy requirements of other metabolic processes, such as substrate cycling, turnover of proteins and other macromolecules, and synthesis and turnover of RNA, that are not included in this network. Moreover, ATP is used in other cellular processes, such as substrate uptake, transport of metabolites across intracellular membranes, and maintenance of transmembrane ion and electrical gradients. Such activities can account for over half of the cellular ATP consumption in microorganisms (Nielsen et al., 2003), and the findings presented here suggest that an even greater proportion of the potential chemical energy released during metabolism in actively growing plant cells could be needed to meet these requirements. Recently, analysis of energy demands in a genome-scale stoichiometric model of Arabidopsis has led to a similar conclusion that only a relatively small proportion of the potential ATP generated by the cell is used for biosynthesis (Poolman et al., 2009). This conclusion has implications for the determination of flux through plant metabolic networks by flux balance analysis. This approach depends on the selection of an appropriate objective function to obtain an optimum solution within the potential flux space defined by the stoichiometry of the biochemical reaction network, plus other constraints, such as the metabolic outputs of the system (Kauffman et al., 2003). Typically, the optimization criterion in microorganisms is maximization of yield of biomass, ATP, or a product of interest, or alternatively, maximization of flux (i.e. maximization of the rate of production of one or more components required for growth; Schuster et al., 2008). A recent study on developing endosperm of barley (Hordeum vulgare) used a refine-ment of this approach in which the objective function was maximization of growth per unit of flux, which is a combination of maximization of biomass yield and minimization of overall flux (Grafahrend-Belau et al., 2009). The logic behind this optimization criterion is that cellular metabolism has evolved to fulfill its functions with the most economic use of available resources (Holzhü tter, 2004), but application of this criterion requires that the major energy-consuming processes have been adequately incorporated into the metabolic model. The demonstration that a large proportion of the potential biochemical energy released through respiration in heterotrophic plant cells is dissipated, or used by processes other than those directly associated with the synthesis of macromolecules and soluble components that accumulate during growth, means that ATP expenditure for such purposes must be accurately assessed if flux balance analysis is to generate meaningful flux maps of the network of central carbon metabolism in plants.

CONCLUSION
Steady-state MFA of multiple labeling experiments in an Arabidopsis cell suspension shows that the experimental observations can be explained by any one of three different models of the subcellular compartmentation of the oxidative PPP. Other biochemical evidence favors the model in which the oxidative steps are duplicated in the cytosol and plastids, and this leads to the conclusion that the bulk of the flux through the oxidative steps is cytosolic. This illustrates the point that interesting features of the metabolic phenotype can be missed with an inappropriate model and emphasizes the importance of basing the model on all the available biochemical evidence.

Cell Culture
Cell suspensions of Arabidopsis (Arabidopsis thaliana ecotype Landsberg erecta) were maintained on an orbital shaker at 150 rpm in 250-mL Erlenmeyer flasks sealed with a double layer of aluminum foil under a 16-h-light/8-h-dark cycle at 22°C (Williams et al., 2008). Every 7 d, 15 mL of cell suspension was subcultured into 100 mL of fresh Murashige and Skoog (MS) medium supplemented with 3% (w/v) Glc, 0.5 mg L 21 1-naphthaleneacetic acid, and 0.05 mg L 21 kinetin. Heterotrophic cell suspensions were produced by subculturing 15 mL of a 7-d-old light-grown cell suspension into 100 mL of fresh medium and incubating in the dark at 22°C.

Incubation of Arabidopsis Cell Culture with 13 C-Labeled Glc
Cells were labeled by subculturing light-grown cells into medium in which the unlabeled Glc was replaced with [1-13 C]Glc, [2-13 C]Glc, or a mixture of 10% [U-13 C 6 ]Glc and 90% unlabeled Glc. Cells were grown in the dark for 7 d at 22°C and then harvested by vacuum filtration, rinsed with 80 mL of MS medium, frozen by immersion in liquid nitrogen, and stored at 280°C. This procedure has been shown to have no discernible effect on the flux distribution in the core metabolic network of Arabidopsis (Kruger et al., 2007a) and to lead to the establishment of an isotopic and metabolic state in the cell culture (Williams et al., 2008).

Extraction Procedures after Stable Isotope Labeling
Soluble metabolites were extracted from frozen cells (approximately 8 g fresh weight) using perchloric acid (Kruger et al., 2007a. Subsequently, NMR spectra were recorded from samples dissolved in 3.2 mL of 2 H 2 O containing 10 mM EDTA, 25 mM 1,4-dioxane, and 10 mM KH 2 PO 4 /K 2 HPO 4 , pH 7.5. Starch and protein were extracted from the insoluble residue using procedures that were similar to those described elsewhere (Williams et al., 2008). Starch was gelatinized by autoclaving a washed aliquot of the residue in 25 mM sodium acetate for 3 h, followed by enzymic digestion with a-amylase and amyloglucosidase. The supernatant, containing Glc, was freeze dried and redissolved in 3.2 mL of 2 H 2 O with 25 mM 1,4-dioxane for NMR analysis. Protein in the perchloric acid extraction pellet was hydrolyzed directly by adding 2 mL of 6 M HCl and heating under vacuum at 110°C for 24 h in a 5-mL Pierce hydrolysis tube (http://www.piercenet.com/). After digestion, the protein hydrolysate was centrifuged and the supernatant was freeze dried. The sample was redissolved in water and freeze dried again after adjusting the pH to 7.5. Finally, the sample was redissolved in 3.2 mL of 2 H 2 O containing 25 mM 1,4-dioxane for NMR analysis.

Quantification of Biosynthetic Outputs
Fluxes to biosynthetic products, and the rate of Glc consumption, were obtained from the redistribution of label following metabolism of [U-14 C]Glc by cell suspensions.

Incubation of Arabidopsis Cell Culture with [U-14 C]Glc
Cells were grown with unlabeled Glc under the conditions used for the stable isotope labeling experiments. After 6 d, 5-mL aliquots of the cell suspension were transferred to sealed 100-mL conical flasks and supplemented with 185 kBq [U-14 C]Glc (specific activity, 10.4 TBq mol 21 ). After 12 or 24 h, the cells were harvested by filtration and washed with two volumes of MS medium, and then the cells and medium were frozen in liquid nitrogen. The production of 14 CO 2 was determined in replicate flasks in which the incubation was terminated at the appropriate time by injection of 2 mL of 12 M formic acid to stop metabolism and release dissolved CO 2 from the medium (Harrison and Kruger, 2008). Released 14 CO 2 was collected in 1.0 mL of 10% (w/v) KOH in a vial suspended within each sealed culture flask.

Extraction of Arabidopsis Cells after 14 C Labeling
Frozen cells were combined with 1 mL of 100% methanol at 220°C and frozen in liquid nitrogen. The sample was thawed on ice for 2 h and then centrifuged at 8,500g for 10 min at 5°C. The supernatant was removed and stored on ice. The pellet was resuspended in 1 mL of cold methanol, and after storage on ice for 5 min, the suspension was centrifuged and the supernatant was combined with the supernatant from the initial extraction. The insoluble cell residue was then extracted twice with 1 mL of 100% methanol for 5 min at 70°C. Samples were centrifuged to separate methanol-soluble material from the insoluble residue, and the supernatants were combined with that from the cold extractions. Finally, the insoluble residue was washed with two further 1-mL aliquots of 100% methanol, and the insoluble material was stored at 280°C.

Ethanol Extraction from the Incubation Medium after 14 C Labeling
[ 14 C]Ethanol was recovered from the cell culture medium by distillation after adding 2 mL of 100% ethanol to the sample to aid recovery.

Fractionation of Methanolic Cell Extracts
Distilled water and chloroform were added to the methanolic fraction to generate a CHCl 3 :CH 3 OH:water mixture (4:4:2, v/v). The mixture was centrifuged at 3,500g for 10 min at 20°C to promote phase separation. The aqueous methanol phase was retained, and the chloroform phase was rinsed with a further two washes of 2 mL of water plus 4 mL of methanol before storage at 220°C. The aqueous methanol fractions from each wash were pooled with the initial aqueous methanol fraction, and this combined aqueous methanol fraction was then dried by rotary evaporation at 25°C and resuspended in 5 mL of water. Total radioactivity in this fraction was measured before further separation into basic, acidic, and neutral fractions by solid-phase extraction on Dowex ion-exchange columns (Kruger et al., 2007a). The aqueous basic fraction, consisting of amino acids, was freeze dried, resuspended in 1 mL of water, and fractionated using an anion-exchange column to separate acidic amino acids from neutral/basic amino acids (Cossins and Beevers, 1963). Each neutral fraction was freeze dried, resuspended in 1 mL of water, and separated into Suc, Glc, and Fru using thin-layer chromatography (Scott and Kruger, 1995).
The chloroform (lipid) fraction was separated into its components by reverse-phase chromatography on an aminopropyl column (LC-NH 2 SPE; Supelco) as described by Kim and Salem (1990).
Total radioactivity in the methanol-insoluble material was determined by incubating a 100-mL aliquot of insoluble residue overnight with 400 mL of Solvable tissue solubilizer (Perkin-Elmer; http://las.perkinelmer.com) followed by liquid scintillation counting. Insoluble material was further analyzed by autoclaving, digestion with amylase, amyloglucosidase, and pronase, and separation into fractions by ion-exchange chromatography (Malone et al., 2006). The resulting aqueous, acidic, and basic fractions represented starch, cell wall material, and protein, respectively. Radioactivity in the remaining undigested residue was determined by washing the residue to remove the supernatant including digested material and enzymes, followed by overnight incubation with tissue solubilizer (Kruger et al., 2007a).

Monitoring Release of 14 CO 2 by Arabidopsis Cell Culture Metabolizing Positionally Labeled [ 14 C]Glc
A 5-d-old Arabidopsis cell suspension was incubated in the dark at 22°C in a 100-mL conical flask in a total volume of 5 mL of culture medium supplemented with 3.7 kBq [1-14 C]-, [2-14 C]-, [3,4-14 C]-, or [6-14 C]Glc. Each flask was then sealed with a rubber bung and aerated on an orbital shaker at 100 rpm. Released 14 CO 2 was collected in 0.5 mL of 10% (w/v) KOH in a vial suspended in the flask. The KOH solution was replaced at 2-h intervals for 12 h and again 24, 36, and 48 h after the beginning of the incubation. Ratios of 14 CO 2 release from positionally labeled Glc were determined as described by Harrison and Kruger (2008).

Determination of Radioactivity
Radioactivity was determined by liquid scintillation counting in a Beckman LS 6500 scintillation counter. An appropriate aliquot of each fraction was mixed with 4 volumes of Optiphase HiSafe scintillation fluid. Counting efficiency was typically greater than 90%.

NMR Spectroscopy
One-dimensional 1 H-decoupled 13 C-NMR spectra were recorded at 20°C and 150.9 MHz on a Varian Unity Inova 600 spectrometer equipped with a 10mm-diameter broadband probe. Waltz 16 1 H decoupling was applied during the acquisition time, and spectra were acquired with either a 6-s relaxation delay and low-power frequency-modulated decoupling to generate the nuclear Overhauser effect or a 29-s relaxation delay without the nuclear Overhauser effect. Large numbers of scans were generally required to generate spectra of sufficient quality for quantitative analysis of the peak intensities, and spectra were accumulated in blocks to check for unwanted spectral changes. For example, a total acquisition time of 30 h was used for the 6-s spectra of the soluble metabolites, and over 40 h was used for the 29-s spectra of protein hydrolysates. Spectra of starch digests were obtained much more quickly, reflecting the high Glc concentration in the samples. Spectra were processed using NUTS software (Acorn NMR), and in some instances, relative peak intensities in 6-s spectra were corrected using relative intensities in the fully relaxed spectra. Spectra were processed after a Lorentzian-Gaussian transformation, with a line broadening of 21 Hz and a Gaussian factor of 0.1, and signal integration was performed after baseline correction. Chemical shifts were measured relative to the 1,4-dioxane signal at 67.30 ppm, and peaks were assigned on the basis of literature values and authentic standards. Positional isotopomers were specified using 0 to denote 12 C, 1 to denote 13 C, and X to denote either 12 C or 13 C.

Metabolic Modeling
Metabolic modeling was performed using 13C-FLUX (version 20050329;Wiechert et al., 2001) according to the protocol outlined in Figures 4 and 5. After defining the metabolic network, the free fluxes were fitted to the biomass measurements and labeling data using the Donlp2 (P. Spelluci, Technische Universität Darmstadt, Germany) routine in 13C-FLUX. This sequential quadratic programming algorithm is superior to the CooolEvoAlpha evolutionary algorithm in 13C-FLUX, in terms of both computational efficiency and avoidance of local minima (Wiechert et al., 2001). The optimizer was used in conjunction with Monte Carlo simulations and bootstrap sampling of the isotopomer abundances based on the measured values and their deviations to generate solutions for the flux map. The bootstrap Monte Carlo sampling method was invoked by the command: Donlp2 "model file name".ftbl -a 0.00001 -m "number of simulations" -bootstrap -logAllF . "output file name".txt&watch "grep Resi "output file name".txt|tail 220". This strategy of combining Donlp2 with bootstrap Monte Carlo sampling was first tested on the network model distributed with 13C-FLUX and on a published model for developing sunflower (Helianthus annuus) embryos (Alonso et al., 2007a) to confirm that it was capable of generating the expected fluxes. Each simulation generated a set of free fluxes after minimizing the sum of squared standardized differences (the residuum) between the experimental measurements and their predicted values. Principal component analysis of the free fluxes obtained from 1,000 simulations was performed using SIMCA-P 11.5 (Umetrics) to identify the set of best fit feasible solutions with low values of the residuum. This set of outputs was used to generate a mean flux solution that was then used as a starting point for a further optimization run with Donlp2. EstimateStat, the statistical analysis component of 13C-FLUX, was used to compare the merits of different modeling schemes and to assign errors to the flux values deduced by the numerical analysis.

Supplemental Data
The following materials are available in the online version of this article.
Supplemental Figure S2. Comparison of relative fractional abundance of isotopomers in carbohydrates derived from cytosolic and plastidic hexose phosphate pools.
Supplemental Figure S3. Distribution of flux estimates for selected steps in the metabolic network from Monte Carlo simulations.
Supplemental Figure S4. Assessment of robustness of global best fit flux solution.
Supplemental Figure S5. Relative expression levels of genes encoding isozymes of the PPP and associated plastid envelope translocators in Arabidopsis cells.
Supplemental Table S1. Calculation of fluxes to major metabolic end products.
Supplemental Table S2. Amino acid composition of protein in cell suspension culture.
Supplemental Table S3. Fluxes to amino acids.
Supplemental Table S4. Summary of biosynthetic output constraints applied in flux determinations using 13C-FLUX.
Supplemental Table S5. Metabolic fluxes in heterotrophic Arabidopsis cells.
Supplemental Table S6. Coenzyme stoichiometries for generation of biosynthetic products from intermediates of the central metabolic network.
Supplemental Table S7. Relative fractional abundance of isotopomers in carbohydrates derived from cytosolic and plastidic hexose phosphate pools.
Supplemental File S1. 13C-FLUX model of the central metabolic network in Arabidopsis cell suspension culture (Microsoft Excel file).