Metabolic Flux Analysis—Linking Isotope Labeling and Metabolic Fluxes

Abstract

Metabolic flux analysis (MFA) is an increasingly important tool to study metabolism quantitatively. Unlike the concentrations of metabolites, the fluxes, which are the rates at which intracellular metabolites interconvert, are not directly measurable. MFA uses stable isotope labeled tracers to reveal information related to the fluxes. The conceptual idea of MFA is that in tracer experiments the isotope labeling patterns of intracellular metabolites are determined by the fluxes, therefore by measuring the labeling patterns we can infer the fluxes in the network. In this review, we will discuss the basic concept of MFA using a simplified upper glycolysis network as an example. We will show how the fluxes are reflected in the isotope labeling patterns. The central idea we wish to deliver is that under metabolic and isotopic steady-state the labeling pattern of a metabolite is the flux-weighted average of the substrates’ labeling patterns. As a result, MFA can tell the relative contributions of converging metabolic pathways only when these pathways make substrates in different labeling patterns for the shared product. This is the fundamental principle guiding the design of isotope labeling experiment for MFA including tracer selection. In addition, we will also discuss the basic biochemical assumptions of MFA, and we will show the flux-solving procedure and result evaluation. Finally, we will highlight the link between isotopically stationary and nonstationary flux analysis.

1. Introduction

Metabolism involves a set of chemical reactions that convert nutrients to energy and macromolecules to sustain cellular survival and growth. Many enzymes are thought to work on production lines: the product coming from the first enzyme is the substrate for the second enzyme. Together, these enzymes form metabolic pathways that fulfill certain metabolic tasks. The catabolic pathways, for example, break down large molecules such as proteins and lipids to make ATP and other small molecule metabolites. The anabolic pathways use biochemical energy stored in ATP to synthesize macromolecules. Interestingly, the catabolism and anabolism are overlapped rather than exclusive. Some enzymes, such as the fructose-bisphosphate aldolase, are shared between the catabolic glycolysis pathway and the anabolic gluconeogenesis pathway. Meanwhile, co-factors such as ATP, acetyl-CoA and NAD(PH) are shared among many enzymatic reactions [1]. These features demonstrate the fact that these biochemical reactions are highly interconnected and essentially form a network. The elucidation of the metabolic network was a huge scientific achievement in the 20th century [2]. Despite the vast knowledge on the structure of the metabolic network, the dynamic partition of metabolic material through different pathways under various physiological and pathological conditions remains incompletely characterized.

In the past two decades, metabolic research has been rapidly expanded due to the huge improvement on technologies for metabolite measurements. These technological advances include nuclear magnetic resonance (NMR) and mass spectrometry (MS) [3,4,5,6,7,8,9]. The large-scale analysis of the metabolite pool sizes, namely metabolomics, has been widely deployed in clinical research and pharmacology [10]. Metabolites can serve as biomarkers for diagnosis. For example, elevated spermine level and decreased citrate level in the prostate fluid has been shown to be a distinct metabolomics profile in prostate cancer patients [11,12]. Although these discoveries are highly valuable for diagnosis, the measurement of metabolite pool size alone reveals little information regarding the biochemical events behind such changes [13,14]. The increased pool size of a metabolite can be the result of either enhanced production or diminished consumption. For example, sildenafil is an inhibitor of phosphodiesterase type 5, which breaks down cyclic guanosine monophosphate (cGMP). Usage of sildenafil increases the pool size of cGMP but the flux through this metabolite is decreased [15]. This example demonstrates that changes in metabolite pool size do not directly translate into flux changes. An analogy for the relationship between the metabolite pool size and metabolic flux is the car traffic. The metabolite pool size represents the density of cars on the road while the flux represents the flow of traffic. Although the density of cars can be measured by taking a static picture of the road, such a measurement alone is insufficient to tell how fast the flow of traffic is.

Unlike the pool size of metabolites that can be directly measured by NMR and MS, the metabolic fluxes are usually unmeasurable. Although the extracellular fluxes such as nutrient uptake rates can be calculated based on the concentration changes over time, the intracellular fluxes cannot be measured this way. Instead, they have to be solved by the approach named flux analysis [16,17,18,19]. Flux analysis is a powerful technique and has been applied in multiple fields of research such as pharmaceutical production, drug design and cancer biology [20,21,22]. For example, about one century ago, Warburg et al. [23] discovered cancer cells generate more lactate than healthy cells under aerobic conditions, a phenomenon now well known as the Warburg effect [24]. However, it remains unclear for many decades how cancer cells could make use of lactate to support its own growth. Recently, with the help of metabolic flux analysis, it became clear that lactate is an important carbon source for citric acid cycle in both healthy and cancer cells [25,26].

Depending on what experimental data are available, the flux analysis can be carried out in different ways. However, all flux analysis methods must resolve one basic problem: the total number of unknown fluxes often exceed the total number of measurement (i.e., the network is undetermined). One way to resolve this issue is by simplifying the flux network and using stoichiometry balance of the intracellular metabolites to reduce the total number of unknown fluxes. This approach is known as stoichiometric flux analysis (SFA) [16,27,28,29]. On the other hand, if a detailed flux network has to be used, one could also resolve this issue by assuming a cellular objective function (e.g., maximum cell growth or ATP production). This approach is known as the flux balance analysis (FBA) [17,30,31,32,33,34]. In addition to satisfying the stoichiometric constraints, the intracellular fluxes must also maximize the objective function.

Despite the popular application of FBA and SFA in system biology and metabolic engineering, these methods have some apparent limitations: (1) it is hard to estimate the fluxes through parallel pathways, cyclic pathways or bidirectional reversible reactions [17,35,36]. (2) The objective function for mammalian cell metabolism is hard to define. Mammalian cells, unlike bacteria in a batch culture, can exhibit different metabolic phenotypes during proliferation, differentiation and senescence [37]. These practical limitations largely come from the limited amount of measurements that are utilized in FBA and SFA. In contrast, metabolic flux analysis (MFA) uses stable isotope tracers and takes labeling patterns of intracellular metabolites as additional input information. While more commonly referred to as ¹³C-MFA, the framework of MFA is not limited to ¹³C labeling. In fact, MFA has been applied to analyze labeling from ¹⁵N, ²H and other tracers [14,38]. As a result, MFA is much more powerful than FBA and SFA in overcoming the two limitations mentioned above [17].

Table 1. Commonly used software tools for metabolic flux analysis (MFA).

Name	Data Source	Main Features	Number of Citations	Reference and Most Cited Applications
13CFLUX2	MS and NMR	Compatible to multi-platform	107	[39,40,41,42]
FiatFlux	GC-MS	13C Glucose tracer, flux ratio	160	[43,44,45,46]
INCA	MS and NMR	Isotopically non-stationary MFA	147	[47,48,49,50,51]
METRAN	MS	Intuitive graphical user interface, confidence interval calculation	183	[52,53,54]
NMR2Flux+	NMR	NMR data, plant network	124	[55,56,57,58]
OpenFLUX	MS	Steady-state 13C MFA, experimental design	154	[59,60,61,62]
OpenMebius	MS	Isotopically non-stationary MFA	47	[63,64,65]

has listed some most commonly used software tools for MFA. In this review, we will demonstrate how MFA is carried out, and we will discuss the principles of using MFA to solve metabolic fluxes.

2. Example on Upper Glycolysis

Without using isotope tracers, we have a very limited insight on how metabolites interconvert in the metabolic network. Taking upper glycolysis pathway as an example. If glucose is consumed at the rate of 100 nmol/h, glyceraldehyde-3-phosphate (GAP) goes to lower glycolysis at the rate of 200 nmol/h (Figure 1A). No more information on the metabolic fluxes can be drawn from this unlabeled experiment. In comparison, when 1,2-¹³C₂-Glucose (Glucose where C-1 and C-2 are labeled with ¹³C) tracer is used (Figure 1B), we can observe the reversibility of upper glycolysis reactions. The fructose 1,6-bisphosphate (FBP) will be labeled in three forms, M+0 (no ¹³C labeling), M+2 (2 carbon atoms are labeled with ¹³C) and M+4 (4 carbon atoms are labeled with ¹³C) (Figure 1B). FBP can be made from the glucose tracer through the intermediates glucose-6-phopshate (Glc6P) and fructose-6-phopshate (Fruc6P). Both intermediates would have the same labeling (M+2) as the glucose tracer since their production only involves phosphorylation and isomerization without breaking the C-C bonds. Therefore, the M+2 FBP represents its synthesis flux mostly directly from Fruc6P (f₁; Figure 1B). The observation of M+0 and M+4 FBP, however, suggests the presence of additional synthesis route for FBP, which involves recombination of the carbon backbones. Indeed, the aldolase reaction that breaks down FBP into dihydroxyacetone phosphate (DHAP) and GAP is reversible (f₂ and f₃; Figure 1B). The observed M+4 FBP has to be generated from M+2 DHAP and M+2 GAP through the reverse reaction of aldolase (f₃; Figure 1B). Similarly, M+0 FBP has to be generated from M+0 DHAP and M+0 GAP. However, since the DHAP derived from Glucose C1–C3 (Glc[1-3]) must be M+2 labeled, there has to be an alternative way to generate M+0 DHAP. In fact, the triose phosphate isomerase (TPI) reaction is also reversible (f₄ and f₅; Figure 1B) converting M+0 DHAP from M+0 GAP. Thus, using isotopic tracers, we have established a qualitative link between the labeling pattern of FBP and the reversibility of the aldolase and TPI reactions. The lower the observed FBP M+2, the higher the reverse flux of aldolase (f₃) is; the higher the observed FBP M+0, the higher the reverse flux of TPI (f₅) is. Such reversibility of aldolase and TPI would not have been measurable without using an isotope labeled glucose tracer.

Although the real metabolic network might be more complicated than we illustrated (e.g., the existence of pentose phosphate pathway), this simple example illustrates the gist of MFA. When we use the suitable isotope labeled tracer, the labeling patterns of intercellular metabolites reflect the fluxes. The goal of MFA is to infer the metabolic fluxes from these isotope labeling patterns. Before we move to the solving process of MFA, we need to examine quantitatively how the labeling patterns are determined by the fluxes and what information is encoded in the labeling patterns.

When the cells are at metabolic steady state, the production rate of an intercellular metabolite equals the consumption rates of it. Furthermore, when the cells are under isotopic steady state, the production rate of a specific labeling fraction equals the consumption rates of it. For FBP, such mass balances on the labeled fractions are described by the following equation.

$$\begin{array}{c}\{\begin{array}{c}{\mathrm{f}}_3\times {\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+0}+{\mathrm{f}}_1\times \mathrm{G}\mathrm{l}{\mathrm{c}}_{\mathrm{M}+0}={\mathrm{f}}_2\times \mathrm{F}\mathrm{B}{\mathrm{P}}_{\mathrm{M}+0}\\ {\mathrm{f}}_3\times {\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+2}+{\mathrm{f}}_1\times \mathrm{G}\mathrm{l}{\mathrm{c}}_{\mathrm{M}+2}={\mathrm{f}}_2\times \mathrm{F}\mathrm{B}{\mathrm{P}}_{\mathrm{M}+2}\\ {\mathrm{f}}_3\times {\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+4}+{\mathrm{f}}_1\times \mathrm{G}\mathrm{l}{\mathrm{c}}_{\mathrm{M}+4}={\mathrm{f}}_2\times \mathrm{F}\mathrm{B}{\mathrm{P}}_{\mathrm{M}+4}\end{array}\end{array}$$

(1)

In Equation (1), we use DHAP⊕GAP to represent the combination of GAP and DHAP. In such a condensation reaction, the substrates are combined independently in terms of their labeling patterns. Therefore, we have the following equations.

$$\begin{array}{c}{\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+0}=\mathrm{D}\mathrm{H}\mathrm{A}{\mathrm{P}}_{\mathrm{M}+0}\times \mathrm{G}\mathrm{A}{\mathrm{P}}_{\mathrm{M}+0}\\ {\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+2}=\mathrm{D}\mathrm{H}\mathrm{A}{\mathrm{P}}_{\mathrm{M}+2}\times \mathrm{G}\mathrm{A}{\mathrm{P}}_{\mathrm{M}+0}+\mathrm{D}\mathrm{H}\mathrm{A}{\mathrm{P}}_{\mathrm{M}+0}\times \mathrm{G}\mathrm{A}{\mathrm{P}}_{\mathrm{M}+2}\\ {\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+4}=\mathrm{D}\mathrm{H}\mathrm{A}{\mathrm{P}}_{\mathrm{M}+2}\times \mathrm{G}\mathrm{A}{\mathrm{P}}_{\mathrm{M}+2}\end{array}$$

(2)

Our first observation is that Equation (1) is a set of linear equations. We can re-write them as

$$\begin{array}{c}{\mathrm{f}}_3\times \left[\begin{array}{c}{\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+0}\\ {\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+2}\\ {\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+4}\end{array}\right]+{\mathrm{f}}_1\times \left[\begin{array}{c}{\mathrm{Glc}}_{\mathrm{M}+0}\\ {\mathrm{Glc}}_{\mathrm{M}+2}\\ {\mathrm{Glc}}_{\mathrm{M}+4}\end{array}\right]={\mathrm{f}}_2\times \left[\begin{array}{c}{\mathrm{FBP}}_{\mathrm{M}+0}\\ {\mathrm{FBP}}_{\mathrm{M}+2}\\ {\mathrm{FBP}}_{\mathrm{M}+4}\end{array}\right]\end{array}$$

(3)

The second observation is, due to the mass balance of FBP, the sum of the production fluxes equals the consumption flux

$$\begin{array}{c}{\mathrm{f}}_1+{\mathrm{f}}_3={\mathrm{f}}_2\end{array}$$

(4)

Therefore, we can replace the consumption flux of FBP with the sum of the production fluxes, and Equation (3) can be written as

$$\begin{array}{c}{\mathrm{f}}_3\times \left[\begin{array}{c}{\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+0}\\ {\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+2}\\ {\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+4}\end{array}\right]+{\mathrm{f}}_1\times \left[\begin{array}{c}{\mathrm{Glc}}_{\mathrm{M}+0}\\ {\mathrm{Glc}}_{\mathrm{M}+2}\\ {\mathrm{Glc}}_{\mathrm{M}+4}\end{array}\right]=\left({\mathrm{f}}_3+{\mathrm{f}}_1\right)\times \left[\begin{array}{c}{\mathrm{FBP}}_{\mathrm{M}+0}\\ {\mathrm{FBP}}_{\mathrm{M}+2}\\ {\mathrm{FBP}}_{\mathrm{M}+4}\end{array}\right]\end{array}$$

(5)

Equation (5) shows that the labeling of FBP is the flux-weighted average of its substrates’ labeling (i.e., the higher the ratio of f₁/f₃, the similar the labeling pattern of FBP to that of Glc). If we move the FBP labeling vector to the left side of the equation, we have:

$$\begin{array}{c}\left[\begin{array}{c}\begin{array}{cc}\begin{array}{c}{\mathrm{Glc}}_{\mathrm{M}+0}-{\mathrm{FBP}}_{\mathrm{M}+0}\\ {\mathrm{Glc}}_{\mathrm{M}+2}-{\mathrm{FBP}}_{\mathrm{M}+2}\\ {\mathrm{Glc}}_{\mathrm{M}+4}-{\mathrm{FBP}}_{\mathrm{M}+4}\end{array}& \begin{array}{c}{\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+0}-{\mathrm{FBP}}_{\mathrm{M}+0}\\ {\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+2}-{\mathrm{FBP}}_{\mathrm{M}+2}\\ {\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+4}-{\mathrm{FBP}}_{\mathrm{M}+4}\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{\mathrm{f}}_1\\ {\mathrm{f}}_3\end{array}\right]=0\end{array}$$

(6)

If we know DHAP⊕GAP, Glc and FBP, we can plug these values into Equation (6) and solve for the ratio between f₁ and f₃. In fact, all three equations of M+0, M+2 and M+4 fractions will give the same f₁/f₃ value. It should be noted that it is not the absolute value of f₁ and f₃ but the ratio of them affects the equation. In fact, f₁ can take any value to fulfill Equation (6) as long as the ratio of f₁/f₃ remains unchanged. This observation brings up the essence of labeling balance equations. The labeling balance equations describe the ratio of converging fluxes.

Similar to the FBP labeling equations, based on the labeling balance of DHAP and GAP we have the following equations

$$\begin{array}{c}\left[\begin{array}{c}\begin{array}{cc}\begin{array}{c}\mathrm{FBP}{\left[1\text{-}3\right]}_{\mathrm{M}+0}-{\mathrm{DHAP}}_{\mathrm{M}+0}\\ \mathrm{FBP}{\left[1\text{-}3\right]}_{\mathrm{M}+2}-{\mathrm{DHAP}}_{\mathrm{M}+2}\end{array}& \begin{array}{c}{\mathrm{GAP}}_{\mathrm{M}+0}-{\mathrm{DHAP}}_{\mathrm{M}+0}\\ {\mathrm{GAP}}_{\mathrm{M}+2}-{\mathrm{DHAP}}_{\mathrm{M}+2}\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{\mathrm{f}}_2\\ {\mathrm{f}}_5\end{array}\right]=0\end{array}$$

(7)

$$\begin{array}{c}\left[\begin{array}{c}\begin{array}{cc}\begin{array}{c}\mathrm{FBP}{\left[4\text{-}6\right]}_{\mathrm{M}+0}-{\mathrm{GAP}}_{\mathrm{M}+0}\\ \mathrm{FBP}{\left[4\text{-}6\right]}_{\mathrm{M}+2}-{\mathrm{GAP}}_{\mathrm{M}+2}\end{array}& \begin{array}{c}{\mathrm{DHAP}}_{\mathrm{M}+0}-{\mathrm{GAP}}_{\mathrm{M}+0}\\ {\mathrm{DHAP}}_{\mathrm{M}+2}-{\mathrm{GAP}}_{\mathrm{M}+2}\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{\mathrm{f}}_2\\ {\mathrm{f}}_4\end{array}\right]=0\end{array}$$

(8)

In the above equations, FBP[1-3] and FBP[4-6] denote the labeling patterns of the first and the last three carbon atoms of FBP, which are converted to DHAP and GAP, respectively, by the aldolase reaction. It is noteworthy that FBP[1-3] and FBP[4-6] are used only to calculate DHAP and GAP. The labeling of FBP[1-6] is not the convolution of FBP[1-3] and FBP[4-6]. The labeling balance of DHAP provides us the ratio between f₂ and f₅, which are the two fluxes making DHAP. The labeling balance of GAP provides us with the ratio between f₂ and f₄, which are the two fluxes making GAP. Equations (6)–(8) can be combined as

$$\begin{array}{c} \left[\begin{array}{c}\begin{array}{ccccc}\begin{array}{c}{\mathrm{Glc}}_{\mathrm{M}+0}-{\mathrm{FBP}}_{\mathrm{M}+0}\\ {\mathrm{Glc}}_{\mathrm{M}+2}-{\mathrm{FBP}}_{\mathrm{M}+2}\\ {\mathrm{Glc}}_{\mathrm{M}+4}-{\mathrm{FBP}}_{\mathrm{M}+4}\end{array}& \begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{c}{\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+0}-{\mathrm{FBP}}_{\mathrm{M}+0}\\ {\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+2}-{\mathrm{FBP}}_{\mathrm{M}+2}\\ {\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+4}-{\mathrm{FBP}}_{\mathrm{M}+4}\end{array}& \begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 0\end{array}\\ \begin{array}{c}0\\ 0\end{array}& \begin{array}{c}\mathrm{FBP}{\left[1\text{-}3\right]}_{\mathrm{M}+0}-{\mathrm{DHAP}}_{\mathrm{M}+0}\\ \mathrm{FBP}{\left[1\text{-}3\right]}_{\mathrm{M}+2}-{\mathrm{DHAP}}_{\mathrm{M}+2}\end{array}& \begin{array}{c}0\\ 0\end{array}& \begin{array}{c}0\\ 0\end{array}& \begin{array}{c}{\mathrm{GAP}}_{\mathrm{M}+0}-{\mathrm{DHAP}}_{\mathrm{M}+0}\\ {\mathrm{GAP}}_{\mathrm{M}+2}-{\mathrm{DHAP}}_{\mathrm{M}+2}\end{array}\\ \begin{array}{c}0\\ 0\end{array}& \begin{array}{c}\mathrm{FBP}{\left[4\text{-}6\right]}_{\mathrm{M}+0}-{\mathrm{GAP}}_{\mathrm{M}+0}\\ \mathrm{FBP}{\left[4\text{-}6\right]}_{\mathrm{M}+2}-{\mathrm{GAP}}_{\mathrm{M}+2}\end{array}& \begin{array}{c}0\\ 0\end{array}& \begin{array}{c}{\mathrm{DHAP}}_{\mathrm{M}+0}-{\mathrm{GAP}}_{\mathrm{M}+0}\\ {\mathrm{DHAP}}_{\mathrm{M}+2}-{\mathrm{GAP}}_{\mathrm{M}+2}\end{array}& \begin{array}{c}0\\ 0\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{\mathrm{f}}_1\\ {\mathrm{f}}_2\\ {\mathrm{f}}_3\\ {\mathrm{f}}_4\\ {\mathrm{f}}_5\end{array}\right]=0\end{array}$$

(9)

Note that Equation (9) is a set of homogeneous linear equations. All fluxes are unitless and can be scaled up or down without violating the equation. Only when the absolute value of f₁ is measured by other experiments, these fluxes can be expressed in absolute values as well.

The general form of the labeling balance equation should be apparent by now. Suppose we have n fluxes making the same metabolite, which has the labeling pattern of $\overline{\mathrm{m}}$. If each substrate i has the labeling pattern of ${\mathrm{m}}_i$ and is contributing to the product through flux f_i, then

$$\begin{array}{c}{\displaystyle {\displaystyle \sum }_{i=1}^n}{f}_i\left({\mathrm{m}}_i-\overline{\mathrm{m}}\right)=0\end{array}$$

(10)

Equation (10) illustrates several fundamental principles of MFA. First, the use of isotope labeling is required to solve fluxes. While this statement seems obvious, it naturally emerges from Equation (10). If no tracer is used, then the labeling pattern of all substrates (${\mathrm{m}}_i$) equals to that of the product ($\overline{\mathrm{m}}$). Therefore, the matrix $({\mathrm{m}}_i-\overline{\mathrm{m}})$ contains only zero vectors and $f_i$ can take any value to satisfy the equation. For the same reason, if the two substrates of one metabolite are identically labeled (e.g., when the U-¹³C glucose tracer is used, all metabolites would be fully labelled with ¹³C), the fluxes from these substrates still cannot be solved. This is the principle guiding the tracer selection in MFA, which will be discussed later in more detail. Second, at the isotope labeling steady state, the labeling pattern of a metabolite is determined only by its production fluxes and is independent of its consumption fluxes. The consumption flux of one metabolite shows up in the labeling balance equations only when it becomes the production flux of another metabolite. Third, metabolite concentration is not in this equation, meaning that under isotopic steady state the fluxes are independent of the metabolite pool size. This is the concept we discussed in the introduction. It is also a natural consequence of the steady-state labeling balance equation. We will show later that under the isotopic non-steady state condition the metabolite pool size is related with fluxes.

3. Generalization of the Labeling Balance Equations

In order to calculate the labeling of DHAP and GAP, we introduced the terms FBP[1-3] and FBP[4-6] in Equations (7) and (8), which are parts of the FBP molecule. These terms are commonly referred to as elementary metabolite units (EMU). EMU is defined as a moiety comprising any distinct subset of a compound’s atoms [66]. The use of EMU is a natural choice when calculating the labeling pattern of metabolites synthesized by combining the carbon backbones. EMU significantly simplifies the flux calculation. Taking FBP as an example, there are totally 2⁶ = 64 different forms in which FBP can be ¹³C-labeled. Using the EMU approach, we are only concerned with three EMUs for FBP: FBP[1-3], FBP[4-6] and FBP[1-6]. Without this simplification process, the calculation of fluxes would become much more difficult. For example, if all the possible EMU should be described, Equation (3) would have to use 128 unknown parameters to describe the labeling of DHAP⊕GAP and FBP, and 64 known parameters to describe the labeling of glucose tracer. The overall workload becomes unnecessarily large.

EMU should also account for the molecular symmetry in the metabolites. In citric acid cycle, for example, succinate is a symmetric molecule. The first and fourth carbon of succinate (Suc [1] and Suc [4]) should always have the same labeling pattern. A sufficiently large “flipping” flux that interconverts Suc [1234] and Suc [4321] can be introduced to the flux model to account for such symmetry. It is also noteworthy that such flipping flux should not be applied to pro-chiral metabolites. Prochirality of a metabolite is often overlooked. In fact, prochirality is easy to identify in the isotope labeling experiment. A molecule is pro-chiral if isotope labeling can convert it from achiral to chiral. Citrate, for example, has two -CH₂COOH groups, one derived from acetyl-CoA and the other one from oxaloacetate. If the acetyl-CoA is ¹³C-labeled, the ¹³C₂-citrate becomes a chiral molecule. This suggests that citrate is a pro-chiral molecule. Therefore, even though the two -CH₂COOH groups are symmetric with respect to a mirror plane, they are not interconvertible.

EMU is the natural choice to describe the labeling pattern of metabolites when they are measured in the forms of mass isotopomer distributions (MIDs) by mass spectrometry. However, if NMR instead of MS is used, the labeling pattern would be described more easily in the form of positional labeling. One such approach is known as cumomers [67]. The cumomer approach is only concerned with the percentage that FBP and other metabolites that are labeled at specific carbon atoms, regardless of the labeling at other positions. In this case, we can rewrite Equations (6)–(8), which describes the labeling balance in EMU terms, using the cumomer approach as the following.

$$\begin{array}{c}\left[\begin{array}{cc}\begin{array}{c}\mathrm{Glc}\left[2\right]-\mathrm{FBP}\left[2\right]\\ \mathrm{Glc}\left[5\right]-\mathrm{FBP}\left[5\right]\end{array}& \begin{array}{c}\mathrm{DHAP}\left[2\right]-\mathrm{FBP}\left[2\right]\\ \mathrm{GAP}\left[2\right]-\mathrm{FBP}\left[5\right]\end{array}\end{array}\right]\left[\begin{array}{c}{\mathrm{f}}_1\\ {\mathrm{f}}_3\end{array}\right]=0\end{array}$$

(11)

$$\left[\begin{array}{cc}\mathrm{FBP}\left[2\right]-\mathrm{DHAP}\left[2\right]& \mathrm{GAP}\left[2\right]-\mathrm{DHAP}\left[2\right]\end{array}\right]\left[\begin{array}{c}{\mathrm{f}}_2\\ {\mathrm{f}}_5\end{array}\right]=0$$

(12)

$$\left[\begin{array}{cc}\mathrm{FBP}\left[5\right]-\mathrm{GAP}\left[2\right]& \mathrm{DHAP}\left[2\right]-\mathrm{GAP}\left[2\right]\end{array}\right]\left[\begin{array}{c}{\mathrm{f}}_2\\ {\mathrm{f}}_4\end{array}\right]=0$$

(13)

In these equations, we are concerned with FBP labeling on C-2 and C-5 position, which is FBP $010000 and FBP $000010 in the cumomer terms. This example demonstrates that EMUs and cumomers follow exactly the same labeling balance equations. They are only different in how the labeling patterns are expressed.

Another way to express the isotope labeling is by enrichment. Enrichment is defined as the percentage of labeling on average atom level. The enrichment of FBP (E_FBP) can be calculated from its labeling pattern.

$$\begin{array}{c}{\mathrm{E}}_{\mathrm{FBP}}={\left[\begin{array}{c}0/6\\ 1/6\\ 2/6\\ 3/6\\ 4/6\\ 5/6\\ 6/6\end{array}\right]}^{\mathrm{T}}\left[\begin{array}{c}{\mathrm{FBP}}_{\mathrm{M}+0}\\ {\mathrm{FBP}}_{\mathrm{M}+1}\\ {\mathrm{FBP}}_{\mathrm{M}+2}\\ {\mathrm{FBP}}_{\mathrm{M}+3}\\ {\mathrm{FBP}}_{\mathrm{M}+4}\\ {\mathrm{FBP}}_{\mathrm{M}+5}\\ {\mathrm{FBP}}_{\mathrm{M}+6}\end{array}\right]\times 100\%\end{array}$$

(14)

The good property of enrichment is that the enrichment of the product in a combination reaction is the carbon number-weighted average of the enrichment of the substrates.

$$\begin{array}{c}{\mathrm{E}}_{\mathrm{DHAP}\oplus \mathrm{GAP}}=\left({\mathrm{E}}_{\mathrm{DHAP}}\times 3+{\mathrm{E}}_{\mathrm{GAP}}\times 3\right)/6\times 100\%\end{array}$$

(15)

Isotopic enrichment is commonly used for radioactive tracers. For stable isotope tracers, this is also a convenient term because the labeling vector is converted to a single number. Nonetheless, the labeling balance follows Equation (10) when using the enrichment to represent the labeling patterns.

4. Basic Assumptions in MFA

Having shown that the labeling balance Equation (10) is the basis for many different versions of MFA, it is time to examine the basic assumptions behind Equation (10) and MFA. The first assumption that we have adopted in this equation is that the fluxes are linear with respect to all the labeled fractions. This means the enzymes and transporters do not discriminate between unlabeled and labeled substrates. In another word, there is no kinetic isotope effect (KIE) [68,69,70,71,72]. As Liuni et al. [69] has shown, ¹²C/¹³C KIE is typically < 1.1 so this assumption is largely valid in carbon labeling experiments. In deuterium (²H)-labeling experiments, however, strong KIE (>2–4) may come in play. The KIE makes the subtraction of a ²H atom slower than the subtraction of a ¹H atom. This rate difference can be measured using isolated enzymes and substrates in a test tube. However, the cellular homeostatic mechanisms may result in a smaller impact on the isotope labelling patterns in cells [73]. One possibility is that the total metabolic fluxes are maintained even when the heavy isotope tracers are used. This is likely the case when the enzyme has remaining capacity at the physiological flux level. Another possibility is that the pathway flux is decreased by the same magnitude that the enzyme reaction in vitro is slowed down due to the heavy-labeled substrate. The actual flux change due to KIE is likely somewhere in between and should be evaluated on a case-by-case basis. To introduce KIE to the labeling balance equation of MFA, the flux coefficients ($f_i$; Equation (10)) should be replaced by row vectors of fluxes specific to each labeled fraction.

The second assumption adopted in Equation (10) is that all the metabolites are well-mixed. We are assuming that each metabolite is labeled homogenously so that its labeling pattern can be expressed in a single vector. This simple assumption has different meanings and implications on different levels in the biological system. On the molecular level, the intermediates in a metabolic pathway can be passed from one enzyme to the next without equilibration with the cellular pool [74]. This process is known as substrate channeling [75,76,77]. The assumption that the metabolites are well-mixed means there is no substrate channeling. On the cellular level, a metabolite can be distributed in different cellular compartments, such as cytosol and mitochondrion [78,79,80]. The well-mixed assumption means either the metabolite is exclusively in one cellular compartment, or there is fast exchange among the compartments so that the labeling patterns are averaged out. On the organ level, we have a mixture of different types of cells. We are assuming these cells share the same labeling patterns of metabolites, or there is fast exchange among cells. When these well-mixed assumptions are violated, separate pools of metabolites should be considered when designing MFA models to address the compartmentalization or substrate channeling issues.

5. Predicting Labeling Patterns and Solving Metabolic Fluxes

Now we are ready to solve the fluxes for our example of upper glycolysis pathway. Since the value of f₁ is measured as 100 nmol/h, our goal is to find a set of fluxes f₂–f₆ that generates a labeling pattern of FBP that best matches the observation. Even though the labeling patterns of DHAP and GAP are not measurable, their labeling patterns can still be predicted along with that of FBP when a set of fluxes are given. Based on the mass balance on each of the labeled fractions of DHAP, GAP and FBP, we have the following equations.

$$\begin{array}{c}\{\begin{array}{c}{\mathrm{f}}_2\times \mathrm{F}\mathrm{B}\mathrm{P}\left[1\text{-}3\right]+{\mathrm{f}}_5\times \mathrm{G}\mathrm{A}\mathrm{P}=\left({\mathrm{f}}_2+{\mathrm{f}}_5\right)\times \mathrm{D}\mathrm{H}\mathrm{A}\mathrm{P}\\ {\mathrm{f}}_2\times \mathrm{F}\mathrm{B}\mathrm{P}\left[4\text{-}6\right]+{\mathrm{f}}_4\times \mathrm{D}\mathrm{H}\mathrm{A}\mathrm{P}=\left({\mathrm{f}}_2+{\mathrm{f}}_4\right)\times \mathrm{G}\mathrm{A}\mathrm{P}\\ {\mathrm{f}}_3\times \mathrm{D}\mathrm{H}\mathrm{A}\mathrm{P}+{\mathrm{f}}_1\times \mathrm{G}\mathrm{l}\mathrm{c}\left[1\text{-}3\right]=({\mathrm{f}}_3+{\mathrm{f}}_1)\times \mathrm{F}\mathrm{B}\mathrm{P}\left[1\text{-}3\right]\\ {\mathrm{f}}_3\times \mathrm{G}\mathrm{A}\mathrm{P}+{\mathrm{f}}_1\times \mathrm{G}\mathrm{l}\mathrm{c}\left[4\text{-}6\right]=({\mathrm{f}}_3+{\mathrm{f}}_1)\times \mathrm{F}\mathrm{B}\mathrm{P}\left[4\text{-}6\right]\end{array}\end{array}$$

(16)

Due to the mass balance of FBP, DHAP and GAP, we have another set of equations.

$$\begin{array}{c}\{\begin{array}{c}{\mathrm{f}}_1+{\mathrm{f}}_3={\mathrm{f}}_2\\ {\mathrm{f}}_2+{\mathrm{f}}_5={\mathrm{f}}_3+{\mathrm{f}}_4\\ {\mathrm{f}}_2+{\mathrm{f}}_4={\mathrm{f}}_3+{\mathrm{f}}_5+{\mathrm{f}}_6\end{array}\end{array}$$

(17)

Since f₁ = 100. There are only two free fluxes f₃ and f₅. All other fluxes can be calculated from these two free fluxes because of the mass balance.

$$\begin{array}{c}\{\begin{array}{c}{\mathrm{f}}_2=100+{\mathrm{f}}_3\\ {\mathrm{f}}_4=100+{\mathrm{f}}_5\\ {\mathrm{f}}_6=200\end{array}\end{array}$$

(18)

Take Equation (18) into Equation (16), we have the following equation.

$$\begin{array}{c}(100+{\mathrm{f}}_3)\times \mathrm{F}\mathrm{B}\mathrm{P}\left[1\text{-}3\right]+{\mathrm{f}}_5\times \mathrm{G}\mathrm{A}\mathrm{P}=\left(100+{\mathrm{f}}_3+{\mathrm{f}}_5\right)\times \mathrm{D}\mathrm{H}\mathrm{A}\mathrm{P}\\ \left(100+{\mathrm{f}}_3\right)\times \mathrm{F}\mathrm{B}\mathrm{P}\left[4\text{-}6\right]+(100+{\mathrm{f}}_5)\times \mathrm{D}\mathrm{H}\mathrm{A}\mathrm{P}=\left(200+{\mathrm{f}}_3+{\mathrm{f}}_5\right)\times \mathrm{G}\mathrm{A}\mathrm{P}\\ {\mathrm{f}}_3\times \mathrm{D}\mathrm{H}\mathrm{A}\mathrm{P}+100\times \mathrm{G}\mathrm{l}\mathrm{c}\left[1\text{-}3\right]=\left(100+{\mathrm{f}}_3\right)\times \mathrm{F}\mathrm{B}\mathrm{P}\left[1\text{-}3\right]\\ {\mathrm{f}}_3\times \mathrm{G}\mathrm{A}\mathrm{P}+100\times \mathrm{G}\mathrm{l}\mathrm{c}\left[4\text{-}6\right]=\left(100+{\mathrm{f}}_3\right)\times \mathrm{F}\mathrm{B}\mathrm{P}\left[4\text{-}6\right]\end{array}$$

(19)

Equation (19) illustrates how we can calculate labeling patterns from the free fluxes. For example, if we take the values of f₃ = 50 and f₅ = 150, the four unknowns of FBP[1-3], FBP[4-6], DHAP and GAP become solvable through Equation (20). Since we used the tracer of 1,2-¹³C₂-glucose, $\mathrm{Glc}{\left[1\text{-}3\right]}_{\mathrm{M}+2}$ and $\mathrm{Glc}{\left[4\text{-}6\right]}_{\mathrm{M}+0}$ are 1. All other forms of Glc[1-3] and Glc[4-6] are 0.

$$\begin{array}{c}\begin{array}{c}\left[\begin{array}{cccc}150& 0& -300& 150\\ 0& 150& 250& -400\\ -150& 0& 50& 0\\ 0& -150& 0& 50\end{array}\right]\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}\mathrm{FBP}\left[1\text{-}3\right]\\ \mathrm{FBP}\left[4\text{-}6\right]\end{array}\\ \mathrm{DHAP}\\ \mathrm{GAP}\end{array}\end{array}\right]=\left[\begin{array}{cc}0& 0\\ 0& 0\\ -100& 0\\ 0& -100\end{array}\right]\left[\begin{array}{c}\begin{array}{c}\mathrm{Glc}\left[1\text{-}3\right]\\ \mathrm{Glc}\left[4\text{-}6\right]\end{array}\end{array}\right]\end{array}\\ \begin{array}{c} =\left[\begin{array}{cc}0& 0\\ 0& 0\\ -100& 0\\ 0& -100\end{array}\right]\left[\begin{array}{c}\begin{array}{cc}\mathrm{Glc}{\left[1\text{-}3\right]}_{\mathrm{M}+0}& \mathrm{Glc}{\left[1\text{-}3\right]}_{\mathrm{M}+2}\\ \mathrm{Glc}{\left[4\text{-}6\right]}_{\mathrm{M}+0}& \mathrm{Glc}{\left[4\text{-}6\right]}_{\mathrm{M}+2}\end{array}\end{array}\right]\end{array}\\ \begin{array}{c} =\left[\begin{array}{cc}0& 0\\ 0& 0\\ -100& 0\\ 0& -100\end{array}\right]\left[\begin{array}{c}\begin{array}{cc}0& 1\\ 1& 0\end{array}\end{array}\right]\end{array}\end{array}$$

(20)

To solve the FBP[1-3], FBP[4-6], DHAP and GAP, we simply need to take the inverse of the left matrix containing the fluxes and multiply it with the other two matrices on the right.

$$\begin{array}{c}\begin{array}{c}\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}\mathrm{FBP}\left[1\text{-}3\right]\\ \mathrm{FBP}\left[4\text{-}6\right]\end{array}\\ \mathrm{DHAP}\\ \mathrm{GAP}\end{array}\end{array}\right]=\left[\begin{array}{c}\begin{array}{cc}\begin{array}{c}\mathrm{FBP}{\left[1\text{-}3\right]}_{\mathrm{M}+0}\\ \mathrm{FBP}{\left[4\text{-}6\right]}_{\mathrm{M}+0}\end{array}& \begin{array}{c}\mathrm{FBP}{\left[1\text{-}3\right]}_{\mathrm{M}+2}\\ \mathrm{FBP}{\left[4\text{-}6\right]}_{\mathrm{M}+2}\end{array}\\ {\mathrm{DHAP}}_{\mathrm{M}+0}& {\mathrm{DHAP}}_{\mathrm{M}+2}\\ {\mathrm{GAP}}_{\mathrm{M}+0}& {\mathrm{GAP}}_{\mathrm{M}+2}\end{array}\end{array}\right]\end{array}\\ \begin{array}{c}= {\left[\begin{array}{cccc}150& 0& -300& 150\\ 0& 150& 250& -400\\ -150& 0& 50& 0\\ 0& -150& 0& 50\end{array}\right]}^{-1}\left[\begin{array}{cc}0& 0\\ 0& 0\\ -100& 0\\ 0& -100\end{array}\right]\left[\begin{array}{c}\begin{array}{cc}0& 1\\ 1& 0\end{array}\end{array}\right]\end{array}\\ \begin{array}{c} =\left[\begin{array}{cc}0.1000& 0.9000\\ 0.8333& 0.1667\\ 0.3000& 0.7000\\ 0.5000& 0.5000\end{array}\right]\end{array}\end{array}$$

(21)

Therefore, we can calculate FBP labeling from Equation (3).

$$\begin{array}{c}\begin{array}{c}\mathrm{FBP}=[{\mathrm{f}}_3\times \left(\mathrm{D}\mathrm{H}\mathrm{A}\mathrm{P}\oplus \mathrm{G}\mathrm{A}\mathrm{P}\right)+{\mathrm{f}}_1\times \mathrm{G}\mathrm{l}\mathrm{c}]/{\mathrm{f}}_2\end{array}\\ \begin{array}{c}={\mathrm{f}}_3/{\mathrm{f}}_2\times \left[\begin{array}{c}{\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+0}\\ {\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+2}\\ {\left(\mathrm{DHAP}\oplus \mathrm{GAP}\right)}_{\mathrm{M}+4}\end{array}\right]+{\mathrm{f}}_1/{\mathrm{f}}_2\times \left[\begin{array}{c}{\mathrm{Glc}}_{\mathrm{M}+0}\\ {\mathrm{Glc}}_{\mathrm{M}+2}\\ {\mathrm{Glc}}_{\mathrm{M}+4}\end{array}\right]\\ \end{array}\\ \begin{array}{c}=\frac{{\mathrm{f}}_3}{{\mathrm{f}}_2}\times \left[\begin{array}{c}0.1500\\ 0.5000\\ 0.3500\end{array}\right]+\frac{{\mathrm{f}}_1}{{\mathrm{f}}_2}\times \left[\begin{array}{c}0.0000\\ 1.0000\\ 0.0000\end{array}\right]=\left[\begin{array}{c}0.0500\\ 0.8333\\ 0.1167\end{array}\right]\end{array}\end{array}$$

(22)

Therefore, by given the values of two free fluxes f₃ = 50 and f₅ = 150, the labeling pattern of FBP can be predicted as 5.00% M+0, 83.33% M+2 and 11.67% M+4. The real MFA process, however, is the reverse of this question (Figure 2A): if the actual labeling pattern of FBP is given, can we find the value of f₃ and f₅ that best fit this observation? Since in most cases there is no direct way to calculate the free fluxes by simply plugging in the numbers, the best-fit fluxes are often solved as an optimization problem. The flux-solving process is essentially a “guessing game”, which is iterative and consists of two steps: (1) calculate the labeling patterns of FBP based on the “guessed flux”, and (2) update the “guessed flux” to reduce the discrepancy between the predicted FBP labeling pattern and the observation. The discrepancy (R) can be defined as the sum of squared residuals (SSR), and it is often inversely weighted by the experimental standard deviations expressed as the measurement covariance matrix with measurement variances located on the diagonal (Σ_Measured) [81].

$$\begin{array}{c}\mathrm{R}={\left({\left[\begin{array}{c}{\mathrm{FBP}}_{\mathrm{M}+0}\\ {\mathrm{FBP}}_{\mathrm{M}+2}\\ {\mathrm{FBP}}_{\mathrm{M}+4}\end{array}\right]}_{\mathrm{Measured}}-{\left[\begin{array}{c}{\mathrm{FBP}}_{\mathrm{M}+0}\\ {\mathrm{FBP}}_{\mathrm{M}+2}\\ {\mathrm{FBP}}_{\mathrm{M}+4}\end{array}\right]}_{\mathrm{Prediction}}\right)}^{\mathrm{T}}{\mathrm{\Sigma }}_{\mathrm{Measured}}^{-1}\left({\left[\begin{array}{c}{\mathrm{FBP}}_{\mathrm{M}+0}\\ {\mathrm{FBP}}_{\mathrm{M}+2}\\ {\mathrm{FBP}}_{\mathrm{M}+4}\end{array}\right]}_{\mathrm{Measured}}-{\left[\begin{array}{c}{\mathrm{FBP}}_{\mathrm{M}+0}\\ {\mathrm{FBP}}_{\mathrm{M}+2}\\ {\mathrm{FBP}}_{\mathrm{M}+4}\end{array}\right]}_{\mathrm{Prediction}}\right)\\ ={\displaystyle {\displaystyle \sum }_{i=0,2,4}}{\left(\frac{{\mathrm{FBP}}_{\mathrm{Measured}}^{\mathrm{M}+\mathrm{i}}-{\mathrm{FBP}}_{\mathrm{Prediction}}^{\mathrm{M}+\mathrm{i}}}{\mathrm{SD}\left({\mathrm{FBP}}_{\mathrm{Measured}}^{\mathrm{M}+\mathrm{i}}\right)}\right)}^2\end{array}$$

(23)

Once the discrepancy is sufficiently small or no better guesses can be made, the optimized flux combination is our best estimation. Figure 2B–D demonstrated such a process. By starting from a random flux combination, the gradient descent algorithm can find the direction to minimize the discrepancy between the labeling measurements and the model prediction. When the discrepancy cannot be further minimized, the optimization completes, and the optimal flux combination is the solution.

From this example, we can see that the calculation of the labeling patterns is the most important part in the flux determination. The labeling calculation gives us the landscape, such as shown in Figure 2B, which describes the difference between the model prediction and the experimental measurements when the flux combination takes different values. The goal of the optimization is to find the path on the landscape to reach the lowest point. There are many generic algorithms that can work for such a purpose, including gradient descent, differential evolution and particle swarm optimization algorithms.

6. Evaluation of MFA Result and Tracer Selection

In the example of upper glycolysis, if we start from the FBP measurement of 5.00% M+0, 83.33% M+2 and 11.67% M+4, we can land on the flux combination of f₃ = 50 and f₅ = 150. This solution is both unique and precise. No other flux combination can generate exactly the same FBP labeling pattern. However, in the real experiments, because of the existence of measurement errors, we may not have a flux combination that perfectly fit the measurements. The optimization algorithms can help us to achieve a locally best solution. However, it is hard to know whether this is the best situation globally. One workaround is to run the optimization algorithm multiple times with a random starting point in each round. If the solutions converge to one flux combination, we have confidence that this is the globally best solution (Figure 2B). If multiple solutions are reported, we can compare the results and choose the best one.

In addition to the best flux combination, we want to know the uncertainty in each flux to fully evaluate the results. The flux uncertainty is commonly represented as the confidence intervals of individual fluxes. When the flux takes a sub-optimal value within the confidence interval, the SSR goes up but is still within the statistically acceptable range. For example, if we assume the measurement of the MID of FBP has 1% error, then the combination of f₃ = 48 and f₅ = 210 would also be an acceptable solution since its SSR (χ² = 1.33) was lower than the threshold assuming 1% error on all labeled fractions (χ²(α = 0.95,df = 1) = 3.84) [81]. In fact, if we plot all the acceptable solutions in the coordinate system of f₃ and f₅, an area of confidence interval is generated (Figure 3A). In this case, the area of confidence is in an oval shape ranging from 43.5 to 57.5 and from 75 to 330 on the x and y axis, respectively. These ranges are also the confidence intervals for f₃ and f₅. A narrow confidence interval on a flux indicates that the SSR is sensitive to the change of this flux. Therefore, this flux can be determined more precisely. On the other hand, a wide confidence interval suggests the SSR is insensitive to the change of this flux and the flux is poorly determined. It should be noted that the range of the confidence interval is determined by the precision of the MID measurements and the sensitivity of the SSR to the fluxes. The flux sensitivity is a local property. When the optimal flux fit changes, the flux confidence interval should be re-evaluated.

By drawing the area of confidence intervals, we can also evaluate the tracer selection for MFA. The choice of the tracer has tremendous impact on the effectiveness of MFA. For example, 1,2-¹³C₂ glucose is an effective tracer to determine the free fluxes f₃ and f₅ in the upper glycolysis network (Figure 3A). In contrast, if we use 50% U-¹³C- glucose + 50% unlabeled glucose, we can only determine f₃ but not f₅ (Figure 3B). This observation motivates us to investigate the principle behind tracer selection to achieve the best effectiveness of MFA.

Essentially, the principle for tracer selection is illustrated in Equation (10). The role of the tracer in MFA is to make the metabolite labeling patterns reflective of the fluxes. More specifically, if one metabolite can be made by two or more reactions, the labeling pattern of this metabolite is the flux-weighted average of the labeling of the substrates. This allows us to determine the ratio of the converging fluxes. However, if the fluxes of two reactions utilize substrates with the same labeling pattern, the product can be only made in the same labeling pattern, and no flux ratio can be determined. In the example of upper glycolysis, when 50% U-¹³C glucose and 50% unlabeled glucose was used as the tracer, the labeling patterns of FBP[1-3] and GAP became identical (50% M+0; 50% M+3). Therefore, no matter how f₅ changes, the labeling pattern of DHAP would remain the same. In another word, the labeling pattern of DHAP can no longer reflect the ratio between f₃ and f₅. Thus, the confidence region of the free fluxes is a vertical belt (Figure 3B). In contrast, the value of f₃ can still be determined. This is because the labeling patterns of the two substrates glucose and DHAP⊕GAP are still different. On the network level, we can assess tracer effectiveness by calculating the sensitivity. The good tracer makes the labeling patterns sensitive to the changes of the flux inputs. Multiple studies have been conducted to evaluate the performance of different tracers on evaluating fluxes in glycolysis [82,83], the pentose phosphate pathway [83,84,85] and the TCA cycle [83,86].

It should be noted that the confidence interval of each flux is not independent from others in the more general cases. In Figure 4, we are demonstrating the flux confidence intervals for the metabolic network presented in the EMU paper [66]. f₂ and f₄ are the two free fluxes in this model (Figure 4A). We are using 2-¹³C₁-A as the tracer, and the MID measurement is made on metabolite F. If we allow 0.3% error on each of the labeled fractions, the confidence intervals of f₂ and f₄ are 16.4-28.3 and 67–405, respectively. The overall confidence region plotted on the f₂–f₄ plane is a belt from lower left to the upper right (Figure 4B). The direction of the belt means when f₂ gets larger, f₄ needs to go larger to minimize the increase of the SSR. The confidence intervals of these two free fluxes are highly correlated. The wide confidence intervals of f₂ and f₄ suggest the fluxes are poorly determined when given the MID measurement of F. This result is seemingly odd because we only have two free fluxes to determine, yet four numbers were measured (M+0–M+3) from the MID of F. The reason that these fluxes are still poorly determined is because the four measurements from the MID of F are not totally independent. Similarly, if we only measure the MID of B in this metabolic network, we also have wide confidence intervals on f₂ and f₄ (Figure 4C). In this case, these two fluxes are anti-correlated. When f₂ gets larger, f₄ needs to go smaller to minimize the increase of the SSR. This observation suggests that we need to measure the MIDs of both B and F to achieve best precision of MFA (Figure 4D). This example demonstrates that the identifiability of flux should be assessed seriously. Simply having more fractional numbers in the MIDs than the number of free fluxes does not guarantee identifiability [87].

Figure 4B,C also show that the confidence region is non-linear and non-convex. Therefore, the calculation of confidence intervals is not quite straightforward, especially in large metabolic network with a large number of free fluxes. The method to accurately calculate confidence intervals has been proposed by Antoniewicz et al. [81]. In brief, the confidence intervals were calculated by (1) moving one target flux away from the best-fit value by a small step; (2) choosing a combination of the other fluxes that minimize the increase of SSR; and (3) calculating the new SSR and repeating step (1) to (3) until the new SSR reached the cutoff for 95% confidence interval. An alternative is the Markov Chain Monte Carlo (MCMC) method. We could sample in the flux space and allow the MID and other measurements to have Gaussian error terms. The distribution of the flux values from the perturbed measurements is used to approximate the flux confidence intervals. The advantage of MCMC method is that it can generate a comprehensive probability distribution of fluxes [88,89].

7. Kinetic Flux Profiling and Isotopically Non-Stationary MFA

So far, all the MFA we have discussed are based on the assumptions of steady states. Essentially, the assumptions are (1) the metabolic steady state is the condition that the metabolite pool sizes and fluxes do not change over time and (2) isotopic steady state is the condition that the labeling patterns of metabolites do not change with time [90]. Experimentally, the isotopic steady state is achieved by incubating or infusing the tracer for sufficiently long time so that the labeling patterns of the metabolites do not change anymore. Under this condition, the ratio of the converging fluxes can be determined if the substrates are labeled in different patterns. Moreover, the size of metabolite pools does not affect the labeling patterns. One limitation of isotopically stationary MFA is that it does not work on pathways where no differential labeling can be observed. For example, CO₂ is the sole carbon source of the photosynthetic product in autotrophic tissues of plant. At isotopic steady state, all the photosynthetic intermediates and products would have exactly the same labeling as CO₂. Consequently, no flux information can be derived [91]. One solution to such a problem is kinetic flux profiling (KFP) [92]. In KFP experiments, we can switch the carbon source from unlabeled to the labeled one and examine the labeling patterns of intermediates and products at different time points. The key idea in KFP is that turnover rate of the metabolite pool is determined by the synthesizing flux and the pool size of metabolite. In fact, KFP is the extension of MFA in isotopic non-steady state case. In the final part of this review paper, we wish to highlight the link between steady state MFA and KFP in isotopic non-steady state [92,93,94,95], and further extend it to the general non-steady state case.

For an intracellular metabolite under isotopic non-steady state, the total pool labeling change with regard to all the labeled fractions can be expressed as

$$\begin{array}{c}\Delta \overline{M}={\displaystyle {\displaystyle \sum }_{i=1}^n}\left({{\displaystyle \int }}_0^T{f}_i{m}_{i}Vdt\right)-{{\displaystyle \int }}_0^T{f}_{out}\overline{m}Vdt\end{array}$$

(24)

$\Delta \overline{M}$ represents the net change in a metabolite from t = 0 to t = T. $\Delta \overline{M}$ can be the amount of a specific labeled fraction, or a vector to describe all labeled fractions. $f_i$ represents the production flux from the i-th substrate. $f_{out}$ represents the sum of the consumption fluxes. These fluxes are normalized to cell volume V. $m_i$ represents the labeling pattern of the i-th substrate and $\overline{m}$ represents the labeling pattern of the product. When expressed in vector forms, the $\Delta \overline{M}$, $m_i$ and $\overline{m}$ describe the same set of labeled fractions.

Under metabolic steady state, where fluxes and the cell volume are constant over time, and the consumption flux equals the sum of the production fluxes, we have

$$\begin{array}{c}\Delta \overline{M}={\displaystyle {\displaystyle \sum }_{i=1}^n}\left(f_{i}V{{\displaystyle \int }}_0^T{m}_{i}dt\right)-f_{out}V{{\displaystyle \int }}_0^T\overline{m}dt\\ ={\displaystyle {\displaystyle \sum }_{i=1}^n}\left(f_{i}V{{\displaystyle \int }}_0^T{m}_{i}dt\right)-\left({\displaystyle {\displaystyle \sum }_{i=1}^n}{f}_i\right)V{{\displaystyle \int }}_0^T\overline{m}dt\\ ={\displaystyle {\displaystyle \sum }_{i=1}^n}\left[f_{i}V{{\displaystyle \int }}_0^T\left(m_i-\overline{m}\right)dt\right]\end{array}$$

(25)

Furthermore, under metabolic steady state, the intracellular metabolite pool size does not change. With a constant metabolite concentration, we also have

$$\begin{array}{c}\Delta \overline{M}=cV\Delta \overline{m}\end{array}$$

(26)

In Equation (26), c is the intracellular concentration of the metabolite and V is the cell volume.

Combining Equations (25) and (26), we have

$$\begin{array}{c}\Delta \overline{m}={\displaystyle {\displaystyle \sum }_{i=1}^n}\left[f_i{{\displaystyle \int }}_0^T\left(m_i-\overline{m}\right)dt\right]/c\end{array}$$

(27)

We can also take derivative of Equation (27), which gives

$$\begin{array}{c}\frac{d\overline{m}}{dt}={\displaystyle {\displaystyle \sum }_{i=1}^n}\left[f_i\left(m_i-\overline{m}\right)\right]/c\end{array}$$

(28)

Equation (28) is the equation describing the labeling change with respect to time for KFP. It is also the isotopic non-steady state generalization of Equation (6). When the isotopic steady state is reached, $d\overline{m}/dt$ = 0, and Equation (28) degenerates to Equation (10). This means the isotopic non-steady state and steady state cases are closely related.

Next, we assume neither isotopic steady state nor metabolic steady state. Starting from Equation (24), considering that the sum of the consumption fluxes can be replaced by the sum of production fluxes minus the net production flux, we can re-write Equation (24) to

$$\begin{array}{c}\begin{array}{c}\Delta \overline{M}={\displaystyle {\displaystyle \sum }_{i=1}^n}\left({{\displaystyle \int }}_0^T{f}_{i}V{m}_{i}dt\right)-{{\displaystyle \int }}_0^T\left({\displaystyle {\displaystyle \sum }_{i=1}^n}{f}_{i}V-f_{Net\_Pro}V\right)\overline{m}dt\end{array}\\ \begin{array}{c}={\displaystyle {\displaystyle \sum }_{i=1}^n}\left({{\displaystyle \int }}_0^T{f}_{i}V{m}_{i}dt\right)-{\displaystyle {\displaystyle \sum }_{i=1}^n}\left({{\displaystyle \int }}_0^T{f}_{i}V\overline{m}dt\right)+{{\displaystyle \int }}_0^T{f}_{Net\_Pro}V\overline{m}dt\end{array}\\ \begin{array}{c}={\displaystyle {\displaystyle \sum }_{i=1}^n}\left[{{\displaystyle \int }}_0^T{f}_{i}V\left(m_i-\overline{m}\right)dt\right]+{{\displaystyle \int }}_0^T{f}_{Net\_Pro}V\overline{m}dt\end{array}\end{array}$$

(29)

Taking the derivative of Equation (29), we have

$$\begin{array}{c}\frac{d\overline{M}}{dt}={\displaystyle {\displaystyle \sum }_{i=1}^n}\left[f_i\left(m_i-\overline{m}\right)V\right]+f_{Net\_Pro}\overline{m}V\end{array}$$

(30)

Without assuming constant metabolite concentration or cell volume, the overall change of the labeling of metabolite pool is

$$\begin{array}{c}\Delta \overline{M}=\Delta \left(cV\overline{m}\right)\end{array}$$

(31)

Taking the derivative of it, we have

$$\begin{array}{c}\frac{d\overline{M}}{dt}=\frac{d\left(cV\overline{m}\right)}{dt}=cV\frac{d\overline{m}}{dt}+\overline{m}\frac{d\left(cV\right)}{dt}\end{array}$$

(32)

The metabolite pool size is only affected by the net production flux

$$\begin{array}{c}\Delta \left(cV\right)={{\displaystyle \int }}_0^T{f}_{Net\_Pro}Vdt\\ \frac{d\left(cV\right)}{dt}=f_{Net\_Pro}V\end{array}$$

(33)

Taking Equation (33) into Equation (32), we have

$$\begin{array}{c}\frac{d\overline{M}}{dt}=cV\frac{d\overline{m}}{dt}+f_{Net\_Pro}\overline{m}V\end{array}$$

(34)

Combining Equation (34) and Equation (30), we have

$$\begin{array}{c}{\displaystyle {\displaystyle \sum }_{i=1}^n}\left[f_i\left(m_i-\overline{m}\right)V\right]=cV\frac{d\overline{m}}{dt}\end{array}$$

(35)

Divide both sides with cV, we have

$$\begin{array}{c}\frac{d\overline{m}}{dt}={\displaystyle {\displaystyle \sum }_{i=1}^n}\left[f_i\left(m_i-\overline{m}\right)\right]/c\end{array}$$

(36)

Equation (36) has the same form as Equation (28). The term c in Equation (36) is a function of time, whereas it is a constant in Equation (28). Equation (36) is the general equation describing the flux and the change in the labeling patterns. It applies to biological systems with changing cell numbers and intracellular metabolite concentrations. Similar to the metabolic and isotopic steady state case in Equation (10), the changes in the labeling patterns are largely affected by the input fluxes. The consumption fluxes are not included in the equation. They can only affect the concentration term c to impact the labeling changes. Comparing Equations (10) and (36), we also see that the metabolite concentration under isotopic steady state is not related to the flux and the labeling patterns. However, the concentration term determines how fast the labeling can change under non-steady state condition. The change of pool size in one metabolite will also result in the changes in the labeling kinetics of other metabolites. For example, when the Glc6P pool expands, the FBP labeling also becomes slower. However, when the isotopic steady state is reached, the pool size of Glc6P does not affect the labeling pattern of FBP (Figure 5).

To solve the fluxes under isotopic non-steady state, we also take the simulation-optimization approach. Given a set of fluxes and the metabolite concentrations, we can predict the labeling patterns at different time points of sampling. The model predictions are compared to the experimental measurements. The flux combinations that generate the labeling kinetics that match the measurements the best are the optimal solution for the isotopic non-steady state MFA.

8. Conclusions

In conclusion, MFA uses stable isotope labeling to solve metabolic fluxes. MFA commonly adopts two assumptions: the metabolites are well-mixed, and the kinetic isotope effect is negligible. Under isotopic steady state, one metabolite’s labeling pattern represents the flux-weighted average of the substrates. Under non-steady state, the rate of change of the labeling pattern is determined by the substrate and product labeling patterns, the production fluxes of this metabolite and the concentration of the metabolite.