Comparison and improvement of algorithms for computing minimal cut sets
 Christian Jungreuthmayer^{1, 2},
 Govind Nair^{1, 2},
 Steffen Klamt^{3} and
 Jürgen Zanghellini^{1, 2}Email author
DOI: 10.1186/1471210514318
© Jungreuthmayer et al.; licensee BioMed Central Ltd. 2013
Received: 18 January 2013
Accepted: 30 October 2013
Published: 6 November 2013
Abstract
Background
Constrained minimal cut sets (cMCSs) have recently been introduced as a framework to enumerate minimal genetic intervention strategies for targeted optimization of metabolic networks. Two different algorithmic schemes (adapted Berge algorithm and binary integer programming) have been proposed to compute cMCSs from elementary modes. However, in their original formulation both algorithms are not fully comparable.
Results
Here we show that by a small extension to the integer program both methods become equivalent. Furthermore, based on wellknown preprocessing procedures for integer programming we present efficient preprocessing steps which can be used for both algorithms. We then benchmark the numerical performance of the algorithms in several realistic mediumscale metabolic models. The benchmark calculations reveal (i) that these preprocessing steps can lead to an enormous speedup under both algorithms, and (ii) that the adapted Berge algorithm outperforms the binary integer approach.
Conclusions
Generally, both of our new implementations are by at least one order of magnitude faster than other currently available implementations.
Keywords
Metabolic network analysis Elementary modes Minimal cut sets Knockout strategies Integer programming Berge’s algorithmBackground
The aim of metabolic engineering is to (re)allocate available cellular resources in order to induce/stimulate cells to produce substances of interest. For instance, by redirecting intracellular carbon fluxes, product yields can be boosted and optimized [1, 2]. However, the identification of engineering targets is not straightforward as cellular metabolism is a highly interconnected and regulated system of reactions. Consequently, naïve interventions sometimes are ineffective or worse, adversely affect other, even quite distant cellular functions. To deal with the complex interactions in cellular metabolism and to identify promising engineering targets several in silico approaches have been developed [39]. Here we are particularly concerned with two algorithms [10, 11], which are based on elementary mode (EM) analysis [12, 13] and eventually compute intervention strategies as minimal cut sets.
EM analysis was successfully used to identify engineering targets for the production of amino acids [14], biofuels [15, 16], and secondary metabolites [17] in various organisms from C. glutamicum[14] to E. coli[15, 16] to S. cerevisiae[18] to A. niger[19]. In fact, EM analysis is ideally suited for metabolic engineering [20, 21] as it allows for an unambiguous and unbiased decomposition of the analyzed network into inseparable, biologically meaningful steadystate pathways. Any intracellular flux distribution can be represented as a properly weighted combination of these EMs. Thus, the full set of EMs describes all possible steadystate functions. Conversely, the cell’s metabolic capabilities can be restricted if EMs are removed from the network. To remove an EM from the network it suffices to delete at least one contributing reaction [12, 13]. However, as each reaction supports more than one EM, other network functionality will be affected, too. Now the question may be phrased as an optimization problem. The task is to find a minimal intervention strategy, which removes all unwanted functionality from the network while, at the same time, keeps desirable network properties.
Recently, Hädicke and Klamt [11] introduced the concept of constrained minimal cut sets (cMCSs) to predict suitable minimal intervention strategies for a given design criterion. They also presented an algorithm (adapted Berge algorithm [11]; see also [22, 23]) by which cMCSs can be computed from EMs. Jungreuthmayer and Zanghellini [10] conceived an alternative method to compute cMCSs by solving a binary integer program (BIP) over the EMs.
By adapting the BIP originally presented in [10] we first demonstrate that both algorithms deliver indeed equivalent results. Inspired by the theory of integer programming, we then develop efficient preprocessing procedures, which allow both methods to handle problems with hundreds of millions of EMs. Finally, by computing intervention strategies in several realistic networks, we benchmark and compare the computational performance of both algorithms.
Methods
EMs are an unbiased way to characterize metabolic networks. An EM is defined by three properties [12, 13]: (i) it is a nontrivial, steady state flux distribution through the network, (ii) it obeys all thermodynamic constraints on reaction reversibilities, and (iii) no subset of an EM exists which also is an admissible flux distribution and obeys (i) and (ii). By this definition, an EM is a minimal, biologically meaningful, steadystate pathway through a network. An EM can be represented as a (flux) vector or by the set of active reactions in the EM. Herein we will mainly use the latter.
In the following we assume that all EMs are known. Several tools to calculate EMs are freely available [2427].
cMCS theory
Here, D and T denote sets of desired and target modes, respectively. The latter contains all EMs, which need to be removed from the network. The former contains all EMs with favorable functionality. An intervention I will be a set of reactions that are deleted (knockedout) in the network. An EM is hit (and becomes inoperable) if at least one reaction of I is part of the EM. The variable n denotes the minimum number of desired EMs, which have to “survive” the intervention. For a given intervention I, we collect all the surviving desired modes in the set D_{ I }.
As each EM represents a unique pathway through a network, removing it from the network means to block that path, which is easily doable by deleting at least one contributing reaction. Thus, to meet condition (2), the task is to find a (minimal) hitting set such that all pathways in T are blocked [see equation (2)]. Mathematically, this problem is also known as dualization of a (hyper)graph, a fundamental problem in discrete mathematics [28]. Several algorithms for calculating hitting sets are available, of which the Berge algorithm [22] has been shown to perform favorably for the problems considered herein [23]. However, minimal hitting sets ensure only that all target modes are hit but do not per se ensure the constraint (3), i.e. the survival of n desired modes.
A simple strategy to fulfill equation (2) in combination with the constraint in equation (3), is to first calculate all possible minimal hitting sets and then, in a second step, to only select those solutions which also obey equation (3). However, the computational performance can be optimized if the constraint equation (3) is checked “on the fly”, leading to the adapted Berge algorithm presented in [11].
A pseudocode of the adapted Berge algorithm can be found in [11], in the following we give a small example to explain basic principles of the Berge algorithm. Consider a hypergraph with hyperedges ε_{1}={a,b} and ε_{2}={a,c} (in our application, ε_{1} and ε_{2} would represent target EMs). The algorithm finds first all minimal hitting sets (cut sets) for the first edge, i.e. γ_{1}={a} and γ_{2}={b}. It then adds the next edge, ε_{2}, and checks whether the already calculated cut sets are also cut sets for the current edge. Since γ_{1} is hitting ε_{2}, γ_{1} is kept unchanged. However, γ_{2} is not a cut set for ε_{2} and, thus, is removed from the list of cut sets. Instead, two new cut sets are created by individually adding each element of ε_{2} to γ_{2}, i.e. γ_{3}={b,a} and γ_{4}={b,c}. To guarantee minimality the algorithm checks if a newly created cut set is a superset of an already existing one. That is, γ_{3} gets removed from the set of cut sets as it is a superset of γ_{1}. Next, a new edge is added to the system and the calculation cycle starts over. Execution stops when all hyperedges are processed. To account for the intervention problem and accelerate the classic algorithm, Hädicke and Klamt suggested to first check if a newly generated cut set is consistent with the constraint (3) and only then check its minimality against all previously calculated cut sets [11]. This modification leads to the adapted Berge algorithm [11] which will be used in the following.
cMCSs can be formulated as a BIP
In a recent paper [10] we showed that if D=n then the intervention I=I(T,D,n) is representable as a BIP. However, even the general problem that at least n out of D modes need to survive the intervention can be formulated as a BIP.
b^{ i } indicates whether reaction i is part of the EM e.
A solution x to equation (1) can be found by solving the following BIP:
Here we used the indices d and t as a reminder that the EM vectors, b_{ i }, are elements of the sets D and T, respectively. The solution vector, x, is the binary representation of a single cMCS, where x^{ i }=0 marks reactions which get deleted, while x^{ i }=1 stands for reactions that remain unaffected by the genetic intervention. The elements of the binary, auxiliary vector, y, indicate whether or not a desirable mode survives the intervention (1 and 0, respectively). Note that our notation uses superscripts to denote coordinates of vectors and subscripts to denote different vectors. Finally, $\left\right\mathit{x}\left\right:=\sum _{i=1}^{m}{x}^{i}$ represents the multilinear norm of x.
Suppose y^{ d }=1, then equation (5c) always holds and can be omitted. Equation (5b) requires that ${x}^{i}\ge {b}_{d}^{i}$, ∀i∈{1,…,m}. Only then the product of ${\mathit{b}}_{d}^{T}$ and x is equal to the norm of b_{ d }. Thus, b_{ d } is included in the final design. In contrast to b_{ d }, b_{ t } will be removed from the network as equation (5d) requires that the product ${\mathit{b}}_{t}^{T}\mathit{x}$ is smaller than the b_{ t }. This is the case only if at least one reaction in b_{ t } is deleted. Except for equation (5e), the systems of equations in this case resembles the BIP problem presented in Jungreuthmayer and Zanghellini [10].
If y^{ d }=0, then equation (5b) is ineffective. Equation (5c) however simplifies into a “kill constraint”, thus eliminating b_{ d } from the surviving modes.
The binary auxiliary variables y=(y^{1},…,y^{D})^{T} were introduced to guarantee that at least n out of D modes survive the intervention. In both cases y counts the number of surviving desired modes, and equation (5e) makes sure that at least n desired modes survive the intervention.
Alternative MCSs may be calculated by excluding existing solutions x_{ j } by adding the following constraints [10] to the set of equations (5a)(5g):
where we used 1 to denote an allonevector. Equation (6a) guarantees that new solutions are found in subsequent steps, whereas equation (6b) prevents the calculation of solutions that are supersets of already existing solutions. Note that the term 1−x_{ j } represents the binary complement of x_{ j }.
The number of constraints added to the BIP can almost be cut in half (in fact, n/2−1) by keeping in mind that the norm of the current solution x_{ j } will never be larger than the previous optimum x_{j−1}. To sequentially calculate all MCSs the full BIP reads
If iteratively applied, the BIP in equation (7) returns all MCSs, x_{ j }, sorted in increasing order of deletions. Note that although the constraint in equation (7e) is redundant, it significantly enhances the computational performance of the BIP solver.
Preprocessing methods
Mathematically, BIPs are classified as NPhard problems. However, extensive research has focused on improving the formulation of BIPs. The basic idea is to use simple logic rules which turn a BIP into a “better” formulation, which is easier to solve [29]. Standard BIP preprocessing rules essentially fix variables, improve bounds or detect inactive constraints [29].
In the following we will be concerned with standard BIP preprocessing methods to reduce the size of the intervention problem in equation (1) but not with the internal structure of the algorithms. These preprocessing procedures will allow to reduce the size of the intervention problem in equation (1), which can then be solved by the Berge algorithm or a BIP. In the following, by “Berge algorithm” we mean the adapted Berge algorithm reported by Hädicke and Klamt [11] which extends the standard Berge algorithm to compute only minimal hitting sets (cut sets) that comply with the constraint (3) on the desired modes [11].
We assume that all EMs are converted to their binary representation according to equation (4). Furthermore, we split the complete set of EMs in three sets, D, N, and T. Here the neutral set, N, contains all (binary) EMs, which are neither element of D nor T.
Step 1. First, we remove all reactions that are simultaneously zero in all EMs of T. These reactions do not support any EM in T. Deleting them will not remove any unwanted mode.
Step 2. Next, essential reactions are identified. If deleting a reaction reduces the number of surviving modes in D to less than n [i.e. violates equation (3)], then this reaction is considered to be essential and cannot be knocked out. A reaction i is essential if D−s^{ i }<n, with ${s}^{i}=\sum _{j=1}^{\left\mathit{D}\right}{x}_{j}^{i}$.
Example of determining essential reactions
R1  R2  R3  R4  R5  R6  R7  R8  

${\mathit{d}}_{1}^{T}$  0  1  0  0  1  0  1  0 
${\mathit{d}}_{2}^{T}$  0  0  0  0  0  1  0  0 
${\mathit{d}}_{3}^{T}$  1  0  1  0  1  0  0  0 
${\mathit{d}}_{4}^{T}$  1  1  0  0  0  0  1  0 
${\mathit{d}}_{5}^{T}$  1  0  1  0  0  0  1  0 
s  3  2  2  0  2  1  3  0 
In general, the more essential reactions we find, the more the system can be reduced. Consequently, it is beneficial if n is large (ideally n=D), as this results in the maximum number of essential reactions. Removing all essential reactions from the system is a critical step that opens the possibility to apply several other preprocessing procedures.
The removal of all essential reactions results in an important change of the system. By definition EMs are nondecomposable, i.e. an EM is not a subset of any other EM. However, if the essential reactions are removed then the residual EMs may become subsets or duplicates of other modes. Hence, the next step is to find all duplicate modes in T and to remove them from the system.
Set of (residual) target modes before and after subsetsuperset elimination
R1  R2  R3  R4  R5  R6  R7  R8  t_{ i }  

${\mathit{t}}_{1}^{T}$  0  1  0  0  0  0  0  0  1  
${\mathit{t}}_{2}^{T}$  0  1  0  0  1  0  0  0  2  * 
${\mathit{t}}_{3}^{T}$  0  0  0  1  0  0  1  0  2  
${\mathit{t}}_{4}^{T}$  1  0  0  1  0  0  1  0  3  * 
${\mathit{t}}_{5}^{T}$  0  0  0  1  0  0  1  1  3  * 
${\mathit{t}}_{6}^{T}$  1  0  1  0  0  0  1  1  3 
The same procedure can be applied to the other modes as well. Mode t_{3} has a norm of 2 and is killed either by knocking out reaction R4 or reaction R7. As both reactions are part of t_{4} and t_{5}, they are certainly removed if mode t_{3} is killed.
Step 4. In a final preprocessing step, we remove duplicate reactions across all EMs in both sets, D and T. Using the illustration in Table 2, this would mean that we remove duplicate columns. Note that this step is most effective after all supersets were removed. For instance, in Table 2 columns R1, R3 and R8 are identical only if t_{2}, t_{4}, and t_{5} are removed. In this step it is not possible to analyze D and T separately. Reactions need to be identical in both sets, D and T, in order to be removed.
Implementation
We implemented the BIP algorithm in C using Gurobi Optimizer 5.0, http://www.gurobi.com/ for solving the BIP problem. The adapted Berge algorithm was implemented in C. The software is available from the authors on request. The simulations were all carried out on an Intel Xeon CPU X5650 @ 2.67GHz under a Linux operating system.
Test cases
Topological properties of the E. coli models used
E2  E1  E0  

Metabolites  63  64  74 
 Internal  49  50  52 
 External  14  14  22 
Reactions  60  64  75 
 Irreversible  26  30  36 
 Reversible  34  34  39 
Elementary modes  55,666  485,169  124,341,216 
Cardinalities for the sets involved in the intervention problem I _{ n } ( T , D , n ) with n∈ { 1 , … , n _{ max } }
E2  E1  

D  5,132  (9%)  46,254  (10%) 
N  18,447  (33%)  217,877  (45%) 
T  32,087  (58%)  221,038  (45%) 
n _{max}  1,120  (2%)  11,436  (2%) 
Results
Berge algorithm outperforms the BIP
It is interesting to observe that over a wide range of values for n the runtime in both methods changes according to a power law (see Figure 3). However, only for the Berge algorithm the exponent remained approximately constant in both test cases.
Preprocessingtimes are essentially independent of n and only scale with the total number of processed EMs. For cases with very few MCSs (see Additional file 2: Figure S1) the Berge algorithm took even less time than the preprocessing.
Preprocessing reduces overall computation time
Runtime of the Berge algorithm with and without preprocessing (PP) for two intervention problems I ( T , D , n ) with differing n
E2  E1  E0  

D  489  (0.88%)  489  (0.10%)  489  (0.00%)  
T  55,177  (99.12%)  484,680  (99.90%)  124,340,727  (100.00%)  
n_{2}=40; n_{2}/D=0.08  
MCSs  4  11  44  
Min. deletions  5  8  16  
Max. deletions  5  8  17  
Without PP  With PP  Without PP  With PP  Without PP  With PP  
System size  60×55,177  13×1  64×484,680  17×3  75×124,340,727  28×6  
Reading EMs (sec)  0.031  (26%)  0.035  (43%)  0.281  (26%)  0.308  (37%)  71.858  (24%)  76.274  (26%) 
Preprocessing (sec)  0.078  (65%)  0.044  (55%)  0.717  (65%)  0.495  (60%)  194.328  (64%)  202.539  (69%) 
Calculate MCSs (sec)  0.011  (9%)  0.001  (2%)  0.100  (9%)  0.021  (3%)  35.866  (12%)  15.294  (5%) 
Total (sec)  0.121  (100%)  0.081  (100%)  1.099  (100%)  0.823  (100%)  302.053  (100%)  294.108  (100%) 
n_{2}=40; n_{2}/D=0.08  
MCSs  72  274  1,720  
Min. deletions  5  6  14  
Max. deletions  7  10  19  
Without PP  With PP  Without PP  With PP  Without PP  With PP  
System size  60×55,177  47×2,295  64×484,680  51×8,664  75×124,340,727  62×321,272  
Reading EMs (sec)  0.031  (10%)  0.033  (25%)  0.277  (7%)  0.308  (20%)  71.084  (3%)  77.771  (5%) 
Preprocessing (sec)  0.078  (26%)  0.083  (65%)  0.715  (17%)  1.162  (76%)  193.805  (7%)  1,493.367  (94%) 
Calculate MCSs (sec)  0.196  (64%)  0.012  (10%)  3.136  (76%)  0.068  (4%)  2,475.110  (90%)  15.909  (1%) 
Total (sec)  0.306  (100%)  0.128  (100%)  4.129  (100%)  1.539  (100%)  2,740.004  (100%)  1,587.053  (100%) 
In the most demanding case (n/D≪1), the Berge algorithm with preprocessing identified 1,720 MCSs in less than 30 minutes in the large E0 model with its 124 million EMs (see Table 5). Only 1% of the computation time was used for the Berge algorithm. Ninetyfour percent of the computation is spent on preprocessing. After preprocessing the initial system of 124 million EMs was reduced to approximately 300,000 modes. In all tested cases with enhanced preprocessing, reading EMs and preprocessing took at least 90% of the total computation time. We repeated the same simulation without preprocessing. While the total runtime with and without preprocessing is comparable if only a few MCSs are found, the runtime savings in MCS calculation more than outweigh the runtime losses due to preprocessing if many MCSs solve an intervention problem. To emphasize this point we show the total runtime as function of the number of MCSs for Figure 3 in the Additional file 2: Figure S1.
Runtime analysis for the Berge algorithm and BIP using several examples from the literature with different design objectives
Runtime (sec)  

Organism  Objective  D  T  # MCS  Min.Δ  Max.Δ  Berge  BIP  CNA 
E. coli[16] (anaerobic)  ethanol  12  4,998  1,048  6  9  0.011  0.287  2.83 
E. coli[16] (aerobic)  ethanol  12  429,264  55,488  11  15  0.883  2.174  547.61 
E. coli[15] (anaerobic)  isobutanol  7  5,615  760  7  10  0.011  0.233  2.69 
E. coli[33] (anaerobic)  nbutanol  7  7,341  2,280  7  10  0.015  0.226  3.43 
Preprocessing strongly reduces the system size
Surprisingly, the computation time does not monotonically increase with the number of solutions [i.e. with the number of additional constraints, see equation (7)] but drops dramatically whenever the norm of the solution decreases. Note that a decreasing norm means an increase in the number of required knockouts. In this model the computation time significantly drops after solution number 11, 53, 116, and 156. At these instances the number of required deletions changes from 11 to 12, to 13, to 14, and to 15, respectively. In all these cases the constraint in equation (7e) decreases too and introduces a tighter bound on the system. This allows the solver’s internal preprocessing to more efficiently compress the system, which in turn brings down the computation time. If, however, the norm of the solution does not change then the computation time scales approximately exponentially with the number of MCSs. This behavior is expected as each solution adds a new constrain to the system, which makes it harder to solve.
Discussion
Recently cMCSs have been introduced to predict optimal intervention strategies in order to achieve an arbitrary metabolic objective [11]. Two algorithmic approaches have been published for their calculation [10, 11]. Here we showed that both methods are equivalent. We addressed the numerical efficiency of both methods in typical design problems and found that in terms of runtime the Berge algorithm is superior compared to BIP.
It may appear as a surprise that the Berge algorithm performs so well even for the large cases presented in this study, especially since the Berge method is known for its unfavorable performance in huge networks [34]. However, here we showed that efficient preprocessing can dramatically reduce the size of the networks. The adapted Berge algorithm could then be run on the reduced systems. Apparently, for small systems the Berge algorithm is effective.
The importance of preprocessing in the calculation of MCSs has been stressed earlier [23]. The preprocessing strategies introduced herein focus especially on the additional constraints posed by cMCSs, whereas [23] dealt only with (unconstrained) MCSs. We were able to show that our implementation outperforms the currently available tool for computing (c)MCSs (see Table 6). The performance gain can be attributed to both the improved preprocessing and the efficient implementation in C. Herein we used standard preprocessing routines, which are frequently applied in BIP [29]. Extensive literature on preprocessing in binary and integer programming is available, see for instance Savelsbergh [35] for a good summary of basic ideas. Since cMCSs can be stated as a BIP, these methods are readily adoptable. However, due to the algorithmic complexity of BIP (solving numerous linear optimizations as part of one BIP, etc.), a full enumeration based on BIP seems not be competitive compared to the Berge algorithm (see Figure 3). Rather the usage of BIP preprocessing rules followed by the Berge algorithm to calculate cMCSs is suggested as an optimal computation strategy.
The efficiency of preprocessing is dependent on the imposed design criterion. In the worst case the set of desired modes is empty (D=∅) and T contains all EMs of a network. This situation corresponds to unconstrained MCSs and thus to a full dualization of the hypergraph spanned by the target modes. Except for step 4, none of our preprocessing routines then provides an advantage and other solvers may be more appropriate [34]. However such cases are not relevant in the context of metabolic engineering, where we want to optimize favorable functionality. To fully utilize the potential of preprocessing the ratio n/D should be close to one. This means that many essential reactions will be removed from the system, and as a result of that many supersets will be detected, too. However, in practice it may suffice if only a few EMs out of the set of desired modes survive, i.e. n/D≪1. Still, preprocessing provides a significant performance gain as indicated in Table 5. The runtime costs of preprocessing will be outweighed by the savings in MCS calculation, if the intervention problem has many solutions. In practice, preprocessing will therefore be favorable, as typical applications have a few thousand solutions (see Table 6).
In our paper [10] we used weights in the objective function of the BIP to take experimental difficulties in the deletion of reactions into account. For instance, some reactions cannot be deleted as they are driven by diffusion, rather than catalyzed by an enzyme. Other reactions, on the other hand, may require the deletion of multiple genes as they are catalyzed by different enzymes in parallel. By an appropriate choice of the weights in the objective function BIP is able to predict the experimentally easiest deletion strategies first [10]. However, in the preprocessing procedures above we did not consider weights in the objective function. Identifying particular solutions in the complete list of MCSs has to be done in a separate postprocessing step (for example by appropriately sorting the output, which can be done quite fast). Thus even with an additional postprocessing step our implementation of Berge’s algorithm will be faster than BIP. Note however that the integration of regulatory information into the cMCS framework is a unique feature of the BIP approach [10].
Both methods, the Berge algorithm as well as BIP, still show room for computational improvements. In the case of the Berge algorithm the computational bottleneck sits in the filtering of potential MCSs to determine if they are, in fact, true MCSs and not supersets of true MCSs [23]. Generating new MCScandidates, however, is very quick. Therefore ways of enhancing the supersetfiltering procedure will be the scope of future work.
One disadvantage of the Berge algorithm is its inability to predict MCSs continuously during the runtime. During execution all MCSs remain candidateMCSs. Only upon termination, when the minimality of all candidate MCSs has been checked against each other, candidateMCSs become MCSs and can be outputted. Thus, even if we were interested in only one solution, the Berge algorithm will  in general  return more than one MCS upon termination. However, other solvers are available [34]. Their adaption for the current situation is the scope of further work.
BIP on the other hand, is able to predict a single solution without the need to enumerate all. In fact, due to the optimization principle only one MCS with a smallest or largest number of deletions can be calculated. In Figure 4 we illustrated the runtime per solution as function of the number of MCSs. The drop in runtime after certain solutions indicates that more advanced preprocessing procedures may further reduce the runtime significantly. In fact, our preprocessing focused on standard procedures like variable fixing. More advanced methods will further reduce the runtime for both the Berge algorithm and the BIP. Additionally we used GUROBI, a commercially available multipurpose optimization toolbox, to solve the BIP. However, a specialized knapsack solver may potentially boost the performance.
Conclusions
We predicted minimal metabolic intervention strategies in typical metabolic engineering problems using two different methods (an adapted Berge algorithm and a BIP). We investigated the numerical performance of these approaches. Both methods significantly profited from the enhanced preprocessing procedures developed here. Under the tested conditions, our implementation of Berge’s algorithm performed best even outperforming other, currently available software.
Abbreviations
 BIP:

Binary integer program
 cMCS:

Constrained minimal cut set
 CNA:

CellNetAnalyzer
 EM:

Elementary mode
 MCS:

Minimal cut set
 PP:

Preprocessing.
Declarations
Acknowledgements
C.J., G.N., and J.Z. gratefully acknowledge the support by the Federal Ministry of Economy, Family and Youth (BMWFJ), the Federal Ministry of Traffic, Innovation and Technology (bmvit), the Styrian Business Promotion Agency SFG, the Standortagentur Tirol and ZIT  Technology Agency of the City of Vienna through the COMETFunding Program managed by the Austrian Research Promotion Agency FFG.
S.K. acknowledges funding by the German Federal Ministry of Education and Research(e:Bio  0316183D) and by the Ministry of Education and Research of SaxonyAnhalt (Research Center “Dynamic Systems: Biosystems Engineering”).
Authors’ Affiliations
References
 Lee JW, Na D, Park JM, Lee J, Choi S, Lee SY: Systems metabolic engineering of microorganisms for natural and nonnatural chemicals. Nat Chem Biol. 2012, 8 (6): 536546. 10.1038/nchembio.970.View ArticlePubMedGoogle Scholar
 Xu P, Ranganathan S, Fowler ZL, Maranas CD, Koffas MA: Genomescale metabolic network modeling results in minimal interventions that cooperatively force carbon flux towards malonylCoA. Metab Eng. 2011, 13 (5): 578587. 10.1016/j.ymben.2011.06.008.View ArticlePubMedGoogle Scholar
 Zomorrodi AR, Suthers PF, Ranganathan S, Maranas CD: Mathematical optimization applications in metabolic networks. Metab Eng. 2012, 14 (6): 672686. 10.1016/j.ymben.2012.09.005.View ArticlePubMedGoogle Scholar
 Kim J, Reed J: OptORF: Optimal metabolic and regulatory perturbations for metabolic engineering of microbial strains. BMC Syst Biol. 2010, 4: 5310.1186/17520509453.PubMed CentralView ArticlePubMedGoogle Scholar
 Choi HS, Lee SY, Kim TY, Woo HM: In Silico identification of gene amplification targets for improvement of lycopene production. Appl Environ Microbiolo. 2010, 76 (10): 30973105. 10.1128/AEM.0011510.View ArticleGoogle Scholar
 Ranganathan S, Suthers PF, Maranas CD: OptForce: An optimization procedure for identifying all genetic manipulations leading to targeted overproductions. PLoS Comput Biol. 2010, 6 (4): e100074410.1371/journal.pcbi.1000744.PubMed CentralView ArticlePubMedGoogle Scholar
 Pharkya P, Burgard AP, Maranas CD: OptStrain: A computational framework for redesign of microbial production systems. Genome Res. 2004, 14 (11): 23672376. 10.1101/gr.2872004.PubMed CentralView ArticlePubMedGoogle Scholar
 Burgard AP, Pharkya P, Maranas CD: Optknock: A bilevel programming framework for identifying gene knockout strategies for microbial strain optimization. Biotechnol Bioeng. 2003, 84 (6): 647657. 10.1002/bit.10803.View ArticlePubMedGoogle Scholar
 Segrè D, Vitkup D, Church GM: Analysis of optimality in natural and perturbed metabolic networks. Proc Natl Acad Sci. 2002, 99 (23): 1511215117. 10.1073/pnas.232349399.PubMed CentralView ArticlePubMedGoogle Scholar
 Jungreuthmayer C, Zanghellini J: Designing optimal cell factories: Integer programing couples elementary mode analysis with regulation. BMC Syst Biol. 2012, 6: 10310.1186/175205096103.PubMed CentralView ArticlePubMedGoogle Scholar
 Hädicke O, Klamt S: Computing complex metabolic intervention strategies using constrained minimal cut sets. Metab Eng. 2011, 13 (2): 204213. 10.1016/j.ymben.2010.12.004. http://www.ncbi.nlm.nih.gov/pubmed/21147248,View ArticlePubMedGoogle Scholar
 Schuster S, Fell DA, Dandekar T: A general definition of metabolic pathways useful for systematic organization and analysis of complex metabolic networks. Nat Biotech. 2000, 18 (3): 326332. 10.1038/73786. http://dx.doi.org/10.1038/73786,View ArticleGoogle Scholar
 Schuster S, Dandekar T, Fell DA: Detection of elementary flux modes in biochemical networks: a promising tool for pathway analysis and metabolic engineering. Trends Biotechnol. 1999, 17 (2): 5360. 10.1016/S01677799(98)012906. http://www.ncbi.nlm.nih.gov/pubmed/10087604,View ArticlePubMedGoogle Scholar
 Becker J, Zelder O, Häfner S, Schröder H, Wittmann C: From zero to herodesignbased systems metabolic engineering of Corynebacterium glutamicum for Llysine production. Metab Eng. 2011, 13 (2): 159168. 10.1016/j.ymben.2011.01.003.View ArticlePubMedGoogle Scholar
 Trinh CT, Li J, Blanch HW, Clark DS: Redesigning Escherichia coli metabolism for Anaerobic production of Isobutanol. Appl Environ Microbiol. 2011, 77 (14): 48944904. 10.1128/AEM.0038211.PubMed CentralView ArticlePubMedGoogle Scholar
 Trinh CT, Unrean P, Srienc F: Minimal Escherichia coli cell for the most efficient production of ethanol from Hexoses and Pentoses. Appl Environ Microbiol. 2008, 74 (12): 36343643. 10.1128/AEM.0270807. http://aem.asm.org/content/74/12/3634.abstract,PubMed CentralView ArticlePubMedGoogle Scholar
 Unrean P, Trinh CT, Srienc F: Rational design and construction of an efficient E. coli for production of diapolycopendioic acid. Metab Eng. 2010, 12 (2): 112122. 10.1016/j.ymben.2009.11.002.PubMed CentralView ArticlePubMedGoogle Scholar
 Xu X, Cao L, Chen X: Elementary flux mode analysis for optimized ethanol yield in anaerobic fermentation of glucose with saccharomyces cerevisiae. Chin J Chem Eng. 2008, 16: 135142. 10.1016/S10049541(08)60052X.View ArticleGoogle Scholar
 Driouch H, Melzer G, Wittmann C: Integration of in vivo and in silico metabolic fluxes for improvement of recombinant protein production. Metab Eng. 2012, 14: 4758. 10.1016/j.ymben.2011.11.002.View ArticlePubMedGoogle Scholar
 Zanghellini J, Ruckerbauer DE, Hanscho M, Jungreuthmayer C: Elementary flux modes in a nutshell: properties, calculation and applications. Biotechnol J. 2013, 8 (9): 10091016. 10.1002/biot.201200269.View ArticlePubMedGoogle Scholar
 Trinh C, Wlaschin A, Srienc F: Elementary mode analysis: a useful metabolic pathway analysis tool for characterizing cellular metabolism. Appl Microbiol Biotechnol. 2009, 81 (5): 813826. 10.1007/s0025300817701.PubMed CentralView ArticlePubMedGoogle Scholar
 Berge C: Hypergraphs: Combinatorics of Finite Sets. 1989, (NorthHolland mathematical library: Volume 45.) Amsterdam: Elsevier Science PublishersGoogle Scholar
 Haus UU, Klamt S, Stephen T: Computing knockout strategies in metabolic networks. J Comput Biol. 2008, 15 (3): 259268. 10.1089/cmb.2007.0229.View ArticlePubMedGoogle Scholar
 Jungreuthmayer C, Ruckerbauer DE, Zanghellini J: regEfmtool: Speeding up elementary flux mode calculation using transcriptional regulatory rules in the form of threestate logic. Biosystems. 2013, 113: 3739. 10.1016/j.biosystems.2013.04.002.View ArticlePubMedGoogle Scholar
 Jevremović D, Trinh CT, Srienc F, Sosa CP, Boley D: Parallelization of Nullspace algorithm for the computation of metabolic pathways. Parallel Comput. 2011, 37 (67): 261278. 10.1016/j.parco.2011.04.002.PubMed CentralView ArticlePubMedGoogle Scholar
 Terzer M, Stelling J: Largescale computation of elementary flux modes with bit pattern trees. Bioinformatics. 2008, 24 (19): 22292235. 10.1093/bioinformatics/btn401.View ArticlePubMedGoogle Scholar
 Kamp Av, Schuster S: Metatool 5.0: fast and flexible elementary modes analysis. Bioinformatics. 2006, 22 (15): 19301931. 10.1093/bioinformatics/btl267.View ArticleGoogle Scholar
 Eiter T, Makino K, Gottlob G: Computational aspects of monotone dualization: A brief survey. Discrete Appl Math. 2008, 156 (11): 20352049. 10.1016/j.dam.2007.04.017.View ArticleGoogle Scholar
 Chen DS, Batson RG, Dang Y: Applied Integer Programming: Modeling and Solution. 2010, Hoboken: John Wiley & Sons, Inc.Google Scholar
 Orth JD, Fleming RMT, Palsson BØ: Reconstruction and use of microbial metabolic networks: the core escherichia coli metabolic model as an educational guide. EcoSalEscherichia coli and Salmonella: Cellular and Molecular Biology. Edited by: Böck A, Curtiss IIIR, Kaper JB, Karp PD, Neidhardt FC, Nyström T, Slauch JM, Squires CL, Ussery D. 2009, Washington: ASM Press, 5699. chapter 10.2.1Google Scholar
 Deutscher J: The mechanisms of carbon catabolite repression in bacteria. Curr Opin Microbiol. 2008, 11 (2): 8793. 10.1016/j.mib.2008.02.007.View ArticlePubMedGoogle Scholar
 Klamt S, SaezRodriguez J, Gilles E: Structural and functional analysis of cellular networks with CellNetAnalyzer. BMC Syst Biol. 2007, 1: 210.1186/1752050912.PubMed CentralView ArticlePubMedGoogle Scholar
 Trinh CT: Elucidating and reprogramming Escherichia coli metabolisms for obligate anaerobic nbutanol and isobutanol production. Appl Microbiol Biotechnol. 2012, 95 (4): 10831094. 10.1007/s0025301241977.View ArticlePubMedGoogle Scholar
 Murakami K, Uno T: Efficient algorithms for dualizing largescale hypergraphs. arXiv:1102.3813. 2011, http://arxiv.org/abs/1102.3813,Google Scholar
 Savelsbergh MWP: Preprocessing and probing techniques for mixed integer programming problems. ORSA J Comput. 1994, 6 (4): 445454. 10.1287/ijoc.6.4.445.View ArticleGoogle Scholar
Copyright
This article is published under license to 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.