Skip to main content

Discovery of Boolean metabolic networks: integer linear programming based approach



Traditional drug discovery methods focused on the efficacy of drugs rather than their toxicity. However, toxicity and/or lack of efficacy are produced when unintended targets are affected in metabolic networks. Thus, identification of biological targets which can be manipulated to produce the desired effect with minimum side-effects has become an important and challenging topic. Efficient computational methods are required to identify the drug targets while incurring minimal side-effects.


In this paper, we propose a graph-based computational damage model that summarizes the impact of enzymes on compounds in metabolic networks. An efficient method based on Integer Linear Programming formalism is then developed to identify the optimal enzyme-combination so as to minimize the side-effects. The identified target enzymes for known successful drugs are then verified by comparing the results with those in the existing literature.


Side-effects reduction plays a crucial role in the study of drug development. A graph-based computational damage model is proposed and the theoretical analysis states the captured problem is NP-completeness. The proposed approaches can therefore contribute to the discovery of drug targets. Our developed software is available at “”.


One of the most important biological processes in organism is metabolism and its dysregulation can contribute to many human diseases. The manipulation of metabolism has been extensively studied in the field of metabolic engineering as many of metabolic process produce commodity and specialty chemicals. In particular, production of industrially valuable products using a microbial host with recombinant technology becomes one of the most successful applications of metabolic engineering [13]. Metabolic network is used to model the behavior of reactions and chemicals [4]. Authors in [5] have proposed a new mathematical model for large metabolic network and demonstrated that the outcome of network expansion is robust against of single or few reactions. Indeed, for the analysis of metabolic networks, many studies have been conducted so as to develop Boolean models. For example, Smart et al. [6] studied the effect of deletion of each enzyme in metabolic network of a Boolean model, and Lemke et al. [7] considered almost the same problem from the viewpoint of the Boolean aspect. Tamura et al. [8] developed an integer linear programming method for Boolean reaction cut (BRC) problem. Furthermore, in [9], the computational complexity of BRC was analyzed. In the Boolean model of metabolic networks, compounds and reactions can be represented by “OR” and “AND” nodes, respectively. Li et al. [10] have developed methods for finding a set of enzymes whose inhibition stops the production of the target compound with a minimum elimination of non-target compounds. Furthermore, Takemoto et al. [11] and Lee et al. [12] estimated the size distribution of the deletion effect of enzymes using a branching process. Enumerating all the possible enzyme combination is infeasible since the number of combination increases exponentially with the number of enzymes. Lu et al. [13] developed an integer linear programming method for designing synthetic metabolic networks by considering minimum reaction insertion in a Boolean model. The proposed method can appropriately solve minimum reaction insertion in a Boolean network.

In pharmaceutics, the development of drugs mainly focus on target identification and lead inhibitor identification [14]. Furthermore, many researchers have done a lot of work on the efficacy of drugs regardless of their side-effects (toxicity). However, toxicity or lack of efficacy may occur when the compounds which are not the intended targets in the metabolic network. Thus how to effectively reduce the side-effect of drugs has become an important and challenging issue. The current works aim to identify the biological targets (could be enzyme or protein) for drugs which can be manipulated so as to produce the desired effect (i.e., curing a disease) with minimum disruptive side-effects [15, 16]. Drug targets, in particular enzymes, are chosen to reduce abnormal metabolites by formulating an optimal combination problem in enzyme combination in metabolic networks [10, 17]. It is known that the function of enzymes is to catalyze reactions and produces metabolites (compounds) in metabolic networks. Human diseases, especially metabolic diseases, may be directly caused by the accumulation of certain compounds due to the malfunctions of enzymes. We denote such compounds as target compounds while the others as non-target compounds. For example, the malfunction of enzyme phenylalanine hydroxylase results in accumulation of the amino acid phenylalanine which causes the disease phenylketonuria [18]. Therefore, there is a great need to identify a set of enzymes which are manipulated by drugs to prevent the buildup of target compounds with minimal damage. Here damage corresponds to the number of non-target compounds whose production are forced to stop by the inhibition of that enzyme (or enzyme set).

Then the problem becomes to identify the optimal enzyme set whose inhibition eliminates the target compounds and at the same time, incurs minimum damage based on the given metabolic network and a set of target compounds. We denote such problem as enzyme combination identification (ECI). Sridhar et al. [19] developed a branch-and-bound algorithm to dynamically explore the search space. Furthermore, two filtering strategies are proposed to prune the search space to guarantee an optimal solution. However, their algorithm is complicated and impractical for us to use it. In this paper, we develop an efficient method based on Integer Linear Programming [2022] which has a wide application in solving NP-hard problems. Furthermore, all instances of the original problem are needed to convert into integer programming formalization so as to apply the existing free solver called CPLEX [23].

This paper has two main contributions. First, we formulate a new biological problem based on Boolean metabolic network to identify the drug target with minimum damage. Second, we propose an effective and efficient method, Integer Linear Programming, (needs O(m+n) variables) to solve the captured problem and it is easy to implement. By integrating the human metabolic networks, we have shown that the proposed approach can accurately identify the target enzyme set for known drugs in the literature. Furthermore, experiments have shown that our proposed method is extremely efficient which can effectively solve the problem in seconds.

The remainder of the paper is organized as follows. In “Problem formulation” section, we formulate the captured problem. “Methods” section presents the algorithm ILP-ECI. In “Results and discussions” section, we discuss the experimental results and give a theoretical illustration. Finally, conclusions are given in the last section.

Problem formulation

A metabolic network can be represented by a directed graph G=(V,A) which captures the interaction between reactions, compound, and enzymes. Let C, R, and E denote the set of compounds, reactions, and enzymes, respectively. Then V C ,V R and V E denote the set of vertices from C, R, and E, where \(\phantom {\dot {i}\!}V_{C}=\{v_{c_{1}},v_{c_{2}},\cdots,v_{c_{m}}\}\) denotes a set of compound nodes, \(\phantom {\dot {i}\!}V_{R}=\{v_{r_{1}},v_{r_{2}},\cdots,v_{r_{n}}\}\) represents a set of reaction nodes and \(\phantom {\dot {i}\!}V_{E}=\{v_{e_{1}},v_{e_{2}},\cdots,v_{e_{l}}\}\) defines a set of enzyme nodes. The subscript m,n,l denote the number of compound nodes, reaction nodes and enzyme nodes, respectively. Here V=V C V R V E and V C V R ={},V C V E ={},V R V E ={} hold. Here {} is the empty set. Let V s V C and V t V C be sets of source nodes and target nodes, respectively. Source nodes are those nodes having no incoming edges and they are the seed compounds. A metabolic network can be defined as follows:

  • For each edge (u,v)A, either (uV C )(vV R ) or (uV R )(vV C |vV E ) holds and “ ” represents “AND” function;

  • Each source node has no incoming edges;

  • For nodes vV s and vV C have at least one incoming edge.

The main problem Enzyme Combination Identification (ECI) in a Boolean network is first described with a simple example and then followed by its mathematical formalization.

A small hypothetical metabolic network is shown in Fig. 1, where circles, rectangles and triangles represent compounds, reactions and enzymes, respectively. For instance, R1 is a reaction, its substrates are C1 and E2 and its product is C2. Here C5 is the target compound in this figure which shall be stopped. To stop the production of C5, reaction R2 must be prevented from taking place. One of the possible solution is to disrupt one of its catalyzing enzymes (E1 for instance). Another is to stop the production of its reactant compounds (i.e., C2 or C3). If C2 is stopped, then two possible ways can achieve this effect (i.e., enzyme E2 or E1). It is noted that C2,R2,C7 and R1 forms a cycle, the inhibition of E1 will result in inactivating reaction R2 which will further makes C2 inproducible. The other way is inhibiting enzyme E2 which makes R1 inactive and further stop the production of compound C2. Therefore, the inhibition of enzyme E1 or E2 can result in stopping the production of the target compound.

Fig. 1
figure 1

A graph constructed for a hypothetical metabolic network

We define the Enzyme Combination Identification as follows.

  • Input: A metabolic network and a set of target compounds T(TC), C is the set of compounds.

  • Output: Find a set of enzymes X(XE) with minimum damage, whose inhibition stops the production of all the compounds in T.

For simplicity, we assume there are no external inputs to all reactions and all the input compounds related to reactions are shown in the network. We note that different compounds and enzymes may have different levels of importance in the metabolic networks. Here we assume that all the compounds and enzymes are of equal importance. We assign binary value (i.e., 0 or 1) to each node V in the Boolean model. “0” means that the corresponding compound is not producible or the corresponding reaction is inactive while “1” means that the corresponding compound is producible and the corresponding reaction is active. For the related enzyme catalyzing the reaction, “1” means that it is active and “0” means that it is inhibited. Let G be such an assignment (G is a function from V to {0,1}). For each node vV, we write v=1 (resp., v=0) if 1 (resp., 0) is assigned to v. Then G can be regarded as a valid assignment if the following conditions are satisfied:

  • (i) For each vV s ,v=1.

  • (ii) For each vV s ,v=1 if and only if there is u such that (u,v)E and u=1.

  • (iii) For each vV r ,v=1 if and only if u=1 holds for all u such that (u,v)E.

Condition (ii) implies that compound nodes correspond to “OR” nodes. Condition (iii) implies that reaction nodes correspond to “AND” nodes, which means that the output is forced to 0 when a node is inactivated.

From Fig. 1, we can see that inhibition of E1 results in the knock out of compounds C4,C7,C8 and C9 in addition to the target compound C5. Then we denote the number of non-target compounds knocked out as the damage, which is caused by the manipulating the enzyme set in the metabolic network. It can be seen that the damage of inhibiting E2 is 2 (i.e., C2 and C4). Compound C7 is still producible because it can be produced by R3 even after the inhibition of E2. The damage effect of inhibition of E1 is 4 (i.e., C4,C7,C8 and C9). Both E1 and E2 are potential drug targets since they can achieve the effect of disrupting the target compound C5. However, E2 is a better drug target than E1 owing to the fact that it causes less damage.


In this section, we introduce integer programming-based methods for ECI. Integer programming, in particular, Integer Linear Programming (ILP) is set to minimize (or maximize) a linear objective function under linear constraints with all the variables taking integer values. In the following, each variable takes including the binary value (i.e., 0 or 1), representing the Boolean values. We apply ILP to ECI since ILP is widely used for solving NP-hard problems.

The ILP formulation for the network in Fig. 1 is as follows: ILP-ECI

$$\begin{array}{*{20}l} {\min} \left\{\sum_{i=1}^{m}TCi\right\} \end{array} $$

subject to

$$\begin{array}{*{20}l} &FC5=1 \end{array} $$
$$\begin{array}{*{20}l} &TR1+FC1+FC7+FE2\geq 1, \\ &FR1+TC1\geq 1, \\ &FR1+TC7\geq 1, \\ &FR1+TE2\geq 1, \end{array} $$
$$\begin{array}{*{20}l} &TR2+FC2+FC3+FE1\geq 1, \\ &FR2+TC2\geq 1,\\ &FR2+TC3\geq 1,\\ &FR2+TE1\geq 1, \end{array} $$
$$\begin{array}{*{20}l} &TR3+FC6+FE1\geq 1, \\ &FR3+TC6\geq 1,\\ &FR3+TE1\geq 1 \end{array} $$
$$\begin{array}{*{20}l} &TC2=TR1 \end{array} $$
$$\begin{array}{*{20}l} &TC4=TR2 \end{array} $$
$$\begin{array}{*{20}l} &TC5=TR2 \end{array} $$
$$\begin{array}{*{20}l} &FC7+TR2+TR3\geq 1,\\ &TC7+FR2\geq 1,\\ &TC7+FR3\geq 1 \end{array} $$
$$\begin{array}{*{20}l} &TC8=TR3 \end{array} $$
$$\begin{array}{*{20}l} &TC9=TR3 \end{array} $$
$$\begin{array}{*{20}l} &TC1=1, \qquad TC3=1, \qquad TC6=1 \end{array} $$
$$\begin{array}{*{20}l} &TX+FX=1 \text{for any X} \end{array} $$

We denote the above formalization as ILP-ECI. Here all variables including the value of reaction compound and enzyme nodes take either 1 or 0. Thus, \(v_{r_{i}}\) can be either 0 or 1, and \(v_{c_{i}}\) and \(v_{e_{i}}\) also take 0 or 1. In this example, \(v_{r_{i}}=0\) (resp. \(v_{r_{i}}=1\)) indicating that the value of reaction i takes 0 is represented by FRi=1 (resp. TRi=1) which implies that the reaction is inactivated (otherwise, the TRi=1 implies the reaction is activated). Therefore, TRi=0 (equivalent to FRi=1) means the corresponding value for true reaction takes 0, which implies the reaction is inactive. And FRi=0 (equivalent to TRi=1) indicates the corresponding value for false reaction takes 0, which implies the reaction is active. Thus, TRi+FRi=1 holds for any node i in the network. Similarly, TCi and FCi are used to represent the values of compound nodes. For instance, TC2=1 means that \(v_{c_{2}}=1\) and in other words, FC2=0 since TCi+FCi=1. Furthermore, \(v_{r_{i}}\) corresponds to “AND” node which implies that if \(v_{e_{i}}=0\) will inactivate \(v_{r_{i}}\).

The objective function (1) means that the damage should be minimized. FCi=1 (or TCi=0) means that a compound \(v_{c_{i}}\) is not producible. Equation (2) means that the target compound \(v_{c_{5}}\) should be 0 after the 0-1 assignment converges. Equation (3) represents the Boolean relation \(v_{r_{1}}=v_{c_{1}}\land v_{c_{7}}\land v_{e_{2}}\). Note that the Boolean relations such as “ ” or “ ” cannot be used in ILP formulation, we need to convert them into linear equations and/or inequations. Actually, “ ” indicates “AND” function and “ ”represents “OR” function. Since x1=x2x3x n can be represented by

$${}\left(x_{1}\!\vee\! \overline{x_{2}} \vee\! \cdots \!\vee \overline {x_{n}}\right)\land\left(\overline{x_{1}} \vee x_{2}\right)\land \left(\overline{x_{1}} \vee x_{3}\right)\land\! \cdots \land \left(\overline{x_{1}} \vee\! x_{n}\right)\,=\,1, $$

the constraint \(v_{r_{1}}=v_{c_{1}}\land v_{c_{7}}\land v_{e_{2}}\) can be converted into

$${}\left(v_{r_{1}}\!\vee \!\overline{v_{c_{1}}}\!\vee\! \overline{v_{c_{7}}}\vee \overline{v_{e_{2}}}\right)\land \left(\overline{v_{r_{1}}}\vee \!v_{c_{1}}\right)\land \left(\overline{v_{r_{1}}}\vee\! v_{c_{7}}\right)\land\left(\overline{v_{r_{1}}}\!\vee \!v_{e_{2}}\right)\,=\,1. $$

Thus Eq. (3) is obtained. Similarly, Eqs. (4)-(5) represent the constraints of \(v_{r_{2}}\) and \(v_{r_{3}}\), respectively.

For a compound node with indegree is 1 which indicates the node has only one incoming edge, the value of the predecessor is just copied. For instance, since \(v_{c_{2}}\) has only one predecessor \(v_{r_{1}}\), \(v_{c_{2}}\) is just copied from \(v_{r_{1}}\) as shown in Eq. (6). Similarly, \(v_{c_{4}}\) is just copied from \(v_{r_{2}}\) which is shown in Eq. (7).

However, for a compound node with indegree more than 1, it is necessary to convert the “ ” relation into linear equation or inequations. Equation (9) represents the Boolean relation \(v_{c_{7}}=v_{r_{2}}\vee v_{r_{3}}\). Since x1=x2x3x n can be represented by

$${}\left(\overline{x_{1}}\!\vee\! {x_{2}} \!\vee\! \cdots \!\vee\! { x_{n}}\right)\land\left({x_{1}} \vee \overline{x_{2}}\right)\land \left({x_{1}} \vee \overline{x_{3}}\right)\land \cdots \land \left(x_{1} \!\vee \overline{x_{n}}\right)\,=\,1, $$

\(v_{c_{7}}=v_{r_{2}}\vee v_{r_{3}}\) can be turned into \(\left (\overline {v_{c_{7}}}\vee v_{r_{2}}\vee v_{r_{3}}\right)\land (v_{c_{7}}\vee \overline {v_{r_{2}}})\land \left (v_{c_{7}}\vee \overline {v_{r_{3}}}\right)\). Thus, Eq. (9) can be obtained.

Equation (6) - Eq. (12) represent the constraints of \(v_{c_{1}},\cdots,v_{c_{9}}\) respectively. Equation (12) means that \(v_{c_{1}}, v_{c_{3}}\) and \(v_{c_{6}}\) are 1 since their indegrees is 0.

Equation (13) means that “T” and “F” correspond to “true (1)” and “false (0)”, respectively, and complement each other. X in Eq. (13) means any component or reaction in the metabolic network.

The above formalization can clearly solve ECI and obtain the correct solution {E2}. Besides, the number of variables is O(m+n) in the above formalization where m and n are the number of compounds and reactions, respectively.

It is noted that solving ILP is NP-complete, however, a problem that can be formalized as ILP is not always NP-complete. Thus in the following, we prove that ECI is NP-complete.


ECI is NP-complete problem with the maximum indegree and outdegree being bounded by 2.


Obviously, the problem is in NP, it suffices to show that it is NP-hardness. The proof is by a polynomial time reduction from minimum edge cover (MEC), which is a problem for a given graph to find the minimum number of edges so that each node is incident to at least one of the selected edges. For instance, E1={e2,e3,e6} is one of optimal solutions of MEC for graph shown in Fig. 2. Let G=(V,E) be an instance of MEC, where V={v1,v2,,v n } and E={e1,e2,,e m }. We then construct the corresponding ECI as below. The metabolite network G=(V c V r V e ,E) is given by

$$\begin{aligned} V_{c}=\left\{c_{1},c_{2},\cdots,c_{m}\right\}\cup{c_{t}}, \quad V_{r}=\left\{r_{1},r_{2},\cdots,r_{n}\right\}, \end{aligned} $$
$$\begin{aligned} &{}E\,=\,\left\{\!\left\{\!c_{i},r_{j}\right\}\!|i\,=\,1,\cdots, m, j\,=\,\!1,\cdots, n, \text{if} \ v_{j} \ \text{is an end point of}\right.\\ &{}\left. \ e_{i}.\right\} \cup \left\{\left\{r_{j},c_{t}\right\}|j=1,\cdots,n\right\} \end{aligned} $$
Fig. 2
figure 2

An instance of minimum edge cover (MEC)

It is noted that the minimum damage is determined uniquely by the inhibition of enzyme set. Furthermore, our objective is to minimize the “damage” (side-effects). Then V e can be regarded as virtual nodes and denoted as an empty set in this case. The ECI problem can be converted into the problem of identifying the minimum set of non-target compounds. Thus the graph for MEC shown in Fig. 2 is converted into ECI shown in Fig. 3. It is clear that this conversion can be done in polynomial time. Then we show that MEC for G has a solution of size z if and only if ECI has a minimum damage of size z. To guarantee that the target compound c t is stopped (i.e., c t =0), it implies that all r j (j=1,,n) takes the value 0. If G has an edge cover of size z, then it follows that the minimum number of c i taking 0 should be z. On the other hand, if the minimum damage for ECI is z, then each r j must be 0 so as to satisfy c t =0, we have at least predecessor of each r j must be included in the minimum damage set. Since there is an edge between c i and r j if and only if v j is incident to e i . Thus {c i |c i minimum damage set} is an edge cover of size z. □

Fig. 3
figure 3

The polynomial time reduction from minimum edge cover (MEC) problem to minimum damage identification problem

Results and discussions

In this section, we verify the biological validity of our proposed method by using known drugs. Besides, the performance of the ILP-ECI algorithm is evaluated by using the execution time which indicates the total time taken by the method. For the experimental data, we extracted the metabolic network information from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database [24]. KEGG is database which provides known drug molecules along the the enzymes they inhibit and their therapeutic category. Then we use drugs at this database as our benchmarks and we report two of them due to the space limitation. And the value in parenthesis that starts with letter “C”, “D”, and “E” (e.g., E1.13.11.34) is the unique identifier assigned to the corresponding compound, drug, and enzyme respectively in KEGG.

The drug we used in this paper is Benoxaprofen (D03080) which inhibits arachidonate 5-lipoxygenase (E1.13.11.34) [25]. This enzyme appears in several networks including arachidonic acid metabolism network (hsa00590). In Pharmacology, the inhibition of 5-lipoxygenase will decrease the biosynthesis of LTB4 (C02165), cysteinylcontaining leukotrienes LTC4 (C02166), LTD4 (C05951) and LTE4 (C05952) (see, for instance, Fig. 4) [26]. According to our graph model, the removal of E1.13.11.34 will eliminate four compounds (LTB4, LTC4, LTD4 and LTE4) in arachidonic acid metabolism network. Furthermore, these four compounds play an important role in the mechanisms of toxic brain damage in acute methanol poisoning in humans [27]. Thus, we chose them as the target compounds. Apart from that, inhibiting this enzyme also eliminates four more compounds (i.e., 5(S)-HPETE(C05356), 5-HETE(C04805), LTA4(C00909), and 20-OH-LTB4(C04853)) [28]. These compounds can be regarded as “damage” in our model. Furthermore, when applying ILP-ECI with LTB4, LTC4, LTD4 and LTE4 as the target compounds, we find LTA4H (E3.3.2.6) and LTC4 synthase (E4.4.1.20) as the optimal enzyme set since the inhibition of these two enzymes eliminates only one non-target compound namely 20-OH-LTB4 (C04853). Thus, it is shown that ILP-ECI potentially finds a better solution in this experiment than the existing drugs since the same target compounds are eliminated by the existing drug in addition to other four compounds. Indeed, recent research validated our model since the anti-inflammatory effect of the levels of LTA4H [29] and LTC4 [30] has been observed. The computational time in this experiment takes only 5.19 s.

Fig. 4
figure 4

Subgraph of arachidonic acid metabolism

Another experiment we conducted is the histidine metabolism network (hsa00340). The enzyme amine oxidase (E. appeared in hsa00340 can be inhibited by drug (Rasagiline (D02562)). According to our graph model, the removal of enzyme E. will result in eliminating two compounds Methylimidazoleacetic acid (C05828) and Methylimidazole acetaldehyde (C05827). It should be noted that the level of pros-methylimidazoleacetic acid is closely related to severity of Parkinson disease in patients [31, 32]. Running ILP-ECI with Methylimidazoleacetic acid and Methylimidazole acetaldehyde as the target compounds finds amine oxidase as the optimal enzyme. It takes only 3.26 s to run ILP-ECI. This experiment verifies that Rasgiline targets the optimal enzyme. An important advantage of our Boolean model is its capability of detecting the lack of substrates where the connectivity-based methods fail to handle this. Another advantage of this model is its capability of handling branches and cycles in a pathway from the source compound to the target compound. However, there are still have some limitations in this method. One of the major drawbacks is the assumption that all the compounds and enzymes are of equal importance. However, different compounds and enzymes may have different levels of importance in the metabolic networks. Our future work will focus on developing other models which include the weight of different nodes.


In this paper, we formulate the optimal enzyme combination identification (ECI) problem as an optimization problem in Boolean metabolic networks. We have proven that ECI in the Boolean model is NP-complete and the target enzyme set is uniquely determined when the target compounds are given. Furthermore, considering that an exhaustive search cannot be used to solve ECI when the network is large, we developed an ILP-based algorithm for ECI. Considering that the computational time of IP-based method is exponential to the number of variables, to improve the scalability of the developed method, it is vital to reduce the number of variables appearing in IP formalization. And our proposed IP-based method needs O(m+n) variables.

The efficiency and effectiveness are validated by the computational experiments in which datasets were downloaded from the KEGG database. The results demonstrate that the proposed model can accurately identify the target enzymes for known successful drugs in the literature. Specifically, ILP-ECI has found a different enzyme set for the target compounds of Benoxaprofen which indicates that our method has a great potential to be better than Benoxaprofen. The reason is that the solution of our algorithm damages only one non-target compound while Benxaprofen damages 4 non-target compounds including the compound damaged by ILP-ECI’s solution. Besides, ILP-ECI has found that the same target enzyme like Rasagiline when its target compounds are given in advance. It is to be noted that our problem needs only O(m+n) variables in the IP formalization. The experiments also show that ILP-ECI can solve the problem in a short time which confirms the efficiency of our algorithm.


  1. Bro C, Regenberg B, Förster J, Nielsen J. In silico aided metabolic engineering of saccharomyces cerevisiae for improved biothanol production. Metab Eng. 2006; 8(2):102–11.

    Article  CAS  PubMed  Google Scholar 

  2. Lee SK, Chou H, Ham TS, Lee TS, Keasling JD. Metabolic engineering of microorganisms for biofuels production: from bugs to synthetic biology to fuels. Curr Opin Biotechnol. 2008; 19:556–63.

    Article  CAS  PubMed  Google Scholar 

  3. Alper H, Jin YS, Moxley JF, Stephanopoulos G. Identifying gene tarets for the metebolic engineering of lycopene biosynthesis in escherichia coli. Metab Eng. 2005; 7:155–64.

    Article  CAS  PubMed  Google Scholar 

  4. Soh KC, Hatzimanikatis V. Dreams of metabolism. Trends Biotechnol. 2010; 28:501–8.

    Article  CAS  PubMed  Google Scholar 

  5. Handorf T, Ebenhöh O, Heinrich R. Expanding metabolic networks: scopes of compounds, robustness, and evolution. J Mol Evol. 2005; 61:498–512.

    Article  CAS  PubMed  Google Scholar 

  6. Smart AG, Amaral LA, Ottino JM. Cascading failure and robustness in metabolic networks. Proc Natl Acad Sci. 2008; 105:13223–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Lemke N, Herédia F, Barcellos CK, Dos Reis AN, Mombach JC. Essentiality and damage in metabolic networks. Bioinformatics. 2004; 20:115–9.

    Article  CAS  PubMed  Google Scholar 

  8. Tamura T, Takemoto K, Akutsu T. Finding minimum reation cuts of metabolic networks under a boolean model using integer programming and feedback vertex sets. Int J Knowl Discov Bioinformatics (IJKDB). 2010; 1:14–31.

    Article  Google Scholar 

  9. Tamura T, Akutsu T. Exact algorithms for fidning a minimum reaction cut under a boolean model of metabolic networks. IEICE Trans Fundam Electron. 2010; 93:1497–1507.

    Article  Google Scholar 

  10. Li Z, Wang RS, Zhang XS, Chen L. Detecting drug targets with minimum side effects in metabolic networks. IET Syst Biol. 2009; 3:523–33.

    Article  CAS  PubMed  Google Scholar 

  11. Takemoto K, Tamura T, Akutsu T. Theoretical estimation of metabolic network robustness against multiple reaction knockouts using branching process approximation. Physica A Stat Mech Appl. 2003; 392:5525–535.

    Article  Google Scholar 

  12. Lee D, Goh KI, Kahng B. Branching process approach for boolean bipartite networks of metabolic reations. Phys Rev E. 2012; 86:027101.

    Article  Google Scholar 

  13. Lu W, Tamura T, Song J, Akutsu T. Integer programming-based method for designing synthetic metabolic networks by minimum reation insertion in a boolean model. Plos ONE. 2014; 9:92637.

    Article  Google Scholar 

  14. Drews J. Drug discovery: A historical perpective. Science. 2000; 287:1960–4.

    Article  CAS  PubMed  Google Scholar 

  15. Smith C. Hitting the target. Nature. 2003; 422:341–7.

    Article  Google Scholar 

  16. Takenaka T. Classical vs reverse pharmacology in drug discovery. BJU Int. 2001; 88:7–10.

    Article  CAS  PubMed  Google Scholar 

  17. Li Z, Wang RS, Zhang XS. Two-stage flux balance analysis of metabolic networks for drug target identification. BMC Syst Biol. 2011; 5:11.

    Article  Google Scholar 

  18. Surtees R, Blau N. The neurochemistry of phenylketonuria. Eur J Pediatr. 2000; 159:109–13.

    Article  Google Scholar 

  19. Sridhar P, Song B, Kahveci T, Ranka S. Mining metabolic networks for optimal drug targets. Pac Symp Biocomput. 2008; 13:291–302.

    Google Scholar 

  20. Schrijver A. Theory of Linear and Integer Programming. New York: John Wiley & Sons, Inc; 1986.

    Google Scholar 

  21. Zhenping L, Zhang S, Wang Y, Zhang XS, Chen L. Alignment of molecular networks by integer quadratic programming. Bioinformatics. 2007; 23:1631–9.

    Article  PubMed  Google Scholar 

  22. Cheng X, Qiu Y, Hou W, Ching WK. Integer programming-based method for observability of singleton attractors in boolean networks. IET Syst Biol. 2017; 11:30–5.

    Article  PubMed  Google Scholar 

  23. IBM 2010, IBM ILOG CPLEX Optimizer. Accessed 2 Mar 2017.

  24. Kanehisa M, Goto S. Kegg: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000; 28:27–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Fogh K, Herlin T, Kragballe K. In vitro inhibition of leukotriene b4 formation by exogeneous 5-lipoxygenase inhibitors is associated with enhanced generation of 15-hydroxy-eicosatetraenoic acid (15-hete) by human neutrophils. Arch Dermatol Res. 1988; 280:430–6.

    Article  CAS  PubMed  Google Scholar 

  26. Lötzer K, Funk CD, Habenicht AJ. The 5-lipoxygenase pathway in arterial wall biology and atherosclerosis. Biochim Biophys Acta. 2005; 1736:30–7.

    PubMed  Google Scholar 

  27. Zakharov S, Kotikova K, Nurieva O, Hlusicka J, Kacer P, Urban P, Vaneckova M, Seidl Z, Diblik P, Kuthan P, Navratil T, Pelclova D. Leukotriene-mediated neuroinflammation, toxic brain damage, and neurodegeneration in acute methanol poisoning. lin Toxicol (Phila). 2017; 55:249–59.

    Article  CAS  Google Scholar 

  28. Sweeney FJ, Eskra JD, Carty TJ. Development of a system for evaluating 5-lipoxygenase inhibitors using human whole blood. Prostaglandins Leukot Med. 1987; 28:73–93.

    Article  CAS  PubMed  Google Scholar 

  29. Rao NL, Dunford PJ, Xue X, Jiang X, Lundeen KA, Coles F, Riley JP, Williams KN, Grice CA, Edwards JP, Karlsson L, Fourie AM. Anti-inflammatory activity of a potent, selective leukotriene a4 hydrolase inhibitor in comparison with the 5-lipoxygenase inhibitor zileuton. J Pharmacol Exp Ther. 2007; 321:1154–60.

    Article  CAS  PubMed  Google Scholar 

  30. Torres-Galván MJ1, Ortega N, Sánchez-García F, Blanco C, Carrillo T, Quiralte J. Ltc4-synthase a-444c polymorphism: lack of association with nsaid-induced isolated periorbital angioedema in a spanish population. Ann Allergy Asthma Immunol. 2001; 87:506–10.

    Article  PubMed  Google Scholar 

  31. Blandina P, Cherici G, Moroni F, Prell GD, Green JP. Release of glutamate from striatum of freely moving rats by pros-methylimidazoleacetic acid. J Neurochem. 1995; 64:788–93.

    Article  CAS  PubMed  Google Scholar 

  32. Prell GD, Khandelwal JK, Burns RS, Blandina P, Morrishow AM, Green JP. Levels of pros-methylimidazoleacetic acid: Correlation with severity of parkinson’s disease in csf of patients and with the depletion of striatal dopamine and its metabolites in mptp-treated mice. J Neural Transm. 1991; 3:1435–63.

    Google Scholar 

Download references


The authors would like to thank anonymous reviewers for your helpful and constructive comments.


This research is supported in part of Natural Science Foundation of SZU (Grant No. 2017058) and National Natural Science Foundation of China NSFC (Grant Nos. 91730301, 11626229 and 11671158). The publication costs are funded by National Natural Science Foundation of China NSFC (Grant Nos. 91730301).

Availability of data and materials

All the data sets are publicly available and can be accessed from the KEGG databases.

About this supplement

This article has been published as part of BMC Systems Biology Volume 12 Supplement 1, 2018: Selected articles from the 16th Asia Pacific Bioinformatics Conference (APBC 2018): systems biology. The full contents of the supplement are available online at

Author information

Authors and Affiliations



YQ designed the research. YQ and WKC proposed the methods and did theoretical analysis. YQ and HJ collected the data. YQ, HJ and XC conducted the experiments and analyze the results. YQ, HJ, WC and XC wrote the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Hao Jiang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Qiu, Y., Jiang, H., Ching, WK. et al. Discovery of Boolean metabolic networks: integer linear programming based approach. BMC Syst Biol 12 (Suppl 1), 7 (2018).

Download citation

  • Published:

  • DOI: