Email updates

Keep up to date with the latest news and content from BMC Bioinformatics and BioMed Central.

Open Access Highly Accessed Methodology article

F2C2: a fast tool for the computation of flux coupling in genome-scale metabolic networks

Abdelhalim Larhlimi12*, Laszlo David345*, Joachim Selbig12 and Alexander Bockmayr34

Author Affiliations

1 Department of Bioinformatics, Institute for Biochemistry and Biology, University of Potsdam, Karl-Liebknecht-Str. 24-25, D-14476 Potsdam, Germany

2 Max-Planck Institute for Molecular Plant Physiology, Am Mühlenberg 1, D-14476 Potsdam, Germany

3 FB Mathematik und Informatik, Freie Universität Berlin, Arnimallee 6, D-14195 Berlin, Germany

4 DFG-Research Center Matheon, Arnimallee 6, D-14195 Berlin, Germany

5 Berlin Mathematical School (BMS), Arnimallee 6, D-14195 Berlin, Germany

For all author emails, please log on.

BMC Bioinformatics 2012, 13:57  doi:10.1186/1471-2105-13-57

The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2105/13/57


Received:16 December 2011
Accepted:23 February 2012
Published:23 April 2012

© 2012 Larhlimi et al.; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

Flux coupling analysis (FCA) has become a useful tool in the constraint-based analysis of genome-scale metabolic networks. FCA allows detecting dependencies between reaction fluxes of metabolic networks at steady-state. On the one hand, this can help in the curation of reconstructed metabolic networks by verifying whether the coupling between reactions is in agreement with the experimental findings. On the other hand, FCA can aid in defining intervention strategies to knock out target reactions.

Results

We present a new method F2C2 for FCA, which is orders of magnitude faster than previous approaches. As a consequence, FCA of genome-scale metabolic networks can now be performed in a routine manner.

Conclusions

We propose F2C2 as a fast tool for the computation of flux coupling in genome-scale metabolic networks. F2C2 is freely available for non-commercial use at https://sourceforge.net/projects/f2c2/files/ webcite.

Background

The huge amount of genomic, transcriptomic and related data has allowed for a fast reconstruction of an increasing number of genome-scale metabolic networks, e.g. [1-7]. In the absence of detailed kinetic information, constraint-based modeling and analysis has recently attracted ample interest due to its ability to analyze genome-scale metabolic networks using very few information [8-10]. Constraint-based analysis is based on the application of a series of constraints that govern the operation of a metabolic network at steady state. This includes the stoichiometric and thermodynamic constraints, which limit the range of possible behaviors of the metabolic network, corresponding to different metabolic phenotypes. Applying these constraints leads to the definition of the solution space, called the steady-state flux cone[11]:

C = { v R n Sv = 0 , v i 0 , for all i Irr } , (1)

where S is the m×n stoichiometric matrix of the network, with m internal metabolites (rows) and n reactions (columns), and the vector v R n gives a flux distribution. Furthermore, Irr⊆{1,…,n} denotes the set of irreversible reactions in the network, and Rev={1,…,n}∖Irr denotes the set of reversible reactions.

The flux cone contains the full range of achievable behaviors of the metabolic network at steady state. Various approaches have been proposed either to search for single optimal behaviors using optimization-based methods [12-16] or to assess the whole capabilities of a metabolic network by means of network-based pathway analysis [11,17-20].

Flux coupling analysis (FCA) is concerned with describing dependencies between reactions [21]. The stoichiometric and thermodynamic constraints not only determine all possible steady-state flux distributions over a network, they also induce coupling relations between the reactions. For instance, some reactions may be unable to carry flux under steady-state conditions. If a non-zero flux through a reaction in steady-state implies a non-zero flux through another reaction, then the two reactions are said to be coupled (see Def. 2 for a formal definition). FCA has been used for exploring various biological questions such as network evolution [22-24], gene essentiality [22], gene regulation [25-27], analysis of experimentally measured fluxes [28,29], or implications of the structure of the human metabolic network for disease co-occurrences [30]. Having a time efficient implementation of FCA is important in such studies.

After introducing the main existing algorithms for flux coupling analysis, we propose in this paper a new algorithm which significantly speeds up the calculation of flux coupling. Our algorithm is based on two main principles. First, we reduce the stoichiometric model as much as possible when parsing the stoichiometric matrix. Second, we use inference rules to minimize the number of linear programming problems that have to be solved. We prove the efficiency of our algorithm by successfully competing with the most recent approach [31]. We show that FCA can now be quickly performed even for very large genome-scale metabolic networks.

Approaches for flux coupling analysis

Several algorithms were developed to calculate flux coupling between reactions. For a comparison among the existing approaches, the reader may refer to [31,32]. In the following, we focus on flux coupling methods based on solving a sequence of linear programming (LP) problems. These methods have proved to be significantly faster than other algorithms.

Definitions

We give a short overview of the important concepts we will use throughout this paper. First, we formally define blocked reactions in a metabolic network.

Definition 1 (Blocked reaction)

Given the steady-state flux cone C, let i∈{1,…,n} be a reaction. If vi=0, for all vC, reaction i is called blocked, otherwise i is unblocked.

In the following, we assume that the flux cone is not trivial, i.e., not all reactions are blocked. Next, we define the (un)coupling relationships between reactions.

Definition 2 (Coupling relations)

Let i,j be two unblocked reactions. The (un)coupling relationships = 0 , = 0 , and ↦ are defined in the following way:

i = 0 j if for all vC, vi=0 implies vj=0.

i = 0 j if for all vC, vi=0 is equivalent to vj=0.

ij if there exists λ≠0 such that for all vC,vj=λvi.

i j if there exists vC such that vi=0 and vj≠0.

Reactions i and j are fully (resp. partially, directionally) coupled if the relation ij (resp. i = 0 j , i = 0 j ) holds. Otherwise, i and j are uncoupled.

Note that ij (resp. i = 0 j ) is equivalent to ji (resp. j = 0 i ). In addition, ij implies i = 0 j , which in turn is equivalent to ( i = 0 j and j = 0 i ) .

As shown in [33] the reversibility type of reactions is a key concept in flux coupling analysis.

Definition 3 (Reversibility types)

A reversible reaction iRev is called fully reversible if there exists a flux vector vC such that vi≠0 and vj=0 for all jIrr. Otherwise, reaction i is called pseudo-irreversible.

Using the reversibility type of reactions, we define the following reaction sets which will be used to determine the cases where flux coupling can occur between reactions.

Frev={ii is fully reversible},

Prev={ii is pseudo-irreversible and there exist v + ,vC such that v i + > 0 , v i < 0 } ,

Irev={iiFrevPrev and vi≠0 for some vC},

Blk={ii is blocked }.

Note that the above reaction sets are disjoint and their union is equal to the set of all reactions, i.e., BlkIrevPrevFrev={1,…,n}.

As the classification of reactions according to their reversibility types has been treated elsewhere [31,33,34], we focus on calculating the flux coupling between reactions given the reaction sets Blk, Irev, Prev and Frev. First, we briefly review the main existing flux coupling methods based on linear programming. Afterwards, we present our new method to speed up the flux coupling analysis.

Flux coupling finder

The Flux Coupling Finder (FCF) algorithm [21] determines blocked reactions as well as coupled reactions by solving a sequence of linear programming (LP) problems. The FCF algorithm requires that each reversible reaction is split into two irreversible reactions, one forward and one backward. This may hamper the application of FCF to determine flux coupling in large genome-scale metabolic networks. Indeed, splitting reversible reactions results in an increase in the number of variables (resp. constraints) by |Rev| (resp. 2|Rev|). Since the FCF algorithm solves four LP problems to identify the maximum and minimum flux ratios for every pair of reactions, the total number of LP problems that have to be solved is very large. Furthermore, the FCF algorithm does not compute directly coupling relationships between reactions. A post-processing step is needed to deduce couplings between reactions (in the original network) from those between reaction directions (in the reconfigured network). More importantly, the FCF algorithm explores exhaustively all possible reaction pairs. This leads to a very big number of LP problems that have to be solved. This strategy may not scale well for genome-scale models of complex microorganisms which involve a large number of reactions.

Reversibility-based flux coupling analysis

To cope with the drawbacks of the FCF algorithm, we used the reversibility type of reactions to develop an improved version (WRP-FCF) of the original FCF method [31,34]. When looking for coupled reactions, WRP-FCF applies linear programming only in those cases where coupling relationships can occur [33]. Namely, given two unblocked reactions i and j, we have exactly three cases where a flux coupling is possible between i and j:

Flux coupling between reversible reactions

If i,jPrev or i,jFrev, i = 0 j , i = 0 j and ij are equivalent. More importantly, only the stoichiometric constraints determine whether i = 0 j holds, independently of the thermodynamic constraints.

Assume that blocked reactions have been identified beforehand and their respective columns in the stoichiometric matrix have been removed. Consider the following LP problem

max { v j : Sv = 0 , v i = 0 , 0 v j 1 } , (2)

and let v be an optimal solution. According to Proposition 6.13 in [34], i = 0 j if and only if v j = 0 .

Flux coupling between (pseudo-)irreversible reactions

If iIrev and jPrev, we only have to check whether i = 0 j . The other coupling relationships cannot occur. Let v1 resp. v2 be optimal solutions of the two LP problems

max { v j : Sv = 0 , v k 0 , k Irr , v i = 0 , v j 1 } , min { v j : Sv = 0 , v k 0 , k Irr , v i = 0 , v j 1 } . (3)

Then i = 0 j if and only if v j 1 = v j 2 = 0 .

Flux coupling between irreversible reactions

If i,jIrev, in analogy with the FCF algorithm, we determine upper and lower bounds Uij and Lij such that 0≤LijvjviUijvj for all vC by solving the LP problems

L ij = min { v i : Sv = 0 , v j = 1 , v k 0 , k Irr } , U ij = max { v i : Sv = 0 , v j = 1 , v k 0 , k Irr } . (4)

Comparison of Lij and Uij allows us to determine whether reactions i and j are coupled using the following rules:

i = 0 j (resp. j = 0 i ) if and only if Lij≠0 (resp. Uij≠ + ),

ji if and only if Lij≠0, Uij≠0 and Lij=Uij.

The improved version WRP-FCF does not make an exhaustive computation for all pairs of reactions and does not require a reconfiguration of the metabolic network. Accordingly, not only is the number of LP problems used by WRP-FCF smaller than the number solved by FCF, but also the LP problems used by WRP-FCF are simpler than those employed by FCF. For mathematical proofs, the reader may refer to [31,34].

Feasibility-based flux coupling analysis

Linear programming can be solved by enumerating a set of adjacent feasible solutions such that at each solution the objective function improves or is unchanged. Accordingly, searching for a feasible solution can be faster than computing an optimal solution [35]. Based on this observation, [31] proposed the FFCA approach for flux coupling analysis transforming the optimization-based LP problems used in the WRP-FCF method into feasibility-based LP problems. The authors compared FFCA with other available flux coupling algorithms, and showed that FFCA is slightly faster than WRP-FCF.

Methods

Before using linear programming to calculate flux coupling between reactions, we preprocess the metabolic network in order to i) reduce the number of variables and constraints of the subsequent LP problems and ii) classify reactions according to their reversibility type. The network reduction is mainly based on the removal of trivially blocked reactions and the merging of the stoichiometric columns corresponding to trivially coupled reactions [36-39]. For this, one can use the kernel of the stoichiometric matrix. Alternatively, we can apply the following reduction rules which require a simple parsing of the stoichiometric matrix and are not time consuming. This strategy allows avoiding potential numerical instabilities related to the computation of a basis of the kernel.

Preprocessing

Certain metabolites, called dead-end metabolites[37], are produced (resp. consumed) by some reactions without being consumed (resp. produced) by other reactions. This concept is illustrated in Figure 1(a) where the dead-end metabolite H is produced by reaction 13 without being consumed by any of the remaining reactions.

thumbnailFigure 1. Exemplary metabolic network MetNet before and after preprocessing. (a) MetNet consists of eight metabolites (A,…,H), and thirteen reactions (1,…,13), whereof six reactions are irreversible. Metabolites are depicted as nodes while reactions are depicted as arrows. Reversible reactions are indicated by double arrowheads. (b) MetNet after preprocessing where dead-end metabolites and blocked reactions were removed and fully coupled reactions were merged iteratively. This resulted in the removal of the blocked reaction 13 and the merging of the pairs of equivalent reactions (1,2), (8,9) and (11,12)).

As stated below, reactions which are consuming or producing dead-end metabolites are blocked.

Observation 1 (Dead-end Metabolite)

Let k∈{1,…,m} be a metabolite. Then, the following hold:

If there exists a reaction i such that Ski≠0 and Skj=0 for all ji, then reaction i is blocked.

If there exists a set of reactions IIrr such that Ski>0 (resp. Ski<0) for all iI and Skj=0 for all jI, then all reactions iI are blocked.

In each of these cases, k is called a dead-end metabolite.

Certain metabolites are involved in exactly two reactions. For instance, in the network MetNet depicted in Figure 1(a), metabolite E is produced/consumed only by reactions 8 and 9. The following observation states that the fluxes through reactions involving such metabolites are proportional to each other.

Observation 2 (Trivial Full Coupling (TFC))

Let i and j be two reactions such that, for some metabolite k∈{1,…,m}, Ski≠0, Skj≠0 and Skl=0 for all li,j. Then, reactions i and j are either blocked or fully coupled.

The identification of dead-end metabolites and their corresponding blocked reactions allows us to reduce the number of metabolites and reactions that matter for identifying coupled reactions. As stated in the following observation, the removal of the rows (resp. columns) in the stoichiometric matrix corresponding to dead-end metabolites (resp. blocked reactions) does not influence the flux coupling between reactions.

Observation 3 (Reduction Rule 1)

Let D be a set of dead-end metabolites and let B be a set of blocked reactions. For convenience, suppose B={s + 1,…,n}. Let S be the submatrix of S formed by the rows Sk with kD and the columns Sl with lB. Let Irr=IrrB. Then, for all pairs of reactions iB and jB,

i = 0 j if and only if v i = 0 implies v j = 0 , for all v R s such that Sv=0 and v p 0 for all p Irr .

ij if and only if there exists λ≠0 such that v j = λ v i , for all v R s with Sv=0 and v p 0 for all p Irr .

The next observation shows that two fully coupled reactions can be represented by only one column in the stoichiometric matrix, without altering the flux coupling between reactions.

Observation 4 (Reduction Rule 2)

Let k,l be two reactions such that for all vC,vl=λvk for some λ≠0. For convenience, suppose l=n and λ>0. Let S be the m×(n−1) matrix defined by S p = S p for all pk,l and S k = S k + λ S l . Let Irr=(Irr∪{k})∖{l} if lIrr, and Irr=Irr otherwise. Then, for all pairs of reactions il and jl,

i = 0 j if and only if v i = 0 implies v j = 0 , for all v R n 1 such that Sv=0 and v p 0 for all p Irr .

ij if and only if there exists λ≠0 such that v j = λ v i , for all v R n 1 with Sv=0 and v p 0 for all p Irr .

Note that when applying the reduction rules of Observations 3 and 4, further metabolites and reactions may fulfill the conditions of Observations 1 and 2. Accordingly, we apply these reduction rules in an iterative way. As an illustration, the reduction of the network MetNet depicted in Figure 1(a) involves two iterations. In the first one, metabolite H and reaction 13 are removed, the pairs of reactions (1,2) and (8,9) are merged and metabolites A and E are removed. In the second iteration, the equivalent reactions (11,12) are merged and metabolite G is removed. The preprocessed network depicted in Figure 1(b) contains only four metabolites and nine reactions.

Certain fully coupled reactions could not be identified using Observation 2. The following lemma proves that all fully coupled reaction pairs can be deduced from the kernel kern ( S ) = { v R n Sv = 0 } of the stoichiometric matrix after the removal of all blocked reactions.

Lemma 1

Let (S,Irr) be a metabolic network with n unblocked reactions. For a pair of reactions (i,j) the following are equivalent:

ij.

there exists λ R { 0 } such that vi=λvj, for all vkern(S).

Proof

⇐" Immediate.

"⇒" Since ij, there exists λ≠0 such that vi=λvj for all vC. Assume by contradiction that there is vkern(S) such that viλvj and let L={lIrrvl<0}. Since every reaction is unblocked, for every lL there exists g(l)C with g l ( l ) = 1 . Let w = v l L v l g ( l ) . Clearly, wkern(S) and wl>0 for all lIrr, thus wC. However, wiλwj, contradicting ij. The required statement follows. □

In analogy with the WRP-FCF and FFCA approaches, we identify the reversibility type of reactions in order to apply linear programming only in those cases where coupling relationships can occur [33]. Here, we use the procedure for reaction classification described in [31,34]. Note that applying the above reduction rules beforehand reduces the number of variables and constraints in the LP problems used for the classification of reactions.

Based on the results above, we propose to apply the preprocessing procedure given in Table 1 before identifying coupled reactions using linear programming. We show later that the preprocessing step turns out to be crucial for obtaining an efficient flux coupling algorithm.

Table 1. Main steps of the preprocessing procedure

Algorithmic improvements

In certain metabolic networks, the conversion of a set of substrates into a set of products can be made by different reactions having the same stoichiometry. A simple example of such reactions are isoenzymes which make the same conversion of substrates into products. This concept is illustrated in Figure 1(a) where both reactions 4 and 5 convert C into D in the same way. This holds also for reaction 7 and the merged equivalent reactions (8,9) in Figure 1(b) showing that the network preprocessing may simplify the identification of such alternative conversions. The flux coupling of such reactions is trivial using the following lemma.

Lemma 2 (Trivial Uncoupling (TUC))

Let i,jIrev and k,lPrevFrev be four reactions. Then, the following holds:

• If Si=αSj for some α>0, then ip and jp for all pBlk.

• If Si=αSj for some α<0, then pi (resp. pj) for all pBlk∪{j} (resp. pBlk∪{i}).

• If Si=αSk for some α≠0, then ip and pi for all pBlk and qk for all qBlk∪{i}.

• If Sk=αSl for some α≠0, then kp and pk for all pBlk∪{l} and lq and ql for all qBlk∪{k}.

Proof

The proofs of the four statements are similar. We only consider the first one. Suppose Si=αSj for some α>0 and let us prove ip. Let pBlk. There exists vC such that vp≠0. Let w R n such that wi=0,wj=αvi + vj and wq=vq for all qi,j. We have wC,wi=0,wp≠0 and so the claim follows. □

The next observation states that metabolites which are involved only in irreversible reactions and consumed or produced by exactly one reaction define trivial directional couplings between these reactions.

Observation 5 (Trivial Directional Coupling (TDC))

Let k be some metabolite such that Skl=0 for all lFrevPrev. Let P={iSki>0} and N={jSkj<0}. If card(P)=1 (resp. card(N)=1), then i = 0 j (resp. j = 0 i ) for all (i,j)∈P×N.

Since directional flux coupling is a transitive relation, the flux (un)coupling between many reactions can simply be deduced from dependencies between reactions whose flux (un)coupling has been determined beforehand. This allows us to significantly reduce the total number of LP problems to be solved. Examples of such inferred flux (un)couplings are given in Figure 1(b). According to the TDC rule, we have ( 1 , 2 ) = 0 ( 8 , 9 ) . By solving the LP problems (4), we obtain 10↦(8,9). We can easily conclude that 10↦(1,2).

Table 2 shows the inferred flux (un)coupling relations we can deduce from known relationships between reactions.

Table 2. Transitivity inferred flux (un)coupling

For some pairs of reactions, we need to solve at least one LP problem. The optimal solution not only determines the flux coupling between the considered pair of reactions, but also allows one to determine many other uncoupled reactions.

Observation 6 (Feasibility Rule)

Let vC be a steady state flux vector and let I={ivi=0} and J={jvj≠0}. Then ij for all (i,j)∈I×J.

In general, most irreversible reactions are uncoupled to each other. Accordingly, the LP problems (4) used to determine coupled irreversible reactions are often unbounded. This limits the use of the feasibility rule, which requires the calculation of a feasible flux vector. In order to optimally use the feasibility rule, instead of solving the LP problems (4) to decide whether two irreversible reactions i,jIrev are coupled, we propose to solve the bounded LP problems

L ij = min { v i : Sv = 0 , v j = 1 , v k 0 , k Irr } , L ji = min { v j : Sv = 0 , v i = 1 , v k 0 , k Irr } , (5)

and to use the following scheme:

i = 0 j (resp. j = 0 i ) if and only if Lij≠0 (resp. Lji≠0),

ji if and only if Lij≠0, Lji≠0 and Lij=1/Lji.

The following observation states that removing a fully reversible reaction does not alter the flux coupling between (pseudo-)irreversible reactions.

Observation 7

Let kFrev be a fully reversible reaction. For convenience, suppose k=n. Let S be the m×(n−1) matrix defined by S p = S p for all pk and let Irr=Irr. Then, for all pairs of reactions iFrev and jFrev,

i = 0 j if and only if v i = 0 implies v j = 0 , for all v R n 1 with Sv=0 and v p 0 for all p Irr .

ij if and only if there exists λ≠0 such that v j = λ v i , for all v R n 1 with Sv=0 and v p 0 for all p Irr .

Let SRev be the submatrix defined by the columns in S corresponding to the reversible reactions and let t be the dimension of kern(SRev). Based on Observation 7, we can remove up to t independent fully reversible reactions without altering the flux coupling between (pseudo-)irreversible reactions. Since certain fully reversible reactions may change their reversibility type after the deletion of a fully reversible reaction, we first remove a randomly chosen reaction iFrev together with the coupled reactions with i. We calculate the impact of this deletion on the dimension of kern(SRev). If this dimension decreases by 1, the deletion is maintained; otherwise the removed reactions are restored to the network. This is repeated until t independent fully reversible reactions and their coupled reactions are removed. We assume that the flux coupling between fully reversible reactions is determined beforehand.

Based on the above results, we propose the Fast Flux Coupling Calculator (F2C2) to determine coupled reactions. The main steps of the F2C2 algorithm are given in Table 3.

Table 3. The main steps of the F2C2 algorithm

Results and discussion

The F2C2 algorithm has been implemented within the MATLAB environment, using CLP (via the Mexclp [40] interface) as the underlying linear programming solver. For benchmarking, we analyzed the following genome-scale metabolic networks: Escherichia coli, iJR904 [1], Saccharomyces cerevisiae, iND750 [2], Methanosarcina barkeri, iAF692 [3], Mycobacterium tuberculosis, iNJ661 [4], Escherichia coli, iAF1260 [5], Homo sapiens, Recon 1 [6] and Escherichia coli, iJO1366 [7]. For the numerically sensitive parts, a tolerance level of 10e-6 was set. All computations were performed using a single Intel Xeon 5160 (3.0 GHz) processor on a 64-bit Debian Linux 6.0 system.

As pointed out in the previous section, part of the performance gain of F2C2 over previous FCA algorithms stems from the fact that the preprocessing steps reduce the network size. This affects the running time on two levels: there are fewer reaction pairs and the LP problems to be solved have reduced size. The dramatic effect of the preprocessing steps on the network size is presented in Table 4.

Table 4. Genome-scale metabolic networks with the number of metabolites (met.) and reactions (reac.) before and after applying the preprocessing steps

The algorithmic improvements further reduce the number of LP problems to be solved. A direct performance comparison between the FFCA and F2C2 algorithms (including the running times and number of LP problems solved) is summarized in Table 5. In all cases, F2C2 outperformed FFCA by several orders of magnitude. In [31] it has been shown that FFCA is more efficient on genome-scale metabolic networks than other flux coupling algorithms.

Table 5. Performance comparison between the FFCA and F2C2 algorithms in terms of the number of LP problems solved (LPs) and their total running times (TRT)

For an intuitive, visual representation of the individual improvements’ impact on the algorithm’s performance, a more in-depth analysis has been performed on the recent metabolic network of E. coli, iJO1366. Four different sets of improvements were cumulatively switched on, and the linear programs solved were plotted for each case (Figure 2). To better highlight the relevant differences, the following modifications were applied to the results. First, 249 (out of 2582) reactions identified as blocked were removed from the images. This is a common preprocessing step in most FCA algorithms. Secondly, the order of reactions was permuted such that the reactions in Irev, Prev and Frev are clustered together. Additionally, in each of these three clusters, the fully coupled reactions were moved towards the end of the segment.

thumbnailFigure 2. Visualization of the LP problems solved using different algorithms. (a) The FFCA algorithm, (b) the FFCA algorithm after applying the preprocessing procedure given in Table 1 without kernel computation (Step 7), (c) the FFCA algorithm after applying the preprocessing procedure and using the kernel of the stoichiometric matrix to identify fully coupled reactions and (d) the F2C2 algorithm given in Table 3. The dashed lines mark the boundary of the Irev, Prev and Frev regions. Colors: Black - an LP problem is solved for the corresponding (ordered) pair of reactions; Gray - the corresponding LP problem is not solved due to one (or both) reactions being eliminated from the network (a preprocessing improvement); White - corresponding LP problem does not get solved due to an algorithmic improvement.

Figure 2(a) plots the LP problems solved in the FFCA algorithm. Applying the simple preprocessing steps without using the kernel (Figure 2(b)), several reactions are found to be fully coupled with others, and as such can be merged together. When Lemma 1 is applied (Figure 2(c)), all fully coupled sets are detected. As a consequence the gray stripes get thicker and the LP problems corresponding to (Prev, Prev) and (Frev, Frev) pairs need not be solved anymore. The use of the algorithmic improvements (Figure 2(d)) filters the pairs in (Irev, Irev) and (Irev, Prev) blocks, greatly reducing the total number of LP problems solved.

Conclusions

We have presented the new flux coupling algorithm F2C2, which outperforms previous methods by orders of magnitude. Flux coupling analysis of genome-scale metabolic networks can now be performed in a routine manner. A software tool F2C2 is freely available for non-commercial use at https://sourceforge.net/projects/f2c2/files/ webcite.

Abbreviations

F2C2: Fast Flux Coupling Calculator; FCA: Flux Coupling Analysis; FFCA: Feasibility-based Flux Coupling Analysis; FCF: Flux Coupling Finder; LP: Linear Program; TDC: Trivial Directional Coupling; TFC: Trivial Full Coupling.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

The F2C2 algorithm was designed by all authors, based on ideas from AL and LD. The implementation and computational experiments were done by LD. The first draft of the manuscript was written by AL, and revised by LD and AB. The final document was approved by all authors.

References

  1. Reed JL, Vo TD, Schilling CH, Palsson BØ: An expanded genome-scale model of Escherichia coli K-12 (iJR904 GSM/GPR).

    Genome Biol 2003, 4:R54. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  2. Duarte NC, Herrgard MJ, Palsson BØ: Reconstruction and validation of Saccharomyces cerevisiae iND750, a fully compartmentalized genome-scale metabolic model.

    Genome Res 2004, 14(7):1298-1309. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  3. Feist AM, Scholten JCM, Palsson BØ, Brockman FJ, Ideker T: Modeling methanogenesis with a genome-scale metabolic reconstruction of Methanosarcina barkeri.

    Mol Syst Biol 2006, 2:2006.0004. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  4. Jamshidi N, Palsson BØ: Investigating the metabolic capabilities of Mycobacterium tuberculosis H37Rv using the in silico strain iNJ661 and proposing alternative drug targets.

    BMC Syst Biol 2007, 1:26. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  5. Feist AM, Henry CS, Reed JL, Krummenacker M, Joyce AR, Karp PD, Broadbelt LJ, Hatzimanikatis V, Palsson BØ: A genome-scale metabolic reconstruction for Escherichia coli K-12 MG1655 that accounts for 1260 ORFs and thermodynamic information.

    Mol Syst Biol 2007, 3:121. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  6. Duarte NC, Becker SA, Jamshidi N, Thiele I, Mo ML, Vo TD, Srivas R, Palsson BØ: Global reconstruction of the human metabolic network based on genomic and bibliomic data.

    Proc Natl Acad Sci U.S.A 2007, 104(6):1777-1782. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  7. Orth JD, Conrad TM, Na J, Lerman JA, Nam H, Feist AM, Palsson BØ: A comprehensive genome-scale reconstruction of Escherichia coli metabolism-2011.

    Mol Syst Biol 2011, 7:535. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  8. Covert MW, Famili I, Palsson BØ: Identifying constraints that govern cell behavior: a key to converting conceptual to computational models in biology?

    Biotechnol Bioeng 2003, 84(7):763-772. PubMed Abstract | Publisher Full Text OpenURL

  9. Price ND, Reed JL, Palsson BØ: Genome-scale models of microbial cells: evaluating the consequences of constraints.

    Nat Rev Microbiol 2004, 2(11):886-897. PubMed Abstract | Publisher Full Text OpenURL

  10. Schellenberger J, Que R, Fleming R, Thiele I, Orth J, Feist A, Zielinski D, Bordbar A, Lewis N, Rahmanian S, Kang J, Hyduke D, Palsson B: Quantitative prediction of cellular metabolism with constraint-based models: the COBRA Toolbox v2.0.

    Nat Protoc 2011, 6:1290-1307. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  11. Clarke BL: Stability of complex reaction networks. In Adv. Chem. Phys. Volume 43. Edited by Prigogine I, Rice SA. John Wiley & Sons, Inc., Hoboken; 1980:1-216. OpenURL

  12. Edwards JS, Palsson BØ: Metabolic flux balance analysis and the In silico analysis of Escherichia coli K-12 gene deletions.

    BMC Bioinf 2000, 1:1. BioMed Central Full Text OpenURL

  13. Segrè D, Vitkup D, Church G: Analysis of optimality in natural and perturbed metabolic networks.

    Proc Natl Acad Sci USA 2002, 99(23):15112-15117. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  14. Kauffman KJ, Prakash P, Edwards JS: Advances in flux balance analysis.

    Curr Opin Biotechnol 2004, 14(5):491-496. OpenURL

  15. Shlomi T, Berkman O, Ruppin E: Regulatory on/off minimization of metabolic flux changes after genetic perturbations.

    Proc Natl Acad Sci USA 2005, 102(21):7695-7700. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  16. Lee JM, Gianchandani EP, Papin JA: Flux balance analysis in the era of metabolomics.

    Brief Bioinf 2006, 7(2):140-150. Publisher Full Text OpenURL

  17. Schilling CH, Letscher D, Palsson BØ: Theory for the systemic definition of metabolic pathways and their use in interpreting metabolic function from a pathway-oriented perspective.

    J Theor Biol 2000, 203(3):229-248. PubMed Abstract | Publisher Full Text OpenURL

  18. Schuster S, Hilgetag C, Woods JH, Fell DA: Reaction routes in biochemical reaction systems: algebraic properties, validated calculation procedure and example from nucleotide metabolism.

    J Math Biol 2002, 45(2):153-181. PubMed Abstract | Publisher Full Text OpenURL

  19. Voss K, Heiner M, Koch I: Steady state analysis of metabolic pathways using Petri nets.

    In Silico Biol 2003, 3(3):367-387. PubMed Abstract | Publisher Full Text OpenURL

  20. Larhlimi A, Bockmayr A: A new constraint-based description of the steady-state flux cone of metabolic networks.

    Discrete Appl Math 2009, 157:2257-2266. Publisher Full Text OpenURL

  21. Burgard AP, Nikolaev EV, Schilling CH, Maranas CD: Flux coupling analysis of genome-scale metabolic network reconstructions.

    Genome Res 2004, 14(2):301-312. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  22. Notebaart RA, Kensche PR, Huynen MA, Dutilh BE: Asymmetric relationships between proteins shape genome evolution.

    Genome Biol 2009, 10:R19. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  23. Pál C, Papp B, Lercher MJ: Adaptive evolution of bacterial metabolic networks by horizontal gene transfer.

    Nat Genet 2005, 37:1372-1375. PubMed Abstract | Publisher Full Text OpenURL

  24. Yizhak K, Tuller T, Papp B, Ruppin E: Metabolic modeling of endosymbiont genome reduction on a temporal scale.

    Mol Syst Biol 2011, 7:479. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  25. Notebaart RA, Teusink B, Siezen RJ, Papp B: Co-regulation of metabolic genes is better explained by flux coupling than by network distance.

    PLoS Comput Biol 2008, 4:e26. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  26. Montagud A, Zelezniak A, Navarro E, de Córdoba P F, Urchueguía JF, Patil KR: Flux coupling and transcriptional regulation within the metabolic network of the photosynthetic bacterium Synechocystis sp. PCC6803.

    Biotechnol J 2011, 6:330-342. PubMed Abstract | Publisher Full Text OpenURL

  27. Szappanos B, Kovács K, Szamecz B, Honti F, Costanzo M, Baryshnikova A, Gelius-Dietrich G, Lercher M, Jelasity M, Myers C, Andrews B, Boone C, Oliver S, Pál C, Papp B: An integrated approach to characterize genetic interaction networks in yeast metabolism.

    Nat Genet 2011, 43(7):656-662. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  28. Suthers PF, Chang YJ, Maranas CD: Improved computational performance of MFA using elementary metabolite units and flux coupling.

    Metab Eng 2010, 12:123-128. PubMed Abstract | Publisher Full Text OpenURL

  29. Bundy JG, Papp B, Harmston R, Browne RA, Clayson EM, Burton N, Reece RJ, Oliver SG, Brindle KM: Evaluation of predicted network modules in yeast metabolism using NMR-based metabolite profiling.

    Genome Res 2007, 17:510-519. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  30. Lee DS, Park J, Kay KA, Christakis NA, Oltvai ZN, Barabasi AL: The implications of human metabolic network topology for disease comorbidity.

    Proc Natl Acad Sci USA 2008, 105(29):9880-9885. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  31. David L, Marashi SA, Larhlimi A, Mieth B, Bockmayr A: FFCA: a feasibility-based method for flux coupling analysis of metabolic networks.

    BMC Bioinf 2011, 12:236. BioMed Central Full Text OpenURL

  32. Xi Y, Chen YPP, Qian C, Wang F: Comparative study of computational methods to detect the correlated reaction sets in biochemical networks.

    Brief Bioinf 2011, 12(2):132-150. Publisher Full Text OpenURL

  33. Larhlimi A, Bockmayr A: A new approach to flux coupling analysis of metabolic networks.

    Computational Life Sciences II , Second International Symposium (CompLife 2006),

    Cambridge, UK, Volume 4216 of Lecture Notes in Computer Science. 2006:205–215

    OpenURL

  34. Larhlimi A: New concepts and tools in constraint-based analysis of metabolic networks. PhD thesis, Freie Universität Berlin 2008

  35. Schrijver A: Theory of Linear and Integer Programming. John Wiley & Sons Inc., NY, USA; 1986. OpenURL

  36. Pfeiffer T, Sánchez-Valdenebro I, Nun̄o JC, Montero F, Schuster S: METATOOL: for studying metabolic networks.

    Bioinformatics 1999, 15:251-257. PubMed Abstract | Publisher Full Text OpenURL

  37. Klamt S, Stelling J, Ginkel M, ED G: FluxAnalyzer: exploring structure, pathways, and flux distributions in metabolic networks on interactive flux maps.

    Bioinformatics 2003, 19(2):261-269. PubMed Abstract | Publisher Full Text OpenURL

  38. Gagneur J, Klamt S: Computation of elementary modes: a unifying framework and the new binary approach.

    BMC Bioinf 2004, 5:175. BioMed Central Full Text OpenURL

  39. Urbanczik R, Wagner C: Functional stoichiometric analysis of metabolic networks.

    Bioinformatics 2005, 21(22):4176-4180. PubMed Abstract | Publisher Full Text OpenURL

  40. Löfberg J: Mexclp. 2006, [http://control.ee.ethz.ch/johanl/clp.php]