Constraint-based analysis of metabolic networks has become increasingly important for describing and predicting the behavior of living organisms [1, 2]. While a growing number of metabolic network reconstructions have become available during the last years, the computational analysis of genome-scale networks with hundreds or thousands of reactions may still be very time-consuming. Therefore, there is a need for more efficient algorithms and tools (see e.g. [3–5]).

Flux coupling analysis (FCA) [6] is a useful method for finding dependencies between fluxes of a metabolic network at steady-state. If a non-zero flux through a reaction in steady-state implies a non-zero flux through another reaction, then the first reaction is said to be coupled to the second reaction. Several studies have used FCA for exploring various biological questions such as network evolution [7–9], gene essentiality [7], analysis of experimentally measured fluxes [10, 11] or gene regulation [12, 13]. Having a time efficient implementation of FCA is important in such studies.

In the rest of this section, we first recall some basic definitions. Then we briefly review the previously proposed FCA methods.

### Mathematical definitions

#### Basic preliminaries

For a metabolic network with *m* (internal) metabolites and *n* reactions, the *stoichiometric matrix S* is an *m* × *n* matrix, where element *S*
_{
ij
} is the stoichiometric coefficient of metabolite *i* in reaction *j*. We denote by *R* = {1, ..., *n*} the set of reactions and by M = {1, ..., *m*} the set of metabolites.

There are two types of reactions in a metabolic network. *Irr* is the set of *irreversible* reactions which are assumed to always have non-negative flux values. *Rev* is the set of *reversible* reactions which are allowed to have positive, negative or zero flux values. If a reversible reaction has a positive (resp. negative) flux, we say that it is working in *forward* (resp. *backwar*d) direction. *Splitting* a reversible reaction *i* means making reaction *i* irreversible and adding one more irreversible reaction *i* + *n* to the network, which works in the backward direction. Without loss of generality, we assume that in the numbering of reactions, the first |*Rev*| reactions are the reversible ones.

If *v* is the vector of metabolic fluxes through all the reactions, the *steady*-*state* condition holds for *v* if *S* · *v* = 0. The steady-state *flux cone* is defined as *C* = {*v* ∈ ℝ^{
n
}| *S* · *v* = 0, *v*
_{
i
}
*≥*0, for all *i* ∈ *Irr*}. The lineality space of *C* is given by *lin*.*space*(*C*) = {*v* ∈ ℝ^{
n
}| *S* · *v* = 0, *v*
_{
i
} = 0, for all *i* ∈ *Irr*}.

If for some reaction *r* ∈ *R*, *v*
_{
r
} = 0 for all *v* ∈ *C*, then *r* is called *blocked*, otherwise *r* is *unblocked*. The set of blocked reactions will be denoted by *Blk*. The set of unblocked reactions can be further partitioned based on the *reversibility type* of reactions [14]. We define *Irev* as the set of all reactions that can work only in one direction, i.e., those reactions that take only non-negative or only non-positive flux values at steady-state. If *i* ∈ *Irev* and *v*
_{
i
} ≤ 0 for all *v* ∈ *C*, one can multiply by -1 the *i*-th column of the stoichiometric matrix, such that in the modified network, we have *v*
_{
i
} ≥ 0 for all *v* ∈ *C*. Without loss of generality, we assume that all reactions in *Irev* are operating in the forward direction.

The set of reactions that can work in both directions at steady-state is further partitioned into two subsets: (*a*) *Prev*, the set of *pseudo*-*irreversible* reactions, and (*b*) *Frev*, the set of *fully reversible* reactions. A reaction *i* is in *Prev* if it can work in both directions, and additionally, for all flux vectors *v* in the lineality space of the flux cone, we have *v*
_{
i
} = 0. Accordingly, we define *Frev* as the set of those reactions, which can work in both directions, and additionally, can have a non-zero flux value if fluxes through all irreversible reactions are set to zero.

#### Flux coupling relations between a pair of reactions [6]

Consider two unblocked reactions *i*, *j* ∈ *R*. If for all *v* ∈ *C*, *v*
_{
i
} ≠ 0 implies *v*
_{
j
} ≠ 0, then *i* is *directionally coupled* to *j*. This will be denoted by *i* → *j*, or equivalently, *j* ← *i*. If for all *v* ∈ *C*, *v*
_{
i
} ≠ 0 implies *v*
_{
j
} ≠ 0 and vice versa, then we say that *i* and j are *partially coupled*, or *i* ↔ *j*. If *i* and *j* are partially coupled, and additionally there exists a constant *k* such that for all *v* ∈ *C*, *v*
_{
i
} ≠ 0 implies *v*
_{
j
}/*v*
_{
i
} = *k*, then *i* and *j* are *fully coupled*, or *i* ⇔ j. If neither *i* → *j* nor *j* → *i* holds, then *i* and *j* are *uncoupled*, or
.

#### Elementary flux patterns

The set *supp*(*v*) = {*i* ∈ *R* | *v*
_{
i
} ≠ 0} of reactions active in *v* is called the *support* of *v*. Suppose *v* ∈ *C* is a flux vector and *A* ⊆ *R* is a subnetwork. The flux pattern of *v* for *A* is defined as *A* ∩ *supp*(*v*), which is the set of those reactions in *A* which have non-zero values in *v*[15]. A flux pattern is called an *elementary flux pattern* (EFP) if it cannot be written as the union of other flux patterns. For studying EFPs, Kaleta et al. [15] assume that the network contains only irreversible reactions. To achieve this, every reversible reaction should be split into two irreversible reactions (forward and backward).

### Approaches to Flux Coupling Analysis

In this subsection, we briefly recall existing approaches to flux coupling analysis. For additional information and the technical details on the implementation of these algorithms, we refer to Additional file 1.

#### Flux Coupling Finder Algorithm (FCF)

The most widely used method for FCA is the Flux Coupling Finder (FCF) algorithm [6]. This approach is based on solving linear programs (LPs), and therefore is an optimality-based method. After finding blocked reactions and splitting reversible reactions, for every pair of unblocked reactions *i* and *j*, two LPs are solved. Depending on the optimal values obtained, the coupling relation between *i* and *j* is determined. There is a post-processing step in FCF. Since the reversible reactions have been split, flux coupling relations for these reactions have to be obtained from the coupling relations of the corresponding irreversible forward and backward reactions.

The FCF algorithm has been successfully used for finding coupling relations in a number of metabolic networks [6–13, 16]. However, this approach is rather time-consuming for genome-scale metabolic networks with thousands of reactions, although it is still one of the fastest FCA methods.

#### FCA based on Minimal Metabolic Behaviors (MMB-FCA)

Larhlimi and Bockmayr [14, 17] have proposed a different strategy for flux coupling analysis. In this approach, a minimal set of generating vectors of the flux cone is computed. Then, the coupling relation for any pair of reactions is inferred based on the co-appearance of non-zero fluxes in the generating vectors.

#### FCA based on elementary flux patterns (EFP-FCA)

Recently, Kaleta et al. [15] introduced the concept of elementary flux patterns (EFPs) for the analysis of minimal active reactions in a "subnetwork", which account for possible steady-state flux distributions in a (much) larger metabolic network. They also presented a method based on mixed-integer linear programming (MILP) to compute EFPs. Kaleta et al. suggested that EFPs can be used for characterizing flux coupling relations (see Supplemental Material in [15]). Consider a subnetwork including two unblocked reactions *i* and *j*. If each of these reactions can have a non-zero flux independently of the other (i.e.
), {*i*} and {*j*} are the only EFPs in this subnetwork. On the other hand, if we assume (without loss of generality) that *i* is directionally coupled to *j*, then the EFPs of this subnetwork are {*i*, *j*} and {*j*}. Finally, if the reactions are partially coupled, we will have only one EFP, which is {*i*, *j*}. With this method, it is not possible to distinguish between partial and full coupling, since flux patterns only contain the information about the activity or inactivity of the fluxes, but not the flux values.

### Further strategies for improving flux coupling analysis

In FCF, and some other FCA approaches, every reversible reaction is split into a forward and a backward reaction. This splitting procedure results in an increase in the number of LPs and also in the size of each LP to be solved. Moreover, a complicated postprocessing step is required to infer the flux coupling relations of the original reversible reactions. For all these reasons, it might be better to perform flux coupling analysis without splitting the reversible reactions. An implementation of FCF without splitting is presented in Additional file 1.

Additionally, Larhlimi and Bockmayr [14] show that depending on the reversibility type of the reactions, only certain flux coupling relations can occur. Two unblocked reactions *i* and *j* can be coupled only if one of the following 4 cases holds (note that initially 3 × 3 = 9 reversibility types for the pair (*i*, *j*) would be possible):

1. *i*, *j* ∈ *Irev*: In this case, *i* and *j* can be directionally, partially or fully coupled.

2. *i* ∈ *Irev* and *j* ∈ *Prev*: The only possibility is *j* → *i*.

3. *i*, *j* ∈ *Prev*: In this case, we can only have *i* ⇔ *j*.

4. *i*, *j* ∈ *Frev*: In this case, we can only have *i* ⇔ *j*.

Therefore, it is enough to determine the reversibility types of *i* and *j*, and then check if the corresponding coupling relation holds. This will be referred to as "Reversibility-Type prunings" (or simply, RT-prunings). RT-prunings were originally used in MMB-FCA as presented in [14].

After RT-prunings, a further improvement can be applied to LP-based FCA methods. Assume that in the stoichiomtric matrix *S*, the columns corresponding to the blocked reactions have been removed. According to Proposition 6.13 in [18], for two reactions *i* and *j* such that *i*, *j* ∈ *Prev* or *i*, *j* ∈ *Frev*, the following two statements are equivalent:

This means that it is sufficient to solve only one LP:

The two reactions *i* and *j* will be fully coupled if and only if in the above LP
. Note that the irreversibility constraints are not included in the above LP. Therefore, not only solving one LP is sufficient, but also the LP becomes smaller. This improvement will be called *Prev*/*Frev*-based improvement (or simply, PF-improvement).

An improved version of FCF, which takes into account all of the above-mentioned improvements, has been suggested in [18]. Here we implemented this improved version of FCF and compared it with the other FCA approaches (see the Results and Discussion section).

### Goals of the present work

FCA is a promising tool for metabolic network analysis. However, to perform FCA most authors seem to use their personal implementation of the FCF algorithm [9, 13]. Only very recently, an implementation of FCF has been included in the latest versions of the FASIMU software [19, 20]. On the other hand, several approaches for FCA have been proposed in the literature. It is not known which of these methods is the fastest in practice. In this paper, we present a novel "feasibility-based" flux coupling analysis method (FFCA) and compare it to previously existing approaches. A corresponding software tool is freely available for non-commercial use.