Metabolic Flux Analysis—Linking Isotope Labeling and Metabolic Fluxes

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.


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

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-13 C 2 -Glucose (Glucose where C-1 and C-2 are labeled with 13 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 13 C labeling), M+2 (2 carbon atoms are labeled with 13 C) and M+4 (4 carbon atoms are labeled with 13 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 1 ; 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 2 and f 3 ; 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 3 ; 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][2][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 4 and f 5 ; 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 3 ) is; the higher the observed FBP M+0, the higher the reverse flux of TPI (f 5 ) is. Such reversibility of aldolase and TPI would not have been measurable without using an isotope labeled glucose tracer.   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.
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.
Our first observation is that Equation (1) is a set of linear equations. We can re-write them as  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.
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.
Our first observation is that Equation (1) is a set of linear equations. We can re-write them as Metabolites 2020, 10, 447

of 21
The second observation is, due to the mass balance of FBP, the sum of the production fluxes equals the consumption flux Therefore, we can replace the consumption flux of FBP with the sum of the production fluxes, and Equation (3) can be written as 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 1 /f 3 , 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: If we know DHAP⊕GAP, Glc and FBP, we can plug these values into Equation (6) and solve for the ratio between f 1 and f 3 . In fact, all three equations of M+0, M+2 and M+4 fractions will give the same f 1 /f 3 value. It should be noted that it is not the absolute value of f 1 and f 3 but the ratio of them affects the equation. In fact, f 1 can take any value to fulfill Equation (6) as long as the ratio of f 1 /f 3 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 In the above equations, FBP [1][2][3] and FBP [4][5][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][2][3] and FBP [4][5][6] are used only to calculate DHAP and GAP. The labeling of FBP [1][2][3][4][5][6] is not the convolution of FBP [1][2][3] and FBP [4][5][6]. The labeling balance of DHAP provides us the ratio between f 2 and f 5 , which are the two fluxes making DHAP. The labeling balance of GAP provides us with the ratio between f 2 and f 4 , which are the two fluxes making GAP. Equations (6)-(8) can be combined as Note that Equation (9)   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 m. If each substrate i has the labeling pattern of m i and is contributing to the product through flux f i , then 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 (m i ) equals to that of the product (m). Therefore, the matrix (m i − 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-13 C glucose tracer is used, all metabolites would be fully labelled with 13 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.

Generalization of the Labeling Balance Equations
In order to calculate the labeling of DHAP and GAP, we introduced the terms FBP [1][2][3] and FBP [4][5][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 6 = 64 different forms in which FBP can be 13 C-labeled. Using the EMU approach, we are only concerned with three EMUs for FBP: FBP [1][2][3], FBP [4][5][6] and FBP [1][2][3][4][5][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 2 COOH groups, one derived from acetyl-CoA and the other one from oxaloacetate. If the acetyl-CoA is 13 C-labeled, the 13 C 2 -citrate becomes a chiral molecule. This suggests that citrate is a pro-chiral molecule. Therefore, even though the two -CH 2 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 Metabolites 2020, 10, 447 7 of 21 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.
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.
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.
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.

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, 12 C/ 13 C KIE is typically < 1.1 so this assumption is largely valid in carbon labeling experiments. In deuterium ( 2 H)-labeling experiments, however, strong KIE (>2-4) may come in play. The KIE makes the subtraction of a 2 H atom slower than the subtraction of a 1 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.

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 1 is measured as 100 nmol/h, our goal is to find a set of fluxes f 2 -f 6 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.
Due to the mass balance of FBP, DHAP and GAP, we have another set of equations.
Since f 1 = 100. There are only two free fluxes f 3 and f 5 . All other fluxes can be calculated from these two free fluxes because of the mass balance.
To solve the FBP [1][2][3], FBP [4][5][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.
Therefore, by given the values of two free fluxes f 3 = 50 and f 5 = 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 3 and f 5 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].
Metabolites 2020, 10, x FOR PEER REVIEW 10 of 21 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. 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.

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 3 = 50 and f 5 = 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 3 = 48 and f 5 = 210 would also be an acceptable solution since its SSR (χ 2 = 1.33) was lower than the threshold assuming 1% error on all labeled fractions (χ 2 (α = 0.95,df = 1) = 3.84) [81]. In fact, if we plot all the acceptable solutions in the coordinate system of f 3 and f 5 , 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 3 and f 5 . 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-13 C 2 glucose is an effective tracer to determine the free fluxes f 3 and f 5 in the upper glycolysis network ( Figure 3A). In contrast, if we use 50% U-13 C-glucose + 50% unlabeled glucose, we can only determine f 3 but not f 5 ( Figure 3B). This observation motivates us to investigate the principle behind tracer selection to achieve the best effectiveness of MFA. 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.  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-13 C glucose and 50% unlabeled glucose was used as the tracer, the labeling patterns of FBP [1][2][3] and GAP became identical (50% M+0; 50% M+3). Therefore, no matter how f 5 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 3 and f 5 . Thus, the confidence region of the free fluxes is a vertical belt ( Figure 3B). In contrast, the value of f 3 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 2 and f 4 are the two free fluxes in this model ( Figure 4A). We are using 2-13 C 1 -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 2 and f 4 are 16.4-28.3 and 67-405, respectively. The overall confidence region plotted on the f 2 -f 4 plane is a belt from lower left to the upper right ( Figure 4B). The direction of the belt means when f 2 gets larger, f 4 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 2 and f 4 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 2 and f 4 ( Figure 4C). In this case, these two fluxes are anti-correlated. When f 2 gets larger, f 4 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].

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, CO2 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 CO2. Consequently, no flux information can be derived [91]. One solution to such a problem is kinetic flux The colored region is showing the 99% confidence region, in order to enlarge the region to be more visible, determined by measuring the labeling patterns of both B and F. The discrepancy between predicted labeling pattern and the measurement is converted to the χ 2 -statistic assuming 0.3% error on the labeled fractions. 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].

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 2 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 2 . 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 ∆M represents the net change in a metabolite from t = 0 to t = T. ∆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 m represents the labeling pattern of the product. When expressed in vector forms, the ∆M, m i and 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 Furthermore, under metabolic steady state, the intracellular metabolite pool size does not change. With a constant metabolite concentration, we also have In Equation (26), c is the intracellular concentration of the metabolite and V is the cell volume.
Combining Equations (25) and (26), we have We can also take derivative of Equation (27), which gives 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, dm/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 Taking the derivative of Equation (29), we have Without assuming constant metabolite concentration or cell volume, the overall change of the labeling of metabolite pool is ∆M = ∆(cVm) Taking the derivative of it, we have The metabolite pool size is only affected by the net production flux Taking Equation (33) Divide both sides with cV, we have 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).

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.
Author Contributions: Y.W. and X.S. wrote the main draft. F.E.W., C.S. and T.Z. revised the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This research is supported, in part, by NIH grants P30CA072720-5923 (XS).

Conflicts of Interest:
The authors declare no conflict of interest.

Code
Availability: All the code used to generate graphs can be accessed at 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.

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.
Author Contributions: Y.W. and X.S. wrote the main draft. F.E.W., C.S. and T.Z. revised the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This research is supported, in part, by NIH grants P30CA072720-5923 (XS).