Chapter 8 Mathematical Modeling of Isotope Labeling Experiments for Metabolic Flux Analysis Shilpa Nargund and Ganesh Sriram Abstract Isotope labeling experiments (ILEs) offer a powerful methodology to perform metabolic flux analysis. However, the task of interpreting data from these experiments to evaluate flux values requires significant mathematical modeling skills. Toward this, this chapter provides background information and examples to enable the reader to (1) model metabolic networks, (2) simulate ILEs, and (3) understand the optimization and statistical methods commonly used for flux evaluation. A compartmentalized model of plant glycolysis and pentose phosphate pathway illustrates the reconstruction of a typical metabolic network, whereas a simpler example network illustrates the underlying metabolite and isotopomer balancing techniques. We also discuss the salient features of commonly used flux estimation software 13CFLUX2, Metran, NMR2Flux+, FiatFlux, and OpenFLUX. Furthermore, we briefly discuss methods to improve flux estimates. A graphical checklist at the end of the chapter provides a reader a quick reference to the mathematical modeling concepts and resources. Key words Metabolic flux analysis, Isotope labeling experiment, Mathematical model, Isotopomer balance, Cumomer, Elementary metabolite unit, Bondomer, 13CFLUX2, Metran, NMR2Flux+, FiatFlux, OpenFLUX

1  Introduction Metabolic fluxes are rates of carbon traffic through metabolic pathways and are important indicators of cell physiology [1, 2]. Quantifying fluxes in large metabolic networks can be a challenging problem [3–5]. Whereas a small number of fluxes (e.g., the exchange of metabolites between cells or tissues and their surrounding medium, or fluxes leading to accumulation of biomass) can be directly measured, a majority of fluxes can only be estimated indirectly. This is especially the case for flux distributions at many intracellular metabolic branchpoints or through cyclic pathways [6]. A powerful technique for indirectly quantifying such intracellular fluxes is isotope-assisted metabolic flux analysis (isotope MFA). This technique involves performing isotope ­

Ganesh Sriram (ed.), Plant Metabolism: Methods and Protocols, Methods in Molecular Biology, vol. 1083, DOI 10.1007/978-1-62703-661-0_8, © Springer Science+Business Media New York 2014

109

110

Shilpa Nargund and Ganesh Sriram

l­abeling experiments (ILEs), which involve feeding nutrients ­containing mixtures of isotopes of certain elements (e.g., 13C and 12 C for carbon, 14N and 15N for nitrogen, 17O and 18O for oxygen) to a cell or tissue culture, and then measuring isotopomer (isotope isomer) abundances of intracellular metabolites and biomass components (e.g., proteinogenic amino acids, lipogenic fatty acids, sugars, nucleotides, and small metabolites) by using nuclear magnetic resonance (NMR) or mass spectrometry (MS) [4, 5]. Various isotopomers of a given metabolite may be generated in an ILE, depending on the atom rearrangements of reactions in metabolic pathways that contributed to the synthesis of the ­ metabolite. The relative abundances of these isotopomers depend on the relative fluxes through these pathways. Therefore, fluxes can be estimated by interpreting the set of measured isotopomer abundances (Iexp) with the aid of a steady-state mathematical model of metabolism that relates fluxes to isotopomer abundances. Mathematically, flux evaluation requires using such a model to simulate a set of isotopomer abundances (Isim) from guessed flux values and successively refining the guessed fluxes until Isim converges with Iexp [4, 6]. In this regard, it is convenient that metabolic pathways and their atom rearrangements are fairly well documented in metabolic databases such as MetaCyc ([7]; also Chapter 10) and KEGG [8]. Despite the availability of this information, the interpretation and analysis of isotope labeling patterns from an ILE is further complicated by (1) reversible reactions and cyclic pathways, which can cause the label to back-mix (e.g., [9]), (2) natural abundances of some isotopes (e.g., 1.1 % of all carbon in nature is 13C), and (3) compartmentalization of metabolism, especially in plant cells [10–13]. These complications cause hundreds or even thousands of isotopomers to arise in even relatively simple metabolic networks. This necessitates computational methods for modeling and simulation of ILEs. This chapter describes the steps involved in (1) building a mathematical model of the metabolic and isotopomer network, (2) simulating an ILE, and (3) using computer programs and optimization algorithms to estimate fluxes from Iexp. Quite a few computer programs and software toward these purposes are now available, including 13CFLUX2 [14], Metran [15, 16], NMR2Flux+ [4, 17], FiatFlux [18], and OpenFLUX [19]. These programs offer a variety of capabilities (Subheading 5.2 provides an overview). This chapter does not focus on the use of a specific flux evaluation program, but presents and discusses the underlying (mathematical) concepts to enable the reader to use any program to estimate fluxes from isotope labeling data. Chapter 9 presents an application of the OpenFLUX platform.

Modeling of Isotope Labeling Experiments

111

2  Materials 1. Computers with operating systems that enable use of C, C++, or MATLAB. See Note 1. 2. Flux evaluation software such as 13CFLUX2 (available from http://www.13cflux.net/13cflux2), Metran (available from http://www.che.udel.edu/mranton/metran.html), NMR2Flux+ (available from the corresponding author on request), FiatFlux (available from http://www.imsb.ethz.ch/ researchgroup/nzamboni/research/Software/fiatflux), and OpenFLUX (available from http://openflux.sourceforge. net). Licenses for some these software may be free of charge for academic use, but may require payment for commercial use. Additionally, their use may require agreements to be signed by end users.

3  Methods 3.1  Steady-State Mathematical Modeling of Plant Metabolism

A metabolic network is an interconnected web of several metabolic reactions and pathways. Each reaction in such a network interconverts metabolites in stoichiometric proportions and according to particular atom rearrangements. Building a model of a metabolic network, especially in the context of an ILE, entails the following steps.

3.1.1  Cataloging Reactions, Atom Rearrangements, and Reversibilities

To model an ILE, a user should first assemble a catalog of all major reactions that are expected to process more than ~1 % of carbon or the element of interest. Plant biochemistry textbooks (e.g., [20, 21]) provide an overview of plant metabolic pathways and the underlying carbon atom rearrangements, whereas review articles (e.g., [13, 22, 23] and several others) shed light on specific pathways and types of reactions. Freely available databases including MetaCyc (which encompasses AraCyc, PoplarCyc, and various others) ([7, 24]; also see Chapter 10) and Kyoto Encyclopedia of Genes and Genomes (KEGG) [8] provide substantial, current information on metabolic reactions in an organism-specific manner as well as molecular structures for metabolites and carbon atom rearrangements for reactions. Reversible reactions are consequential in the modeling of ILEs, as they can cause the label to back-mix and form unique isotopic patterns. Although textbooks, reviews, and pathway databases provide some information about whether a reaction is reversible, this information should strictly be obtained from the change in Gibbs energy (ΔG) accompanying the reaction at physiological conditions. Irreversible (spontaneous) reactions

112

Shilpa Nargund and Ganesh Sriram

have s­ignificantly negative ΔG values (≤~−3.0 kcal/mol) at physiological conditions. In comparison, reversible reactions have ΔG values relatively closer to 0 kcal/mol. Databases including MetaCyc provide standard Gibbs energies of formation (ΔGf) of compounds, from which the standard Gibbs energy changes (ΔG°) can be calculated (e.g., [25]). See Note 2. Additionally, the catalog of reactions should include extracellular fluxes including nutrient uptake from and product secretion into extracellular space as well as the effluxes of metabolites toward the biomass components including proteinogenic amino acids, lipogenic fatty acids and glycerol, nucleotides, starch, cell wall, and other carbohydrates (e.g., [26]). 3.1.2  Modeling Compartmentalization of Reactions in Subcellular Organelles

A unique feature of plant metabolism is its extensive compartmentalization, with significant shuttling of metabolites between compartments [10, 11, 23]. Not only are metabolic reactions and pathways distributed between several subcellular compartments, but there is also some degree of duplication between the organelles. For example, plant cells operate glycolysis and the pentose phosphate pathway (PPP) both in the cytosol and the plastid [13]; there is also evidence of a mitochondrial glycolysis [27]. The biosynthesis pathways of some amino acids can be exclusively plastidic (e.g., histidine, branched chain amino acids, and aromatic amino acids), whereas other amino acids can be synthesized in the plastid, cytosol, and mitochondrion [21]. Isoprenoid synthesis can occur via a cytosolic mevalonate pathway or a plastidic methyl erythrose pathway, both of which begin with intermediates of glycolysis in the respective compartments [28, 29]. Evidence on pathway compartmentalization is available, although not comprehensively, from the literature and metabolic pathway databases. Therefore, a hybrid approach consisting of scanning the literature, querying metabolic databases and using bioinformatic tools is necessary to obtain an accurate compartmentalized model. Chapter 12 discusses this approach in greater detail. Additionally, it may be necessary to reassess the compartmentalization in an initial model if it cannot satisfactorily account for measured isotopomer abundances. Pools of the same metabolite in different subcellular compartments exchange by means of transport fluxes that are often mediated by intercompartmental transporter proteins [23, 30, 31]. Such pools may or may not be in isotopic equilibrium with each other. Therefore, the metabolic model should incorporate not only the fluxes through the replicated pathways in different compartments but also the fluxes corresponding to transport of metabolites across organellar membranes. Naturally, the model should include separate compartmental pools of the metabolites of a replicated pathway (see Subheading 4 for an example).

Modeling of Isotope Labeling Experiments

113

3.1.3  Simplifying the Model

Although plant metabolism typically includes thousands of ­reactions and there are some excellent kinetic studies of secondary metabolism [32, 33], current isotope MFA techniques are able to focus on central and to some extent, intermediary carbon metabolism. This is mainly because (1) branched and cyclic pathways, which isotope MFA can elucidate, are more common in central carbon metabolism; (2) isotopomers of the products of central carbon metabolism are easily measured (see Chapters 6 and 7); and (3) large isotopomer models that include primary and secondary metabolism may involve substantial computational load (demonstrated in Subheading 3.2), although articles describing such attempts for bacterial metabolism have begun appearing [34]. Therefore, users may wish to, at least initially, simplify their models so that they include only central or intermediary carbon metabolism. Additionally, for steady-state isotope MFA, models can be simplified without loss of information by lumping (consolidating) together sequential reactions that lie in an unbranched pathway and that do not involve carbon atom rearrangements (see Subheading 4 for an example and Note 3).

3.1.4  Writing the Mathematical Model

On the most basic level, a mathematical model of an ILE consists of mass balance equations for metabolites and isotopomers. For any metabolite:

{v

in

}

+ v generation + {vout + vconsumption } = vaccumulation



(1)

where vin is the sum of all fluxes that transport the metabolite into the network, vgeneration is the sum of all fluxes that generate the metabolite through chemical reaction, vout is the sum of all fluxes that transport the metabolite out of the network, and vconsumption is the sum of all fluxes that consume the metabolite through chemical reaction. At metabolic steady state, there is no accumulation of intracellular metabolites in the network. Therefore

{v

in

} {

}

+ v generation + vout + vconsumption = 0



(2)

Figure  1 depicts an illustrative metabolic subnetwork. This subnetwork has two extracellular fluxes vin (entry of metabolite A) and vout (exit of metabolite E) that can potentially be measured directly, as well as five intracellular fluxes v1, v2, v3, vb, and vd that cannot be measured directly. The metabolite balance equations for this subnetwork are (3) v − v = 0  (a) in

1

(b)

vb − v1 = 0 

(4)

(c)

v1 + v3 − v 2 = 0 

(5)

114

Shilpa Nargund and Ganesh Sriram

...

cell

vb

B vin

A

v1

C

v2

E

vout

v3 D vd

…

vin 1 0 0 0 1 0

1 1

0 0 0 0 0 1

1 0

0 0 0

0

0

0

0

vb

0

0

0

vd

0 0

1 0 1

1 1 0

0

v1

0

0

v2

0

1 v3

0

vout

Fig. 1 Metabolic fluxes in a simple subnetwork. Arrows represent fluxes; uppercase letters represent metabolites. Substrate A enters the cell and E leaves the cell; the corresponding extracellular fluxes vin and vout can usually be measured directly. The ellipses and dashed lines indicate that B and D are downstream of an unknown part of the network, whose fluxes toward B and D are lumped into two fluxes vb and vd. The fluxes v1, v2, v3, vb, and vd are intracellular fluxes. In networks involving parallel intracellular pathways, loops or cycles, such fluxes can be estimated by fitting isotopomer data to a mathematical model of the metabolic network

(d)

v d − v3 = 0 

(e)

v 2 − vout = 0 

(6) (7)

These equations can be represented by a matrix equation: 1 0  0  0 0

vinn  0 0 −1 0 0 0  vb  1 0 −1 0 0 0  v d    0 0 1 −1 1 0  v1  =  0 1 0 0 −1 0  v 2    0 0 0 1 0 −1 v3  v   out 

0  0    0    0  0

(8)



or

S⋅v = 0

(9)

Modeling of Isotope Labeling Experiments

115

where S is a “stoichiometric” matrix and v is a vector that consolidates all fluxes. 3.1.5  Choosing Free Fluxes

The number of fluxes in this subnetwork (seven: vin, vout, v1, v2, v3, vb, vd) is two less than the number of metabolite balances (five: A, B, C, D, E). This difference is typically equal to the number of “degrees of freedom” of the network. For a network having two degrees of freedom, such as this one, two fluxes will need to be measured or chosen to determine the values of all fluxes in the network. After considering the measured (extracellular) fluxes, the intracellular fluxes that will need to be chosen are called “free” fluxes. Once certain fluxes are measured and the correct number of free fluxes is chosen, the remaining fluxes in the network (dependent fluxes) can be expressed as linear combinations of the measured and free fluxes. Let us suppose that in the subnetwork of Fig. 1, either vin or vout, but not both, can be directly measured. This leaves one degree of freedom, i.e., one flux needs to be chosen to completely determine the values of all fluxes. In such simple networks with a small number of branches, the free flux choices are obvious. For example, if vin but not vout is measured, then v3, vd, or v2 are valid free flux choices, but v1 is not (the measurement of v1 does not provide additional information since the stoichiometry makes v1 identical to vin). Conversely, if vout but not vin is measured, then v1, vb, v2, or vd are valid free flux choices, but v2 is not (since the stoichiometry makes v2 identical to vout). However, in larger metabolic networks, the free flux choices are often not so trivial. Free flux choice then needs to be performed by using a systematic, mathematical method. For a metabolic network with m metabolites and n net fluxes, the number of free fluxes f is given by f = n − rank(S). This is as long as no metabolite balances are linearly dependent. This condition is usually satisfied in medium-sized or larger networks. Alternatively, f = n − m − p, where p is the number of measured fluxes. Several, often an enormous number of, sets of free fluxes may exist. A valid set of free fluxes can be determined by expressing Eq. 9 as

Sm ⋅ v m + Sx ⋅ v x = 0

(10)

where v = [vm|vx] is a partitioning of v based on whether fluxes are measured/free (vm) or unknown (vx), and S = [Sm|Sx] is a corresponding partitioning of S. This gives

( )

v x = S x−1 ⋅ S m ⋅ v m



(11)

For a valid v = [vm|vx] partitioning, the matrix Sx in Eq. 11 has to be invertible. Therefore, one method of determining free fluxes in a metabolic network is

116

Shilpa Nargund and Ganesh Sriram

1. Combinatorially generate various free flux sets. 2. Examine each set to see if it results in an invertible Sx. There are several alternative methods to determine free flux sets (e.g., [19]), some of which rely on the technique of singular value decomposition [35]. 3.2  Simulating an ILE in the Metabolic Network Model

For many realistic metabolic networks, the only method to estimate the values of the free fluxes (and ultimately all the fluxes) is to guess and iteratively refine a set of free fluxes that best account for the isotopomer abundances Iexp obtained by conducting an ILE [6]. Iteration is necessary because the isotopomer balance equations are implicit in the unknown fluxes; therefore, they cannot be inverted trivially to express the unknown fluxes as explicit functions of the measured isotopomer abundances. However, the methods described below always afford a “forward” simulation of the ILE, which involves calculation of a set of isotopomer abundances Isim from a known or guessed set of fluxes. Flux estimation algorithms therefore use a feasible, initial flux guess and employ global optimization to successively refine previous flux guesses until Iexp and Isim converge. Simulating the ILE, a critical step of the flux estimation process, requires both the stoichiometries and the carbon atom rearrangements of reactions in the model. These can be combined into isotopomer balances. Figure 2a depicts carbon atom rearrangements and Fig. 2b illustrates isotopomer balances for metabolite C in the example network of Fig. 1. Solving these balances by using given or guessed values of the free fluxes provides Isim of the metabolite C (and also of metabolite E, whose steady-state isotopomer distribution is identical to that of metabolite C). Given a measured isotopomer distribution Iexp for E, the free fluxes can be guessed or adjusted within the stoichiometrically feasible flux space to obtain an Isim that best matches Iexp. A metabolite with n carbon atoms has 2n isotopomers and correspondingly 2n isotopomer balance equations. These equations are bilinear for reactions involving condensation of two reactant metabolites into a product. For instance, in Fig. 2b, the first term involving the flux v1 renders the equations bilinear. In steady-state scenarios, lumping of metabolites can reduce the number of isotopomer balance equations. Researchers have developed several methods to solve the bilinear isotopomer balances analytically, including cumomer [36, 37] (see Note 6), bondomer [9, 38] (see Note 7), and elementary metabolite unit (EMU) [15] (see Note 8) balancing. All these methods transform the bilinear balance equations into cascades of linear equations that can be analytically solved in sequence. The cumomer and EMU balancing methods keep track of 13C and 12C atoms in metabolites, just as isotopomers do. These two methods are applicable to all ILEs and provide results identical to those obtained with isotopomer balancing. However, cumomer ­balancing

Modeling of Isotope Labeling Experiments

a

117

v1

A+B

v1

C

v2

E

v3

1 2

A

+

1 2 3 C

3

B

v2

1 2 3 E

v3

D3 2 1

D Flux balance on metabolite C v1 + v3 = v2

b

Isotopomer balances on metabolite C

C

: v1 (

A B

D

) = v2 (

C

.

) + v3 (

)

: v1 (

.

) + v3 (

) = v2 (

)

: v1 (

.

) + v3 (

) = v2 (

)

: v1 (

.

) + v3 (

) = v2 (

)

: v1 (

.

) + v3 (

) = v2 (

)

: v1 (

.

) + v3 (

) = v2 (

)

: v1 (

.

) + v3 (

) = v2 (

)

: v1 (

.

) + v3 (

) = v2 (

)

Fig. 2 Isotopomer balancing in the mathematical simulation of an ILE. The intracellular fluxes v1, v2, and v3 in the example network of Fig. 1 can be estimated by finding flux values that best account for the experimentally measured isotopomer abundances. For this, the ILE has to be mathematically simulated through isotopomer balancing. (a) Flux balance equation and carbon atom rearrangements of the reactions involving C. Carbon atom rearrangements that occur in the three reaction v1, v2, and v3 are depicted with respect to the metabolite C, with each circle representing a carbon atom. For instance, the first carbon atom of C can originate from the first carbon atom of A or third carbon atom of D. (b) Balances for the eight possible isotopomers of the metabolite C. (C is a three-­carbon metabolite and therefore has 23 isotopomers and 23 isotopomer balance equations.) Empty circles represent 12C; filled circles represent 13C. The isotopomer which is being balanced is listed at the left before its equation. Clearly, the isotopomer equations are nonlinear because of the condensation reaction v1

could be significantly faster and mathematically convenient (as it provides analytical instead of numerical solutions) than isotopomer balancing. EMU balancing, in turn, could be significantly faster than cumomer balancing for most metabolic networks. Bondomer balancing keeps track of intact and biosynthetic c­ arbon– carbon bonds in metabolites and therefore provides a different perspective. However, it is only applicable to ILEs employing a single, uniformly labeled carbon source with at least one carbon–carbon bond. Additionally, a recently reported method called fluxomer balancing decouples bilinear terms by using a single variable that consolidates both the fluxes and isotopomers [39].

118

Shilpa Nargund and Ganesh Sriram

3.3  Optimization and Statistical Methods for Flux Estimation

Because the isotopomer balances are implicit in the unknown fluxes, the fluxes have to be evaluated from the measured isotopomer abundances iteratively. Therefore, the fluxes are unknown parameters that minimize some error criterion, and the flux evaluation becomes a parameter estimation problem. Global optimization algorithms are essential for solving this problem by intelligently and efficiently searching the flux space iteratively to identify the flux sets that best account for the measured isotopomers Iexp. These algorithms usually minimize a χ2 metric that represents the difference between Iexp and Isim corresponding to the most recent flux guess. This metric is defined as:



N

Χ =∑ 2

j =1

(lexp

j

− lsim j σ 2j

)

2

(12)

where Iexpj and Isimj are the jth (of N) elements of Iexp and Isim, respectively, and σj is the measurement standard deviation of the jth element of Iexp. Global optimization algorithms are desirable in flux evaluation, and suitable algorithms include simulated annealing (used in [4]), sequential quadratic programming [6] and genetic algorithms. To capture all possible flux solutions that could account for the measured Iexp, it is necessary to perform the flux evaluation several hundred times, by beginning different iterations with flux guesses spread over the allowable flux space. Furthermore, it is necessary to obtain standard deviations or statistical confidence intervals for fluxes. This is a complex procedure due to the nonlinear relationship between fluxes and isotopomer abundances. Thus far, this has been performed by employing one of the methods described below. 3.3.1  Bootstrap Monte Carlo Simulation

This method is widely used in parameter estimation. In isotope MFA, it has been employed by Schmidt et al. [40] as well as by Shanks et al. ([4], see Supplementary Material IV). In this method (Fig. 3), a random number generator is used to first perturb the original measurement dataset (consisting of isotopomer abundances and standard deviations) to generate a large number ­(typically K ~500 or more) of synthetic datasets normally distributed around the original dataset. Then the flux evaluation routine uses optimization to independently evaluate fluxes from each synthetic dataset, thus resulting in K synthetic flux distributions that will typically be normally distributed around the flux distribution obtained from the original dataset. These synthetic flux distributions can be used to generate standard deviations or confidence intervals for the fluxes, as described in Press et al. [35]. The advantages of the Monte Carlo method are that it captures the nonlinearities in the flux–measurement relationship well,

Modeling of Isotope Labeling Experiments Iterative flux evaluation procedure; solves the inverse problem of M = f(v)

Random number generator

M = f(v)

measurements1

Iterative flux evaluation

fluxes1

measurements2

Iterative flux evaluation

fluxes2

…

…

measurementsK

Iterative flux evaluation

(evaluated) fluxes, v

fluxesK

{

{

(original) measurements, M

119

distribution of synthetic measurements

distribution of evaluated fluxes

statistics (SDs, conf. intervals) for fluxes

Fig. 3 Evaluation of statistical confidence intervals for fluxes

and it is precise and requires minimal familiarity with analytical statistics [35]. A major disadvantage is that it can be very time consuming, as the flux evaluation has to be repeated for each of the hundreds of synthetic datasets. 3.3.2  Linearization Method

This method was first used by Wiechert et al. [41]. Here, the function mapping isotopomers and fluxes is linearized, and the resulting linear relationship is used to transform the measurement standard deviations into flux standard deviations.

3.3.3  Method Suited to Metabolic Networks

A method introduced by Antoniewicz et al. [42] is particularly suitable for metabolic networks. In this method, an optimization routine first finds the minimum of a χ2 metric corresponding to the difference between measured and simulated isotopomer abundances (Eq.  12). Then, the confidence interval for a particular flux is determined by exploring outward from the χ2 minimum increasing the value of the flux of interest to a small amount, keeping this flux constant and varying all other fluxes until a new minimum is obtained. The flux of interest is then increased in a stepwise manner until (1) an optimum equal to the χ2 threshold for one degree of freedom is obtained (this threshold is 2.71 for a 90 % confidence interval and 3.84 for a 95 % confidence interval), or (2) the flux reaches its stoichiometric lower or upper bound.

120

Shilpa Nargund and Ganesh Sriram

4  E xample: Compartmentalized Model of Plant Glycolysis and Pentose Phosphate Pathways Figure  4 illustrates a compartmentalized model of the glycolysis and PPP network in a heterotrophic plant cell that uses glucose (Glc) as its sole carbon source. This model has 14 metabolites in two subcellular compartments—the cytosol and the plastid. In addition to cytosolic glucose (Glc) and CO2 (not assigned to a compartment since it is reasonable to assume that it readily diffuses between the cytosol and the plastid), the cytosolic metabolites are glucose-6-phosphate (G6P), fructose-6-phosphate (F6P), lumped triose phosphates (T3P), lumped pentose-5-phosphates (P5P), sedoheptulose-7-phosphate (S7P), and erythrose-4-phosphate (E4P). Their plastidic counterparts are G6PP, F6PP, T3PP, P5PP, S7PP, and E4PP. This model has 23 net fluxes: vInp, vpgif, vpfk, vgap, vg6pdh, vcellwall, vtktAf, vtalf, vtktBf, vCO2x, vg6pt, vp5pt, vt3pt, vg6pdhp, vpgifp, vpfkp, vtktAfp, vtalfp, vtktBfp, vgapp, vstarch, vhis, vphe. Some important steps of this network reconstruction are as follows. 1. Because this network features no carbon atom rearrangements between the three-carbon glycolytic metabolites (glyceraldehyde-3-­phosphate, dihydroxyacetone phosphate, 1,3-­bisphosphoglycerate, 3-phosphoglycerate, phosphoenolpyruvate, CYTOSOL

vg6pt

PLASTID

vp5pt

vphe

vCO2x

Glc

PPP

CO2

S7P

vInp

G6P vpgif

P5P

vg6pdh vcellwall

vtktAf

T3P

F6P

T3P vgap

vpgifp

PPP S7Pp

P5Pp

vg6pdhp vstarch

vtktBf

vt3pt

vpfkp Glycolysis

E4Pp

vtalfp

vtktAfp

T3Pp

F6Pp

F6Pp

F6P vpfk Glycolysis

G6Pp

vtalf

vhis

CO2

E4P

vtktBfp

T3Pp vgapp

Fig. 4 Model of compartmentalized glycolysis and pentose phosphate pathway (PPP) in plant cells. This figure is adapted from another publication authored by us [12], which appeared in print before this chapter. Arrows represent fluxes. Glycolysis and PPP duplicated in cytosol and plastid contain 13 metabolites and 23 net fluxes. The plastidic compartment is enclosed by a dashed box. Dashed arrows indicate metabolites leaving the system. Flux balances allow dependent fluxes (shaded in gray) to be expressed as linear combinations of measured and free fluxes. vInp and vgapp are measured fluxes. The fluxes vg6pdh, vtktAf, vg6pdhp, and vg6pt are free fluxes. Subscript “p” indicates metabolites in the plastidic compartment. E4P erythrose-4-phophate, F6P fructose-6-­phosphate, G6P glucose-6-phosphate, Glc glucose, P5P pentose phosphates, S7P sedoheptulose7-phosphate, T3P (lumped) triose phosphates

Modeling of Isotope Labeling Experiments

121

and pyruvate), they are lumped into a single metabolite T3P to simplify glycolysis. For the same reason, the rapidly equilibrating five-carbon PPP metabolites ribose 5-phosphate, ribulose 5-phosphate, and xylulose 5-phoshpate are lumped into a metabolite called pentose phosphate (P5P). 2. To capture the duplication of glycolysis and the PPP in the cytosol and the plastid, this network uses separate cytosolic and plastidic pools of the metabolites of these pathways (e.g., G6P and G6PP) as well as transport fluxes vg6pt, vt3pt, and vp5pt that enable exchange of these metabolite pools across the plastidic membrane. 3. Two fluxes in this network, vg6pdh and vg6pdhp—that correspond to decarboxylation reactions—are designated irreversible because of their Gibbs-free energies range between −0.79 and −61.92 kJ/mol, depending on pH. 4. The sink fluxes of this network are the evolution of CO2 (vCO2x), conversion of T3P (vgap), T3Pp (vgapp), P5Pp (vhis), and E4Pp (vphe) into amino acids or nucleotides, G6P into cell wall (vcellwall) and G6Pp into starch (vstarch). After writing the flux balance equations for the 12 intracellular metabolites that do not enter or leave the network (i.e., all metabolites except Glc and CO2), this network has 23 − 12 = 11 degrees of freedom. Seven extracellular fluxes including the glucose uptake flux (vInp) as well as the effluxes of T3P (vgap), T3Pp (vgapp), G6P (vcellwall), G6Pp (vstarch), P5Pp (vhis), and E4Pp (vphe) can be measured, leaving four degrees of freedom. Therefore, we designated the four fluxes vg6pdh, vtktAf, vg6pdhp, and vg6pt as the free fluxes to be estimated by conducting ILEs.

5  Discussion 5.1  Software Available for MFA

Because the mathematical operations described above can become quite complex for realistic networks, it is imperative to automate part or all of them, so as to obtain accurate results in reasonable time. Over the years, researchers have developed various software to simulate ILEs on metabolic networks, estimate fluxes using optimization algorithms, and find flux confidence intervals using statistical methods. These software include 13CFLUX (now replaced by its most recent version 13CFLUX2), FiatFlux, NMR2Flux, Metran, and OpenFLUX (see references in Subheading  1), although this is a non-exhaustive list. These software differ mainly in their ILE simulation technique (see discussion below), the optimization algorithm(s) available and the statistical method used for estimating flux standard deviations. Additionally, depending on the capabilities of a specific software, ancillary or home-built computational tools may be required to determine free fluxes and to process MS or NMR data.

122

Shilpa Nargund and Ganesh Sriram

FiatFlux, developed by Zamboni et al. [18] on MATLAB, ­ rovides flux ratios at metabolic branchpoints by solving probabilp ity equations that relate mass isotopomer distributions to relative fluxes at the branchpoints. FiatFlux is an open source software and is freely available for academic use. Metabolic networks of certain microorganisms are preconfigured in the publicly available version of the software, thus requiring no stoichiometric reaction or atom rearrangement information to be input by the user. This version works with select labels of glucose, e.g., 1-13C and uniformly (U)-13C and can only accept isotopomer information obtained from gas chromatography (GC)-MS analysis of TBDMS-­ derivatized amino acids. The software is capable of directly processing this raw GC-MS data, thereby demanding minimum user intervention. Additionally, it has a graphical user interface that makes it user-friendly. NMR2Flux+, developed by Sriram et al. [4, 17] and written in C, simulates ILEs by using cumomer or bondomer balancing. This program is freely available to academic researchers upon request. The user needs to input reaction stoichiometries and carbon atom rearrangements through a comma separated value (CSV) file, from which the program assembles a metabolic model. The user also needs to determine and specify the free flux choices a priori. NMR2Flux+ accepts GC-MS data, which should be preprocessed (MS data corrected for natural abundances of elements, NMR peaks integrated) and supplied in the form of CSV files. The program is capable of conducting sensitivity analyses, thus allowing users to design ILEs. It also facilitates simultaneous flux estimation from multiple parallel ILEs (see Subheading 5.3). To evaluate fluxes, NMR2Flux+ uses simulated annealing as its optimization algorithm and bootstrap Monte Carlo simulations to compute standard deviations of fluxes. Metran was developed by Antoniewicz et al. [15, 16] on MATLAB and is freely available for academic use. It employs EMU balancing to simulate ILEs. Users can describe reactions and atom rearrangements via a graphical user interface; Metran determines free fluxes based on this input. The software accepts raw MS and MS–MS data and corrects it before fitting it to the simulated data. Metran is capable of conducting sensitivity analysis and thus design of experiments. It is also capable of estimating fluxes simultaneously form multiple parallel ILEs. Its optimization algorithm is based on sequential quadratic programming. To estimate the standard deviations of fluxes, Metran uses an algorithm developed by Antoniewicz et al. [42] (Subheading 3.3.3). OpenFLUX is an open source software developed by Quek et al. [19] on MATLAB. It also employs EMU balancing to simulate ILEs. Metabolic reactions and atom rearrangements can be defined in a spreadsheet. The user can specify a subset of the free fluxes, and software determines the remaining free fluxes. OpenFLUX accepts

Modeling of Isotope Labeling Experiments

123

corrected MS data in spreadsheets. The FMINCON in MATLAB’s optimization toolbox serves as its optimization algorithm. To estimate standard deviations of fluxes, it uses the algorithm used in Metran [42] (Subheading 3.3.3). Wiechert and coworkers pioneered the development of MFA software with the C++-based program 13CFLUX [43]. As this book went to press, this program was replaced with the more advanced 13CFLUX2 [14], also developed in C++ with state-of-­ the-art capabilities and a graphical user interface. 13CFLUX2 is freely available for academic purposes. It simulates ILEs with cumomer or EMU balancing techniques. Reactions and carbon atom rearrangements can be specified in extensible markup language (XML) files, and the software uses this information for generating the stoichiometric matrix and cumomer or EMU balances. The user can either specify all or a few of the free flux choices, leaving the program to ascertain the rest. 13CFLUX2 is capable of handling NMR, MS, and MS–MS isotopomer measurements as well as generic measurements all of which have to be processed as described in Subheading 5.1. The user can choose between multiple optimization algorithms including sequential quadratic programming, nonlinear programming or certain commercially available optimization libraries. The program allows a priori design of ILEs and simultaneous flux estimation from multiple parallel ILEs. Additionally, it is integrated with a flux visualization program Omix and with MATLAB. It is also compatible for use on cluster computing. 5.2  Preparing MS and NMR Data for Flux Evaluation

Raw MS and NMR data need to be processed to obtain isotopomer abundances. Particularly, MS data has to be corrected for the presence of naturally abundant isotopes other than metabolic carbon in the metabolites and the derivatizing reagents [44]. NMR peaks and peaklets need to be integrated to obtain isotopomer abundances [4, 45]. This is discussed in detail in another chapter by us in this series [12]. Often, both MS and NMR are capable of measuring only linear combinations of isotopomer abundances of any given metabolite (Fig. 5), although it is sometimes possible to obtain all 2n isotopomer abundances for certain metabolites. For instance, the two-carbon amino acid Glycine (Gly) derivatized with t-butyldimethylsilyl (TBDMS) tags [46] forms two distinct fragments on a single-quadrupole MS: 1. Gly[2], which contains carbon atom C-2 and has the mass distribution [m + 0, m + 1]; 2. Gly[12], which contains carbon atoms 1 and 2 and has the mass distribution [m + 0, m + 1, m + 2]; Together, these fragments provide enough information to unambiguously determine the abundances of all 22 = 4 isotopomers of Gly. Contrastingly, the three-carbon amino acid alanine (Ala)

124

Shilpa Nargund and Ganesh Sriram

a

NMR HSQC peak for Asp a

b

Mass isotopomers of Ala 2 3

singlet =

+

m+0 =

+

doublet 1 =

+

m+1 =

+

doublet 2 =

+

m+2 =

+

double doublet =

+

c

+

+

Mass isotopomers of Ala 1 2 3 m+0 = m+1 =

+

+

m+2 =

+

+

m+3 =

Fig. 5 Linear combinations of isotopomer abundances are measureable by NMR and MS. Empty circles ­represent 12C; filled circles represent 13C. (a) The 2-D NMR peak corresponding to the Asp α carbon atom can split into four peaklets (a singlet, two doublets and a double doublet), the relative areas under which are proportional to the relative abundances of the linear combinations of the four illustrated isotopomers. A singlet is obtained when the α carbon atom is 13C and its immediate neighbors are both 12C. Doublet 1 and doublet 2 are obtained when the α carbon atom and only one of its immediate neighbors are 13C. The double doublet is obtained when the α carbon and both its neighbors are 13C. The isotope label on the γ position does not affect the peak splitting because of the absence of long-range coupling between it and the α carbon atom. (b, c) Mass spectrometry fragments metabolites and measures their masses. If m is the molecular weight of a fragment with n carbon atoms then the mass of its isotopomers can range from m + 0 to m + n depending on the isotope labeling of the carbon atoms. These mass isotopomers represent linear combinations of positional isotopomers

derivatized with t-butyldimethylsilyl (TBDMS) tags also forms two fragments (Ala[23] and Ala[123]); however, these fragments do not provide enough information to unambiguously determine the abundances of all 23 = 8 Ala isotopomers. Nevertheless, ILE simulations can make use of the partial information in molecules such as Ala. Newer MS techniques such as tandem MS–MS [47] hold promise in providing information about substantial numbers of isotopomers of metabolites (see Note 5). 5.3  Improving Flux Estimates from MFA

Several methods, including those listed below, enable an experimenter to improve the identifiability of fluxes obtainable from an ILE. 1. Measure isotopomer abundances of several biomass components. Proteinogenic amino acids, the most widely measured biomass components, are synthesized via several metabolic pathways that collectively span most of primary carbon metabolism. However, measuring isotopomer abundances of additional biomass components such as carbohydrates including sugars

Modeling of Isotope Labeling Experiments

125

and starch, soluble metabolites including organic acids, lipids, and nucleotides will unveil flux information on several more pathways thus affording a network-wide enhancement of the accuracy and confidence of flux estimation. 2. Iteratively refine metabolic network model to improve fit. Often, an initially chosen network model has to either be simplified or augmented on the basis of preliminary (dis)agreements between measured isotopomer abundances and those simulated assuming the network. For example, Shanks and coworkers [48] observed that in Catharanthus roseus hairy roots fed with U-13C sucrose, the isotopomer abundances of leucine (Leu) were not accounted for (fitted) by a model that included conventional Leu synthesis from its precursors acetyl-CoA and pyruvate. The augmentation of a Leu degradation and cycling pathway improves the fit. In such case, the labeling patterns in the isotopomers themselves hold clues on the refinements required in the network. In the example above, the isotopomer abundances around Leu γ and δ1 were least accounted for by the initially assumed model, suggesting a Leu degradation and cycling pathway that separate these particular atoms of Leu. 3. Design isotopic labels to improve flux identifiability. A priori design of ILEs is possible using established statistical methods [36, 49, 50]. Such design identifies isotope labels and also particular isotopomer abundance measurements that are crucial in improving flux identifiabilities. 4. Simultaneous flux estimation from multiple parallel labeling experiments. This technique, wherein the labeling data from two or more ILEs with different isotopic labels is fit to a model simultaneously, improves accuracy of flux estimates by increasing the available labeling information (see Note 4). A graphical summary of the elements required for mathematically modeling ILEs (Fig. 6) provides a checklist covering information and tools required for flux evaluation.

6  Notes 1. Each flux evaluation program may have different hardware and operating system requirements. Users should consult the supplier for more details. 2. The Gibbs energy change of a reaction at physiological conditions (ΔG) is often different in value, and sometimes even in sign, from the Gibbs energy change of the same reaction ΔG° at standard conditions (25 °C, 1 atm, water at 55.6 M, pH 7.0, all other reactants and products at 1.0 M). For the illustrative reaction:

126

Shilpa Nargund and Ganesh Sriram

Information needed

Visual representation

Reaction network, reversibilities, transport and biomass accumulation reactions

v1

A+B

v2

C

Mathematical equation

Where to find information KEGG, MetaCyc, PlantCyc, AraCyc, PoplarCyc, literature

S.v = 0

E

v3

D Atom rearrangements

KEGG, MetaCyc, PlantCyc, AraCyc, PoplarCyc, literature

v1 1 2

A

+

3

1 2 3

B

C

v2

1 2 3

E

v3 D 3 2 1

Analytical solution to isotopomer balance equations

Linear combinations of isotopomer abundances of metabolites measured by MS or NMR

:

v1 (

.

) + v3 (

) = v2 (

)

:

v1 (

.

) + v3 (

) = v2 (

)

:

v1 (

.

) + v3 (

) = v2 (

)

References: cumomer [37, 51], bondomer [9, 38], EMU [15]

Several mathematical techniques

References: [4, 12, 44, 45]

Mass isotopomers of Ala 2 3

m+0 =

+

m+1 =

+

m+2 =

+

+

+

Fig. 6 Graphical overview of mathematical modeling concepts

kb



ν A A + νBB  νCC + νD D kf



ΔG and ΔG° are related as:



ìï C vC C vD DG = DG° + RT ln í Cv Dv ïî CAA CBB

üï ý ïþ

(13)

Thus, knowledge of cellular conditions, including metabolite concentrations and organellar pH (either measured directly or estimated from the literature) can enable calculation of ΔG from ΔG°. 3. While modeling and simulating ILEs toward estimating carbon fluxes, it is advisable to minimize the number of reactions in which carbon atom rearrangements do not occur. An example of a reaction in this category is transamination, which shuffles keto and amine groups between reactants but leave the carbon backbone intact. Thus, such a reaction will not influence carbon labeling patterns of metabolites.

Modeling of Isotope Labeling Experiments

127

4. While processing isotopomer data and extracellular measurements from parallel ILEs that only differ in their isotope label, users may find that the extracellular flux measurements may exhibit errors that are significantly larger than errors in isotopomer abundance measurements from replicate ILEs. To account for these, correspondingly large standard deviations on these measurements should be included. 5. Isotopomer abundances, especially those measured by MS, tend to have very small measurement errors. We frequently see errors of 0.005 or less, for isotopomer abundances in the range [0, 1], even across biological replicates. When using χ2 minimization to evaluate fluxes, isotopomer abundances with very small errors may be unduly favored due to the nature of the χ2 metric (see Eq. 12), and the resultant flux distribution could be biased toward these measurements. Furthermore, an error of 0.005 is close to the least count of a (GC-)MS. Therefore, it is advisable to artificially increase measurement errors below the least count of the instrument to the least count. 6. Cumomer (cumulative isotopomer) balancing, introduced by Wiechert et al. [37, 51], is a methodology wherein the bilinear isotopomer balance equations are transformed to cascades of linear equations. Cumomers are defined as sums of certain isotopomers. For example, the cumomers of a three-carbon metabolite C are represented as cxxx, c1xx, cx1x, cxx1, c11x, cx11, c1x1, and c111, where “1” indicates an atom that is definitively 13 C-labeled and “x” indicates an atom whose labeling state is undetermined. Contrastingly, the isotopomers of the same metabolite are represented as c000, c100, c010, c001, c110, c011, c101, and c111, where “1” indicates a definitively 13C-labeled atom and “0” indicates a definitively 12C-labeled atom. Thus, a ­metabolite with n carbon atoms has 2n cumomers. Cumomer abundances can be calculated from isotopomer abundances using [51]:



c xxx =

1

∑c

i , j ,k = 0

ijk

=1

(14)

Cumomer balances can be written analogously to isotopomer balances. The weight of a cumomer is defined to be equal to the number of definitively labeled carbon atoms in the cumomer. The bilinearities that occur in isotopomer balances due to ­condensation reactions are addressed by solving the cumomer balances in order of weight [51], thereby reducing the coupling between the isotopomer abundances. Given a set of fluxes in a metabolic network, all 2n cumomer abundances (and thereby isotopomer abundances) of all metabolites can be simulated analytically via cumomer balancing.

128

Shilpa Nargund and Ganesh Sriram

7. Van Winden et al. [38] and Sriram et al. [9] introduced and developed bondomer balancing, which is useful for modeling ILEs that use only U-13C and naturally abundant versions of a single carbon source. Bondomers (bond isomers) are defined as molecules of a metabolite that differ in the origin of the bonds between their carbon atoms. Bonds are distinguished as intact (the atoms connected by the bond were derived from the same molecule of the carbon source and remained connected through metabolism) or biosynthetic (the atoms connected by the bond were obtained from different molecules of the carbon sources and were biosynthetically reassembled). Probability equations (derived from, e.g., [52]) can be used to interconvert bondomers and isotopomers. Bondomer balances factor the natural abundance of 13C. Therefore, networks to which bondomer balancing can be applied have considerably fewer bondomer balances than isotopomer balances [9]. 8. Elementary metabolite unit (EMU) balancing, developed by Antoniewicz et al. [15], overcomes bilinearities in isotopomer balances through a new definition, the EMU. An EMU is defined as a fragment consisting of a subset of the atoms of a metabolite [15]. The atoms in an EMU do not have to be contiguous or be connected by bonds. A metabolite with n atoms has 2n − 1 EMUs, only a fraction of which are required in simulating an ILE. EMU balancing involves a bottom-up approach wherein only the EMUs necessary for simulating the abundances of given metabolites are invoked. The size of an EMU is defined as the number of atoms present in the EMU. The nonlinearities in isotopomer balances are decoupled by solving EMU balances in the order of their sizes. The EMU technique is especially advantageous when several isotopes are being tracked (13C, 15N, and 17O).

Author Contributions SN and GS conceived the chapter. SN wrote a draft of the manuscript and prepared most of the figures. GS wrote sections of the manuscript, prepared Fig. 3, critically edited the manuscript and prepared the final version. Both authors approved the final version of the manuscript.

Acknowledgments  This work was funded by the National Science Foundation (award number IOS 0922650).

Modeling of Isotope Labeling Experiments

129

References 1. Stephanopoulos G (2002) Metabolic engi- 13. Kruger NJ, von Schaewen A (2003) The oxidative pentose phosphate pathway: structure neering: perspective of a chemical engineer. and organisation. Curr Opin Plant Biol AIChE J 48:920–926 6:236–246 2. Stephanopoulos G (1999) Metabolic fluxes and metabolic engineering. Metab Eng 14. Weitzel M, Nöh K, Dalman T, Niedenführ S, Stute B, Wiechert W (2012) 13CFLUX2 – 1:1–11 high-performance software suite for 3. Giegé P, Heazlewood JL, Roessner-Tunali U, 13C-metabolic flux analysis. Bioinformatics Millar AH, Fernie AR, Leaver CJ, Sweetlove 29:143–145. doi:10.1093/bioinformatics/ LJ (2003) Enzymes of glycolysis are functionbts646 ally associated with the mitochondrion in 15. Antoniewicz MR, Kelleher JK, Stephanopoulos Arabidopsis cells. Plant Cell 15:2140–2151 G (2007) Elementary metabolite units 4. Sriram G, Fulton DB, Iyer VV, Peterson JM, (EMU): a novel framework for modeling isoZhou R, Westgate ME, Spalding MH, Shanks topic distributions. Metab Eng 9:68–86 JV (2004) Quantification of compartmented metabolic fluxes in developing soybean 16. Yoo H, Antoniewicz MR, Stephanopoulos G, Kelleher JK (2008) Quantifying reductive carembryos by employing biosynthetically boxylation flux of glutamine to lipid in a brown directed fractional 13C labeling, two-­ adipocyte cell line. J Biol Chem 283: dimensional [13C, 1H] nuclear magnetic reso20621–20627 nance, and comprehensive isotopomer balancing. Plant Physiol 136:3043–3057 17. Sriram G, Rahib L, He J-S, Campos AE, Parr LS, Liao JC, Dipple KM (2008) Global meta 5. Ratcliffe RG, Shachar-Hill Y (2006) Measuring bolic effects of glycerol kinase overexpression multiple fluxes through plant metabolic netin rat hepatoma cells. Mol Genet Metab works. Plant J 45:490–511 93:145–159 6. Wiechert W (2001) 13C metabolic flux analy 18. Zamboni N, Fischer E, Sauer U (2005) sis. Metab Eng 3:195–206 FiatFlux – a software for metabolic flux analysis 7. Caspi R, Altman T, Dreher K, Fulcher CA, from 13C-glucose experiments. BMC Subhraveti P, Keseler IM, Kothari A, Bioinformatics 6:209 Krummenacker M, Latendresse M, Mueller LA et al (2012) The MetaCyc database of met- 19. Quek L-E, Wittmann C, Nielsen L, Kromer J (2009) OpenFLUX: efficient modelling softabolic pathways and enzymes and the BioCyc ware for 13C-based metabolic flux analysis. collection of pathway/genome databases. Microb Cell Factor 8:25 Nucleic Acids Res 40:D742–D753. doi:10.1093/nar/gkr1014 20. Heldt H-W, Piechulla B (2010) Plant biochemistry. Academic, New York, NY 8. Kanehisa M, Goto S, Kawashima S, Nakaya A (2002) The KEGG databases at GenomeNet. 21. Singh BK (1998) Plant amino acids (books in Nucleic Acids Res 30:42–46 soils, plants, & the environment). CRC, Boca Raton, FL 9. Sriram G, Shanks JV (2004) Improvements in metabolic flux analysis using carbon bond 22. Herrmann KM, Weaver LM (1999) The shikilabeling experiments: bondomer balancing mate pathway. Annu Rev Plant Physiol Plant and Boolean function mapping. Metab Eng Mol Biol 50:473–503 6:116–132 23. Linka N, Weber APM (2010) Intracellular 10. Lunn JE (2007) Compartmentation in plant metabolite transporters in plants. Mol Plant metabolism. J Exp Bot 58:35–47 3:21–53 11. Shanks JV (2005) Phytochemical engineering: 24. Zhang P, Foerster H, Tissier CP, Mueller L, combining chemical reaction engineering with Paley S, Karp PD, Rhee SY (2005) MetaCyc plant science. AIChE J 51:2–7 and AraCyc. Metabolic pathway databases for plant research. Plant Physiol 138:27–37 12. Nargund S, Sriram G (2013) Designer labels for plant metabolism: statistical design of iso- 25. Smith JM, Ness HCV, Abbott MM (2005) tope labeling experiments for improved quanIntroduction to chemical engineering thermotification of flux in complex plant metabolic dynamics, 7th edn. McGraw Hill Higher networks. Mol Biosyst 9:99–112 Education, New York, NY

130

Shilpa Nargund and Ganesh Sriram

26. Sriram G, Gonzalez-Rivera O, Shanks JV (2006) Determination of biomass composition of Catharanthus roseus hairy roots for metabolic flux analysis. Biotechnol Prog 22: 1659–1663 27. Graham JWA, Williams TCR, Morgan M, Fernie AR, Ratcliffe RG, Sweetlove LJ (2007) Glycolytic enzymes associate dynamically with mitochondria in response to respiratory demand and support substrate channeling. Plant Cell 19:3723–3738 28. Lohr M, Schwender J, Polle JEW (2012) Isoprenoid biosynthesis in eukaryotic phototrophs: a spotlight on algae. Plant Sci 185–186:9–22 29. Schwender J, Seemann M, Lichtenthaler H, Rohmer M (1996) Biosynthesis of isoprenoids (carotenoids, sterols, prenyl side-chains of chlorophylls and plastoquinone) via a novel pyruvate/glyceraldehyde 3-phosphate non-­ mevalonate pathway in the green alga Scenedesmus obliquus. Biochem J 316:73–80 30. Knappe S, Flugge UI, Fischer K (2003) Analysis of the plastidic phosphate translocator gene family in Arabidopsis and identification of new phosphate translocator-homologous transporters, classified by their putative substrate-­ binding site. Plant Physiol 131: 1178–1190 31. Weber APM, Linka N (2011) Connecting the plastid: transporters of the plastid envelope and their role in linking plastidial with cytosolic metabolism. Annu Rev Plant Biol 62: 53–77 32. Boatright J, Negre F, Chen X, Kish CM, Wood B, Peel G, Orlova I, Gang D, Rhodes D, Dudareva N (2004) Understanding in vivo benzenoid metabolism in petunia petal tissue. Plant Physiol 135:1993–2011 33. Colón AM, Sengupta N, Rhodes D, Dudareva N, Morgan J (2010) A kinetic model describes metabolic response to perturbations and distribution of flux control in the benzenoid network of Petunia hybrida. Plant J 62:64–76 34. Ravikirthi P, Suthers PF, Maranas CD (2011) Construction of an E. coli genome-scale atom mapping model for MFA calculations. Biotechnol Bioeng 108:1372–1382 35. Press WH, Teukolsky SA, Vetterling WT, Flannery BP (2007) Numerical recipes. The art of scientific computing, 3rd edn. Cambridge University Press, Cambridge 36. Möllney M, Wiechert W, Kownatzki D, de Graaf AA (1999) Bidirectional reaction steps in metabolic networks: IV. Optimal design of isotopomer labeling experiments. Biotechnol Bioeng 66:86–103

37. Wiechert W, Wurzel M (2001) Metabolic isotopomer labeling systems: Part I: global dynamic behavior. Math Biosci 169:173–205 38. Van Winden WA, Heijnen JJ, Verheijen PJT (2002) Cumulative bondomers: a new concept in flux analysis from 2D [13C, 1H] COSY NMR data. Biotechnol Bioeng 80:731–745 39. Srour O, Young J, Eldar Y (2011) Fluxomers: a new approach for 13C metabolic flux analysis. BMC Syst Biol 5:129 40. Schmidt K, Nørregaard LC, Pedersen B, Meissner A, Duus JO, Nielsen JO, Villadsen J (1999) Quantification of intracellular metabolic fluxes from fractional enrichment and 13C-13C coupling constraints on the isotopomer distribution in labeled biomass components. Metab Eng 1:166–179 41. Wiechert W, de Graaf AA (1997) Bidirectional reaction steps in metabolic networks: I. Modeling and simulation of carbon isotope labeling experiments. Biotechnol Bioeng 55: 101–117 42. Antoniewicz MR, Kelleher JK, Stephanopoulos G (2006) Determination of confidence intervals of metabolic fluxes estimated from stable isotope measurements. Metab Eng 8: 324–337 43. Wiechert W, Mollney M, Petersen S, de Graaf AA (2001) A universal framework for 13C metabolic flux analysis. Metab Eng 3: 265–283 44. Winden WAV, Wittmann C, Heinzle E, Heijnen JJ (2002) Correcting mass isotopomer distributions for naturally occurring isotopes. Biotechnol Bioeng 80:477–479 45. Van Winden W, Schipper D, Verheijen P, Heijnen J (2001) Innovations in generation and analysis of 2D [13C, 1H] COSY NMR spectra for metabolic flux analysis purposes. Metab Eng 3:322–343 46. Klapa MI, Aon JC, Stephanopoulos G (2003) Ion-trap mass spectrometry used in combination with gas chromatography for high-­ resolution metabolic flux determination. Biotechniques 34:832–836, 838, 840 passim 47. Choi J, Antoniewicz MR (2011) Tandem mass spectrometry: a novel approach for metabolic flux analysis. Metab Eng 13:225–233 48. Sriram G, Fulton DB, Shanks JV (2007) Flux quantification in central carbon metabolism of Catharanthus roseus hairy roots by 13C labeling and comprehensive bondomer balancing. Phytochemistry 68:2243–2257 49. Libourel IGL, Gehan JP, Shachar-Hill Y (2007) Design of substrate label for steady state flux measurements in plant systems using

Modeling of Isotope Labeling Experiments the metabolic network of Brassica napus embryos. Phytochemistry 68:2211–2221 50. Metallo CM, Walther JL, Stephanopoulos G (2009) Evaluation of 13C isotopic tracers for metabolic flux analysis in mammalian cells. J Biotechnol 144:167–174 51. Wiechert W, Möllney M, Isermann N, Wurzel M, de Graaf AA (1999) Bidirectional reaction

131

steps in metabolic networks: III. Explicit solution and analysis of isotopomer labeling systems. Biotechnol Bioeng 66:69–85 52. Szyperski T (1995) Biosynthetically directed fractional 13C-labeling of proteinogenic amino acids. An efficient analytical tool to investigate intermediary metabolism. Eur J Biochem 232:433–448