Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC1994654

Formats

Article sections

Authors

Related links

Metab Eng. Author manuscript; available in PMC 2007 September 26.

Published in final edited form as:

Published online 2006 September 17. doi: 10.1016/j.ymben.2006.09.001

PMCID: PMC1994654

NIHMSID: NIHMS16238

Department of Chemical Engineering, Bioinformatics and Metabolic Engineering Laboratory, Massachusetts Institute of Technology, Cambridge MA 02139, USA

*corresponding author: Gregory Stephanopoulos, Department of Chemical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, Tel.: 617-253-4583, Fax.: 617-253-3122 ; Email: ude.tim@petsgerg

Metabolic Flux Analysis (MFA) has emerged as a tool of great significance for metabolic engineering and mammalian physiology. An important limitation of MFA, as carried out via stable isotope labeling and GC/MS and NMR measurements, is the large number of isotopomer or cumomer equations that need to be solved, especially when multiple isotopic tracers are used for the labeling of the system. This restriction reduces the ability of MFA to fully utilize the power of multiple isotopic tracers in elucidating the physiology of realistic situations comprising complex bioreaction networks. Here, we present a novel framework for the modeling of isotopic labeling systems that significantly reduces the number of system variables without any loss of information. The elementary metabolite unit (EMU) framework is based on a highly efficient decomposition method that identifies the minimum amount of information needed to simulate isotopic labeling within a reaction network using the knowledge of atomic transitions occurring in the network reactions. The functional units generated by the decomposition algorithm, called elementary metabolite units, form the new basis for generating system equations that describe the relationship between fluxes and stable isotope measurements. Isotopomer abundances simulated using the EMU framework are identical to those obtained using the isotopomer and cumomer methods, however, require significantly less computation time. For a typical ^{13}C-labeling system the total number of equations that needs to be solved is reduced by one order-of-magnitude (100s EMUs vs. 1000s isotopomers). As such, the EMU framework is most efficient for the analysis of labeling by multiple isotopic tracers. For example, analysis of the gluconeogenesis pathway with ^{2}H, ^{13}C, and ^{18}O tracers requires only 354 EMUs, compared to more than 2 million isotopomers.

Accurate flux determination is of great importance for the analysis of cell physiology in fields ranging from metabolic engineering to the study of human metabolic disease (Brunengraber et al., 1997; Hellerstein, 2003; Stephanopoulos, 1999). Initially, metabolic flux analysis (MFA) relied solely on balancing fluxes around metabolites within an assumed network stoichiometry. However, stoichiometric constraints and external flux measurements often did not provide enough information to estimate all fluxes of interest. A more powerful method for flux determination in complex biological systems is based on the use of stable isotopes (Wiechert, 2001). Metabolic conversion of isotopically labeled substrates generates molecules with distinct labeling patterns, i.e. isotope isomers (isotopomers), that can be detected by mass spectrometry (MS) and nuclear magnetic resonance (NMR) spectroscopy (Szyperski, 1995; Szyperski, 1998). Isotope measurements provide many additional independent constraints for MFA. It has been shown that at metabolic and isotopic steady state the isotopomer composition of metabolic intermediates is fully and uniquely determined by the cell’s flux state and the administered isotopic label. Quantitative interpretation of isotopomer data requires the use of mathematical models that describe the relationship between metabolic fluxes and the observed isotopomer abundances (Antoniewicz et al., 2006). Similar to metabolite balancing, balances can be set up around all isotopomers of a particular metabolite. Schmidt et al. (Schmidt et al., 1997) described an elegant method for automatically generating the complete set of isotopomer balances using a matrix based method. More recently, Wiechert et al. (Wiechert et al., 1999) introduced the concept of cumulative isotopomers (cumomers) and provided an efficient procedure for solving isotopomer models. In the forward calculation, isotopomer models are used to simulate a unique profile of isotopomer measurements for given steady-state fluxes. In MFA we are concerned with the more challenging inverse problem, i.e. to determine the flux state of the cell from measurements of isotopomer distributions. Analytical solutions to the inverse problem are only available for very simple systems. Thus, in general, fluxes in complex biological systems are determined by iterative least-squares fitting procedures, where at each iteration the forward problem is solved for an updated set of fluxes.

Isotopomers are defined as isomers of a metabolite that differ only in the labeling state of their individual atoms, for example, ^{13}C vs. ^{12}C in carbon-labeling studies, and ^{2}H vs. ^{1}H in hydrogen-labeling studies. For a metabolite comprising *N* atoms that may be in one of two (labeled or unlabeled) states, 2* ^{N}* isotopomers are possible. Consequently, the number of isotopomers can increase dramatically when multiple tracers are applied. Consider for example glucose (C

Here, we present a novel framework for modeling isotopic tracer systems that significantly reduces the number of system variables without any loss of information. The elementary metabolite units (EMU) framework is a bottom-up modeling approach that is based on a highly efficient decomposition algorithm that identifies the minimum amount of information needed to simulate isotopic labeling within a reaction network. The functional units generated by the decomposition algorithm, called elementary metabolite units (EMUs), form the new basis for generating system equations that describe the relationship between fluxes and isotopomer abundances. Isotopomer abundances simulated using the EMU framework are identical to those obtained using the isotopomer and cumomer methods, however, require significantly fewer variables. We will show that for a typical carbon-13 labeling system the total number of variables and equations that needs to be solved is reduced by one order-of-magnitude (100s EMUs vs. 1000s isotopomers/cumomers). As such, the EMU framework is most efficient for the analysis of labeling by multiple isotopic tracers. We have applied the EMU method to study the gluconeogenesis pathway probed with multiple isotopic tracers, which required only 354 EMUs compared to more than 2×10^{6} isotopomers/cumomers.

We define an elementary metabolite unit of a compound to be a moiety comprising any distinct subset of the compound’s atoms. Consider, for example, metabolite A consisting of 3 atoms. An EMU is a subset of any number of these 3 atoms. The size of an EMU is the number of atoms that are included in the EMU. There are 7 possible EMUs for metabolite A: 3 EMUs of size 1 (A_{1}, A_{2}, A_{3}), 3 EMUs of size 2 (A_{12}, A_{13}, A_{23}), and 1 EMU of size 3 (A_{123}), where the subscript denotes the atoms that are included in the EMU (Figure 1). Note that atoms in an EMU are not necessarily connected by chemical bonds, for example consider EMU A_{13}. In general, for a metabolite comprising *N* atoms 2* ^{N}* -1 EMUs are possible.

Elementary metabolite units (EMU) are distinct subsets of the compound’s atoms. There are 7 EMUs for a 3-atom metabolite A. The subscript in the first column and the shaded areas in the second column denote the atoms that are included in the EMU. **...**

In this paper, we will illustrate that the EMU framework can be used for the simulation of isotopic labeling within a reaction network using the minimum number of variables, all of which are expressed in terms of EMUs. In most cases, only a very small fraction of all possible EMUs is required to simulate the isotopic labeling. This approach is fundamentally different from the isotopomer and cumomer methods, where the simulation model always uses the complete set of all possible isotopomer/cumomers. In the next sections, we will first illustrate the EMU approach for simulating MS measurements, and in section 2.8 we will show how NMR measurements can be simulated using EMUs.

First, we need to introduce the concept of an EMU reaction. Figure 2 shows three types of biochemical reactions that we distinguish: a condensation reaction, a cleavage reaction, and a unimolecular reaction. For each reaction type in Figure 2 we would like to determine the minimum amount of information that is needed to determine the mass isotopomer distribution (MID) of product C, i.e. EMU C_{123}. For the condensation reaction, MID of C_{123} is fully determined by MIDs of EMUs A_{12} and B_{1}. For example, the M+0 abundance of C_{123} is equal to the product of M+0 abundances of A_{12} and B_{1}, i.e. C_{123,M+0}=A_{12,M+0}·B_{1,M+0}. The full MID of C_{123} is obtained from the convolution (or Cauchy product, denoted by ‘×’) of MIDs of A_{12} and B_{1}, i.e. C_{123}=A_{12}×B_{1} (Figure 2). For the cleavage and unimolecular reactions, the MID of C_{123} is equal to the MID of A_{123}. Note that for the cleavage reaction atoms of A that are not transferred to C_{123} are not considered in the EMU reaction, i.e. their labeling doesn’t affect the labeling state of C. Note also that EMU reactions are always size balanced, i.e. the EMU product is formed either from EMUs of the same size, or by condensation of smaller EMUs such that the total size of substrate EMUs equals the size of the EMU product. Thus, there can only be two types of EMU reactions: condensation and unimolecular EMU reactions. Table 1 shows the EMU reactions corresponding to the biochemical reactions in Figure 2.

We will now describe an algorithm that systematically decomposes any biochemical reaction network into EMU reactions using the knowledge of atomic transitions occurring in the network reactions. These EMU reactions will then form the basis for generating model equations for isotopic simulations (next section). Consider the example network shown in Figure 3 that will be used to illustrate the theory behind EMU decomposition. In this network, metabolite A is the sole substrate and metabolites E and F are two network products. The intermediary metabolites B, C and D are assumed to be at metabolic and isotopic steady state. The stoichiometry and atom transitions for the reactions are given in Table 2. The structural input that is required for EMU decomposition is threefold:

- The assumed reaction network stoichiometry
- Atom transitions for all reactions in the network
- List of metabolites/metabolite fragments that need to be simulated

In this example, we would like to set up the simplest possible model to simulate the MID of metabolite F, i.e. EMU F_{123}. The following algorithm systematically identifies all EMU reactions that are needed for this simulation model. First, we identify in the network model all EMU reactions that form EMU F_{123}. In this case, F_{123} is formed only in reaction 6 from EMU D_{123}. We record this EMU reaction and the newly identified EMU(s) and repeat this process for all new EMUs, starting with the largest EMU size, i.e. D_{123}. Here, D_{123} is formed in two reactions, i.e. in reaction 2 from B_{123}, and in reaction 5 from B_{23}+C_{1}. Next, B_{123} is formed in reactions 1 and 3 from A_{123} and D_{123}, respectively. Note that A_{123} is a network substrate, i.e. it is not produced by any reaction in the network, and D_{123} was already considered in the previous step. Thus, we have identified all EMU reactions of size 3 that need to be considered. Next, we proceed with EMUs of size 2 that were previously identified; here, B_{23} is formed in reaction 1 and 3 from A_{23}, and D_{23}, respectively, etc. We complete this process for all new EMUs of size 2 and 1, until the EMUs have been traced back to EMUs of network substrates or EMUs that were already accounted for. Table 3 shows the complete EMU decomposition for this system. In this case, 18 EMU reactions were identified connecting 14 EMUs. Of these 14 EMUs, 10 EMUs correspond to intermediary metabolites whose labeling is unknown, and 4 EMUs are fully defined by the choice of substrate labeling of metabolite A. It should be clear that the described decomposition algorithm is exhaustive, unsupervised, and always identifies the minimal set of EMU reactions that need to be considered in the simulation model. Furthermore, this algorithm is easily implemented and is computationally efficient, i.e. the decomposition completes within seconds even for the largest network model that we have considered. The main advantage of the EMU decomposition is that metabolites are never broken into smaller pieces than is strictly required to describe the labeling state of the selected metabolites. In contrast, the isotopomer and cumomer frameworks always use all 2* ^{N}* isotopomers per metabolite to simulate the system. In this example, the complete set of 36 isotopomers describe the system, with 28 unknown isotopomers and 8 substrate isotopomers. Thus, in this example the number of system variables was reduced by more than 50% using the EMU framework.

The EMU reactions obtained from network decomposition form the new basis for generating system equations. EMU reaction networks are much less connected than isotopomer systems, and we can group EMU reactions into decoupled reaction networks based on: (i) the EMU size, and (ii) network connectivity (see section 3.4 for an example). For the example network we obtain three decoupled reaction networks of EMU size 1, 2, and 3, respectively (Figure 4). Similar to metabolite and isotopomer balancing, we can set up balances around all EMUs. In general, EMU balances can be written as a set of equations that are linear in the unknown EMU variables:

Decoupled EMU reaction networks generated for the simulation of metabolite F. Note that the EMU reaction networks contain only these EMUs that are strictly required to simulate F.

$${\text{A}}_{1,\text{k}}(\text{v})\xb7{\text{X}}_{1,\text{k}}={\text{B}}_{1,\text{k}}(\text{v})\xb7{\text{Y}}_{1,\text{k}}({{\text{y}}_{1}}^{\text{in}})$$

(1a)

$${\text{A}}_{2,\text{k}}(\text{v})\xb7{\text{X}}_{2,\text{k}}={\text{B}}_{2,\text{k}}(\text{v})\xb7{\text{Y}}_{2,\text{k}}({{\text{y}}_{2}}^{\text{in}},{\text{X}}_{1})$$

(1b)

...

$${\text{A}}_{\text{n},\text{k}}(\text{v})\xb7{\text{X}}_{\text{n},\text{k}}={\text{B}}_{\text{n},\text{k}}(\text{v})\xb7{\text{Y}}_{\text{n},\text{k}}({{\text{y}}_{\text{n}}}^{\text{in}},{\text{X}}_{\text{n}-1},\dots ,{\text{X}}_{1})$$

(1c)

Matrices X_{i,k} and Y_{i,k} represent to the unknown (balanced), and known EMU variables, respectively, and y_{i}^{in} are EMUs of network substrates. The first subscript denotes the EMU size and the second the subnetwork number, i.e. in case there are multiple decoupled EMU networks of the same EMU size. For brevity will omit the second subscript if there is only one subnetwork of a certain EMU size. Each row in matrices Y_{i,k} and X_{i,k} contains the MID for the corresponding EMU. The EMU balances written in matrix form corresponding to the EMU networks in Figure 4 are shown below.

$$\left[\begin{array}{ccccc}-{\text{v}}_{4}& {\text{v}}_{4}& 0& 0& 0\\ 0& -{\text{v}}_{1}-{\text{v}}_{3}& {\text{v}}_{3}& 0& 0\\ 0& {\text{v}}_{2}& -{\text{v}}_{2}-{\text{v}}_{5}& {\text{v}}_{5}& 0\\ 0& 0& 0& -{\text{v}}_{1}-{\text{v}}_{3}& {\text{v}}_{3}\\ {\text{v}}_{5}& 0& 0& {\text{v}}_{2}& -{\text{v}}_{2}-{\text{v}}_{5}\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{C}}_{1}\\ {\text{B}}_{2}\\ {\text{D}}_{2}\\ {\text{B}}_{3}\\ {\text{D}}_{3}\end{array}\right]\hspace{0.17em}=\hspace{0.17em}\left[\begin{array}{cc}0& 0\\ -{\text{v}}_{1}& 0\\ 0& 0\\ 0& -{\text{v}}_{1}\\ 0& 0\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{A}}_{2}\\ {\text{A}}_{3}\end{array}\right]$$

(2a)

$$\left[\begin{array}{cc}-{\text{v}}_{5}-{\text{v}}_{2}& {\text{v}}_{2}\\ {\text{v}}_{3}& -{\text{v}}_{1}-{\text{v}}_{3}\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{D}}_{23}\\ {\text{B}}_{23}\end{array}\right]=\left[\begin{array}{cc}-{\text{v}}_{5}& 0\\ 0& -{\text{v}}_{1}\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{B}}_{3}\times {\text{C}}_{1}\\ {\text{A}}_{23}\end{array}\right]$$

(2b)

$$\left[\begin{array}{ccc}-{\text{v}}_{6}& {\text{v}}_{6}& 0\\ 0& -{\text{v}}_{5}-{\text{v}}_{2}& {\text{v}}_{2}\\ 0& {\text{v}}_{3}& -{\text{v}}_{1}-{\text{v}}_{3}\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{F}}_{123}\\ {\text{D}}_{123}\\ {\text{B}}_{123}\end{array}\right]=\left[\begin{array}{cc}0& 0\\ -{\text{v}}_{5}& 0\\ 0& -{\text{v}}_{1}\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{B}}_{23}\times {\text{C}}_{1}\\ {\text{A}}_{123}\end{array}\right]$$

(2c)

The product term B_{3}×C_{1} in Eq. 2b represents the convolution of MIDs of B_{3} and C_{1} (see section 2.2). Matrices X_{2} and Y_{2} written out in full for the example network are shown below.

$$\left[\begin{array}{c}{\text{D}}_{23}\\ {\text{B}}_{23}\end{array}\right]=\left[\begin{array}{ccc}{\text{D}}_{23,\text{M}+0}& {\text{D}}_{23,\text{M}+1}& {\text{D}}_{23,\text{M}+2}\\ {\text{B}}_{23,\text{M}+0}& {\text{B}}_{23,\text{M}+1}& {\text{B}}_{23,\text{M}+2}\end{array}\right]$$

(3)

$$\left[\begin{array}{c}{\text{B}}_{3}\times {\text{C}}_{1}\\ {\text{A}}_{23}\end{array}\right]=\left[\begin{array}{ccc}{\text{B}}_{3,\text{M}+0}\xb7{\text{C}}_{1,\text{M}+0}& ({\text{B}}_{3,\text{M}+0}\xb7{\text{C}}_{1,\text{M}+1}+{\text{B}}_{3,\text{M}+1}\xb7{\text{C}}_{1,\text{M}+0})& {\text{B}}_{3,\text{M}+1}\xb7{\text{C}}_{1,\text{M}+1}\\ {\text{A}}_{23,\text{M}+0}& {\text{A}}_{23,\text{M}+1}& {\text{A}}_{23,\text{M}+2}\end{array}\right]$$

(4)

The size of matrix A_{i} depends on the number of balanced EMUs, and the size of matrix B_{i} depends on the number of input EMUs in that EMU network. Thus, for *n* balanced and *m* input EMUs, A_{i} is an *n* × *n* matrix and B_{i} is an *n* × *m* matrix. Note that matrices A_{i} and B_{i} are always strictly linear functions of fluxes, i.e. dA_{i}/dv_{j} and dB_{i}/dv_{j} are constant for any metabolic system. To simulate the isotopic labeling distribution for the selected metabolites, EMU balances are solved sequentially starting with the smallest EMU size network. Since the matrix Y_{1} is known, i.e. it is fully determined by EMUs of network substrate A, we can easily calculate X_{1} using standard linear algebra techniques:

$${\text{X}}_{\text{i}}={{\text{A}}_{\text{i}}}^{-1}\xb7{\text{B}}_{\text{i}}\xb7{\text{Y}}_{\text{i}}$$

(5)

For subsequent EMU networks, matrices Y_{i} may depend on previously determined EMUs of smaller size. Thus, for larger EMU sizes we must first update matrix Y_{i} and then calculate X_{i}. After all EMUs have been computed we can simply read out the simulated MIDs for the selected metabolites of interest from the rows of matrices X_{i}. For example, the MID of F_{123} is found in the first row of matrix X_{3}.

Flux estimation algorithms and algorithms used for calculating flux confidence intervals require repeated calculation of first order derivatives, i.e. sensitivities of simulated measurements with respect to fluxes (Antoniewicz et al., 2006). Here, we present analytical expressions for the calculation of these sensitivities from EMU balances. In general, EMU balances are expressed in the following matrix form:

$${\text{A}}_{\text{i}}\xb7{\text{X}}_{\text{i}}={\text{B}}_{\text{i}}\xb7{\text{Y}}_{\text{i}}$$

(6)

where matrices A_{i} and B_{i} are strictly linear functions of fluxes and matrices X_{i} and Y_{i} are nonlinear functions of fluxes. Implicit differentiation of Eq. 6 yields:

$$\frac{{\text{dA}}_{\text{i}}}{\text{dv}}\xb7{\text{X}}_{\text{i}}+{\text{A}}_{\text{i}}\xb7\frac{{\text{dX}}_{\text{i}}}{\text{dv}}=\frac{{\text{dB}}_{\text{i}}}{\text{dv}}\xb7{\text{Y}}_{\text{i}}+{\text{B}}_{\text{i}}\xb7\frac{{\text{dY}}_{\text{i}}}{\text{dv}}$$

(7)

Note that matrices dA_{i}/dv_{j} and dB_{i}/dv_{j} are constant for a given metabolic system. After rearrangement of Eq. 7 we obtain the following general expression for dX_{i}/dv:

$$\frac{{\text{dX}}_{\text{i}}}{\text{dv}}={{\text{A}}_{\text{i}}}^{-1}\xb7\left(\frac{{\text{dB}}_{\text{i}}}{\text{dv}}\xb7{\text{Y}}_{\text{i}}+{\text{B}}_{\text{i}}\xb7\frac{{\text{dY}}_{\text{i}}}{\text{dv}}-\frac{{\text{dA}}_{\text{i}}}{\text{dv}}\xb7{\text{X}}_{\text{i}}\right)$$

(8)

The above expression may be simplified for the smallest EMU size network, where all terms in matrix Y are given by network substrates and thus are constant, i.e. dY/dv=0. For EMU networks of larger EMU size, matrices Y_{i} may depend on previously determined EMUs of smaller size. In the example network Y_{2} was given by:

$${\text{Y}}_{2}=\left[\begin{array}{c}{\text{B}}_{3}\times {\text{C}}_{1}\\ {\text{A}}_{23}\end{array}\right]$$

(9)

Where, B_{3} and C_{1} are EMUs of size 1 that are calculated from EMU size 1 balances. Application of the product rule yields the following expression for the first order derivative of the convolution of B_{3} and C_{1} (i.e. first row of Y_{2}):

$$\frac{\text{d}}{\text{dv}}\left({\text{B}}_{3}\times {\text{C}}_{1}\right)=\frac{{\text{dB}}_{\text{3}}}{\text{dv}}\times {\text{C}}_{1}+{\text{B}}_{3}\times \frac{{\text{dC}}_{\text{1}}}{\text{dv}}$$

(10)

where matrices dC_{1}/dv and dB_{3}/dv are obtained from the previously determined matrix dX_{1}/dv. Note that the second row of Y_{2} contains A_{23} which is a network substrate EMU, i.e. it is considered constant (dA_{23}/dv=0). Thus, we obtain the following expression for dY_{2}/dv:

$$\frac{{\text{dY}}_{\text{2}}}{\text{dv}}=\left[\begin{array}{c}\frac{{\text{dB}}_{\text{3}}}{\text{dv}}\times {\text{C}}_{1}+{\text{B}}_{3}\times \frac{{\text{dC}}_{\text{1}}}{\text{dv}}\\ 0\end{array}\right]$$

(11)

where (dB_{3}/dv)×C_{1} denotes the 2D-convolution of matrix (dB_{3}/dv) and vector C_{1}. Thus, first order derivatives for any EMU size can be obtained sequentially from the EMU balances. Figure 5 summarizes the general procedure for simulating labeling distributions and calculating first order derivatives of simulated measurements with respect to fluxes using the EMU framework.

For non-zero fluxes matrices A_{i} are compartmental matrices, and from compartmental matrix theory we know that therefore matrices A_{i} are always invertible (Anderson, 1983). In other words, the EMU approach will always compute a unique and stable solution for all EMUs. The most time-consuming operation in solving EMU balances and calculating sensitivities is the inversion of matrices A_{i}, or rather matrix factorization of A_{i}. In general, the computational time for LU decomposition increases with the size of the matrix, i.e. the number of EMU nodes for each subnetwork. We found empirically that for sparsely connected EMU networks, such as the ones shown in Figure 4, the computation time increased linearly with the number of EMUs, i.e. *O*(*n*). For more highly connected networks, for example, the EMU networks corresponding to central carbon metabolism of *E. coli* (section 3.4), the computational time increased to the third power of the number of EMUs, i.e. *O*(*n ^{3}*). Thus, it will often be worthwhile to reduce the number of EMUs by eliminating EMU nodes that have a single influx. This will be illustrated in detail in section 3.3.

There is a number of considerations regarding the chemical structure of metabolites that need to be taken into account for EMU simulation models. Here, we will discuss in detail how to account for chiral, prochiral and rotationally symmetric metabolites, and how biochemically equivalent hydrogen, oxygen and other atoms should be modeled within the EMU framework. These considerations are equally important for the construction of isotopomer models, however, until now they have not received proper attention.

Tetrahedral carbon atoms with four different ligands are called chiral, whereas the term prochiral applies to carbon atoms that hold two stereoheterotopic groups (Moss, 1996). Many biological metabolites contain one or more chiral and/or prochiral carbon atoms. It is well known that biochemical reactions are highly stereospecific, i.e. enzymes differentiate between prochiral atoms and prochiral atom groups. Therefore, it is important to keep track of individual prochiral atoms in a network model and assign stereospecific atom transitions to all biochemical reactions. Consider for example the enzymatic reaction catalyzed by aconitase that converts citrate to isocitrate (Figure 6). Three of the six carbon atoms of citrate are prochiral, i.e. C2, C3 and C4. The enzyme aconitase stereospecifically transfers the pro-R hydrogen from the pro-R arm (i.e. C1-C2) of citrate to C3 of isocitrate, and produces only one of four possible stereoisomers of isocitrate, i.e. (2R,3S)-isocitrate. Note also that the prochirality of C4 is not altered by aconitase. The absolute stereochemistry for many bioreactions has been worked out in detail and can be found in biochemistry textbooks and other general literature.

Biological molecules often contain groups of atoms that are biochemically indistinguishable. If we assume that there is no metabolic channeling, then this results in rapid equilibration of labeling at these positions. Consider for example the chemical structure of pyruvate shown in Figure 7. The three hydrogen atoms at C3 are biochemically equivalent, i.e. enzymes cannot distinguish between these atoms. Furthermore, the two oxygen atoms at C1 are biochemically equivalent due to resonance stabilization. Thus, not all of pyruvate’s EMUs are independent. For example, there are six equivalent EMUs of pyruvate that contain three carbon atoms, two of the three hydrogen atoms at C3, and one of two oxygen atom at C1. Figure 7 shows the six equivalent EMUs. We can predict the number of equivalent EMUs for any EMU as follows:

Six equivalent EMUs of pyruvate. Three hydrogen atoms of pyruvate at C3 are biochemically equivalent, and two oxygen atoms at C1 are equivalent due to resonance stabilization. As such, there are six equivalent EMUs of pyruvate containing all three carbon **...**

$$\text{No}.\hspace{0.17em}\text{of\hspace{0.17em}equivalent\hspace{0.17em}EMUs}=\prod _{\begin{array}{l}\text{for\hspace{0.17em}each\hspace{0.17em}group\hspace{0.17em}of}\\ \text{equivalent\hspace{0.17em}atoms}\end{array}}\left(\begin{array}{c}\text{No}.\hspace{0.17em}\text{of\hspace{0.17em}equivalent\hspace{0.17em}atoms\hspace{0.17em}in\hspace{0.17em}the}\hspace{0.17em}\text{group}\\ \text{No}.\hspace{0.17em}\text{of\hspace{0.17em}these\hspace{0.17em}atoms\hspace{0.17em}in\hspace{0.17em}the\hspace{0.17em}EMU}\end{array}\right)$$

(12)

When setting up EMU balances, we need to account for the presence of equivalent EMUs. We propose the following method. First, whenever a new EMU is generated during EMU network decomposition, we identify all equivalent EMUs for that EMU. Then, for each equivalent EMU we find the EMU reaction(s) that produce that EMU, and divide the contribution from each reaction by the total number of equivalent EMUs. Note that this way we introduce only one EMU variable for each group of equivalent EMUs. For example, consider the enzymatic reaction catalyzed by malic enzyme that converts malate to pyruvate shown in Figure 8 (arbitrary numbering of atoms). The six equivalent EMUs of pyruvate from Figure 7 are produced by the following six EMU reactions:

Malic enzyme converts malate to pyruvate. Note that one of the three hydrogen atoms at C_{#6} of pyruvate is derived from the solvent. The two prochiral hydrogen atoms of malate at C_{#7}, which are biochemically distinct, become indistinguishable after malate **...**

$${\text{Mal}}_{13478}+\text{H}\to {\text{Pyr}}_{134678}$$

(13a)

$${\text{Mal}}_{134789}\to {\text{Pyr}}_{134679}$$

(13b)

$${\text{Mal}}_{13479}+\text{H}\to {\text{Pyr}}_{134689}$$

(13c)

$${\text{Mal}}_{12478}+\text{H}\to {\text{Pyr}}_{124678}$$

(13d)

$${\text{Mal}}_{124789}\to {\text{Pyr}}_{124679}$$

(13e)

$${\text{Mal}}_{12479}+\text{H}\to {\text{Pyr}}_{124689}$$

(13f)

We can represent these 6 equivalent EMUs of pyruvate by a single EMU variable {Pyr_{134678}}, where the brackets denote the presence of multiple equivalent EMUs. Next, we also identify that Mal_{13478} and Mal_{12478}; Mal_{13479} and Mal_{12479}; and Mal_{134789} and Mal_{124789} are equivalent EMUs. Taken together, we obtain the following overall EMU reaction for the production of {Pyr_{134678}}:

$$\left\{{\text{Pyr}}_{134678}\right\}=\frac{1}{3}\left\{{\text{Mal}}_{13478}\right\}\times \text{H}+\frac{1}{3}\left\{{\text{Mal}}_{13479}\right\}\times \text{H}+\frac{1}{3}\left\{{\text{Mal}}_{134789}\right\}$$

(14)

Thus, {Pyr_{134678}} is formed from three EMUs of malate, i.e. {Mal_{13478}}, {Mal_{13479}}, and {Mal_{134789}}; two groups are of EMU size 5 that combine with H^{+}, and one is of EMU size 6. Note also that the two prochiral hydrogen atoms of malate that are initially biochemically distinct become indistinguishable after malate is converted to pyruvate.

A number of metabolites of central carbon metabolism are rotationally symmetric, i.e. they are superposable on themselves by rotation. In isotopic labeling studies these molecules also cause scrambling of isotopic labeling. It is important to note that there is a clear difference between molecules with a rotation axis and molecules with a center of inversion (Moss, 1996). The two types of symmetry have different characteristics. Only molecules with a rotation axis are superposable on themselves. Figure 9 shows the structures of (2S,3R)-butane-1,2,3,4-tetraol (i.e. erythritol) which has a center of inversion, and (2R,3R)-butane-1,2,3,4-tetraol which has a rotation axis. Carbon atoms C1 and C4, and C2 and C3 of erythriol are chemically equivalent, i.e. they react identically in chemical reactions and have the same chemical properties, however, in enzymatic reactions these groups are biochemically distinct. In contrast, carbon atoms C1 and C4, and C2 and C3 of (2R,3R)-butane-1,2,3,4-tetraol are both chemically and biochemically equivalent, i.e. they are not distinguished by enzymes, and thus result in label scrambling.

Differences between molecules with a rotation axis and center of inversion. (2S,3R)-Butane-1,2,3,4-tetraol (i.e. erythritol) has a center of inversion and is not superposable on itself. Thus, carbon atoms C1 and C4, and C2 and C3 of erythritol are biochemically **...**

The most important examples of rotationally symmetric molecules in central metabolism are fumarate, succinate, and D-mannitol. Figure 10 shows the structure of fumarate (with arbitrary numbering of atoms). It should be clear that fumarate has a rotation axis, i.e. after rotating 180° fumarate superposes on itself. Furthermore, the two oxygen atoms at the first carbon and the two oxygen atoms at the last carbon atom are biochemically equivalent due to resonance stabilization. These characteristics of fumarate need to be considered when we identify equivalent EMUs of fumarate. It should also be clear that the number of equivalent EMUs may increase due to rotational symmetry. EMUs of rotationally symmetric molecules may have twice as many equivalent EMUs as nonsymmetric molecules, assuming a rotational symmetry of 180 degrees. The multiplier will increase if the angle is less, e.g. in benzene where the rotational symmetry is 60 degrees the multiplier would be 6 (=360/6)

Equivalent EMUs of fumarate, a rotationally symmetric molecule. The following four EMUs are equivalent: Fum_{12467}, Fum_{13467}, Fum_{45689}, and Fum_{4568,10} (numbering of fumarate atoms is arbitrary). Shaded areas indicate atoms included in the EMUs.

$$\text{No}.\hspace{0.17em}\text{of\hspace{0.17em}equivalent\hspace{0.17em}EMUs}\le \text{2}\times \prod _{\begin{array}{l}\text{for\hspace{0.17em}each\hspace{0.17em}group\hspace{0.17em}of}\\ \text{equivalent\hspace{0.17em}atoms}\end{array}}\left(\begin{array}{c}\text{No}.\hspace{0.17em}\text{of}\hspace{0.17em}\text{equivalent\hspace{0.17em}atoms\hspace{0.17em}in\hspace{0.17em}the\hspace{0.17em}group}\\ \text{No}.\hspace{0.17em}\text{of\hspace{0.17em}these}\hspace{0.17em}\text{atoms\hspace{0.17em}in\hspace{0.17em}the\hspace{0.17em}EMU}\end{array}\right)$$

(15)

Figure 10 shows, for example, the four equivalent EMUs of fumarate Fum_{12467}. Note also, for example, that fumarate EMU F_{1468} has no equivalent EMUs, because its rotational equivalent is F_{1468} itself. Hence, when we enumerate equivalent EMUs it is very important that we separate the effects of rotational symmetry, which is a global characteristic of a molecule, from the effects of equivalent hydrogen and oxygen atoms, which are local characteristics. To better illustrate this, consider hydrogen atoms #5 and #7 of fumarate in Figure 10. These atoms cannot be treated as equivalent atoms, because that would incorrectly identify Fum_{12465} as being equivalent to Fum_{12467}. Taking these considerations into account, Figure 11 shows the complete algorithm for decomposing metabolic networks into EMU reaction networks.

Thus far, we have shown how MS measurements can be simulated using the EMU approach. We will now illustrate the method for simulating NMR measurements, in particular we will show how fractional enrichments and NMR fine spectra can be simulated.

Fractional enrichments are a measure of the fractional abundance of ^{13}C-atoms at specific carbon positions in a molecule. For example, Wiechert et al. (Wiechert et al., 1997) measured fractional enrichments of 25 carbon atoms of amino acids. In the EMU framework fractional enrichments are modeled by EMUs of size 1, i.e. containing a single carbon atom. Thus, to simulate fractional enrichments we can simply solve EMU balances of size 1. Network decomposition is accomplished with the same algorithm as was described before (section 2.3). In the EMU balances (Eq. 1), X_{i} and Y_{i} can now be represented by vectors that contain the fractional labeling of carbon atoms. Note that the EMU method for simulating fractional enrichments is very similar to the atom mapping matrix method that was originally proposed by Zupke and Stephanopoulos (Zupke and Stephanopoulos, 1994), and the weight-1 cumomer method as proposed by Wiechert et al. (Wiechert et al., 1999).

Data obtained from 2D [^{13}C,^{1}H] COSY spectra, also known as NMR fine spectra, provide information on the relative amount of ^{13}C-^{13}C and ^{13}C-^{12}C carbons at specific carbon positions, where the observed carbon atom is always ^{13}C-labeled and the adjacent carbon atoms are either labeled or unlabeled (Szyperski, 1995). For a secondary carbon atom, NMR fine spectra can be expressed as ratios of four isotopomer fractions, i.e. A_{010}, A_{011}, A_{110}, and A_{111}:

$$\left[\begin{array}{c}\text{singlet}\\ \text{doublet\hspace{0.17em}}3\\ \text{doublet\hspace{0.17em}}1\\ \text{doublet\hspace{0.17em}doublet}\end{array}\right]=\left[\begin{array}{c}{\text{A}}_{010}\\ {\text{A}}_{011}\\ {\text{A}}_{110}\\ {\text{A}}_{111}\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}{({\text{A}}_{010}+{\text{A}}_{011}+{\text{A}}_{110}+{\text{A}}_{111})}^{-1}$$

(16)

We can obtain these four isotopomer fractions from cumomer fractions A_{x1x}, A_{x11}, A_{11x}, and A_{111}, as was shown previously by van Winden et al. (van Winden et al., 2002). According to the definition by Wiechert et al. (Wiechert et al., 1999), A_{x1x} denotes the weight-1 cumomer fraction for which the second atom is ^{13}C-labeled and the other two atoms are labeled or unlabeled, i.e. x = *‘0 or 1’*. Alternatively, we have derived that the same isotopomer fractions can also be obtained from the following four EMUs: A_{2}, A_{23}, A_{12}, and A_{123}. We can easily show that cumomer fraction A_{x1x} is equal to the M+1 abundance of EMU A_{2}, i.e. the fraction of fully labeled EMU A_{2}; weight-2 cumomer fractions A_{11x} and A_{x11} are equal to the M+2 abundances of EMUs A_{23} and A_{12}, respectively (i.e. fully labeled EMUs); and weight-3 cumomer fraction A_{111} is equal to the M+3 abundance of EMU A_{123} (i.e. fully labeled EMUs).

$$\left[\begin{array}{c}{\text{A}}_{010}\\ {\text{A}}_{110}\\ {\text{A}}_{011}\\ {\text{A}}_{111}\end{array}\right]={\left[\begin{array}{cccc}1& 1& 1& 1\\ 0& 1& 0& 1\\ 0& 0& 1& 1\\ 0& 0& 0& 1\end{array}\right]}^{-1}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{A}}_{\text{x}1\text{x}}\\ {\text{A}}_{\text{x}11}\\ {\text{A}}_{11\text{x}}\\ {\text{A}}_{111}\end{array}\right]={\left[\begin{array}{cccc}1& 1& 1& 1\\ 0& 1& 0& 1\\ 0& 0& 1& 1\\ 0& 0& 0& 1\end{array}\right]}^{-1}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{A}}_{2}\\ {\text{A}}_{23}\\ {\text{A}}_{12}\\ {\text{A}}_{123}\end{array}\right]$$

(17)

Thus, we can simulate NMR fine spectra either by solving weight-1, 2, and 3 cumomer balances (van Winden et al., 2002), or alternatively by solving EMU balances for the EMUs A_{2}, A_{23}, A_{12}, and A_{123}. For convenience, X_{i} and Y_{i} can be represented as vectors (not matrices) that contain the fractional abundances of the fully labeled EMUs. It should be clear that the number of EMUs generated for the EMUs of size 1, 2, and 3 will always be smaller, or equal to the number of weight-1, 2, and 3 cumomers. Therefore, it is always more efficient to simulate NMR fine spectra using the EMU framework.

The EMU framework can be extended even further to capture NMR fine spectra for tertiary carbon atoms, and long-range ^{13}C-^{13}C couplings. In that case, the following eight EMUs need to be simulated: A_{2}, A_{12}, A_{23}, A_{24}, A_{123}, A_{234}, A_{124}, and A_{1234}. After simulation of these EMUs, we can convert the calculated fractional abundances of fully labeled EMUs to isotopomer fractions using Eq. 18, and then to NMR signal intensities using Eq. 19.

$$\left[\begin{array}{c}{\text{A}}_{0100}\\ {\text{A}}_{0101}\\ {\text{A}}_{0110}\\ {\text{A}}_{0111}\\ {\text{A}}_{1100}\\ {\text{A}}_{1101}\\ {\text{A}}_{1110}\\ {\text{A}}_{1111}\end{array}\right]={\left[\begin{array}{cccccccc}1& 1& 1& 1& 1& 1& 1& 1\\ 0& 1& 0& 1& 0& 1& 0& 1\\ 0& 0& 1& 1& 0& 0& 1& 1\\ 0& 0& 0& 1& 0& 0& 0& 1\\ 0& 0& 0& 0& 1& 1& 1& 1\\ 0& 0& 0& 0& 0& 1& 0& 1\\ 0& 0& 0& 0& 0& 0& 1& 1\\ 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]}^{-1}\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{A}}_{2}\\ {\text{A}}_{24}\\ {\text{A}}_{23}\\ {\text{A}}_{234}\\ {\text{A}}_{12}\\ {\text{A}}_{124}\\ {\text{A}}_{123}\\ {\text{A}}_{1234}\end{array}\right]$$

(18)

$$\left[\begin{array}{c}\text{singlet}\\ \text{doublet\hspace{0.17em}}4\\ \text{doublet\hspace{0.17em}}3\\ \text{double\hspace{0.17em}doublet\hspace{0.17em}34}\\ \text{doublet\hspace{0.17em}1}\\ \text{double\hspace{0.17em}doublet\hspace{0.17em}14}\\ \text{double\hspace{0.17em}doublet\hspace{0.17em}13}\\ \text{quadruple\hspace{0.17em}doublet}\end{array}\right]=\left[\begin{array}{c}{\text{A}}_{0100}\\ {\text{A}}_{0101}\\ {\text{A}}_{0110}\\ {\text{A}}_{0111}\\ {\text{A}}_{1100}\\ {\text{A}}_{1101}\\ {\text{A}}_{1110}\\ {\text{A}}_{1111}\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}{\left({\text{A}}_{0100}+{\text{A}}_{0101}+{\text{A}}_{0110}+{\text{A}}_{0111}+{\text{A}}_{1100}+{\text{A}}_{1101}+{\text{A}}_{1110}+{\text{A}}_{1111}\right)}^{-1}$$

(19)

Note that in Eqs. 18 and 19 we have assumed that the observed carbon atom is C2.

In this example we will compare the EMU framework for simulating mass isotopomer distributions to the isotopomer and cumomer frameworks. Consider the simple network model that was introduced in section 2.3 (Figure 3). The assumed steady-state fluxes and labeling of substrate A are shown in Figure 3. The solution to the EMU balances from Eq. 2 are shown below.

Solution of EMU balances for reaction network of EMU size 1

$$\left[\begin{array}{c}{\text{C}}_{1}\\ {\text{B}}_{2}\\ {\text{D}}_{2}\\ {\text{B}}_{3}\\ {\text{D}}_{3}\end{array}\right]={\left[\begin{array}{ccccc}-20& 20& 0& 0& 0\\ 0& -150& 50& 0& 0\\ 0& 110& -130& 20& 0\\ 0& 0& 0& -150& 50\\ 20& 0& 0& 110& -130\end{array}\right]}^{-1}\xb7\hspace{0.17em}\left[\begin{array}{cc}0& 0\\ -100& 0\\ 0& 0\\ 0& -100\\ 0& 0\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]=\left[\begin{array}{cc}0.0667& 0.9333\\ 0.0667& 0.9333\\ 0.2000& 0.8000\\ 0.9333& 0.0667\\ 0.8000& 0.2000\end{array}\right]$$

Solution of EMU balances for reaction network of EMU size 2

$$\left[\begin{array}{c}{\text{D}}_{23}\\ {\text{B}}_{23}\end{array}\right]={\left[\begin{array}{cc}-130& 110\\ 50& -150\end{array}\right]}^{-1}\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{cc}-20& 0\\ 0& -100\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{ccc}0.0622& 0.8756& 0.0622\\ 0& 1& 0\end{array}\right]=\left[\begin{array}{ccc}0.0133& 0.9733& 0.0133\\ 0.0044& 0.9911& 0.0044\end{array}\right]$$

Solution of EMU balances for reaction network of EMU size 3

$$\begin{array}{c}\left[\begin{array}{c}{\text{F}}_{123}\\ {\text{D}}_{123}\\ {\text{B}}_{123}\end{array}\right]={\left[\begin{array}{ccc}-80& 80& 0\\ 0& -130& 110\\ 0& 50& -150\end{array}\right]}^{-1}\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{cc}0& 0\\ -20& 0\\ 0& -100\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{cccc}0.0003& 0.0702& 0.9253& 0.0041\\ 0& 1& 0& 0\end{array}\right]\\ =\left[\begin{array}{cccc}0.0001& 0.8008& 0.1983& 0.0009\\ 0.0001& 0.8008& 0.1983& 0.0009\\ 0.0000& 0.9336& 0.0661& 0.0003\end{array}\right]\end{array}$$

Thus, we find that the simulated MID of F is 0.0 mol% (M+0), 80.1 mol% (M+1), 19.8 mol% (M+2), and 0.1 mol% (M+3). These simulated abundances are identical to those obtained using the isotopomer and cumomer methods. The main difference between the methods is the number of equations that was solved to simulate the labeling of metabolite F. For the isotopomer method, 28 nonlinear isotopomer balances were solved using Newton’s iterative method. With the cumomer method, 4 linear problems of size 4, 11, 10, and 3, respectively, were solved using standard linear algebra techniques. Note that the total number of cumomers was the same as the number of isotopomers. In contrast, with the EMU model 3 linear problems of size 5, 2, and 3, respectively, were solved. Table 4 summarizes the comparison between the three modeling methods for this simple example.

The second example that we consider is the simplified model of the tricarboxylic acid (TCA) cycle shown in Figure 12. The stoichiometry and atom transitions for the eight reactions are given in Table 5. In this network, acetyl coenzyme A and aspartate are the two substrates, and glutamate and carbon dioxide are the products. Here, we simulated the steady-state labeling distribution of glutamate assuming a mixture of 25% [2-^{13}C]AcCoA and 25% [1,2-^{13}C]AcCoA as the tracer input. The assumed flux distribution is shown in Figure 12. In this example, we only considered the labeling of carbon atoms and for simplicity ignored natural isotope enrichments. The algorithm described in section 2.3 decomposed the TCA cycle network into 4 decoupled EMU reaction networks that are shown in Figure 13. The total number of balanced EMUs was 24, i.e. 8 EMUs in the first network (EMU size 1), 5 EMUs in the second (EMU size 2), 8 EMUs in the third (EMU size 3), and 3 EMUs in the fourth network (EMU size 5). The EMU balances for the four decoupled networks are shown below.

Simplified model of the tricarboxylic acid cycle. Abbreviations of metabolites: OAC, oxaloacetate; Asp, aspartate; AcCoA, acetyl coenzyme A; Cit, citrate; AKG, α-ketoglutarate; Glu, glutamate; Suc, succinate. The assumed fluxes have arbitrary **...**

EMU reaction networks generated for glutamate from EMU network decomposition of the TCA cycle. The complete molecule of glutamate corresponds to EMU Glu_{12345}. Subscripts denote carbon atoms that are included in the EMUs. Abbreviations of metabolites are **...**

Stoichiometry and atom transitions for reactions of the TCA cycle. This network model was used to simulate the steady-state mass isotopomer distribution of glutamate.

$$\left[\begin{array}{cccccccc}-{\text{v}}_{6}-{\text{v}}_{8}& {\text{v}}_{6}& .& .& .& .& .& .\\ {\scriptstyle \frac{1}{2}}{\text{v}}_{7}& -{\text{v}}_{5}-{\text{v}}_{7}& {\scriptstyle \frac{1}{2}}{\text{v}}_{7}& {\text{v}}_{5}& .& .& .& .\\ .& {\text{v}}_{6}& -{\text{v}}_{6}-{\text{v}}_{8}& .& .& .& .& .\\ .& .& .& -{\text{v}}_{4}& {\scriptstyle \frac{1}{2}}{\text{v}}_{4}& {\scriptstyle \frac{1}{2}}{\text{v}}_{4}& .& .\\ .& .& .& .& -{\text{v}}_{2}& .& {\text{v}}_{2}& .\\ .& .& .& .& .& -{\text{v}}_{2}& .& {\text{v}}_{2}\\ {\text{v}}_{1}& .& .& .& .& .& -{\text{v}}_{1}& .\\ .& .& .& .& .& .& .& -{\text{v}}_{1}\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{OAC}}_{2}\\ {\text{Fum}}_{2}\\ {\text{OAC}}_{3}\\ {\text{Suc}}_{2}\\ {\text{AKG}}_{3}\\ {\text{AKG}}_{4}\\ {\text{Cit}}_{3}\\ {\text{Cit}}_{4}\end{array}\right]=\left[\begin{array}{ccc}-{\text{v}}_{8}& .& .\\ .& .& .\\ .& -{\text{v}}_{8}& .\\ .& .& .\\ .& .& .\\ .& .& .\\ .& .& .\\ .& .& -{\text{v}}_{1}\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{Asp}}_{2}\\ {\text{Asp}}_{3}\\ {\text{AcCoA}}_{2}\end{array}\right]$$

$$\left[\begin{array}{ccccc}-{\text{v}}_{8}-{\text{v}}_{6}& {\text{v}}_{6}& .& .& .\\ {\text{v}}_{7}& -{\text{v}}_{7}-{\text{v}}_{5}& {\text{v}}_{5}& .& .\\ .& .& -{\text{v}}_{4}& {\text{v}}_{4}& .\\ .& .& .& -{\text{v}}_{2}& {\text{v}}_{2}\\ .& .& .& .& -{\text{v}}_{1}\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{OAC}}_{23}\\ {\text{Fum}}_{23}\\ {\text{Suc}}_{23}\\ {\text{AKG}}_{34}\\ {\text{Cit}}_{34}\end{array}\right]=\left[\begin{array}{cc}-{\text{v}}_{8}& .\\ .& .\\ .& .\\ .& .\\ .& -{\text{v}}_{1}\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{Asp}}_{23}\\ {\text{OAC}}_{2}\times {\text{AcCoA}}_{2}\end{array}\right]$$

$$\begin{array}{l}\left[\begin{array}{cccccccc}-{\text{v}}_{6}-{\text{v}}_{8}& {\text{v}}_{6}& .& .& .& .& .& .\\ {\scriptstyle \frac{1}{2}}{\text{v}}_{7}& -{\text{v}}_{5}-{\text{v}}_{7}& {\scriptstyle \frac{1}{2}}{\text{v}}_{7}& {\text{v}}_{5}& .& .& .& .\\ .& {\text{v}}_{6}& -{\text{v}}_{6}-{\text{v}}_{8}& .& .& .& .& .\\ .& .& .& -{\text{v}}_{4}& {\scriptstyle \frac{1}{2}}{\text{v}}_{4}& {\scriptstyle \frac{1}{2}}{\text{v}}_{4}& .& .\\ .& .& .& .& -{\text{v}}_{2}& .& {\text{v}}_{2}& .\\ .& .& .& .& .& -{\text{v}}_{2}& .& {\text{v}}_{2}\\ .& .& .& .& .& .& -{\text{v}}_{1}& .\\ .& .& .& .& .& .& .& -{\text{v}}_{1}\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{OAC}}_{123}\\ {\text{Fum}}_{123}\\ {\text{OAC}}_{234}\\ {\text{Suc}}_{123}\\ {\text{AKG}}_{234}\\ {\text{AKG}}_{345}\\ {\text{Cit}}_{234}\\ {\text{Cit}}_{345}\end{array}\right]\\ =\left[\begin{array}{cccc}-{\text{v}}_{8}& .& .& .\\ .& .& .& .\\ .& -{\text{v}}_{8}& .& .\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ .& .& -{\text{v}}_{1}& .\\ .& .& .& -{\text{v}}_{1}\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{Asp}}_{123}\\ {\text{Asp}}_{234}\\ {\text{OAC}}_{23}\times {\text{AcCoA}}_{2}\\ {\text{OAC}}_{2}\times {\text{AcCoA}}_{12}\end{array}\right]\\ \left[\begin{array}{ccc}-{\text{v}}_{3}& {\text{v}}_{3}& .\\ .& -{\text{v}}_{2}& {\text{v}}_{2}\\ .& .& -{\text{v}}_{1}\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{Glu}}_{12345}\\ {\text{AKG}}_{12345}\\ {\text{Cit}}_{12345}\end{array}\right]\\ =\left[\begin{array}{c}.\\ .\\ -{\text{v}}_{1}\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[{\text{OAC}}_{123}\times {\text{AcCoA}}_{12}\right]\end{array}$$

Rotational symmetric molecules fumarate and succinate were modeled as described in section 2.7, and the following groups of EMUs were identified as equivalent: Fum_{2} and Fum_{3}; Suc_{2} and Suc_{3}; Fum_{123} and Fum_{234}; Suc_{123} and Suc_{234}. Note that the EMU model only contains EMUs of size 5 and smaller, even though citrate contains 6 carbon atoms. This is direct consequence of the EMU decomposition, and clearly illustrates the advantage of the EMU framework. The decomposition will never generate EMUs of larger size than the simulated metabolites, i.e. in this case 5 (= size of Glu_{12345}). In contrast, the isotopomer and cumomer methods used 2^{6}=64 variables for citrate. The 24 EMUs constituted a significant reduction from the complete set of 176 isotopomers/cumomers that were required to describe this system (a reduction of 86%). The cumomer model consisted of seven subproblems of size 6, 28, 53, 52, 28, 8, and 1, respectively. As expected, all three modeling methods, i.e. EMU, isotopomer, and cumomer, simulated identical mass isotopomer abundances for glutamate: 34.64 mol% (M+0), 26.95 mol% (M+1), 27.08 mol% (M+2), 8.07 mol% (M+3), 2.86 mol% (M+4), and 0.39 mol% (M+5).

In section 2.6 we noted that the computational effort for solving EMU balances depends on the number of EMU nodes. It is often possible to reduce the number of EMU variables in the EMU networks by eliminating EMU nodes with a single influx. Note that no information is lost in this process. We applied this strategy to simplify the EMU networks for the TCA cycle example. The reduced EMU networks are shown in Figure 14. In this case, the number of EMUs was reduced from 24 to 9 EMUs, i.e. 95% reduction compared to the complete set of 176 isotopomers. The corresponding EMU balances are shown below.

$$\left[\begin{array}{ccc}-{\text{v}}_{8}-{\text{v}}_{6}& {\text{v}}_{6}& 0\\ {\scriptstyle \frac{1}{2}}{\text{v}}_{7}+{\scriptstyle \frac{1}{2}}{\text{v}}_{5}& -{\text{v}}_{5}-{\text{v}}_{7}& {\scriptstyle \frac{1}{2}}{\text{v}}_{7}\\ 0& {\text{v}}_{6}& -{\text{v}}_{8}-{\text{v}}_{6}\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{OAC}}_{2}\\ {\text{Fum}}_{2}\\ {\text{OAC}}_{3}\end{array}\right]=\left[\begin{array}{ccc}-{\text{v}}_{8}& 0& 0\\ 0& 0& -{\scriptstyle \frac{1}{2}}{\text{v}}_{5}\\ 0& -{\text{v}}_{8}& 0\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{Asp}}_{2}\\ {\text{Asp}}_{3}\\ {\text{AcCoA}}_{2}\end{array}\right]$$

$$\left[\begin{array}{cc}-{\text{v}}_{8}-{\text{v}}_{6}& {\text{v}}_{6}\\ {\text{v}}_{7}& -{\text{v}}_{5}-{\text{v}}_{7}\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{OAC}}_{23}\\ {\text{Fum}}_{23}\end{array}\right]=\left[\begin{array}{cc}-{\text{v}}_{8}& 0\\ 0& -{\text{v}}_{5}\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{Asp}}_{23}\\ {\text{OAC}}_{2}\times {\text{AcCoA}}_{2}\end{array}\right]$$

$$\left[\begin{array}{ccc}-{\text{v}}_{8}-{\text{v}}_{6}& {\text{v}}_{6}& 0\\ {\scriptstyle \frac{1}{2}}{\text{v}}_{7}& -{\text{v}}_{5}-{\text{v}}_{7}& {\scriptstyle \frac{1}{2}}{\text{v}}_{7}\\ 0& {\text{v}}_{6}& -{\text{v}}_{8}-{\text{v}}_{6}\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{OAC}}_{123}\\ {\text{Fum}}_{123}\\ {\text{OAC}}_{234}\end{array}\right]=\left[\begin{array}{cccc}-{\text{v}}_{8}& 0& 0& 0\\ 0& 0& -{\scriptstyle \frac{1}{2}}{\text{v}}_{5}& -{\scriptstyle \frac{1}{2}}{\text{v}}_{5}\\ 0& -{\text{v}}_{8}& 0& 0\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{c}{\text{Asp}}_{123}\\ {\text{Asp}}_{234}\\ {\text{OAC}}_{23}\times {\text{AcCoA}}_{2}\\ {\text{OAC}}_{3}\times {\text{AcCoA}}_{12}\end{array}\right]$$

$${\text{Glu}}_{12345}={\text{OAC}}_{123}\times {\text{AcCoA}}_{12}$$

Note that with the reduced EMU model we can simulate the labeling of glutamate in this system for any steady-state fluxes and any substrate labeling by solving four very simple linear problems of size 3, 2, 3, and 1, respectively. The solutions to the EMU balances for the assumed steady-state fluxes and labeling of acetyl-CoA are shown below. As expected, the simulated MID of glutamate was identical to that obtained with the full EMU model and the isotopomer and cumomer models. Table 6 summarizes the comparison of modeling methods for the TCA cycle example.

Solution of EMU balances for reaction network of EMU size 1

$$\left[\begin{array}{c}{\text{OAC}}_{2}\\ {\text{Fum}}_{2}\\ {\text{OAC}}_{3}\end{array}\right]={\left[\begin{array}{ccc}-175& 125& 0\\ 62{\scriptstyle \frac{1}{2}}& -125& 37{\scriptstyle \frac{1}{2}}\\ 0& 125& -175\end{array}\right]}^{-1}\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{ccc}-50& 0& 0\\ 0& 0& -25\\ 0& -50& 0\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{cc}1& 0\\ 1& 0\\ 0.5& 0.5\end{array}\right]=\left[\begin{array}{cc}0.8333& 0.1667\\ 0.7667& 0.2333\\ 0.8333& 0.1667\end{array}\right]$$

Solution of EMU balances for reaction network of EMU size 2

$$\left[\begin{array}{c}{\text{OAC}}_{23}\\ {\text{Fum}}_{23}\end{array}\right]={\left[\begin{array}{cc}-175& 125\\ 75& -125\end{array}\right]}^{-1}\xb7\hspace{0.17em}\left[\begin{array}{cc}-50& 0\\ 0& -50\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{ccc}1& 0& 0\\ 0.4167& 0.5000& 0.0833\end{array}\right]=\left[\begin{array}{ccc}0.7083& 0.2500& 0.0417\\ 0.5917& 0.3500& 0.0583\end{array}\right]$$

Solution of EMU balances for reaction network of EMU size 3

$$\begin{array}{c}\left[\begin{array}{c}{\text{OAC}}_{123}\\ {\text{Fum}}_{123}\\ {\text{OAC}}_{234}\end{array}\right]={\left[\begin{array}{ccc}-175& 125& 0\\ 37{\scriptstyle \frac{1}{2}}& -125& 37{\scriptstyle \frac{1}{2}}\\ 0& 125& -175\end{array}\right]}^{-1}\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{cccc}-50& 0& 0& 0\\ 0& 0& -25& -25\\ 0& -50& 0& 0\end{array}\right]\hspace{0.17em}\xb7\hspace{0.17em}\left[\begin{array}{cccc}1& 0& 0& 0\\ 1& 0& 0& 0\\ 0.3542& 0.4792& 0.1458& 0.0208\\ 0.4167& 0.2917& 0.2500& 0.0417\end{array}\right]\\ =\left[\begin{array}{cccc}0.6927& 0.1927& 0.0990& 0.0156\\ 0.5698& 0.2698& 0.1385& 0.0219\\ 0.6927& 0.1927& 0.0990& 0.0156\end{array}\right]\end{array}$$

Solution of EMU balance for reaction network of EMU size 5

$$[{\text{Glu}}_{12345}]=[\begin{array}{cccccc}0.3464& 0.2695& 0.2708& 0.0807& 0.0286& 0.0039\end{array}]$$

In this example we applied the EMU framework to a realistic metabolic network model of *E. coli* central carbon metabolism. The network was comprised of 73 reactions (with corresponding carbon atom transitions) of which 17 were assumed reversible, utilizing 76 metabolites with 5 substrates, 65 balanced intracellular metabolites, and 6 products. The network model included reactions for the glycolysis, pentose phosphate pathway, Entner-Doudoroff pathway, TCA cycle, product formation, amphibolic reactions, one-carbon metabolism, and amino acid biosynthesis reactions. For this network we simulated the mass isotopomer distributions of 26 amino acid fragments that can be measured experimentally by GC/MS. Table 7 provides an overview of the simulated amino acid fragments. In this case, the network model was decomposed into 14 decoupled EMU reaction networks of EMU size 1 to 9. Table 8 summarizes the EMU decomposition. The largest EMU subnetwork was EMU size 1 subnetwork (i.e. one-carbon EMUs), which contained 141 EMUs. The total number of EMUs was 307, which was reduced to 223 EMUs by eliminating EMU nodes with a single influx. In comparison, there were 4,612 isotopomers/cumomers required to simulate the same amino acid fragments, i.e. a reduction of 93–95%. It is interesting to note that there were 241 carbon atoms in this network model, but only 141 EMUs of size 1 were required for the simulation model. Thus, clearly not all individual carbon atoms needed to be simulated in this network. In contrast, the cumomer method required balancing of all 241 weight-1 cumomers. The simulated mass isotopomer distributions from all three methods (i.e. EMU, isotopomer, cumomer) were identical (results not shown).

Selected ion fragments of TBDMS derivatized amino acids that were simulated using the EMU and isotopomer methods.

In this final example we consider the pathway of gluconeogenesis shown in Figure 15. We constructed a detailed biochemical network model for this pathway, where we considered all carbon, hydrogen, and oxygen atom transitions. This pathway is suitable for probing with multiple isotopic tracers, i.e. ^{13}C, ^{2}H, and ^{18}O. In this example, the goal was to simulate the mass isotopomer distribution of glucose. Glucose is the main product of gluconeogenesis and is easily analyzed by GC/MS. The gluconeogenesis network model was comprised of 24 reactions utilizing 21 metabolites, with 5 substrates (oxaloacetate, glycerol, glycogen, NADH, and water), 14 balanced intracellular metabolites, and 2 products (glucose, and CO_{2}). Table 9 shows the number of carbon, hydrogen, and oxygen atoms for each of the 21 metabolites in this system. Note that we only considered stable, i.e. carbon-bound, hydrogen atoms for each metabolite (see section 1.2). Simulation of this system using the isotopomer and cumomer methods is impossible, because that would require 2,637,120 variables. With the EMU approach, however, the network was decomposed into 60 decoupled EMU reaction networks of EMU size 1 to 19, with 493 total EMUs, which was further reduced to 354 EMUs by eliminating EMU nodes with a single influx. Table 10 shows the complete list of EMU networks generated from EMU decomposition for this example. Note that the largest EMU network contained only 12 EMUs (9 in the reduced EMU model). With this model, simulation of the mass isotopomer distribution of glucose for any fluxes and labeling input took less than 0.1 sec (on a Pentium III 1.6 GHz computer). Thus, we have reduced the computational complexity of simulating the gluconeogenesis pathway from an impossible problem to solve to a set of linear subproblems that are trivial to solve. In Table 11 we compare the EMU method vs. the isotopomer/cumomer methods, considering alternative labeling strategies with combinations of ^{13}C, ^{2}H, and ^{18}O tracers. In all cases the EMU method was superior compared to the isotopomer/cumomer methods.

Metabolites in the gluconeogenesis pathway. For each metabolite we considered solely stable, i.e. carbon-bound, hydrogen atoms.

Complete list of EMUs generated for the simulation model of glucose from EMU network decomposition of the gluconeogenesis pathway. The mass isotopomer distribution of glucose, including all carbon, hydrogen, and oxygen atoms was simulated. This required **...**

Metabolic flux analysis (MFA) provides key parameters for quantifying physiology in fields ranging from metabolic engineering to the analysis of human metabolic disease. The mathematical centerpiece of MFA is the simulation model that describes labeling dynamics under transient or stationary conditions. An important limitation of MFA as carried out via stable-isotope labeling and stable-isotope measurements is the large number of isotopomer equations that need to be solved, especially when multiple isotopic tracers are used for the labeling of the system. The isotopomer/cumomer modeling framework is a generic top-down modeling strategy. It provides the most detailed description of the labeling state of a system given by the isotopomer fractions of all metabolites. It has been the assumption that this description of a system’s labeling state is required to interpret isotope data (Wiechert and Wurzel, 2001). However, alternative modeling methods have been developed for specific isotope measurements and for specific input of labeling. For example, it is well known that fractional enrichments of carbon atoms can be simulated efficiently using atom mapping matrices, the method that was originally proposed by Zukpe and Stephanopoulos (Zupke and Stephanopoulos, 1994). More recently, Van Winden et al. (van Winden et al., 2002) developed the concept of bondomers that allows efficient simulation of NMR fine spectra and MS data without the use of isotopomers. However, the bondomer method is valid only for ^{13}C-labeling experiments where a single uniformly ^{13}C-labeled substrate is applied. If multiple carbon sources are present, then all substrates need to be uniformly ^{13}C-labeled with the same enrichment. This requirement significantly limits the applicability of the bondomer method.

In this contribution we have developed a novel approach for modeling isotopomer distributions that significantly reduces the size of the computational problem without any loss of information. This approach is valid for any stable-isotope measurement and any labeling input. The elementary metabolite unit (EMU) framework is a bottom-up modeling approach and is based on a highly efficient decomposition algorithm that identifies the minimum amount of information that is required for isotopic simulations. We have shown that for realistic metabolic networks EMU models are orders-of-magnitude smaller than corresponding isotopomer and cumomer models, especially when multiple isotopic tracers are applied. As such, the EMU approach is the preferred method for simulating MS and NMR measurements for metabolic flux and network analysis. Important applications of the EMU method include the analysis of metabolic pathways using a combination of multiple isotopic tracers, e.g. gluconeogenesis pathway in section 3.5, and the analysis of nonstationary systems. Due to the significant model-size reduction the EMU framework allows efficient analysis of dynamic systems. For example, the time required for dynamic simulations of the *E. coli* network is on the order of seconds for the EMU model, compared to an hour for the corresponding cumomer model. These and other applications of the EMU approach will be explored in more detail in subsequent studies.

We acknowledge the support of the DuPont-MIT Alliance and the NIH Metabolomics Roadmap Initiative DK070291.

**Publisher's Disclaimer: **This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

- Anderson DH. Compartmental modeling and tracer kinetics. Springer-Verlag; New York: 1983.
- Antoniewicz MR, Kelleher JK, Stephanopoulos G. Determination of confidence intervals of metabolic fluxes estimated from stable isotope measurements. Metab Eng. 2006;8:324–337. [PubMed]
- Brunengraber H, Kelleher JK, Des RC. Applications of mass isotopomer analysis to nutrition research. Annu Rev Nutr. 1997;17:559–596. [PubMed]
- Hellerstein MK. In vivo measurement of fluxes through metabolic pathways: the missing link in functional genomics and pharmaceutical research. Annu Rev Nutr. 2003;23:379–402. [PubMed]
- Moss GP. Basic terminology of stereochemistry. Pure Appl Chem. 1996;68:2193–2222.
- Schmidt K, Carlsen KM, Nielsen J, Villadsen J. Modeling isotopomer distributions in biochemical networks using isotopomer mapping matrices. Biotechnol Bioeng. 1997;55:831–840. [PubMed]
- Stephanopoulos G. Metabolic fluxes and metabolic engineering. Metab Eng. 1999;1:1–11. [PubMed]
- Szyperski T. Biosynthetically directed fractional 13C-labeling of proteinogenic amino acids. An efficient analytical tool to investigate intermediary metabolism. Eur J Biochem. 1995;232:433–448. [PubMed]
- Szyperski T. 13C-NMR, MS and metabolic flux balancing in biotechnology research. Q Rev Biophys. 1998;31:41–106. [PubMed]
- van Winden WA, Heijnen JJ, Verheijen PJ. Cumulative bondomers: a new concept in flux analysis from 2D [13C,1H] COSY NMR data. Biotechnol Bioeng. 2002;80:731–745. [PubMed]
- Wiechert W. 13C metabolic flux analysis. Metab Eng. 2001;3:195–206. [PubMed]
- Wiechert W, Mollney M, Isermann N, Wurzel M, de Graaf AA. Bidirectional reaction steps in metabolic networks: III. Explicit solution and analysis of isotopomer labeling systems. Biotechnol Bioeng. 1999;66:69–85. [PubMed]
- Wiechert W, Siefke C, deGraaf AA, Marx A. Bidirectional reaction steps in metabolic networks: II. Flux estimation and statistical analysis. Biotechnology and Bioengineering. 1997;55:118–135. [PubMed]
- Wiechert W, Wurzel M. Metabolic isotopomer labeling systems. Part I: global dynamic behavior. Math Biosci. 2001;169:173–205. [PubMed]
- Zupke C, Stephanopoulos G. Modeling of isotope distributions and intracellular fluxes in metabolic networks using atom mapping matrices. Biotechnology Progress. 1994;10:489–498.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |