Fast Computation of Minimal Cut Sets in Metabolic Networks with a Berge Algorithm That Utilizes Binary Bit Pattern Trees

Christian
Marie
Jürgen

Pages: pp. 1329-1333

Abstract—Minimal cut sets are a valuable tool for analyzing metabolic networks and for identifying optimal gene intervention strategies by eliminating unwanted metabolic functions and keeping desired functionality. Minimal cut sets rely on the concept of elementary flux modes, which are sets of indivisible metabolic pathways under steady-state condition. However, the computation of minimal cut sets is nontrivial, as even medium-sized metabolic networks with just 100 reactions easily have several hundred million elementary flux modes. We developed a minimal cut set tool that implements the well-known Berge algorithm and utilizes a novel approach to significantly reduce the program runtime by using binary bit pattern trees. By using the introduced tree approach, the size of metabolic models that can be analyzed and optimized by minimal cut sets is pushed to new and considerably higher limits.

Keywords—Elementary mode analysis; minimal cut sets; gene knockout; bit pattern; tree code

Introduction

Elementary flux mode (EFM) analysis is a well-established method to unbiasedly decompose (metabolic) networks into unique, indecomposable steady-state pathways [ 1], [ 2]. The indecomposability of EFMs implies that the deletion of any single reaction contributing to an EFM completely disables it. A cut set (CS) is defined as a set of reaction deletions that eliminates all target EFMs to be killed. A CS is called a minimal cut set (MCS) if none of its subsets is a CS [ 3]. MCSs are of special interest for metabolic engineering, as they require the least effort when biologically implemented. It has been shown that CSs are especially useful for two types of tasks:

1. The calculation of all MCSs of an unconstrained system, which means that every determined MCS kills all EFMs of the investigated systems. Numerous applications of this method have been suggested, for example, evaluating structural robustness and fragility, identifying targets in rational drug design, and predicting phenotypes [ 4].
2. EFMs can be utilized to split the functional units of a metabolic system into two groups: a) desired functions and b) unwanted functions [ 5]. A constrained cut set (cCS) is defined as a set of reaction deletions that eliminates all unwanted EFMs and keeps—at least some of—the desired EFMs. Such a strategy can be used to optimize microbiological production hosts to efficiently produce substances of interest, for example, ethanol [ 6]. This is achieved by assigning all EFMs that are involved in efficient ethanol production to the set of desired EFMs. All other EFMs are assigned to the set of unwanted EFMs. Consequently, the biological implementation of the computed MCSs of such a system results in an optimized microorganism. Therefore, MCSs are a valuable tool to identify and realize optimal gene intervention strategies, as has been shown by Trinh et al. [ 7].

The MCS computation of small metabolic networks is simple and computationally not demanding. However, the computation of MCSs suffers from a major disadvantage. The number of EFMs of metabolic networks grows combinatorially with the systems size [ 8]. Even medium-sized networks with approximately 100 reactions can easily have several hundred million EFMs. The huge number of EFMs makes the computation of MCSs a challenge and results in very long runtimes that do not allow the use of MCSs for medium-scale or large-scale metabolic models. Recently, several methods have been developed to compute efficiently MCSs using hitting set algorithms [ 4], [ 9], [ 3], [ 5], [ 10], [ 11] and binary linear programming [ 12]. The performance evaluation of MCS computation algorithms is nontrivial and currently no methods are known that are output-polynomial [ 10]. In this study, we introduce an improved version of the hitting set approach that was originally reported by Haus et al. [ 3]. The approach of Haus et al. significantly increased the performance of MCS programs by utilizing Berge's algorithm [ 13].

The Berge algorithm mainly processes three data sets:

1. the set of all EFMs to be killed,
2. the set of preminimal cut sets (preMCS), which kill all EFMs that have already been processed, and
3. the new preMCS candidates.

The Berge algorithm executes two nested loops (see Fig. 1). The outer loop iterates through all EFMs that have to be killed. The inner loop runs over all existing preMCSs. Inside the inner loop, it is tested if the current preMCS kills the EFM that is being analyzed. If the EFM is killed, then the preMCS is valid and the next preMCS is tested. However, if the EFM is not eliminated by the tested preMCS, then the preMCS is invalid. In this case, the preMCS is used to generate new preMCS candidates and, then, it is removed from the set of preMCSs. New candidates are simply created by taking the invalid preMCS and by adding a single element for each reaction that the processed EFM carries a flux. For example, if an insufficient preMCS contains two reactions ({R1, R7}) and the EFM that caused the preMCS to fail the test has four flux-carrying reactions ({R4, R6, R9, R11}), four new candidates are created with three reactions each ({R1, R7, R4}, {R1, R7, R6}, {R1, R7, R9}, {R1, R7, R11}). After the last EFM has been processed, all remaining preMCSs are MCSs that kill all EFMs that have to be deleted. However, the Berge algorithm suffers from a bottleneck, which is caused by the need to remove preMCS candidates that are supersets of already computed preMCS [ 3]. Consequently, a superset test is required to guarantee that all CSs are MCSs. Using a linear search approach to find subsets in the set of existing preMCS results in a search time of $O$(n^2) . Various other subset search approaches that scale better than the linear approach have been reported, for example, in [ 14] and [ 15] methods were reported that can identify subsets in time O$(n^2/\log n)$ . To reduce the runtime spent in the subset search of Berge's algorithm, we implemented a binary bit pattern tree approach [ 16] that can be run in parallel on multiple CPU cores. The concept of binary bit pattern trees is not new, but we are not aware that they have been used in the context of MCS computation. By using our tree approach, the size of metabolic models that can be analyzed and genetically optimized by MCSs is pushed to new and significantly higher limits. Besides the importance of MCSs in metabolic engineering, the computation of MCSs is also an essential problem in discrete mathematics with many applications in computer science, artificial intelligence, and game theory [ 17]. These fields may also benefit from the idea of using binary bit pattern trees to speed up the computation of MCSs. In this study, we explain the main concept of our tree method and compare the performance of the tree method with a linear search algorithm.

Figure    Fig. 1. Flowchart of the essential parts of the Berge algorithm.

Methods

Our MCS calculator utilizes the Berge algorithm [ 13] as demonstrated by Haus et al. [ 3]. In Berge's algorithm, the majority of the computation time is spent on the procedure that checks whether or not a new MCS candidate is a superset of an already determined preMCS [ 3] (see Table 2). If this is the case, then the new preMCS candidate is dismissed. Otherwise, the candidate is added to the list of preMCS.

A simple approach to perform this superset test is to sequentially access each preMCS of the complete set of existing preMCS and to compare it with the new preMCS candidate. As soon as a subset is found for a new candidate, the search procedure is stopped and the next candidate is tested.

However, this approach scales badly with an increasing number of preMCS. In particular, this is true if there do not exist any subsets for a candidate, in which case the complete set of existing preMCSs is examined. Note that the number of preMCSs can be huge and, in general, is much larger than the final number of computed MCSs (see Fig. 2, which shows the number of preMCSs as a function of processed EFMs for the presented benchmark system).

Figure    Fig. 2. Number of preMCS over number of processed EFMs for the system used throughout this paper.

Typically, CSs are expressed in form of bit patterns, where “1” stands for a reaction that is knocked out by the CS and “0” means that the corresponding reactions is not affected. A new candidate is a superset of an existing preMCS, if the result of a Boolean AND-operation of the candidate and the preMCS is equal to the preMCS itself (see Table 1). This implies that a candidate is not a superset if the bit pattern of a preMCS contains a “1” at a bit position, where the candidate is “0.” This fact can be used to speed up the superset check of the Berge algorithm by organizing the set of existing preMCSs in a favorable way as a binary tree and by traversing through this tree in an intelligent fashion. Thereby, a significant amount of superset tests can be avoided that would have to be performed in the case of a linear search.

Table 1. Implementation of the Superset Test Utilizing a Boolean AND-Operation

As shown in Fig. 3, our implementation of the Berge algorithm uses a binary tree to structure the preMCSs. A leaf node of the tree is always the bit pattern of a preMCS and is represented by rectangles with rounded corners. Whether a preMCS is attached left or right to its parent node is determined by the bit value of the preMCS at a certain bit position. In Fig. 3 a, preMCS is attached to the left side if the bit at the inspected bit position is set and to the right side if it is not set. For example, the bit position of the root node is 3. Consequently, all preMCSs with a “1” at position 3 are attached to the subtree at the left side of the root node, whereas all preMCSs on the right side of the root node are “0” at bit position 3. The bit position that is used at each level of the tree (e.g., in Fig. 3 at level 0 bit number 3, at level 1 bit number 4, at level 2 bit number 2, and so on) is determined during a preprocessing step before the Berge algorithm is started and has a significant influence on the runtime of the program. The bit pattern value of a nonleaf node is created by Boolean AND-operations of the bit patterns of the node's subtrees. This AND-operation has the effect that “1”s that are set in all subnodes are propagated from the bottom of the tree (leaf nodes) to the top (root node). If a candidate is checked, the superset test is started on the root node. The left subtree of a node is always investigated first. If no subset could be found in the left subtree, the right subtree is explored. However, the left side of a subtree is only entered if the tested candidate does not contain a “0,” where all subnodes and, hence, all preMCS have a “1.” Such a case is shown in Fig. 3, where the candidate “100101” is tested (blue-colored dash-dotted line). As the bit pattern of the root nodes (“000000”) does not contain any “1”s, the left subtree of the root node is entered. The bit pattern of the left child node (“001000”) indicates that in all nodes and, hence, in all preMCSs attached to it the bit at position 3 is set. As bit number 3 of the tested candidate “100101” is not set, none of the attached preMCSs can be a subset of the candidate and this subtree is not further investigated. Therefore, the right subtree of the root node is entered next. In this subtree, all leaf nodes are accessed to find a subset. As no subset is found, the candidate is accepted and can be added to the set of preMCSs. A different situation is illustrated by the candidate “011101” (red-colored dashed line). Since the bit at position 3 is set, the left subtree of the nonleaf node “001000” must also be checked. However, this check does not return a valid subset. Hence, the right subtree is investigated. In this subtree, a subset is detected (“011001”). Consequently, the check procedure is immediately terminated and the candidate is dismissed.

Figure    Fig. 3. Binary bit pattern tree of all existing preMCS. Leaf nodes containing bit patterns of preMCS are represented by rectangles with rounded corners. The bit patterns of the nonleaf nodes are built by Boolean AND-operations of the node's subtrees.

Before the Berge algorithm is started, extensive preprocessing of the set of provided EFMs is performed. The preprocessing comprises the following steps:

1. eliminating reactions that must not be knocked out, as this would result in the deletion of (too many) wanted modes,
2. removing duplicate modes,
3. removing modes that are supersets of other modes, and
4. combining duplicate reactions.

The preprocessing step 4 is illustrated in Fig. 1 as the step “compress EFMs,” which is the reverse operation of the step “expand MCS.” For details, we refer to [ 18].

To benchmark the performance of our bit pattern tree algorithm, we used a system with 114,614 EFMs. The EFMs were calculated with the open source program regEfmtool [ 19], [ 20]. The used network has been described in [ 21]. It models the central carbon metabolism of generic plant cells and consists of 78 reactions (33 reversible and 45 irreversible) and 55 internal metabolites. The network describes the following main pathways: TCA cycle, glycolysis, pentose phosphate, respiration, sucrose, and starch synthesis. All of them belong to one of these four compartments: vacuole, mitochondria, plastid, and cytosol. Transporters have been added for intercompartmental exchange. For the benchmark, we computed MCSs of the unconstrained system, which means that all EFMs of the system have to be killed. An unconstrained system has the effect that the preprocessing steps 1 to 3 mentioned earlier cannot be applied and, hence, the size of the system is only marginally reduced during preprocessing. Consequently, using an unconstrained system results in a high computational workload. The total number of MCSs of the chosen system is 2,815,375, where the minimum and maximum number of knockouts is 4 and 18, respectively (see Fig. 5). However, because of preprocessing step 4, which combines duplicate reactions, the Berge algorithm only needs to compute 93,009 MCSs. In a postprocessing step these 93,009 MCSs are expanded to the final number of 2,815,375. The concept of CS expanding is illustrated in Fig. 4, which shows two networks. Network (a) depicts the original toy system that has two EFMs and seven MCSs that remove both EFMs: $\{R1\}$ , $\{R4\}$ , $\{R5\}$ , $\{R2,R6\}$ , $\{R3,R6\}$ , $\{R2,R7\}$ , and $\{R3,R7\}$ . Network (b) is the toy system after combining duplicate reactions, which only has three MCS: $\{R1\}$ , $\{R4\}$ , and $\{R23,R67\}$ . However, in a postprocessing step the three MCSs can be expanded to obtain the full set of MCSs of the original network. Using the compressed network instead of the full network results in a strong reduction of the program's total execution time [ 18].

Figure    Fig. 4. Network (a) is the original system and network (b) the system after combining duplicate reactions (R2 and R3 R23; R4 and R5 R45; R6 and R7 R67). Both networks have two EFMs. However, network (a) has seven MCSs, whereas network (b) has only three MCSs. In a postprocessing step, the three MCSs of the compressed network can be expanded to the seven MCSs of the original network. Using the compressed network to do the MCS computation results in strong reduction of the total runtime.

Figure    Fig. 5. Number of MCSs as a function of number of knockouts for the metabolic system used in the presented benchmark.

Results

The results of the benchmark are listed in Table 2. The table clearly shows that the tree approach is superior to the linear search algorithm, for example, in single-thread mode the tree algorithm is approximately 30 times faster than the linear search algorithm. The performance gain is even higher if the MCSs of even bigger systems are computed (see Table 3). The benchmark also shows that most of the time is spent on the superset test (at least 84 percent), which is consistent with the observations of others [ 3]. The time lost managing the bit pattern tree accumulated to approximately 40 seconds. A part of these 40 seconds is spent in the candidate generation procedure, where invalid preMCSs are eliminated and, hence, must be removed from the tree. The other task of the tree management is to add valid candidates to the preMCS tree.

Table 2. Comparison of MCS Computation Runs with (a) Linear Search, (b) Tree Search with Random Bitorder, and (c) Tree Search with a Bit Order Derived from EFM Properties for a System with 114,614 EFMs

As shown in Table 2, the bit position used to create the tree (see Fig. 3) has a strong influence on the total runtime (compare test case (b) and test case (c)). In our benchmark, the best performance was achieved when the order of the bit position was determined by the frequency of active reactions in the set of EFMs to be killed. The more often a certain reaction was carrying a flux in the EFMs, the higher up in the tree the reaction was used to split the set of preMCSs. For example, given in Fig. 3, this means that the reaction represented by bit position 2 occurs most frequently in the EFMs to be killed. The results obtained by using this EFM statistics approach to determine the bit order are shown in test case (c)  tree—EFM. Test case (b) tree—random used a random bit order. For test case (b), five runs with varying bit orders were used and the results averaged. To obtain comparable results, the same random bit order was used for each set of runs with varying numbers of threads. Using the bit order derived from the reaction frequency resulted in a performance gain of approximately a factor 3 compared to the random bit order.

The achieved performance improvement can be explained by the number of comparison events that have to be performed (see Table 2). In the linear case $20.4\cdot 10^{12}$ subset tests were done in total, whereas only $30.9\cdot 10^{9}$ were required for the tree approach shown in test case (c). Although, in test case (a) 660 times more comparison events occurred than in test case (c), the runtime improvement was much lower than the factor 660. This is mainly caused by the additional complexity of the tree code and by unavoidable effects such as an increased number of cache mismatches, which occur more frequently if memory is accessed randomly instead of linearly.

Note that Table 2 also illustrates that the linear search scales much better with an increasing number of threads than the tree approach, as the total runtime for the linear search was reduced by a factor of 7.1 when 10 threads were used, whereas the tree implementation only gained a factor of 1.7. Analyzing and improving the multithreaded tree code implementation will be scope for future work.

Results of MCS computations for five metabolic networks are listed in Table 3. The table clearly shows that the performance gain grows with the number of preMCSs that have to be determined during the computation. The number of EFMs is not necessarily indicative of the number of MCSs. Our P. pastoris model has 15 million EFMs and roughly 1.8 million MCS, while our A. thaliana contains only 1.7 million EFMs but has 3.7 billion MCSs. Note that we terminated the linear search runs of the two large systems, as the execution time was extraordinarily high. Based on the average time spent to search for a single candidate in the set of preMCSs, an extrapolated runtime was estimated. Table 3 shows that the tree approach outperforms the linear search algorithm. Moreover, the tree approach allows to study metabolic systems that could not be analyzed with a linear search.

Table 3. Comparison of MCS Computation with Linear Search and with Tree Search for Different Metabolic Networks Using 10 Parallel Threads

Conclusion

We presented a program that computes MCSs by utilizing a hitting set algorithm originally invented by Berge [ 3]. The bottleneck of the Berge algorithm is the search procedure that is required to guarantee that all computed MCSs are truly minimal. This is done by performing a test that verifies that a new MCS candidate is not a superset of any of the already determined preMCS [ 5]. A simple but slow approach to test a new candidate against all already existing preMCSs is to linearly iterate through the complete set of existing preMCSs and check if this set contains a subset of the candidate. To speed up the superset test, we stored the existing preMCSs in a binary tree structure. This tree structure can be used to skip a significant amount of the necessary comparison events between a candidate and the existing preMCSs. Consequently, the total runtime for the computation of the complete set of MCSs is strongly reduced, and hence, numerical analyses and gene intervention strategies of much bigger systems can be calculated by our approach. Note that the determination of the MCSs requires the calculation of the EFMs first. Even though significant progress has been made in recent years on the efficient computation of EFMs [ 16], [ 22], [ 23], in many situations the EFM computation is the limiting factor of MCS analyses. However, recently alternative approaches have been reported which are based on the duality of EFMs and MCSs [ 3], [ 24]. These approaches allow the computation of MCSs directly from the stoichiometric matrix without having to calculate the EFMs beforehand.

The Berge algorithm is well suited to compute MCSs of systems with a moderate number of columns/reactions (up to several hundred) and a high number of rows/modes (several million). However, systems with other characteristics (e.g., a high number of columns and a low number of rows) cannot be investigated easily with the Berge algorithm—even if a tree approach is used—as the number of preMCSs grows dramatically with each processed row/mode. The huge number of preMCSs results in a tremendous memory consumption and a very slow execution that make the Berge algorithm unsuitable for these type of systems.

Besides the importance of MCSs in metabolic engineering, the calculation of MCSs is also an essential problem in discrete mathematics with numerous applications in computer science, artificial intelligence, and game theory [ 17]. These fields might also benefit from the idea of using binary bit pattern trees to speed up the computation of MCSs.

Acknowledgments

This work was supported 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 COMET-Funding Program managed by the Austrian Research Promotion Agency FFG.

References

• 1. S. Schuster, D.A. Fell, and T. Dandekar, “A General Definition of Metabolic Pathways Useful for Systematic Organization and Analysis of Complex Metabolic Networks,” Nature Biotechnology, vol. 18, no. 3, pp. 326-332, Mar. 2000.
• 2. S. Schuster, T. Dandekar, and D.A. Fell, “Detection of Elementary Flux Modes in Biochemical Networks: A Promising Tool for Pathway Analysis and Metabolic Engineering,” Trends in Biotechnology, vol. 17, no. 2, pp. 53-60, Feb. 1999.
• 3. U.-U. Haus, S. Klamt, and T. Stephen, “Computing Knock-Out Strategies in Metabolic Networks,” J. Computational Biology, vol. 15, no. 3, pp. 259-268, Apr. 2008.
• 4. S. Klamt, and E.D. Gilles, “Minimal Cut Sets in Biochemical Reaction Networks,” Bioinformatics, vol. 20, no. 2, pp. 226 -234, Jan. 2004.
• 5. O. Hädicke, and S. Klamt, “Computing Complex Metabolic Intervention Strategies Using Constrained Minimal Cut Sets,” Metabolic Eng., vol. 13, no. 2, pp. 204-213, Mar. 2011.
• 6. C.T. Trinh, P. Unrean, and F. Srienc, “Minimal Escherichia coli Cell for the Most Efficient Production of Ethanol from Hexoses and Pentoses,” Applied and Environmental Microbiology, vol. 74, no. 12, pp. 3634-3643, http://aem.asm.org/content/74/12/3634.abstract, June 2008.
• 7. C.T. Trinh, R. Carlson, A. Wlaschin, and F. Srienc, “Design, Construction and Performance of the Most Efficient Biomass Producing E. coli Bacterium,” Metabolic Eng., vol. 8, no. 6, pp. 628-638, Nov. 2006.
• 8. S. Klamt, and J. Stelling, “Combinatorial Complexity of Pathway Analysis in Metabolic Networks,” Molecular Biology Reports, vol. 29, no. 1, pp. 233-236, Mar. 2002.
• 9. S. Klamt, “Generalized Concept of Minimal Cut Sets in Biochemical Networks,” Biosystems, vol. 83, nos. 2/3, pp. 233-247, Feb. 2006.
• 10. M. Hagen, “Lower Bounds for Three Algorithms for Transversal Hypergraph Generation,” Discrete Applied Math., vol. 157, pp. 1460-1468, 2009.
• 11. K. Murakami, and T. Uno, “Efficient Algorithms for Dualizing Large-Scale Hypergraphs,” Proc. Meeting on Algorithm Eng. and Experiments (ALENEX '13), pp. 1-13, 2013.
• 12. C. Jungreuthmayer, and J. Zanghellini, “Designing Optimal Cell Factories: Integer Programing Couples Elementary Mode Analysis with Regulation,” BMC Systems Biology, vol. 6, article 103, 2012.
• 13. C. Berge, Hypergraphs: Combinatorics of Finite Sets, vol. 45, first ed. North Holland, Aug. 1989.
• 14. H. Sheni, and D.J. Evans, “Fast Sequential and Parallel Algorithms for Finding Extremal Sets,” Int'l J. Computer Math., vol. 61, pp. 195-211, 1996.
• 15. P. Pritchard, “On Computing the Subset Graph of a Cellection of Sets,” J. Algorithms, vol. 33, pp. 187-203, 1999.
• 16. M. Terzer, and J. Stelling, “Large-Scale Computation of Elementary Flux Modes with Bit Pattern Trees,” Bioinformatics, vol. 24, no. 19, pp. 2229-2235, Oct. 2008.
• 17. T. Eiter, K. Makino, and G. Gottlob, “Computational Aspects of Monotone Dualization: A Brief Survey,” Discrete Applied Math., vol. 156, no. 11, pp. 2035-2049, June 2008.
• 18. C. Jungreuthmayer, G. Nair, S. Klamt, and J. Zanghellini, “Comparison and Improvement of Algorithms for Computing Minimal Cut Sets,” to be published in BMC Bioinformatics, 2013.
• 19. C. Jungreuthmayer, D. Ruckerbauer, and J. Zanghellini, “regEfmtool: Speeding up Elementary Flux Mode Calculation Using Transcriptional Regulatory Rules in the Form of Three-State Logic,” BioSystems, vol. 113, pp. 37-39, July 2013.
• 20. C. Jungreuthmayer, and J. Zanghellini, “Regulatory Elementary Flux Mode Tool, regEfmtool,” http://www.biotec.boku.ac.at/ regulatoryelementaryfluxmode.html, 2012.
• 21. M. Beurton-Aimar, B. Beauvoit, A. Monier, F. Vallée, M. Dieuaide-Noubhani, and S. Colombié, “Comparison between Elementary Flux Modes Analysis and 13C-Metabolic Fluxes Measured in Bacterial and Plant Cells,” BMC Systems Biology, vol. 5, article 95, 2011.
• 22. D. Jevremovi , C.T. Trinh, F. Srienc, C.P. Sosa, and D. Boley, “Parallelization of Nullspace Algorithm for the Computation of Metabolic Pathways,” Parallel Computing, vol. 37, nos. 6/7, pp. 261-278, 2011.
• 23. C. Jungreuthmayer, D. Ruckerbauer, and J. Zanghellini, “Utilizing Gene Regulatory Information to Speed Up the Calculation of Elementary Flux Modes,” arXiv, http://arxiv.org/abs/1208.1853, p. 1208.1853v1, 2012.
• 24. K. Ballerstein, A. von Kamp, S. Klamt, and U.-U. Haus, “Minimal Cut Sets in a Metabolic Network Are Elementary Modes in a Dual Network,” Bioinformatics, vol. 28, no. 3, pp. 381-387, Dec. 2011.