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
Hädicke and Klamt [11] defined cMCSs as solutions I of an intervention problem
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}.
A proper solution I of equation (1) is a set of reactions obeying two conditions: First, the removal of the reactions in I will delete the complete target set, T, from the network
and no subset of I will do so. This is the MCS property. To be a constrained MCS the intervention I will keep at least n desirable EMs unaffected, i.e.
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.
Let e be an EM of a metabolic network with m reactions, fulfilling the steady state condition, and b=b(e) its binary representation,
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, represents the multilinear norm of x.
Suppose y^{d}=1, then equation (5c) always holds and can be omitted. Equation (5b) requires that , ∀i∈{1,…,m}. Only then the product of 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 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 .
Consider the example in Table 1 with D=5 and n=3. R1 and R7 are essential reactions, as for them D−s^{i}=5−3=2<3=n, which indicates that knocking out R1 or R7 will kill more desirable modes than allowed. We note that if D=n, all active reactions are essential.
Table 1. Example of determining essential reactions
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.
Step 3. Next, we screen T to find and remove residual EMs that are supersets of other residual EMs in T. Consider the target set (of residual EMs), T, shown in Table 2. The modes are sorted in order of ascending norm. The example illustrates that mode t_{1} can be removed by knocking out reaction R2. However, knocking out reaction R2 also kills t_{2}, as t_{2} is a superset of t_{1}.
Table 2. Set of (residual) target modes before and after subsetsuperset elimination
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/ webcite 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
We used the E. coli core model, E0, [30] and two smaller models, E1 and E2, which were derived from the E0 model by removing several reactions. Compared to the E0 model, glucose was considered as the only carbon substrate for E1 and exchange of αketoglutarate, acetaldehyde, acetate, formate, lactate, and pyruvate was not allowed. In addition to these modifications we also removed the glyoxylate shunt and the (NAD and NADP dependent) malic enzymes to obtain the E2 model from E1. All three models are illustrated in Figure 1. Their main topological properties are summarized in Table 3. A list of reactions for these models may be found in the (Additional file 1).
Figure 1. Overview of the different E. coli models. For simplicity we only show pathways. Cofactors etc. are suppressed. Metabolites contributing to the biomass are depicted in yellow. Pathways included in the E2 model are indicated in red. E1 contains the red and blue pathways only, while E0 [30] incloses all reactions, including the noncolored pathways. A detailed listing of all models may be found in the Additional file 1.
Table 3. Topological properties of the E. coli models used
Additional file 1. tar gzipped archive of the SBML files for the models E0, E1, and E2.
Format: ZIP Size: 14KB Download file
To test the numerical efficiency of the implemented MCS algorithms we set up the following intervention problems: We first identify the most efficient EMs in all models. Efficiency is defined as the product between growth rate and ethanol secretion. Next, we classify all EMs to be desirable, whose ethanol secretion is larger or equal than the excretion of the most efficient EMs. Targets are all other modes that do not utilize ethanol. Modes which take up ethanol (negative secretion) are considered neutral, as ethanol uptake is repressed in the presence of glucose in the growth medium [31]. Therefore these modes do not need to be targeted. In Figure 2 we illustrate the intervention problem and the choice of D, N, and T for the E2 model. The major properties of the design criteria for the different E. coli models are listed in Table 4.
Figure 2. Phenotypic plot of all EMs in E2. Flux values are normalized to the total carbon influx (C flux_{in}). EMs are color coded with respect to the modes’ efficiency, defined as the product between normalized growth and normalized ethanol secretion (etoh_{out}). The most efficient EM is marked by an arrow. The areas D, N, and T indicate the corresponding sets of EMs for the intervention problem I_{n}(T,D,n) with n ∈ {1,…,n_{max}}.
Table 4. Cardinalities for the sets involved in the intervention problem I_{n}(T,D,n) with n∈ {1,…,n_{max}}
Results
Berge algorithm outperforms the BIP
Figure 3 shows the computation time to calculate all MCSs using either method as a function of the minimally required number n of surviving desired EMs. We used the design criteria outlined above. At n=2 we found 81,168 and 441,095 MCSs in E2 and E1, respectively. (The number of MCSs as function of n may be found in Additional file 2: Figure S1.) In all tested situations the adapted Berge algorithm clearly outperforms the BIP. Even in the most demanding case (n=2), the Berge algorithm calculated all MCSs in E1 in less than 10 min. On the other hand, it already took the BIP 22 hours to calculate all 331 MCSs for n=85 in the smaller E2 model. In the same situation the Berge implementation finished in 0.4 sec.
Figure 3. Runtime for the Berge algorithm (left panel) and the BIP (right panel) for the models E2 (red) and E1 (blue), respectively. The runtimes for preprocessing (circles) and calculating MCSs are plotted as function of n. n is the lower bound of desired EMs, which must survive the intervention. Table 4 gives the size of the intervention problem for both models. As the trend is clearly visible, we did not evaluate the BIP for n < 80.
Additional file 2. Total runtime with and without preprocessing for the Berge algorithm.
Format: PDF Size: 70KB Download file
This file can be viewed with: Adobe Acrobat Reader
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
To test the impact of our preprocessing procedures we set up identical intervention problems for all models. That is, we solved I_{0}=I(T_{0},D,n), I_{1}=I(T_{1},D,n), and I_{2}=I(T_{2},D,n), where we used the indices 0, 1, and 2 to denote the dependence on the models E0, E1, and E2, respectively. We used identical sets of desired EMs in all models, i.e. D_{0}=D_{1}=D_{2}=D and n_{0}=n_{1}=n_{2}=n. T_{i},i∈{0,1,2} consisted of all EMs not contained in D. Values for D, n, T_{i} and the runtimes for the Berge algorithm in two different cases (n/D≈1 and n/D≪1) may be found in Table 5.
Table 5. Runtime of the Berge algorithm with and without preprocessing (PP) for two intervention problems I(T,D,n) with differing n
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.
Finally in Table 6 we show several examples from the literature, which can be easily and efficiently solved by either method. As a comparison we have also listed runtimes using the current version (version 2012.1) of CellNetAnalyzer (CNA) [32]. CNA uses a MATLAB implementation of the Berge algorithm. However, its preprocessing capabilities are less developed. That is why both programs, our Bergealgorithm and BIP, outperform CNA in all instances by at least one order of magnitude. Note however that CNA uses a MATLAB script, while our programs are implemented in C. A significant part of the performance difference may therefore be attributed to the slower performance of MATLAB compared to native executables written in C.
Table 6. Runtime analysis for the Berge algorithm and BIP using several examples from the literature with different design objectives
Preprocessing strongly reduces the system size
In Figure 4 we used a BIP and show the computation time as a function of the number of MCSs for the aerobic E. coli[33] model of Trinh et al.[16] (see line number 2 in Table 6 for model details). Note that although Table 6 lists 55,488 different MCSs, the BIP (and our Berge algorithm for that matter) only needs to calculate 164 solutions. Due to preprocessing the original network is reduced from 71 reactions and 429,275 EMs to an equivalent system with only 23 columns and 28 rows. This smaller system has 164 MCSs, which are then reconstructed to the full set of MCSs by expanding duplicated columns. A similar observation may also be made in Table 5. In these examples the system size is at least reduced by a factor of 30 (case E2, n_{2}=40).
Figure 4. Computation time of the BIP as function of the number of MCSs for the aerobicE. coli[33]model of Trinhet al.[16](line number 2 in Table6). The accumulated computation time for the respective models is plotted in the lower panel. The top panel shows the computation time for each solution. Note that it is impossible to calculate the runtime per solution for the Berge algorithm as the algorithm does not allow to continuously output the solution.
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.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
CJ implemented the algorithms and participated in the design of the study. GN coded the models and helped in running the analysis. SK participated in the design of the study. JZ conceived of the study, and participated in its design and coordination. CJ, SK and JZ wrote the manuscript. All authors read and approved the final manuscript.
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”).
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. PubMed Abstract  Publisher Full Text

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. PubMed Abstract  Publisher Full Text

Zomorrodi AR, Suthers PF, Ranganathan S, Maranas CD: Mathematical optimization applications in metabolic networks.
Metab Eng 2012, 14(6):672686. PubMed Abstract  Publisher Full Text

Kim J, Reed J: OptORF: Optimal metabolic and regulatory perturbations for metabolic engineering of microbial strains.
BMC Syst Biol 2010, 4:53. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

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. Publisher Full Text

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):e1000744. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Pharkya P, Burgard AP, Maranas CD: OptStrain: A computational framework for redesign of microbial production systems.
Genome Res 2004, 14(11):23672376. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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. PubMed Abstract  Publisher Full Text

Segrè D, Vitkup D, Church GM: Analysis of optimality in natural and perturbed metabolic networks.
Proc Natl Acad Sci 2002, 99(23):1511215117. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Jungreuthmayer C, Zanghellini J: Designing optimal cell factories: Integer programing couples elementary mode analysis with regulation.
BMC Syst Biol 2012, 6:103. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Hädicke O, Klamt S: Computing complex metabolic intervention strategies using constrained minimal cut sets.
Metab Eng 2011, 13(2):204213.
http://www.ncbi.nlm.nih.gov/pubmed/21147248 webcite
PubMed Abstract  Publisher Full Text 
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.
http://dx.doi.org/10.1038/73786 webcite
Publisher Full Text 
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.
http://www.ncbi.nlm.nih.gov/pubmed/10087604 webcite
PubMed Abstract  Publisher Full Text 
Becker J, Zelder O, Häfner S, Schröder H, Wittmann C: From zero to hero–designbased systems metabolic engineering of Corynebacterium glutamicum for Llysine production.
Metab Eng 2011, 13(2):159168. PubMed Abstract  Publisher Full Text

Trinh CT, Li J, Blanch HW, Clark DS: Redesigning Escherichia coli metabolism for Anaerobic production of Isobutanol.
Appl Environ Microbiol 2011, 77(14):48944904. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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.
http://aem.asm.org/content/74/12/3634.abstract webcite
PubMed Abstract  Publisher Full Text  PubMed Central Full Text 
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. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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. Publisher Full Text

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. PubMed Abstract  Publisher Full Text

Zanghellini J, Ruckerbauer DE, Hanscho M, Jungreuthmayer C: Elementary flux modes in a nutshell: properties, calculation and applications.
Biotechnol J 2013, 8(9):10091016. PubMed Abstract  Publisher Full Text

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. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Berge C: Hypergraphs: Combinatorics of Finite Sets. (NorthHolland mathematical library: Volume 45.) Amsterdam: Elsevier Science Publishers; 1989.

Haus UU, Klamt S, Stephen T: Computing knockout strategies in metabolic networks.
J Comput Biol 2008, 15(3):259268. PubMed Abstract  Publisher Full Text

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. PubMed Abstract  Publisher Full Text

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. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Terzer M, Stelling J: Largescale computation of elementary flux modes with bit pattern trees.
Bioinformatics 2008, 24(19):22292235. PubMed Abstract  Publisher Full Text

Kamp Av, Schuster S: Metatool 5.0: fast and flexible elementary modes analysis.
Bioinformatics 2006, 22(15):19301931. PubMed Abstract  Publisher Full Text

Eiter T, Makino K, Gottlob G: Computational aspects of monotone dualization: A brief survey.
Discrete Appl Math 2008, 156(11):20352049. Publisher Full Text

Chen DS, Batson RG, Dang Y: Applied Integer Programming: Modeling and Solution. Hoboken: John Wiley & Sons, Inc.; 2010.

Orth JD, Fleming RMT, Palsson BØ: Reconstruction and use of microbial metabolic networks: the core escherichia coli metabolic model as an educational guide. In 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. Washington: ASM Press; 2009:5699.
chapter 10.2.1

Deutscher J: The mechanisms of carbon catabolite repression in bacteria.
Curr Opin Microbiol 2008, 11(2):8793. PubMed Abstract  Publisher Full Text

Klamt S, SaezRodriguez J, Gilles E: Structural and functional analysis of cellular networks with CellNetAnalyzer.
BMC Syst Biol 2007, 1:2. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Trinh CT: Elucidating and reprogramming Escherichia coli metabolisms for obligate anaerobic nbutanol and isobutanol production.
Appl Microbiol Biotechnol 2012, 95(4):10831094. PubMed Abstract  Publisher Full Text

Murakami K, Uno T: Efficient algorithms for dualizing largescale hypergraphs.
arXiv:1102.3813 2011.
http://arxiv.org/abs/1102.3813 webcite

Savelsbergh MWP: Preprocessing and probing techniques for mixed integer programming problems.
ORSA J Comput 1994, 6(4):445454. Publisher Full Text