Discovery of Boolean metabolic networks: integer linear programming based approach

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 [1–3]. 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 [20–22] 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.

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, R₁ is a reaction, its substrates are C₁ and E₂ and its product is C₂. Here C₅ is the target compound in this figure which shall be stopped. To stop the production of C₅, reaction R₂ must be prevented from taking place. One of the possible solution is to disrupt one of its catalyzing enzymes (E₁ for instance). Another is to stop the production of its reactant compounds (i.e., C₂ or C₃). If C₂ is stopped, then two possible ways can achieve this effect (i.e., enzyme E₂ or E₁). It is noted that C₂,R₂,C₇ and R₁ forms a cycle, the inhibition of E₁ will result in inactivating reaction R₂ which will further makes C₂ inproducible. The other way is inhibiting enzyme E₂ which makes R₁ inactive and further stop the production of compound C₂. Therefore, the inhibition of enzyme E₁ or E₂ can result in stopping the production of the target compound.

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 E₁ results in the knock out of compounds C₄,C7,C₈ and C₉ in addition to the target compound C₅. 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 E₂ is 2 (i.e., C₂ and C₄). Compound C₇ is still producible because it can be produced by R₃ even after the inhibition of E₂. The damage effect of inhibition of E₁ is 4 (i.e., C₄,C7,C₈ and C₉). Both E₁ and E₂ are potential drug targets since they can achieve the effect of disrupting the target compound C₅. However, E₂ is a better drug target than E₁ 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.

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}}\).

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).

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 {E₂}. 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.

Proof

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, E₁={e₂,e₃,e₆} is one of optimal solutions of MEC for graph shown in Fig. 2. Let G=(V,E) be an instance of MEC, where V={v₁,v₂,⋯,v _n } and E={e₁,e₂,⋯,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} $$

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. □

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.

Another experiment we conducted is the histidine metabolism network (hsa00340). The enzyme amine oxidase (E.1.4.3.4) appeared in hsa00340 can be inhibited by drug (Rasagiline (D02562)). According to our graph model, the removal of enzyme E.1.4.3.4 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.