Which sets of elementary flux modes form thermodynamically feasible flux distributions?

Elementary flux modes (EFMs) are non‐decomposable steady‐state fluxes through metabolic networks. Every possible flux through a network can be described as a superposition of EFMs. The definition of EFMs is based on the stoichiometry of the network, and it has been shown previously that not all EFMs are thermodynamically feasible. These infeasible EFMs cannot contribute to a biologically meaningful flux distribution. In this work, we show that a set of thermodynamically feasible EFMs need not be thermodynamically consistent. We use first principles of thermodynamics to define the feasibility of a flux distribution and present a method to compute the largest thermodynamically consistent sets (LTCSs) of EFMs. An LTCS contains the maximum number of EFMs that can be combined to form a thermodynamically feasible flux distribution. As a case study we analyze all LTCSs found in Escherichia coli when grown on glucose and show that only one LTCS shows the required phenotypical properties. Using our method, we find that in our E. coli model < 10% of all EFMs are thermodynamically relevant.


Introduction
Elementary flux mode analysis (EFMA) is a key concept in constraint-based modeling, which allows a metabolic network to be decomposed into irreducible functional building blocks, called elementary flux modes (EFMs). An EFM corresponds to a minimal set of reactions that can operate at steady state, thereby using all irreversible reactions in the appropriate direction [1]. Here minimal means that no reaction can be removed from the set without losing the ability to form a non-zero steadystate flux. EFMs represent functional units in a metabolic network. In fact, every steady-state flux can be represented as a nonnegative linear combination of EFMs without cancelations [S. M€ uller and G. Regensburger, unpublished results;2,3]; it is always possible to find such a 'conformal sum' [4,5]. The no-cancelation rule guarantees that every reaction proceeds in the same direction in all contributing EFMs and accounts for the fact that, for given metabolite concentrations, only one direction of a reversible reaction is thermodynamically feasible (TF).
From a systems-biology point of view, the advantage of an unbiased decomposition of a metabolic network into EFMs lies in the ability to fully characterize the metabolic capabilities of an organism. This can be used, for instance, to analyze cellular robustness [6] or in metabolic engineering to turn wild-type organisms into so-called 'networks of minimal functionality' [7]. These (mutant) networks are typically made up of very few, desired EFMs while all the unwanted (wild-type) functionality is eliminated by appropriately selected gene knockouts [8][9][10][11][12]. Therefore, EFMA is an ideal tool for metabolic engineering and synthetic biology to rationally design optimal cell factories [13].
However, a complete EFMA is currently limited to medium-scale metabolic models since the number of EFMs explodes with the size of the metabolic network [14]. Using a massively parallelized approach, the largest complete EFMA reported to date found almost two billion EFMs in a metabolic reconstruction of Phaeodactylum tricornutum with 318 reactions [15]. Alternative methods do not aim to find all EFMs, but limit their scope to find subsets. Subsets are selected randomly [16,17] or based on support information [18] or subject to additional constraints [19][20][21]. The question remains, however, whether or not these subsets are biologically relevant. To find relevant EFMs, Rezola et al. [22] used gene expression data and identified small subsets of EFMs that successfully characterized and described key metabolic features in different tissues. Similarly, Jol et al. [23] and our group [24,25] used experimentally determined metabolomes to identify all EFMs that are TF. However, even if two EFMs are TF, it does not necessarily mean that their convex superposition is TF as well [23]. This raises the question which TF EFMs can be combined in a convex superposition such that the resulting flux distribution will again be TF.
Here, we expand on our earlier work on the thermodynamics of EFMs [24,25] and present a mixed integer linear program (MILP) that identifies the largest thermodynamically consistent sets (LTCSs) of EFMs. These LTCSs are characterized by the fact that every nonnegative linear combination of its EFMs results in a TF flux distribution. Moreover, we show that physico-chemical constraints alone already severely limit the metabolic capabilities of an organism since only a small fraction of all EFMs are required to represent TF flux distributions. This confirms the hypothesis that, under given conditions, only a few EFMs are actually biologically relevant and accessible to an organism [26].

Notation and model assumptions
We consider a metabolic core model of Escherichia coli [25], referred to as M-glc, to study growth on minimal medium (containing ammonia, oxygen, phosphate, protons and water) with glucose as the sole carbon source. The model is characterized by its stoichiometric matrix S 2 R m Â r with m = 76 internal metabolites and r = 101 reactions of which 48 are reversible. The network contains n tot = 169 916 EFMs of which n = 32 374 EFMs are TF. Reactions and intracellular metabolites are thermodynamically characterized by their Gibbs free energy of reaction D r G and their standard transformed Gibbs energy of formation D f G 00 , respectively. The latter is estimated using a pH of 7 and an ionic strength I = 0.15 M at a temperature T = 310.15 K (37°C), according to Alberty [27].

Thermodynamic EFMA
Thermodynamic constraints are often utilized to augment classical constraint-based approaches [28]. For instance, Hoppe et al. [29] developed a metabolomicsintegrated flux-balance analysis. Similarly, we developed thermodynamic EFMA (tEFMA) [25], a computational tool that calculates all TF EFMs in a metabolic reconstruction. tEFMA exploits the fact that, according to the second law of thermodynamics, an EFM e i is TF if and only if all reactions j which support e i proceed in the direction of negative Gibbs free energy [30], that is, D r G j < 0 for all reactions j with e i j [ 0. Based on this fundamental property, tEFMA avoids the calculation of thermodynamically infeasible EFMs, which drastically reduces the computational burden of an EFMA and makes the analysis of large scale metabolic networks feasible without losing any biologically relevant information [24]. Thus, given an experimentally measured cellular metabolome, the standard Gibbs free energy of formation for as many cellular metabolites as possible, and a metabolic reconstruction, tEFMA returns the complete set of TF EFMs consistent with the measurements.
Largest thermodynamically consistent sets tEFMA computes all n TF EFMs of a network. However, not every set of TF EFMs is necessarily thermodynamically consistent. For example, two EFMs that utilize the same reversible reaction in different directions cannot be active simultaneously. This is illustrated in the simple example network in Fig. 1. Suppose that all four EFMs are TF. Yet, EFM2 and EFM3 cannot be active at the same time since D r G cannot be smaller than zero for both directions of the reversible reaction R3. In other words, thermodynamics implies the no-cancelation rule mentioned in the introduction. In this example, the sets {EFM1, EFM2, EFM4} and {EFM1, EFM3, EFM4} are the LTCSs, and its elements can contribute to a TF steady-state flux.
For instance, consider the TF steady-state flux v T = (v 1 ,. . .,v 5 ) = (1, 2, À1, 2, 1). Obviously, the flux can be decomposed as v = EFM2 + 2 9 EFM3 (see   2 for an illustration). Although EFM2 and EFM3 are TF individually, they cannot be active simultaneously (no-cancelation rule). Hence, this decomposition is not thermodynamically consistent. In contrast, the representation v = EFM1 + EFM3 + EFM4 is thermodynamically consistent. Indeed, these three EFMs form one of the LTCSs above, and every TF steady-state flux can be represented by elements of one LTCS. The example raises the question of how all LTCSs can be computed systematically, given a set of TF EFMs.
Definition: LTCS. A set of TF EFMs is called thermodynamically consistent if every nonnegative linear combination of its elements is TF. Moreover, a set of TF EFMs is called an LTCS.
(i) if the set is thermodynamically consistent, and (ii) if no other TF EFM can be added to the set without losing thermodynamic consistency.
We first determine an LTCS L 1 of maximal cardinality. Alternative LTCSs of maximal cardinality or LTCSs of lower cardinality, L l with l > 1, can be found by successively excluding already existing LTCSs; see constraint (2) below An LTCS L 1 of maximal cardinality is an optimal solution to the MILP max where lnðc min k =c 0 Þ lnðc k =c 0 Þ lnðc max k =c 0 Þ: ð1gÞ We use the superscript in k 1 to denote its association with the LTCS L 1 . Briefly, k 1 indicates the presence (k 1 i ¼ 1) or absence (k 1 i ¼ 0) of EFM e i in the LTCS L 1 , and we maximize P n i ¼ 1 k 1 i , that is, the cardinality of L 1 , by varying the contributing EFMs and  the (logarithms of the) metabolite concentrations c k , as stated in Eqn (1a). Most importantly, d ij indicates if EFM e i is supported by reaction j, and q j counts the number of EFMs supported by reaction j, as stated in Eqns (1d) and (1c). If at least one EFM is supported by reaction j, that is, q j ≥ 1, then this reaction must be TF, according to the main constraint (Eqn 1b). Equivalently, if reaction j is infeasible, then q j = 0 and hence k 1 i is forced to 0 for all EFMs i supported by reaction j.
Finally, Eqns (1e) and (1f) determine the Gibbs free energy of reaction j, given the (logarithms of the) metabolite concentrations. Thereby, S kj denotes the elements of the stoichiometry matrix. The inequalities (Eqn 1g) constrain the metabolite concentrations.
Alternative optima and suboptimal solutions k l with l > 1 (and the corresponding LTCS L l ) can be found by successively excluding already existing solutions k j with j 2 {1,. . .,lÀ1} from the MILP [9]. This is achieved by successively adding the constraint The process terminates when the MILP becomes infeasible, and no further solutions are found.
In the following, we computed all LTCSs for E. coli grown on minimal medium with glucose as the sole carbon source. Subsequently, we narrowed down the number of LTCSs to one by successively considering additional yield, expression and flux data.
LTCSs are much smaller than the set of TF EFMs As a matter of fact, there are 40 LTCSs for E. coli grown on minimal medium with glucose as the sole carbon source. The largest LTCS contains 15 560 EFMs. This corresponds to only 9% of all EFMs or 47% of the TF EFMs (see Fig. 3).
Moreover, we found that in general the relative frequency of biomass producing EFMs is larger in LTCSs of smaller cardinality. Still, the average number Biomass producing EFMs Not biomass producing EFMs of biomass producing EFMs per LTCS (1984 AE 69%) appeared to be more stable than the average number of EFMs per LTCS (4316 AE 86%).

Yield data identifies biologically relevant LTCSs
An LTCS represents the metabolic capabilities of E. coli, under the conditions specified in the model. To characterize these capabilities, we used maximum yield parameters Y xi=glc . Y xi=glc was defined as the maximum of the yields of all EFMs in an LTCS for a specific product x i . To identify the biologically relevant LTCSs, we used typical growth parameters obtained by Andersen and Meyenburg [31]. Figure 4 shows different maximal yields for each LTCS in comparison with measured data. Note that those maximal yields were only achieved by a few EFMs within an LTCS. Most of the EFMs in any given LTCS had a smaller or even zero yield. Thus, every yield between the maximum and zero can be achieved by a suitable combination of a maximum yield EFM and a zero yield EFM. In particular, if the measured yield is below the maximum, then it can be achieved by an appropriate combination of EFMs. Conversely, if the achievable maximum yield of an LTCS is below the measured value, then no combination of EFMs can result in the observed yield, and those LTCSs can be excluded from further analysis. We found that only 12 out of 40 LTCSs were consistent with the measured yields (see Fig. 5). These 12 sets can be calculated directly if the measured yields are used as additional constraints in the MILP (Eqn 1). Then the modified MILP reads where and Here r u counts the number of EFMs that have a certain minimal yield of metabolite u, cf. Eqn (3f), and ɛ iu indicates if EFM i has the required yield for metabolite u, cf. Eqn (3g). The constraint (Eqn 3c) ensures that at least one EFM has the required yield. All other subequations are also found in the original MILP (Eqn 1).

Expression data further reduces the number of relevant LTCSs
We further analyzed the remaining 12 LTCSs, using expression data. Six LTCSs, L 8 ; L 12 ; L 15 ; L 16 ; L 21 and L 38 , had an active fumarate reductase (FrdABCD), but had an inactive succinate dehydrogenase (SdhCDAB) in all their EFMs. This is in contrast to experimental findings since, under aerobic   Each LTCS is denoted by its set index and printed in the same color as in Fig. 7. The letters A to L denote different segments in the diagram, along with the number of TF EFMs in these segments. Only L 1 (full line) was found to be consistent with yield, expression, and flux data. All other LTCSs (dashed lines) were eliminated, which turned out to be consistent with independent 13 C flux data.

1787
The conditions, sdhCDAB is optimally expressed [32] while the frdABCD operon is repressed [33]. Thus only six LTCSs, L 1 ; L 3 ; L 5 ; L 6 ; L 7 and L 18 , were found to be consistent with the data. The Venn diagram in Fig. 6 singled out L 3 and L 18 due to their lack of overlap with the other LTCSs (see the next section for a mechanistic characterization of these sets). In contrast to L 3 and L 18 , the four remaining LTCSs share some core functionality, represented by a large fraction of common EFMs (segment A). We investigated if these EFMs were characterized by their supports and observed that on average shared EFMs (contained in several LTCSs) were shorter than EFMs unique to an LTCS (see Fig. 7). Functionally, EFMs in segment A do not invoke the pentose phosphate pathway and do not produce biomass. However, some produce maintenance energy (see Table 1). In all other segments, biomass production is feasible. For a complete listing see Table 2.

Flux data pinpoint a single relevant LTCS
The metabolic capabilities of E. coli when grown aerobically on minimal medium under glucose limited conditions are fully described by the six LTCSs in Fig. 6. To further narrow down the number of LTCSs, we analyzed the different segments in Fig. 6, using flux data.
All EFMs in the segments D, E and J (see L 3 and L 18 in Fig. 6) were characterized by a reverse flux across glucose-6-phosphate isomerase (Pgi), directed towards glucose 6-phosphate (see Table 1). Under the standard growth conditions investigated here, Pgi is forward active [34]. Thus we were able to eliminate L 3 and L 18 from the set of relevant LTCSs. (Note that segment G, which is the largest subset of L 3 , was not removed since EFMs in G have zero flux across Pgi.).
In L 1 (as in every other LTCS) all reversible reactions have a fixed direction (due to the no-cancelation rule). The predicted directions were fully consistent with independent 13 C flux data [35].
To summarize, we found that from 40 LTCSs only L 1 was consistent with all data. L 1 contains 15 559 TF EFMs, which represent only 9% of all EFMs. More specifically, 4486 EFMs produce biomass, 2024 produce maintenance energy, and 54 produce both. In fact, L 1 is composed of several segments: all EFMs in segment A do not invoke the pentose phosphate pathway and do not produce biomass, but are responsible for the production of maintenance energy; all EFMs in segment B produce biomass without invoking TktA or TktB; all EFMs in segment G do not carry flux across Pgi, but are able to produce biomass and/or maintenance energy. Finally, EFMs in segment H are not characterized by a single common property.

General remarks on LTCSs
Using E. coli as an example, we outlined a procedure that narrowed down the feasible solution space and eventually identified a single LTCS. The success of such an analysis is dependent on the quality of the measured metabolome. However, the general concept of LTCSs is not affected by the metabolome's quality. In fact, we can find LTCSs even if the metabolome is unknown. As soon as a network contains at least two EFMs that are supported by a reversible reaction carrying fluxes in opposite directions, different LTCSs exist. Moreover, the cardinalities of these LTCSs are always smaller than the total number of EFMs in a system (see for instance Fig. 1, where 4 ¼ n tot [ jL 1 j ¼ jL 2 j ¼ 3). Thus even in the absence of a measured metabolome it is useful to look at LTCSs, as only then is a thermodynamically consistent understanding guaranteed. Although an LTCS is less complex than the complete set of EFMs, one now has to analyze multiple LTCSs. In general the number of LTCSs scales combinatorially with the number of reversible reactions in a network. Practically, that is why an accurately measured metabolome is essential.

Discussion
Every intercellular flux distribution is TF and can be decomposed into TF EFMs. However, the reverse is not necessarily true. That is, the conformal superposition of any two TF EFMs is itself not necessarily TF.
Here we developed a method that identifies the largest sets of TF EFMs that are thermodynamically consistent (LTCSs). Within an LTCS every nonnegative linear combination of its elements results in a TF flux distribution. A necessary condition for an LTCS is that all reactions supporting the EFMs of an LTCS operate in the same direction. This is known as the no-cancelation rule [2,3].
Geometrically, an LTCS spans a TF subcone of the flux cone (see Fig. 8). In fact, thermodynamic constraints segment the flux cone into LTCSs.
Although subcones corresponding to LTCSs may overlap (see Fig. 6), each LTCS has unique metabolic capabilities.
We found that shorter TF EFMs are more likely to be elements in multiple LTCSs than longer TF EFMs. More specifically, we found that segment A in Fig. 6 contains only EFMs producing maintenance energy, but no EFM producing biomass. This can be understood considering the number of reactions involved. Whereas biomass requires the production of many precursors which involves long pathways, maintenance energy requires only the production of adenosine triphosphate (ATP), which can be achieved by short routes. Since every reaction has to comply with the second law of thermodynamics, the likelihood of thermodynamical feasibility decreases with increasing number of contributing reactions. Every metabolic phenotype can be described by an LTCS. Conversely, if an EFM is not part of an LTCS, it is biologically irrelevant, since it does not contribute to a thermodynamically consistent decomposition of a TF flux distribution. In general, decompositions into EFMs are not unique, and several methods have been proposed [36][37][38][39][40]. However, none of these methods takes thermodynamics into account, which may lead to inconsistent decompositions. Therefore, it is even more important to identify those LTCSs that consistently describe a phenotype.
We outlined a systematic procedure to identify biologically relevant LTCSs. Based on the integration of additional (omics) data, we successively narrowed down the number of LTCSs. In fact, most LTCSs were found to be inconsistent with commonly available growth parameters. Further consideration of expression and flux data eventually identified a single LTCS that characterizes the phenotype. The additional information could have been used to adapt the network first. Doing so would have reduced the number of LTCSs from 40 to four (L 1 ; L 2 ; L 31 ; and L 36 ), and a comparison with growth parameters would have identified L 1 , the same LTCS as before. However, for less studied organisms detailed data may not be available. In this case, an analysis of phenotypical properties like in Tables 1 and 2 will identify the most valuable piece of information to narrow down the number of LTCSs.
Our method is able to compute LTCSs if all TF EFMs are known, and we showed that the set of TF EFMs characterizing a phenotype is smaller than the set of all TF EFMs. However, currently our method does not allow computation of LTCSs directly. It would be desirable to enumerate only the biologically relevant EFMs, which would facilitate an unbiased analysis of metabolic systems even on a genome-scale level. Recent progress enabled the selective calculation of subsets of EFMs [21] and the identification of relevant regulated EFMs [41,42]. Combining these ideas with our current approach may lead to promising lines of future research.

Methods
We used the software package tEFMA [24] together with published metabolite concentration data [43] to calculate all TF EFMs in a core metabolic model of E. coli [25] growing on a glucose limited minimal medium. For all unmeasured metabolites in this model we used conservative default concentration ranges between c min k ¼ 10 À7 m and c max k ¼ 1 m. D f G 0 values were obtained from the online version of EQUILIBRATOR [44].
For two metabolites (ubiquinol-8 and biomass) no D f G values were available. Reactions to which those metabolites contributed were not checked for thermodynamic feasibility to avoid false conclusions [25].
The set of Eqns (1) and (2) were solved with the IBM ILOG CPLEX Optimization Studio, version 12.5. A Perl-script that sets up the systems equations and invokes the CPLEX LP solver, and the metabolic model and all data used in this study are available at https://github.com/mpgerstl/ltcsCalculator. Note that CPLEX is a commercial software product although free academic licenses are available.