Modeling isotopic fractionation
A simple network is given in Fig. 1a as an example, with the corresponding carbon atom transitions shown in Fig. 1b. The labeling state of a metabolite containing n carbon atoms is denoted by the coefficients (a, b, …, n), which take a value of 0 if the atom in the corresponding carbon position is 12C or a value of 1 if the atom is 13C. For instance, B
0,1,1 represents the isotopomer of the metabolite B which is unlabeled in position 1 and labeled in positions 2 and 3. In the absence of KIEs, each isotopomer of a given metabolite reacts with a probability equal to its abundance relative to the total metabolite pool. Each isotopomer A
a,b,c
is thus produced at the rate v
1 weighted by \( {\overline{S}}_{a,b,c} \), which denotes the concentration of the isotopomer S
a,b,c
relative to the total S pool (\( {\displaystyle \sum_{a,b,c}}{S}_{a,b,c} \)). Similarly, A
a,b,c
is consumed at the rate v
2 weighted by Ā
a,b,c
, which denotes the concentration of the isotopomer A
a,b,c
relative to the total A pool. Therefore, the balance around the isotopomer A
a,b,c
is:
$$ \frac{d{A}_{a,b,c}}{dt}={v}_1.{\overline{S}}_{a,b,c}-{v}_2.{\overline{A}}_{a,b,c} $$
(4)
where v
1 and v
2 are functions of the concentrations of reactants and effectors, and of the kinetic parameters of the enzymes catalyzing these reactions. This equation can be rewritten using the fluxomer variable introduced by [29]. A fluxomer is the rate at which a metabolic reaction transforms one or several substrate(s) isotopomers into product(s) isotopomers. Fluxomers are usually expressed as rate fractions [28, 29], which allows a direct mapping between fluxomers and isotopomer abundances. However, this is not relevant to calculate isotopomers derivatives. Therefore, fluxomers were expressed as absolute rates rather than rate fractions. In the absence of KIEs, the fluxomer f
a,b,c1
which consumes S
a,b,c
is defined by:
$$ {f}_1^{S_{a,b,c}}={v}_1.{\overline{S}}_{a,b,c} $$
(5)
and the mass balance around A
a,b,c
can be rewritten as:
$$ \frac{d{A}_{a,b,c}}{dt}={f}_1^{S_{a,b,c}}-{f}_2^{A_{a,b,c}} $$
(6)
It must be noted that the term "fluxomer" is somewhat misleading in the context of this study since it does not define a flux (which is a term used for the rate when the reaction is part of a whole system) but rather a local rate (i.e. the rate when considering the reaction isolated from the rest of the system). Bearing that in mind, we keep this terminology for sake of simplicity, but note the important distinction between these terms.
KIEs affect the rate constant of reactions depending on the isotopic composition of their substrate(s). Consequently, KIEs were included by weighting each fluxomer with a coefficient α which represents the KIE of a given enzyme with respect to its substrate(s) isotopomers. For instance, each coefficient \( {\alpha}_{v_2}^{A_{a,b,c}} \) is defined as:
$$ {\alpha}_{v_2}^{A_{a,b,c}}=\frac{v_2^{A_{a,b,c}}}{v_2^{A_{0,0,0}}} $$
(7)
where \( {v}_2^{A_{a,b,c}} \) (\( {v}_2^{A_{0,0,0}} \)) represents the rate of reaction 2 when isotopomer A
a,b,c
(A
0,0,0) is its unique substrate. The mass balance around A
a,b,c
finally becomes:
$$ \frac{d{A}_{a,b,c}}{dt}={\alpha}_{v_1}^{S_{a,b,c}}.{f}_1^{S_{a,b,c}}-{\alpha}_{v_2}^{A_{a,b,c}}.{f}_2^{A_{a,b,c}} $$
(8)
To calculate the mass balances for the isotopomers of pool B, two additional points must be taken into account. First, it is worth recalling that bidirectional isotope exchange that arise from reversibility significantly impacts the distribution of isotopes through the network [31]. Forward and backward reaction rates must be considered separately in the case of reversible reactions, like v
3 in the model considered here. Second, in the case of reactions with more than one substrate, the impact of KIEs on reaction rate may depend on the isotopic composition of each substrate. Therefore, one fluxomer must be defined for each combination of substrates isotopomers. Reaction v
3 was first decomposed into its forward (v
forw3
) and backward (v
back3
) components, each of them being then decomposed into fluxomers:
$$ {f}_{3,\ forw}^{B_{a,b,c}}={v}_3^{forw}.{\overline{B}}_{a,b,c} $$
(9)
$$ {f}_{3,\ back}^{C_a,{D}_{b,c}}={v}_3^{back}.{\overline{C}}_a.{\overline{D}}_{b,c} $$
(10)
and the mass balance around B
a,c,b
is:
$$ \frac{d{B}_{a,c,b}}{dt}={\alpha}_{v_2}^{A_{a,b,c}}.{f}_2^{A_{a,b,c}}+{\alpha}_{v_3}^{C_a}.{\alpha}_{v_3}^{D_{c,b}}.{f}_{3,\ back}^{C_a,{D}_{c,b}}-{\alpha}_{v_3}^{B_{a,c,b}}.{f}_{3,\ forw}^{B_{a,c,b}}-{\alpha}_{v_4}^{B_{a,c,b}}.{f}_4^{B_{a,c,b}} $$
(11)
Using the same formalism, the mass balances around the isotopomers of C, D and E are:
$$ \frac{d{C}_a}{dt}={\displaystyle \sum_{b,c}}\left({\alpha}_{v_3}^{B_{a,b,c}}.{f}_{3,\ forw}^{B_{a,b,c}}\right)-{\displaystyle \sum_{b,c}}\left({\alpha}_{v_3}^{C_a}.{\alpha}_{v_3}^{D_{b,c}}.{f}_{3,\ back}^{C_a,{D}_{b,c}}\right)-{\alpha}_{v_5}^{C_a}.{f}_5^{C_a} $$
(12)
$$ \frac{d{D}_{b,c}}{dt}=\sum_a\left({\alpha}_{v_3}^{B_{a,b,c}}.{f}_{3, forw}^{B_{a,b,c}}\right)-\sum_a\left({\alpha}_{v_3}^{C_a}.{\alpha}_{v_3}^{D_{b,c}}.{f}_{3, back}^{C_a,{D}_{b,c}}\right)-{\alpha}_{v_6}^{D_{b,c}}.{f}_6^{D_{b,c}} $$
(13)
$$ \frac{d{E}_{a,b,c}}{dt}={\alpha}_{v_4}^{B_{a,b,c}}.{f}_4^{B_{a,b,c}}-{\alpha}_{v_7}^{E_{a,b,c}}.{f}_7^{E_{a,b,c}} $$
(14)
The system of ordinary differential equations (ODEs) defined by Eqs. 8, 11, 12, 13, 14 can be transformed into matrix notation by defining a matrix U where rows correspond to isotopomers and columns correspond to fluxomers. This matrix contains the stoichiometric coefficients of each substrate(s) and product(s) isotopomer(s) for each fluxomer f
p
r
(where p represents a particular (set of) substrate(s) isotopomer(s) consumed through the reaction r). A vector f is then defined that contains all the fluxomers (calculated from the kinetic rate laws of each reaction, the concentrations of reactants and effectors, and the relative abundances of substrate(s) isotopomers), and a vector a that contains the coefficients α
p
r
for each fluxomer f
p
r
. The isotopomer balances can be calculated by:
$$ \frac{dI}{dt}=U.\left(f(t)\circ a\right) $$
(15)
where I is the vector of isotopomer concentrations and ∘ is the element-wise product operator, i.e. (f(t) ∘ a)
x
= f
x
(t). a
x
. The dynamics (and steady-states) of all the isotopomers can be simulated by solving this system of ODEs. Finally, fluxes and metabolite concentrations can be calculated by summing fluxomers and isotopomer concentrations, respectively. The numerical procedure and its implementation are summarized in Fig. 2.
Impact of KIEs on the operation of E. coli central metabolism in 13C-labeling experiments
We applied this framework to evaluate the impact of KIEs on the operation of the glycolytic and pentose phosphate pathways of the model bacterium Escherichia coli grown on 13C-labeled glucose. The kinetic model published by [30] was used as a scaffold to construct the isotopic model, with the carbon transition network taken from [6] (Fig. 3). KIEs measured on five enzymes were included (glucose-6-phosphate dehydrogenase, G6PDH; 6-phosphogluconate dehydrogenase, GND; ribulose-5-phosphate epimerase, RPE; fructose-1,6-bisphosphate aldolase, ALD; and pyruvate dehydrogenase, PDH), with the corresponding α
p
r
coefficients listed in Table 1. Steady-states were simulated for glucose at natural abundance and six different glucose labelings: U-13C-glucose, a mixture containing 80 % of 1-13C- and 20 % of U-13C-glucose, an equimolar mixture of 12C- and U-13C-glucose, 2,3,4,5,6-13C-glucose, 1,2-13C-glucose and 3,4-13C-glucose. These label inputs are representative of those commonly used in 13C-ILEs [6, 32–35].
For the metabolic state represented by the model, the relative changes of fluxes and metabolite concentrations caused by KIEs for each glucose input are shown in Fig. 4a and b, respectively. Impact on these variables is low (<2 %), with maximal changes observed when cells are grown on U-13C-glucose. Since in this situation all the metabolites are fully labeled, impact of KIEs on reaction rates is also maximal. For instance, assuming that KIEs are cumulative at the level of a single reaction [28], the reaction rate of the ribulose-5-phosphate epimerase (RPE) is reduced by 4.14 % when RB5P is fully labeled. However, the predicted change of flux (−1.7 %) is lower than that, which can be explained by the distribution of flux control over several reactions of the network. Moreover, the signs of the changes of fluxes and metabolite concentrations depend on the label input. This indicates that the impact of KIEs on the system cannot be intuitively inferred from their impact at the reaction level, even qualitatively. The impact of KIEs must be investigated at the system-level and not at the level of an isolated reaction or metabolic node.
The predicted changes of metabolite concentrations are significantly lower than the precision commonly reached in metabolomics (circa 10 %) [18, 36], hence they cannot be detected with current methods. The same situation is observed for fluxes, with changes caused by KIEs slightly lower than the flux precision commonly obtained in 13C-metabolic flux analysis (13C-MFA) [5, 6, 12]. However, the predicted flux changes caused by KIEs are one order of magnitude higher than the flux precision estimated using parallel labeling experiments with different label inputs [37]. Therefore, fluxes obtained in high-resolution 13C-MFA studies might be biased, and KIEs should be considered during flux estimation and sensitivity analysis when such high degree of precision can be reached.
Impact of KIEs on isotopic patterns of central metabolites
Isotopic data collected in 13C-ILEs can be exploited quantitatively to infer information on the structure and the operation of metabolic systems. Since KIEs may impact isotopic data and ultimately jeopardize the validity of the biological interpretations, their impact must be assessed rigorously. This was performed for the various isotopic data that can be acquired with modern analytical platforms: isotopomer distributions [11, 38, 39], isotopologue distributions [7, 40] and molecular enrichments [40, 41]. The corresponding (steady-state) isotopic datasets were simulated in the absence or presence of KIEs, for three different label inputs: glucose at natural abundance, a mixture containing 80 % of 1-13C-glucose and 20 % of U-13C-glucose, and an equimolar mixture of 12C- and U-13C-glucose. The distributions of errors caused by KIEs are shown in Fig. 5 for each isotopic dataset and label input.
The sensitivity of isotopic data to KIEs strongly depends on the isotopic information considered and on the label input. As a general trend, the impact of KIEs is maximal for the equimolar mixture of 12C- and U-13C-glucose, and minimal for glucose at natural abundance. The predicted errors caused by KIEs are surprisingly low – below 0.001 – for virtually all metabolites and isotopic data. The only exception is pyruvate, with predicted errors up to 0.0045 for its isotopologues, 0.0047 for its isotopomers, and 0.0045 for its molecular enrichment. These errors originate exclusively from KIEs of pyruvate dehydrogenase (PDH). Here again, the system-level impact of KIEs on the isotopic distributions strongly differs from the impact that can be expected for isolated reactions. Moreover, impacts of individual KIEs do not add up at the system level, contrary to what was previously hypothesized [28].
This example illustrates how the framework developed here can be used to evaluate the sensitivity of each isotopic data to KIEs. The less reliable measurements with respect to KIEs can be identified and discarded from the datasets exploited in quantitative approaches such as 13C-MFA. In the situations investigated here, the biases caused by KIEs on the distribution of isotopes into the glycolytic and pentose phosphate pathway metabolic intermediates are significantly lower than the precision and accuracy obtained by MS (or even MS/MS) [40, 42] and NMR [10], which is around 1 %, and cannot be detected. Many other central metabolic enzymes are likely subject to KIEs, which may significantly increase their predicted impact on isotopic patterns. However, the present results support the assumption that KIEs can be neglected when studies focus on the glycolytic and pentose phosphate pathways and are based on isotopic information obtained from metabolic intermediates.
Isotopic data are robust to KIEs under a broad range of metabolic states
The above results indicate that isotopic data are robust to KIEs locally, i.e. near the given metabolic state investigated. However, metabolic systems are highly non-linear and the impact of KIEs might strongly differ depending on the metabolic state considered [28]. To grasp the global impact of KIEs on isotopic distributions, we first generated 10,000 sets of random enzyme levels for which we calculated the steady-states of the system with 80 % 1-13C- glucose + 20 % U-13C-glucose as label input. Although it is clear that these sets are not representative of those expressed in vivo, they are expected to result in different metabolic states. This was confirmed by analyzing the distribution of some systemic steady-state variables, e.g. the glucose uptake flux varied between 0.02 and 1.6 mmol/gDW/s, and the carbon channeled into glycolysis varied between 0 and 98 % of the glucose uptake flux (Fig. 6). For each metabolic state, errors caused by KIEs on the three isotopic datasets were calculated and are summarized in Fig. 7. Errors remain very low for the vast majority of isotopic data, with 99.7 % of absolute isotopomer errors below 0.002, 99.0 % of absolute isotopologue errors below 0.003, and 97 % of enrichment errors below 0.003. Here again errors strongly depend on the metabolite considered. For all isotopic datasets, the highest errors are related to pyruvate, consistently with the above results. Still, errors for this metabolite errors are minor in most cases (<0.003 for 89 % of its isotopomers and 77 % of its isotopologues). A more detailed analysis revealed that only the proportions of the unlabeled and fully labeled pyruvate isotopomers (and isotopologues) are significantly impacted by KIEs (Fig. 8), errors regarding other isotopomers (and isotopologues) remain very low (<0.002) under all the metabolic states. Errors related to other metabolites are minor (<0.003) for all the datasets, for instance the highest absolute error for glyceraldehyde-3-phosphate (GAP) isotopologues is 0.002. These results highlight a global robustness of the isotopic data measured on glycolytic and pentose phosphate pathway intermediates, and soften the conclusions of [28].
Reversibility is a key determinant of the robustness of isotopic patterns
In a linear pathway with an irreversible step (i.e. not subject to feedback inhibition by its product(s)), reactions downstream the irreversible step do not exert any flux control [43]. Reversibility therefore participates in the robustness of the metabolic operation to KIEs by contributing to the distribution of (flux and concentration) control via feedback regulation. In ILEs, another consequence of reversibility is the bidirectional exchange of isotopes that occurs between substrates and products of reversible reactions. To test if this exchange also plays a role in the robustness of isotope distribution to KIEs, a model of the same network was generated without explicitly considering isotope exchange in the isotopomer balances, i.e. without decoupling forward and backward reaction rates when constructing the equation system. It is important to mention that all the systemic properties that emerge from metabolic regulation (in particular the distributed control) are identical in both models since the rate laws remain unchanged. With the two models, the impact of KIEs on isotopologue abundances was simulated with and without isotope exchange (Fig. 9). Errors caused by KIEs were significantly higher for all metabolites when exchange was not explicitly considered, with a 10-fold increase of the mean of absolute errors (0.0014 vs 0.00015). Moreover, errors increased for all metabolites. For instance, while errors on fructose-6-phosphate (F6P) and dihydroxyacetone phosphate (DHAP) isotopologues are negligible with exchange (<0.0002 for both metabolites), they increase up to 0.006 for F6P and 0.008 for DHAP without exchange. Hence, reversibility appears as a key determinant of the robustness of isotopic data to KIEs, both by contributing to the distribution of control across the network and by enabling bidirectional exchange of isotopes between metabolites.
This is consistent with the 13C-depletion observed in amino acids [44] since they are produced by less reversible biosynthetic pathways and are therefore expected to be more sensitive to KIEs. Conversely, no major isotopic fractionation was detected on central metabolites in a recent study aiming to produce standard samples that contain metabolites with controlled and predictable isotopic contents [40]. Our simulations support these conclusions and strengthen the reliability of the proposed strategy. These results stress the necessity to take into account reversibility, or at least to evaluate the potential consequences of neglecting this important property of metabolic systems, when integrating isotopic data into quantitative models.
